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 "mhd-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::MHDSolutions::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);
 
 
      //This should be the stationary case assembleResidual case!
 
      void assembleResidualSystem(
                                  const HArDCore2D::MHDSolutions::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;
     
 
      virtual double convective_stabilization_parameter(size_t iT,
                                                     const HArDCore2D::MHDSolutions::IExactSolution * isolution,
                                               const Eigen::VectorXd & uT_old,
                                               const Eigen::VectorXd & bT_old,
                                                     const size_t cell_degree,
                                                     const double viscosity,
                                                     const double magnetic_diffusivity,
                                                     const double current_time) const;  
 
      virtual double magnetic_convective_stabilization_parameter(size_t iT,
                                                     const HArDCore2D::MHDSolutions::IExactSolution * isolution,
                                               const Eigen::VectorXd & uT_old,
                                               const Eigen::VectorXd & bT_old,
                                                     const size_t cell_degree,
                                                     const double viscosity,
                                                     const double magnetic_diffusivity,
                                                     const double current_time) const;  
                                               
                                                                                  
    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::MHDSolutions::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::MHDSolutions::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::MHDSolutions::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 Eigen::VectorXd & bT_old,
                                      const HArDCore2D::MHDSolutions::IExactSolution * isolution,
                                      const double current_time);
 
      virtual std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _magnetic_convective_stabilization_term(
                                              size_t iT, 
                                              const Eigen::VectorXd & uT_old,
                                              const Eigen::VectorXd & bT_old,
                                              const HArDCore2D::MHDSolutions::IExactSolution * isolution,
                                              const double current_time);                                                                                                            
      
      std::pair<Eigen::MatrixXd, Eigen::VectorXd>
      _compute_local_contribution(
                                  size_t iT,
                                  const HArDCore2D::MHDSolutions::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::MHDSolutions::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