hho-nsa.hpp

Go to the documentation of this file.

// Implementation of the HHO scheme, with full gradient, for the Navier-Stokes variable density
 
 
#ifndef HHONSA_HPP
#define HHONSA_HPP
 
#include <iostream>
 
#include <boost/math/constants/constants.hpp>
 
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
 
#include <mesh.hpp>
// #include <mesh_builder.hpp>
#include <local_static_condensation.hpp>
#include <linearsolver.hpp>
 
#include <vhhospace.hpp>
#include <BoundaryConditions/BoundaryConditions.hpp>
#include <integralweight.hpp>
 
 
namespace HArDCore2D
{
  //------------------------------------------------------------------------------
  // Navier-Stokes model
  //------------------------------------------------------------------------------
 
 
  struct VarDensityNavierStokes
  {
    typedef Eigen::SparseMatrix<double> systemMatrixType;
    
    typedef std::function<VectorRd(const VectorRd &)> MomentumForcingTermType;
    typedef std::function<double(const VectorRd &)> CompressibilityForcingTermType;
    typedef std::function<double(const VectorRd &)> DensityForcingTermType;
    
    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;
    typedef std::function<double(const VectorRd &)> VolumicFractionType;
    
    typedef IntegralWeight ViscosityType;
 
    VarDensityNavierStokes(
           const VHHOSpace & vhho_space,      
           const GlobalDOFSpace & p_space,    
           const GlobalDOFSpace & phi_space,    
           const BoundaryConditions & BC,     
           bool use_threads,                  
           bool stokes = false,               
           std::ostream & output = std::cout  
           );
 
    void assembleLinearHHOSystem(
      const MomentumForcingTermType & f,        
      const CompressibilityForcingTermType & g, 
      const VolumicFractionType & rhog,         
      const ViscosityType & mu,                 
      const Eigen::VectorXd & u0,               
      const Eigen::VectorXd & uold,             
      const double & time_step                  
      );
      
    void assembleLinearDGSystem(
                                const CompressibilityForcingTermType & h,
                                const Eigen::VectorXd & uprho,
                                const VolumicFractionType & rho,
                                const double & time_step
                                );
                    
    inline size_t nloc_sc_u() const
    {
      return m_nloc_sc_u;
    }
 
    inline size_t nloc_sc_p() const
    {
      return m_nloc_sc_p;
    }
 
    inline size_t numSCDOFs_u() const
    {
      return m_vhhospace.mesh().n_cells() * m_nloc_sc_u;
    }
 
    inline size_t numSCDOFs_p() const
    {
      return m_vhhospace.mesh().n_cells() * m_nloc_sc_p;
    }
 
    inline size_t numSCDOFs() const
    {
      return numSCDOFs_u() + numSCDOFs_p();
    }
 
    inline size_t numDirDOFs() const
    {
      return m_BC.n_dir_edges()*m_vhhospace.numLocalDofsEdge();
    }
    
    inline size_t dimVelocity() const
    {
      return m_vhhospace.dimension();
    }
 
    inline size_t dimPressure() const
    {
      return m_pspace.dimension();
    }
    
    inline size_t dimVolumicFraction() const
    {
      return m_phispace.dimension();
    }
    
    inline size_t numNonSCDOFs() const
    {
      return dimVelocity() + dimPressure() - numSCDOFs() + 1;
    }
 
    inline size_t sizeSystemHHO() const
    {
      return numNonSCDOFs() - numDirDOFs();
    }
    
    inline size_t sizeSystemDG() const
    {
      return dimVolumicFraction();
    }
 
    inline const VHHOSpace & vhhospace() const
    {
      return m_vhhospace;
    }
    
    inline const GlobalDOFSpace & pspace() const
    {
      return m_pspace;
    }
    
    inline const GlobalDOFSpace & phispace() const
    {
      return m_phispace;
    }
 
    inline const Mesh & mesh() const
    {
      return m_vhhospace.mesh();
    }
    
    inline const systemMatrixType & systemMatrixHHO() const {
      return m_A_hho;
    }
 
    inline systemMatrixType & systemMatrixHHO() {
      return m_A_hho;
    }
    
    inline const Eigen::VectorXd & systemVector() const {
      return m_b_hho;
    }
 
    inline Eigen::VectorXd & systemVector() {
      return m_b_hho;
    }
 
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    inline const systemMatrixType & scMatrix() const {
      return m_sc_A_hho;
    }
 
    inline Eigen::VectorXd & scVector() {
      return m_sc_b_hho;
    }
    
    inline const systemMatrixType & systemMatrixDG() const {
      return m_A_dg;
    }
      
    inline systemMatrixType & systemMatrixDG() {
      return m_A_dg;
    }
     
    inline const Eigen::VectorXd & systemVectorDG() const {
      return m_b_dg;
    }
      
    inline Eigen::VectorXd & systemVectorDG() {
      return m_b_dg;
    }
    
    inline bool & isStokes() { 
      return m_stokes; 
    }
 
    Eigen::VectorXd interpolate(
                                const VelocityType & u, 
                                const PressureType & p, 
                                const VolumicFractionType & phi,
                                const int doe_cell = -1, 
                                const int doe_face = -1 
                                ) const;    
 
