hypre.hpp

Go to the documentation of this file.

#ifndef HYPRE_HPP
#define HYPRE_HPP
 
#include <ostream>
#include <Eigen/Sparse>
 
#include <BoundaryConditions/BoundaryConditions.hpp>
 
#include <local_static_condensation.hpp>
 
#include <unordered_set>
 
 
#include "discrete-space.hpp"
#include "local-space-dimensions.hpp"
#include "hho-interpolate.hpp"
#include "ns-solutions.hpp"
 
namespace HArDCore2D {  
  
  namespace DSL {
 
    struct HYPREParameters {
      bool static_condensation = false; // true;
      bool use_threads = true;
      bool upwinding = true;
      double stabilization_parameter= 1.0;
      double viscosity = 1.;
      double magnetic_diffusivity = 1.;
      std::ostream & output = std::cout;
    };
 
    class HYPRE {
    public:
      typedef Eigen::SparseMatrix<double> SystemMatrixType;
 
      typedef std::function<VectorRd(const VectorRd &)> MomentumForcingTermType;
      typedef std::function<VectorRd(const VectorRd &)> MagneticForcingTermType;
      typedef std::function<double(const VectorRd &)> CompressibilityForcingTermType;
 
      typedef std::function<VectorRd(const VectorRd &)> VelocityType;
      typedef std::function<MatrixRd(const VectorRd &)> VelocityGradientType;
 
      typedef std::function<double(const VectorRd &)> PressureType;
      typedef std::function<VectorRd(const VectorRd &)> PressureGradientType;
 
 
      HYPRE(
            const Mesh & mesh,
            size_t degree,
            const BoundaryConditions & bc,
            const HYPREParameters & parameters = {}
            );
      virtual ~HYPRE() {}
      
      virtual Eigen::VectorXd interpolate(
                      const VelocityType & u, 
                      const PressureType & p  
                      ) const = 0;
  
      template<typename F>
      Eigen::VectorXd pkpoPkPkPkInterpolate(
                                            const F & v,
                                            const InterpolateParameters & parameters = {}
                                            ) const;      
 
      size_t degree() const
      {
        return m_degree;
      }
 
      double stabilizationParameter() const
      {
        return m_stabilization_parameter;
      }
 
      const std::vector<Eigen::MatrixXd> potential() const
      {
        return m_potential;
      }
 
      const std::vector<Eigen::MatrixXd> stabilization() const
      {
        return m_stabilization;
      }
 
      const Eigen::MatrixXd & potential(size_t iT) const
      {
        return m_potential[iT];
      }
 
      size_t numLocalStaticallyCondensedDofsVelocity() const
      {
        return m_nloc_sc_u;
      }
 
      size_t numLocalStaticallyCondensedDofsPressure() const
      {
        return m_nloc_sc_p;
      }
 
      size_t numStaticallyCondensedDofsVelocity() const
      {
        return m_velocity_space->mesh().n_cells() * m_nloc_sc_u;
      }
 
      size_t numStaticallyCondensedDofsPressure() const
      {
        return m_velocity_space->mesh().n_cells() * m_nloc_sc_p;
      }
    
      size_t numStaticallyCondensedDofs() const
      {
        return 2*(numStaticallyCondensedDofsVelocity() +numStaticallyCondensedDofsPressure());
      }
 
      size_t numDirichletDofs() const
      {
        return 2*(m_bc.n_dir_edges() * m_velocity_space->numberOfLocalEdgeDofs());
      }
 
      size_t dimVelocity() const
      {
        return m_velocity_space->dimension();
      }
 
      size_t dimPressure() const
      {
        return m_pressure_space->dimension();
      }
    
      size_t numNonStaticallyCondensedDofs() const
      {
        return  2*dimVelocity() + 2*dimPressure() + 2 - numStaticallyCondensedDofs();
      }
 
      size_t sizeSystem() const
      {
        return numNonStaticallyCondensedDofs() - numDirichletDofs();
      }
 
      size_t numNonStaticallyCondensedDofsVelocity() const
      {
        return dimVelocity() - numStaticallyCondensedDofsVelocity();
      }
 
      size_t numNonStaticallyCondensedDofsPressure() const
      {
        return dimPressure() - numStaticallyCondensedDofsPressure();
      }
 
 
      inline const DiscreteSpace & velocitySpace() const
      {
        return *m_velocity_space;
      }
 
