sddr-navier-stokes.hpp

Go to the documentation of this file.

// Implementation of the discrete de Rham sequence for the Navier-Stokes problem in curl-curl formulation.
//
// Author: Jerome Droniou (jerome.droniou@monash.edu)
//
 
 
#ifndef SDDR_STOKES_HPP
#define SDDR_STOKES_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 <BoundaryConditions/BoundaryConditions.hpp>
#include <BoundaryConditions/BChandlers.hpp>   
 
#include <sxgrad.hpp>
#include <sxcurl.hpp>
#include <sxdiv.hpp>
 
namespace HArDCore3D
{
 
  struct StokesNorms
  {
    StokesNorms(double norm_u, double norm_curl_u, double norm_p, double norm_grad_p):
      u(norm_u),
      curl_u(norm_curl_u),
      p(norm_p),
      grad_p(norm_grad_p),
      hcurl_u( std::sqrt( std::pow(norm_u, 2) + std::pow(norm_curl_u, 2) ) ),
      hgrad_p( std::sqrt( std::pow(norm_p, 2) + std::pow(norm_grad_p, 2) ) )
      {
        // do nothing
      };
 
    // Constructor without arguments
    StokesNorms() : StokesNorms(0., 0., 0., 0.) {};
    
    double u; 
    double curl_u; 
    double p; 
    double grad_p; 
    double hcurl_u; 
    double hgrad_p; 
  
  };
 
  struct NavierStokes
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> ForcingTermType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> VelocityType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> VorticityType;
    typedef std::function<double(const Eigen::Vector3d &)> PressureType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> PressureGradientType;
    typedef IntegralWeight ViscosityType;
 
    NavierStokes(
           const DDRCore & ddrcore,           
           const BoundaryConditions & BC,     
           bool use_threads,                  
           std::ostream & output = std::cout  
         );
 
    //--------------------------------------------------
    //    Methods to assemble and solve
    //--------------------------------------------------    
 
 
    Eigen::VectorXd vectorDirBC(
                                const VelocityType & u,    
                                const PressureType & p     
                               );
 
    void computeLinearComponents(
                                const ForcingTermType & f,      
                                const VelocityType & u,   
                                const VorticityType & omega,  
                                const double & reynolds     
                                );
 
    double assembleLinearSystem(
                              const Eigen::VectorXd & Xn,  
                              const double & reac         
                              );
 
    double normGlobalRHS(
                        const std::vector<Eigen::VectorXd> & locRHS  
                        );
 
    Eigen::VectorXd newtonIncrement( 
                                   LinearSolver<SystemMatrixType> & solver,           
                                   const double & navier_scaling,   
                                   const double & relax,            
                                   double & norm_dX                 
                                   );
 
    std::vector<StokesNorms> computeStokesNorms(
                const std::vector<Eigen::VectorXd> & list_dofs 
                ) const;
 
    std::pair<StokesNorms, StokesNorms> computeContinuousErrorsNorms(
                const Eigen::VectorXd & v,            
                const VelocityType & u,               
                const VorticityType & curl_u,         
                const PressureType & p,               
                const PressureGradientType & grad_p   
                ) const;
 
    double computeFlux(
                  const Eigen::VectorXd & u,            
                  const Hull & surface                  
                  ) const;
 
    //------------------------------------------------------------------
    //    Getters: dimension of space and numbers of unknowns of various types
    //------------------------------------------------------------------
    
    inline size_t nDOFs_up() const
    {
      return m_sxcurl.dimension() + m_sxgrad.dimension();
    }
 
    // ------ Dirichlet DOFs ------
    inline size_t nDOFs_dir_edge_u() const
    {
      return m_nDOFs_dir_edge_u;
    }
 
    inline size_t nDOFs_dir_face_u() const
    {
      return m_nDOFs_dir_face_u;
    }
 
    inline size_t nDOFs_dir_u() const
    {
      return m_nDOFs_dir_edge_u + m_nDOFs_dir_face_u;
    }
 
    inline size_t nDOFs_dir_vertex_p() const
    {
      return m_nDOFs_dir_vertex_p;
    }
 
    inline size_t nDOFs_dir_edge_p() const
    {
      return m_nDOFs_dir_edge_p;
    }
 
