ns-solutions.hpp

Go to the documentation of this file.

#ifndef NS_SOLUTIONS_HPP
#define NS_SOLUTIONS_HPP
 
#include <Eigen/Dense>
 
#include "MeshBuilder2D.hpp"
 
namespace HArDCore2D
{
  namespace NavierStokesSolutions
  {   
    //------------------------------------------------------------------------------
    // Interface for exact solutions
 
    struct IExactSolution {
      IExactSolution(std::string name, bool is_stokes, double viscosity);
      virtual ~IExactSolution();
      
      inline std::string name() const {
        return m_name;
      }
 
      inline bool isStokes() const {
        return m_is_stokes;
      }
 
      inline double viscosity() const {
        return m_viscosity;
      }
    
      virtual Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const = 0;
      virtual double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const = 0;
      virtual Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const = 0;
      virtual double pressure(const Eigen::Vector2d & x, const double & t = 0) const = 0;
      virtual Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const = 0;
      virtual Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const = 0;
 
    protected:
      std::string m_name;
      bool m_is_stokes;
      double m_viscosity;
    };
  
    //------------------------------------------------------------------------------
    // Linear velocity
 
    struct LinearVelocity : public IExactSolution {
      LinearVelocity();
      
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
    
    //------------------------------------------------------------------------------
    // Linear pressure
 
    struct LinearPressure : public IExactSolution {
      LinearPressure();
      
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // Quadratic velocity
 
    struct QuadraticVelocity : public IExactSolution {
      QuadraticVelocity();
      
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // Trigonometric
 
    struct Trigonometric : public IExactSolution {
      Trigonometric();
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // A viscosity-dependent solution of the Stokes problem to test
    // pressure-robustness
 
    struct PressureRobustStokes : public IExactSolution {
      PressureRobustStokes(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // A viscosity-dependent solution of the Navier-Stokes problem to test
    // pressure-robustness
 
    struct PressureRobustNavierStokes : public IExactSolution {
      PressureRobustNavierStokes(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
    
    //------------------------------------------------------------------------------
    // A solution with zero velocity and non-zero pressure to test
    // pressure-robustness
 
    struct OnlyPressure : public IExactSolution {
      OnlyPressure(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // Kovasznay
 
    struct Kovasznay : public IExactSolution {
      Kovasznay(const double & viscosity, std::ostream & output = std::cout);
 
      static inline Mesh2D::MeshBuilder::TransformationType domain_transformation = [](const std::array<double, 2> & x) {
        return std::array<double, 2>{2. * x[0] - 0.5, 2. * x[1]};
      };
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
 
    private:
      double m_Re;
      double m_lambda;
    };
 
    //------------------------------------------------------------------------------
    // A time-dependent solution of the Navier-Stokes problem to test
    // time-stepping routines
 
    struct TimeDependentNavierStokes : public IExactSolution {
      TimeDependentNavierStokes(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // A time-dependent solution of the Navier-Stokes problem to test
    // time-stepping routines
 
    struct TimeDependentStokes : public IExactSolution {
      TimeDependentStokes(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
    //------------------------------------------------------------------------------
    // A divergence-free sinusoidal solution with time dependence to test timestepping
    // implementations
 
    struct TrigTimeDep : public IExactSolution {
      TrigTimeDep(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
 
    private:
      double m_param;
    };
 
    //------------------------------------------------------------------------------
    // A solution with only time-dependence to check time-stepping routine
 
    struct OnlyTimeDep : public IExactSolution {
      OnlyTimeDep(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
 
    private:
      double m_param;
    };
 
    //------------------------------------------------------------------------------
    // A solution with linear dependence in space and time
 
    struct LinearTemporal : public IExactSolution {
      LinearTemporal(const double & viscosity);
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
 
    private:
      double m_param;
    };
 
    //------------------------------------------------------------------------------
    // A lid-driven cavity test problem
    // Note: an exact solution does not exist - this test exists only to set initial
    // and boundary conditions
 
    struct LidDrivenCavity : public IExactSolution {
      LidDrivenCavity();
    
      Eigen::Vector2d momentumForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;
      double compressibilityForcingTerm(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d velocity(const Eigen::Vector2d & x, const double & t = 0) const;   
      double pressure(const Eigen::Vector2d & x, const double & t = 0) const;
      Eigen::Matrix2d velocityGradient(const Eigen::Vector2d & x, const double & t = 0) const;    
      Eigen::Vector2d pressureGradient(const Eigen::Vector2d & x, const double & t = 0) const;
    };
 
  } // namespace NavierStokesSolutions
 
} // namespace HArDCore2D
#endif