      inline const DiscreteSpace & pressureSpace() const
      {
        return *m_pressure_space;
      }
 
      inline const Mesh & mesh() const
      {
        return m_velocity_space->mesh();
      }
 
      inline const SystemMatrixType & systemMatrix() const {
        return m_A;
      }
 
      inline SystemMatrixType & systemMatrix() {
        return m_A;
      }
 
      inline const Eigen::VectorXd & systemVector() const {
        return m_b;
      }
 
      inline Eigen::VectorXd & systemVector() {
        return m_b;
      }
    
      inline const SystemMatrixType & staticallyCondensedMatrix() const {
        return m_sc_A;
      }
 
      inline const Eigen::VectorXd & staticallyCondensedVector() const {
        return m_sc_b;
      }
 
      inline const BoundaryConditions & boundaryConditions() const {
        return m_bc;
      }
 
      void assembleResidualSystem(
                                  const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                  const Eigen::VectorXd & uph_old,
                                  const Eigen::VectorXd & uph_prev,
                                  const Eigen::VectorXd & brh_old,
                                  const Eigen::VectorXd & brh_prev,
                                  const double time_step,
                                  const double current_time);
 
 
      //[Vito] This should be the stationary case assembleResidual case!
 
      void assembleResidualSystem(
                                  const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                  const Eigen::VectorXd & uph_old,
                                  const Eigen::VectorXd & brh_old) {
                                    Eigen::VectorXd uph_prev = Eigen::VectorXd::Zero(uph_old.size());
                                    Eigen::VectorXd brh_prev = Eigen::VectorXd::Zero(brh_old.size());
                                    assembleResidualSystem(isolution,
                                                           uph_old,
                                                           uph_prev,
                                                           brh_old,
                                                           brh_prev,
                                                           std::numeric_limits<double>::infinity(),
                                                           0.);
                                  };
 
      Eigen::VectorXd reconstructSolution(const Eigen::VectorXd & uhp_solsystem) const;
 