    std::vector<double> computeEnergyNorms(
                                           const std::vector<Eigen::VectorXd> & list_dofs   
                                           ) const;
                                           
    std::vector<double> computePressureL2Norm(
                                           const std::vector<Eigen::VectorXd> & list_dofs   
                                           ) const;
    
    std::vector<std::tuple<double,double,double>> computeDGNorms(const std::vector<Eigen::VectorXd> & list_dofs, const VelocityType & u) const;
    
    // Add to hho-nsa.hpp in the public section:
std::vector<double> computePhysicalEnergies(
                       const Eigen::VectorXd & uprho,           
                       const MomentumForcingTermType & f,       
                       const ViscosityType & mu                 
                       ) const;
                                           
    Eigen::VectorXd pressureVertexValues(
                       const Eigen::VectorXd & p    
                       ) const;
                       
    Eigen::VectorXd volumicFractionVertexValues(
                                                const Eigen::VectorXd & phi 
                                                ) const;
    
    void writeToVtuFile(
                           const Eigen::VectorXd & upphi,
                           const Eigen::VectorXd & upphiI,
                           const std::unique_ptr<Mesh> & mesh_ptr,
                           const std::string & solution_name,
                           const std::string & mesh_name,
                           const size_t & degree,
                           const int & time_it,
                           const double & reynolds,
                           const std::vector<std::string> & plot_variables,
                           const std::string & output_path
                           );
                           
    // Write data to file 
    void writeTimeStepData(
                           const std::string & filename,
                           const int & time_it,
                           const double & time,
                           const double & dt,
                           const double & L2_u_err,
                           const double & H1_u_err,
                           const double & Energy_err,
                           const double & L2_rho_err,
                           const double & cf_rho_err,
                           const double & upw_rho_err,
                           const double & E_kinetic,
                           const double & E_potential,
                           const double & E_dissipation
                           ) const;
                           
    void writeFileHeader(
                         const std::string & filename,
                         const std::string & solution_name,
                         const int & solution_number,
                         const std::string & mesh_name,
                         const size_t & degree,
                         const std::unique_ptr<Mesh> & mesh_ptr,
                         const double & viscosity
                         ) const;
    
    Eigen::VectorXd imposingBC(
                                const VelocityType & u,
                                const Eigen::VectorXd & uph
                                ) const;
    double getUpwindDensity(
                            const size_t & iE,
                            const size_t & iT,
                            const Eigen::VectorXd & uprho
                            ) const;
    Eigen::VectorXd newtonRaphson(
                                   LinearSolver<VarDensityNavierStokes::systemMatrixType> & solver,
                                   const std::unique_ptr<Mesh> & mesh_ptr,
                                   const size_t & newton_maxit,
                                   const double & newton_tol,
                                   const Eigen::VectorXd & u_newt,
                                   const Eigen::VectorXd & un,
                                   const MomentumForcingTermType & f,
                                   const CompressibilityForcingTermType & g,
                                   const VolumicFractionType & rhog,
                                   const ViscosityType & mu,
                                   const double & time_step
                                   );
    
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                            
                                const MomentumForcingTermType & f,      
                                const CompressibilityForcingTermType & g,    
                                const VolumicFractionType & rhog,
                                const ViscosityType & mu,              
                                const Eigen::VectorXd & u0,
                                const Eigen::VectorXd & uold,
                                const double & time_step
                                );
 
    LocalStaticCondensation _compute_static_condensation(const size_t & iT) const;
 
    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                                     
                                      );
                                      
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _compute_local_cell_dG_contribution(
           size_t iT,
           const DensityForcingTermType & h,
           const Eigen::VectorXd & uprho,
           const double & time_step
           );
      
      Eigen::MatrixXd _compute_local_iface_dG_contribution(
           size_t iF,     // local edge index
           size_t iT1,    // global element index
           size_t iT2,
           const Eigen::VectorXd & uprho
           );
      
      std::pair<Eigen::MatrixXd, Eigen::VectorXd> _compute_local_bface_dG_contribution(
           size_t iF,
           size_t iTb,
           const Eigen::VectorXd & uprho,
           const VolumicFractionType & rho
           );
    
    const VHHOSpace & m_vhhospace;
    const GlobalDOFSpace & m_pspace; 
    const GlobalDOFSpace & m_phispace;
    const BoundaryConditions & m_BC;
    bool m_use_threads;
    bool m_stokes;
    std::ostream & m_output;
    
    const size_t m_nloc_sc_u; // Number of statically condensed velocity DOFs in each cell
    const size_t m_nloc_sc_p; // Number of statically condensed pressure DOFs in each cell
    
    systemMatrixType m_A_hho;     // Matrix and RHS of statically condensed system
    Eigen::VectorXd m_b_hho;
    systemMatrixType m_sc_A_hho;  // Static condensation operator and RHS (to recover statically condensed DOFs)
    Eigen::VectorXd m_sc_b_hho;
    
    systemMatrixType m_A_dg;     // Matrix and RHS of dG system
    Eigen::VectorXd m_b_dg;
      
    double m_stab_par;
  };
 
} // end of namespace HArDCore3D
 
#endif