    inline size_t nDOFs_dir_face_p() const
    {
      return m_nDOFs_dir_face_p;
    }
 
    inline size_t nDOFs_dir_p() const
    {
      return m_nDOFs_dir_vertex_p + m_nDOFs_dir_edge_p + m_nDOFs_dir_face_p;
    }
 
    // ------ statically condensed and global UKN -----
    inline size_t nSCUKN_u() const
    {
      return m_nloc_sc_u.sum();
    }
 
    inline size_t nSCUKN_p() const
    {
      return m_nloc_sc_p.sum();
    }
 
    inline size_t nSCUKN() const
    {
      return nSCUKN_u() + nSCUKN_p();
    }
    
    inline size_t nUKN_u() const
    {
      return m_sxcurl.dimension() - nDOFs_dir_u();
    }
    
    inline size_t nUKN_p() const
    {
      return m_sxgrad.dimension() - nDOFs_dir_p();
    }
 
    inline size_t nUKN_lambda() const
    {
      return m_nUKN_lambda;
    }
 
    inline size_t nUKN_uplambda() const
    {
      return nUKN_u() + nUKN_p() + nUKN_lambda();
    }
 
    inline size_t nGLUKN_u() const
    {
      return nUKN_u() - nSCUKN_u();
    }
    
    inline size_t nGLUKN_p() const
    {
      return nUKN_p() - nSCUKN_p();
    }
 
    inline size_t nGLUKN() const
    {
      return nUKN_uplambda() - nSCUKN();
    }
 
    //--------------------------------------------------------
    //  Getters: underlying SDDR spaces
    //--------------------------------------------------------
    
    inline const SXGrad & sxGrad() const
    {
      return m_sxgrad;
    }
 
    inline const SXCurl & sxCurl() const
    {
      return m_sxcurl;
    }
 
    inline const SXDiv & sxDiv() const
    {
      return m_sxdiv;
    }
 
    //-------------------------------------------------
    //    Matrices and RHS
    //-------------------------------------------------
    
    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 & scMatrix() const {
      return m_sc_A;
    }
 
    inline Eigen::VectorXd & scVector() {
      return m_sc_b;
    }
 
    inline std::vector<Eigen::VectorXd> & rhsSourceNaturalBC() {
      return m_rhs_source_naturalBC;
    }
 
    //-------------------------------------------
    //  Getters: various parameters
    //-------------------------------------------
    
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    inline const bool & iterativeSolver() const {
      return m_itsolver;
    }
 
    inline bool & iterativeSolver() {
      return m_itsolver;
    }
 
    inline double & Reynolds() {
      return m_reynolds;
    }
 
  private:
    /* First matrix: (CT ., CT .)_div
       Second matrix: (., GT .)_curl
       RHS: of size velocity + pressure + nb lagrange multiplier (0 or 1), and includes natural BCs
    */
    std::tuple<Eigen::MatrixXd, Eigen::MatrixXd, Eigen::VectorXd, Eigen::MatrixXd>
    _compute_linear_contributions(
                                const size_t iT,                      
                                const Eigen::VectorXd & interp_f,     
                                const ForcingTermType & f             
                                );
 