     /* double convective_stabilization_parameter(size_t iT,
                                                const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                const size_t cell_degree,
                                                const double viscosity,
                                                const double magnetic_diffusivity,
                                                const double current_time) const {
        // Construct quadrature scheme on element
        const Cell & T = *m_velocity_space->mesh().cell(iT);
        QuadratureRule quad = generate_quadrature_rule(T, 3 * cell_degree);
 
 
        // Local convective stabilization bilinear form assumed to be given by
        // j_T(w_T, v_T) = Beta_T Sum_{F in F_T} Int_{F} (w_F - w_T).(v_F - v_T)
        // where Beta_T = max(c_s, ||u_T||_{L^infty(T)}). We therefore need an
        // estimate of ||u_T||_{L^infty(T)} which we obtain by taking the largest
        // l^infty norm of u_T on the quadrature points on T.
        double beta_T = 0;                                  
 
        if (m_upwinding){
          double c_s =1.e-4;                                    // Safeguard parameter
          double diff = std::max(viscosity, magnetic_diffusivity);
 
          for (size_t iqn = 0; iqn < quad.size(); iqn++) {
            // Stabilisation based on max-norm of discrete solution
            // VectorRd uT_old_iqn = VectorRd::Zero();
            // for (size_t i = 0; i < cell_dofs->dimension(); i++) {
            //   double uT_old_i = uT_old(m_velocity_space->localOffset(T) + i);
            //   uT_old_iqn += uT_old_i * cell_dofs_quad[i][iqn];
            // } // for i
            // Stabilisation based on max-norm of analytic solution
          
 
            Eigen::VectorXd uT_old_iqn = isolution->velocity(quad[iqn].vector(), current_time);
            Eigen::VectorXd bT_old_iqn = isolution->magnetic_field(quad[iqn].vector(), current_time);
 
            double local_norm =    std::max(uT_old_iqn.lpNorm<Eigen::Infinity>(), bT_old_iqn.lpNorm<Eigen::Infinity>());
            if (local_norm * T.diam() / diff > 1.0){
                beta_T = std::max(beta_T, std::max(local_norm, c_s));
            }
 
          } // for iqn
        }  
 
        // ALTERNATIVE APPROACH (based on maximum absolute value of u.n on the boundary of T)
        // double beta_T = 1.e-4;                                    // Safeguard parameter
        // for (size_t iE = 0; iE < T.n_edges(); iE++) {
        //   const Edge & E = *T.edge(iE);
        //   auto nTE = T.edge_normal(iE);
        //   QuadratureRule quad_E = generate_quadrature_rule(E, 3 * cell_degree);
        //   for (size_t iqn = 0; iqn < quad_E.size(); iqn++) {
        //     Eigen::VectorXd uE_old_iqn = isolution->velocity(quad_E[iqn].vector());
        //     beta_T = std::max(beta_T, std::abs(uE_old_iqn.dot(nTE)));
        //   }
        // }
 
        return beta_T;
      };*/
 
      
            double convective_stabilization_parameter(size_t iT,
                                                     const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                     const size_t cell_degree,
                                                     const double viscosity,
                                                     const double magnetic_diffusivity,
                                                     const double current_time = 0.) const {
                // Construct quadrature scheme on element
                const Cell & T = *m_velocity_space->mesh().cell(iT);
                QuadratureRule quad = generate_quadrature_rule(T, 3 * cell_degree);
 
                // Compute element diameter and maximum diffusion coefficient
                //const double h_T = T.diam();
                //const double diff = std::max(viscosity, magnetic_diffusivity);
 
                double beta_T = 0.0;
 
                if (m_upwinding) {
                    for (size_t iqn = 0; iqn < quad.size(); ++iqn) {
                        // Evaluate analytic velocity at quadrature point
            Eigen::VectorXd uT_old_iqn = isolution->velocity(quad[iqn].vector(), current_time);
            Eigen::VectorXd bT_old_iqn = isolution->magnetic_field(quad[iqn].vector(), current_time);
            double local_norm = std::max(uT_old_iqn.lpNorm<Eigen::Infinity>(), bT_old_iqn.lpNorm<Eigen::Infinity>());
 
                        // Local Péclet number
                        //double pe_T = (local_norm * h_T) / (diff + 1e-12);
 
                        // Smooth tapering: 0 for Pe_T << 1, ~1 for Pe_T >> 1
                        //double taper = std::max(0.0, std::tanh( 0.5* (pe_T - 1.0)));
 
                        // Apply tapered stabilization
                        beta_T = std::max(beta_T, local_norm);
                    }
                }
 
        //std::cout<<" the value of beta_T is:"<< beta_T <<std::endl;
 
                return beta_T;
            };
 
  
    
    
      double magnetic_convective_stabilization_parameter(size_t iT,
                                                const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                const size_t cell_degree,
                                                const double viscosity,
                                                const double magnetic_diffusivity,
                                                const double current_time = 0. ) const {
                // Construct quadrature scheme on element
                const Cell & T = *m_velocity_space->mesh().cell(iT);
                QuadratureRule quad = generate_quadrature_rule(T, 3 * cell_degree);
 
                // Compute element diameter and maximum diffusion coefficient
                //const double h_T = T.diam();
                //const double diff = std::max(viscosity, magnetic_diffusivity);
 
                double gamma_T = 0.0;
 
                if (m_upwinding) {
                    for (size_t iqn = 0; iqn < quad.size(); ++iqn) {
                        // Evaluate analytic velocity at quadrature point
            Eigen::VectorXd uT_old_iqn = isolution->velocity(quad[iqn].vector(), current_time);
            Eigen::VectorXd bT_old_iqn = isolution->magnetic_field(quad[iqn].vector(), current_time);
            double local_norm = std::max(uT_old_iqn.lpNorm<Eigen::Infinity>(), bT_old_iqn.lpNorm<Eigen::Infinity>());
 
                        // Local Péclet number
                        //double pe_T = (local_norm * h_T) / (diff + 1e-12);
 
                        // Smooth tapering: 0 for Pe_T << 1, ~1 for Pe_T >> 1
                        //double taper = std::max(0.0, std::tanh(0.5 * (pe_T - 1.0)));
 
 
                        // Apply tapered stabilization
                        gamma_T = std::max(gamma_T, local_norm);
                    }
                }
 
                return gamma_T;
            };
     
 
 
 
 
 
    protected:
      void _construct_local_operators(size_t iT);
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _time_deriv_term(size_t iT,
                          const Eigen::VectorXd & uT_old,
                          const Eigen::VectorXd & uT_prev,
                          const double time_step); //  = std::numeric_limits<double>::infinity()
 
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _viscous_term(size_t iT, const Eigen::VectorXd & uT_old);
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _magnetic_diffusive_term(size_t iT, const Eigen::VectorXd & bT_old);
 
      
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd> _velocity_pressure_coupling(
                                              size_t iT,
                                              const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                              const Eigen::VectorXd & uT_old,
                                              const Eigen::VectorXd & pT_old,
                          const double & current_time = 0.
                                              );
 
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd> _magnetic_field_magnetic_pressure_coupling(
                                              size_t iT,
                                              const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                              const Eigen::VectorXd & bT_old,
                                              const Eigen::VectorXd & rT_old,
                          const double & current_time = 0.
                                              );                    
 
 
      virtual Eigen::MatrixXd _forcing_terms(
                         size_t iT,
                         const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
               const double current_time = 0.
                         );
 
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _linearized_convective_term(size_t iT, const Eigen::VectorXd & uT_old);
 
