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 "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 = true;
      bool use_threads = true;
      double stabilization_parameter = 1.;
      double viscosity = 1.;
      double permeability = 0.;
      std::ostream & output = std::cout;
    };
 
    class HYPRE {
    public:
      typedef Eigen::SparseMatrix<double> SystemMatrixType;
 
      typedef std::function<VectorRd(const VectorRd &)> MomentumForcingTermType;
      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 numStaticallyCondensedDofsVelocity() + numStaticallyCondensedDofsPressure();
      }
 
      size_t numDirichletDofs() const
      {
        return 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 dimVelocity() + dimPressure() - numStaticallyCondensedDofs() + 1;
      }
 
      size_t sizeSystem() const
      {
        return numNonStaticallyCondensedDofs() - numDirichletDofs();
      }
 
      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 double time_step,
                                  const double current_time);
 
      void assembleResidualSystem(
                                  const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                  const Eigen::VectorXd & uph_old) {
    Eigen::VectorXd uph_prev = Eigen::VectorXd::Zero(uph_old.size());
    assembleResidualSystem(isolution,
                   uph_old,
                   uph_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 Eigen::VectorXd uT_old,
      //                                           const double current_time = 0.) const {
      //   // Construct quadrature scheme on element
      //   const Cell & T = *m_velocity_space->mesh().cell(iT);
 
      //   // 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.
      //   QuadratureRule quad = generate_quadrature_rule(T, 3 * cell_degree);
      //   const auto & cell_dofs = m_velocity_space->get<PolynCellType>("CellDOFs")[iT];
      //   auto cell_dofs_quad = evaluate_quad<Function>::compute(*cell_dofs, quad);
      //   double beta_T = 1.e-4;                                    // Safeguard parameter
      //   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);
      //     beta_T = std::max(beta_T, uT_old_iqn.lpNorm<Eigen::Infinity>());
      //   } // 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;
      // };
      
      virtual double convective_stabilization_parameter(size_t iT,
                                                        const HArDCore2D::NavierStokesSolutions::IExactSolution * isolution,
                                                        const size_t cell_degree,
                                                        const double viscosity,
                                                        const Eigen::VectorXd uT_old,
                                                        const double current_time = 0.) 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);
      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> _darcy_term(size_t iT, const Eigen::VectorXd & uT_old);
 
      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::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.);
      
 
      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 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) {
    Eigen::VectorXd uph_prev = uph_old;
    return _compute_local_contribution(iT,
                       isolution,
                       uph_old,
                       uph_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;
      double m_stabilization_parameter;
      double m_viscosity;
      double m_permeability;
      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