    Eigen::MatrixXd _compute_Stokes_matrix(
                                        const size_t iT           
                                        );
 
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_jacobian_rhs(
                          const size_t iT,   
                          const Eigen::VectorXd & Xn,     
                          const double & reac       
                          );
 
    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                                     
                                      );
    
    // Constructor parameter
    const DDRCore & m_ddrcore;
    bool m_use_threads;
    std::ostream & m_output;
    SerendipityProblem m_ser_pro;
    const SXGrad m_sxgrad;
    const SXCurl m_sxcurl;
    const SXDiv m_sxdiv;
    const Eigen::VectorXi m_nloc_sc_u; // Nb of statically condensed DOFs for velocity in each cell (cell unknowns)
    const Eigen::VectorXi m_nloc_sc_p; // Nb of statically condensed DOFs for pressure in each cell (cell unknowns)
 
    // Boundary conditions
    BoundaryConditions m_BC;
    size_t m_nDOFs_dir_edge_u;
    size_t m_nDOFs_dir_face_u;
    size_t m_nDOFs_dir_vertex_p;
    size_t m_nDOFs_dir_edge_p;
    size_t m_nDOFs_dir_face_p;
    size_t m_nUKN_lambda;
    std::vector<std::pair<size_t,size_t>> m_locSCUKN;    // Location, among the DOFs, of the statically condensed unknowns (removing Dirichlet and non-SC unknowns); the unknowns are between m_locSCUCK[i].first and m_locSCUCK[i].first+m_locSCUKN[i].second for all i
    std::vector<std::pair<size_t,size_t>> m_locGLUKN;    // Location, among the DOFs, of the globally coupled unknowns (removing Dirichlet and non-global unknowns); the unknowns are between m_locGLUCK[i].first and m_locGLUCK[i].first+m_locGLUKN[i].second for all i
    Eigen::VectorXi m_DOFtoSCUKN;     // Maps a dof i to its SC unknown in the recovery operator, or to -1 if the DOF is not a SCUKN
    Eigen::VectorXi m_DOFtoGLUKN;     // Maps a dof i to its GL unknown in the global system, or to -1 if the DOF is a not a GLUKN
 
    // Scheme and solver parameters    
    double m_stab_par;
    bool m_itsolver;    // true if we use the iterative solver, false uses a direct solver
 
    double m_reynolds = 1.;   // Initialised at 1 by default, has to be adjusted after loading the value
 
    // Local matrices and RHS (including BCs) for the linear part of the model. The boolean to flag that these quantities have been computed
    std::vector<Eigen::MatrixXd> m_CTuCTu;    // curl - curl matrix
    std::vector<Eigen::MatrixXd> m_PTuGTp;    // velocity - grad pressure matrix
    std::vector<Eigen::VectorXd> m_rhs_source_naturalBC;   // RHS from the source and natural BCs
    std::vector<Eigen::MatrixXd> m_massT;     // mass-matrix for pseudo-temporal iterations
    bool m_linear_computed;
    // Global matrix and RHS of statically condensed system, at each Newton iteration
    SystemMatrixType m_A;   
    Eigen::VectorXd m_b;
    // Global static condensation operator and RHS (to recover statically condensed DOFs), at each Newton iteration
    SystemMatrixType m_sc_A; 
    Eigen::VectorXd m_sc_b;
 
  };
 
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
  using std::pow;
 
  // To scale the pressure, the viscosity and the nonlinear term (curl u) x u
  double pressure_scaling = 1.;
  double reynolds = 1.;
  double navier_scaling = 1.;
  
  //------------------------------------------------------------------------------
  // Trigonometric solution
  //------------------------------------------------------------------------------
 
  static NavierStokes::VelocityType
  trigonometric_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d(
                                             0.5 * sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                             0.5 * cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                             -cos(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2))
                                             );
                    };
 
  static NavierStokes::VorticityType
  trigonometric_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return 3. * PI * Eigen::Vector3d(
                                                      cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      -sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      0.
                                                      );
                         };
 
  static NavierStokes::PressureType
  trigonometric_p = [](const Eigen::Vector3d & x) -> double {
                      return pressure_scaling * sin(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2));
                    };
 
  static NavierStokes::PressureGradientType
  trigonometric_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return pressure_scaling * 2. * PI * Eigen::Vector3d(
                                                      cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      sin(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2))
                                                      );
                         };
 
  static NavierStokes::ForcingTermType
  trigonometric_curl_u_cross_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d (3.0*PI*sin(2*PI*x.x())*pow(sin(2*PI*x.z()), 2)*cos(2*PI*x.x())*pow(cos(2*PI*x.y()), 2),
                                              3.0*PI*sin(2*PI*x.y())*pow(sin(2*PI*x.z()), 2)*pow(cos(2*PI*x.x()), 2)*cos(2*PI*x.y()),
                                              1.5*PI*pow(sin(2*PI*x.x()), 2)*sin(2*PI*x.z())*pow(cos(2*PI*x.y()), 2)*cos(2*PI*x.z()) + 1.5*PI*pow(sin(2*PI*x.y()), 2)*sin(2*PI*x.z())*pow(cos(2*PI*x.x()), 2)*cos(2*PI*x.z())
                                              );
                    };
 
  static NavierStokes::ForcingTermType
  trigonometric_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return (1./reynolds) * 6. * std::pow(PI, 2) * Eigen::Vector3d(
                                                                    sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                                                    cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                                                    -2. * cos(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2))
                                                                    )
                        + navier_scaling * trigonometric_curl_u_cross_u(x)
                        + trigonometric_grad_p(x);
                    };
 
  //------------------------------------------------------------------------------
  // Constant solution
  //------------------------------------------------------------------------------
 
  static NavierStokes::VelocityType
  constant_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return Eigen::Vector3d(1., 1., 1.);
             };
 
  static NavierStokes::VorticityType
  constant_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static NavierStokes::ForcingTermType
  constant_curl_u_cross_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static NavierStokes::PressureType
  constant_p = [](const Eigen::Vector3d & x) -> double {
               return 0.;
             };
 
  static NavierStokes::PressureGradientType
  constant_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static NavierStokes::ForcingTermType
  constant_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return Eigen::Vector3d::Zero();
             };
 
  //------------------------------------------------------------------------------
  // Linear solution
  //------------------------------------------------------------------------------
 
  static NavierStokes::VelocityType
  linear_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return 1e3*Eigen::Vector3d(x(0), -x(1), 0.);
             };
 
  static NavierStokes::VorticityType
  linear_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static NavierStokes::ForcingTermType
  linear_curl_u_cross_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static NavierStokes::PressureType
  linear_p = [](const Eigen::Vector3d & x) -> double {
               return pressure_scaling * (x(0) - x(1));
             };
 
  static NavierStokes::PressureGradientType
  linear_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return pressure_scaling * Eigen::Vector3d(1., -1., 0.);
                  };
 
  static NavierStokes::ForcingTermType
  linear_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return linear_grad_p(x);
             };
 
  //------------------------------------------------------------------------------
  // Vertical flow
  //------------------------------------------------------------------------------
 
  constexpr double scal_u = 1.;
  static NavierStokes::VelocityType
  vertical_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return scal_u * Eigen::Vector3d( cos(PI*x.y()), cos(PI*x.x()), 1.);
                    };
 
  static NavierStokes::VorticityType
  vertical_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return scal_u * Eigen::Vector3d(0, 0, -PI*sin(PI*x.x()) + PI*sin(PI*x.y()) );
                         };
 
  static NavierStokes::PressureType
  vertical_p = trigonometric_p;
 
  static NavierStokes::PressureGradientType
  vertical_grad_p = trigonometric_grad_p;
 
  static NavierStokes::ForcingTermType
  vertical_curl_u_cross_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return std::pow(scal_u, 2) 
                           * Eigen::Vector3d (
                                             -(-PI*sin(PI*x.x()) + PI*sin(PI*x.y()))*cos(PI*x.x()),
                                             (-PI*sin(PI*x.x()) + PI*sin(PI*x.y()))*cos(PI*x.y()), 
                                             0.
                                             );
                    };
 
  static NavierStokes::ForcingTermType
  vertical_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return (1./reynolds) * scal_u * Eigen::Vector3d(pow(PI, 2)*cos(PI*x.y()), pow(PI, 2)*cos(PI*x.x()), 0.)
                        + navier_scaling * vertical_curl_u_cross_u(x)
                        + vertical_grad_p(x);
//                        return Eigen::Vector3d::Zero();
                    };
 
 
  //------------------------------------------------------------------------------
  // Impose BC: p=1 at entrance x=0, u.n=1 at exit x=1 on lower bottom corner
  //------------------------------------------------------------------------------
  static NavierStokes::VelocityType
  pressflux_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return ( lower_bottom_corner_x_one.is_in(x) ? VectorRd(1., 0., 0.) : VectorRd::Zero());
                    };
 
  static NavierStokes::VorticityType
  pressflux_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                          return Eigen::Vector3d::Zero();
                        };
 
  static NavierStokes::PressureType
  pressflux_p = [](const Eigen::Vector3d & x) -> double {
                          // return (x.x()<1e-8 ? 1. : 0.);
                          return  pressure_scaling * (-x.z());
                        };
 
  static NavierStokes::PressureGradientType
  pressflux_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                         return Eigen::Vector3d::Zero();
                      };
 
  static NavierStokes::ForcingTermType
  pressflux_curl_u_cross_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d::Zero();
                    };
 
  static NavierStokes::ForcingTermType
  pressflux_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d::Zero();
                    };
 
 
} // end of namespace HArDCore3D
 
#endif