      virtual std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::VectorXd>
      _linearized_convective_mixed_term(size_t iT, const Eigen::VectorXd & wT_old, const Eigen::VectorXd & vT_old);
 
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _convective_stabilization_term(size_t iT,
                                     const Eigen::VectorXd & uT_old,
                                     const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                     const double current_time = 0.);
 
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _magnetic_convective_stabilization_term(size_t iT,
                                     const Eigen::VectorXd & bT_old,
                                     const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                     const double current_time = 0.);                               
      
      std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _compute_local_contribution(
                                  size_t iT,
                                  const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                  const Eigen::VectorXd & uph_old,
                                  const Eigen::VectorXd & uph_prev,
                                  const Eigen::VectorXd & brh_old,
                                  const Eigen::VectorXd & brh_prev,
                                  const double time_step,
                                  const double current_time);
 
      std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _compute_local_contribution(
                                  size_t iT,
                                  const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                  const Eigen::VectorXd & uph_old,
                                  const Eigen::VectorXd & brh_old) {
                                    Eigen::VectorXd uph_prev = uph_old;
                                    Eigen::VectorXd brh_prev = brh_old;
 
                                    return _compute_local_contribution(iT,
                                                                       isolution,
                                                                       uph_old,
                                                                       uph_prev,
                                                                       brh_old,
                                                                       brh_prev,
                                                                       std::numeric_limits<double>::infinity(),
                                                                       0.);
 
      }
      
      LocalStaticCondensation _compute_static_condensation(const size_t & iT);
      
      void _assemble_local_contribution(
                                        size_t iT,
                                        const std::pair<Eigen::MatrixXd, Eigen::VectorXd> & lsT,
                                        std::list<Eigen::Triplet<double> > & A1,
                                        Eigen::VectorXd & b1,
                                        std::list<Eigen::Triplet<double> > & A2,
                                        Eigen::VectorXd & b2
                                        );
      
      size_t  m_degree;
      const BoundaryConditions & m_bc;
      bool m_static_condensation;
      bool m_use_threads;
      bool m_upwinding;
 
      double m_stabilization_parameter;
      double m_viscosity;
      double m_magnetic_diffusivity;
 
      std::ostream & m_output;
 
      size_t m_nloc_sc_u;
      size_t m_nloc_sc_p;
      
      std::vector<Eigen::MatrixXd> m_discrete_l2;
      std::vector<Eigen::MatrixXd> m_gradient;
      std::vector<Eigen::MatrixXd> m_gradient_rhs;
      std::vector<Eigen::MatrixXd> m_potential;
      std::vector<Eigen::MatrixXd> m_stabilization;
 
      std::unique_ptr<DiscreteSpace> m_velocity_space;
      std::unique_ptr<DiscreteSpace> m_pressure_space;
 
 
      SystemMatrixType m_A;
      Eigen::VectorXd m_b;
 
      SystemMatrixType m_sc_A;
      Eigen::VectorXd m_sc_b;
      
      std::vector<int> m_nsc_to_dof_map;
      std::vector<int> m_sc_to_dof_map;
 
    };
 
  } // namespace DSL
  
} // namespace HArDCore2D
 
#endif