navier_stokes.h

#ifndef __NAVIER_STOKES__
#define __NAVIER_STOKES__

#include "euler.h"
#include "parameters/parameters_navier_stokes.h"

namespace PHiLiP {
namespace Physics {

template <int dim, int nstate, typename real>
class NavierStokes : public Euler <dim, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using PhysicsBase<dim,nstate,real>::dissipative_flux;
    using PhysicsBase<dim,nstate,real>::source_term;
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    NavierStokes( 
        const Parameters::AllParameters *const                    parameters_input,
        const double                                              ref_length,
        const double                                              gamma_gas,
        const double                                              mach_inf,
        const double                                              angle_of_attack,
        const double                                              side_slip_angle,
        const double                                              prandtl_number,
        const double                                              reynolds_number_inf,
        const bool                                                use_constant_viscosity,
        const double                                              constant_viscosity,
        const double                                              temperature_inf = 273.15,
        const double                                              isothermal_wall_temperature = 1.0,
        const thermal_boundary_condition_enum                     thermal_boundary_condition_type = thermal_boundary_condition_enum::adiabatic,
        std::shared_ptr< ManufacturedSolutionFunction<dim,real> > manufactured_solution_function = nullptr,
        const two_point_num_flux_enum                             two_point_num_flux_type = two_point_num_flux_enum::KG);

    const double viscosity_coefficient_inf;
    const bool use_constant_viscosity;
    const double constant_viscosity;
    const double prandtl_number;
    const double reynolds_number_inf;
    const double isothermal_wall_temperature;
    const thermal_boundary_condition_enum thermal_boundary_condition_type;

protected:    

    const double sutherlands_temperature; 
    const double freestream_temperature; 
    const double temperature_ratio; 


public:

    template<typename real2>
    std::array<dealii::Tensor<1,dim,real2>,nstate> 
    convert_conservative_gradient_to_primitive_gradient (
        const std::array<real2,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &conservative_soln_gradient) const;

    template<typename real2>
    dealii::Tensor<1,dim,real2> compute_temperature_gradient (
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient) const;

    template<typename real2>
    real2 compute_viscosity_coefficient (const std::array<real2,nstate> &primitive_soln) const;

    template<typename real2>
    real2 compute_viscosity_coefficient_sutherlands_law (const std::array<real2,nstate> &primitive_soln) const;

    template<typename real2>
    real2 scale_viscosity_coefficient (const real2 viscosity_coefficient) const;

    template<typename real2>
    real2 compute_scaled_viscosity_coefficient (const std::array<real2,nstate> &primitive_soln) const;

    template<typename real2>
    real2 compute_scaled_heat_conductivity_given_scaled_viscosity_coefficient_and_prandtl_number (
        const real2 scaled_viscosity_coefficient, 
        const double prandtl_number_input) const;

    template<typename real2>
    real2 compute_scaled_heat_conductivity (const std::array<real2,nstate> &primitive_soln) const;

    template<typename real2>
    dealii::Tensor<1,dim,real2> compute_heat_flux (
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient) const;

    template<typename real2>
    dealii::Tensor<1,dim,real2> compute_heat_flux_given_scaled_heat_conductivity_and_temperature_gradient (
        const real2 scaled_heat_conductivity,
        const dealii::Tensor<1,dim,real2> &temperature_gradient) const;

    template<typename real2>
    dealii::Tensor<1,3,real2> compute_vorticity (
        const std::array<real2,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &conservative_soln_gradient) const;

    real compute_vorticity_magnitude_sqr (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_vorticity_magnitude (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_enstrophy (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_vorticity_based_dissipation_rate_from_integrated_enstrophy (
        const real integrated_enstrophy) const;

    real compute_pressure_dilatation (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    dealii::Tensor<2,dim,real> compute_deviatoric_strain_rate_tensor (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_deviatoric_strain_rate_tensor_magnitude_sqr (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_deviatoric_strain_rate_tensor_based_dissipation_rate_from_integrated_deviatoric_strain_rate_tensor_magnitude_sqr (
        const real integrated_deviatoric_strain_rate_tensor_magnitude_sqr) const;

    template<typename real2>
    dealii::Tensor<2,dim,real2> 
    extract_velocities_gradient_from_primitive_solution_gradient (
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient) const;

    template<typename real2>
    dealii::Tensor<2,dim,real2> 
    compute_strain_rate_tensor (
        const dealii::Tensor<2,dim,real2> &vel_gradient) const;

    real compute_strain_rate_tensor_magnitude_sqr (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_strain_rate_tensor_based_dissipation_rate_from_integrated_strain_rate_tensor_magnitude_sqr (
        const real integrated_strain_rate_tensor_magnitude_sqr) const;


    template<typename real2>
    dealii::Tensor<2,dim,real2> 
    compute_viscous_stress_tensor_via_scaled_viscosity_and_strain_rate_tensor (
        const real2 scaled_viscosity_coefficient,
        const dealii::Tensor<2,dim,real2> &strain_rate_tensor) const;

    template<typename real2>
    dealii::Tensor<2,dim,real2>
    compute_viscous_stress_tensor (
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient) const;

    std::array<dealii::Tensor<1,dim,real>,nstate> 
    dissipative_flux (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const override;

    dealii::Tensor<1,dim,real> compute_scaled_viscosity_gradient (
        const std::array<real,nstate> &primitive_soln,
        const dealii::Tensor<1,dim,real> &temperature_gradient) const;

    dealii::Tensor<2,nstate,real> dissipative_flux_directional_jacobian (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const dealii::Tensor<1,dim,real> &normal) const;

    dealii::Tensor<2,nstate,real> dissipative_flux_directional_jacobian_wrt_gradient_component (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const dealii::Tensor<1,dim,real> &normal,
        const int d_gradient) const;

    std::array<real,nstate> dissipative_source_term (
        const dealii::Point<dim,real> &pos) const;

    std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_soln,
        const real current_time) const override;

    dealii::Tensor<2,nstate,real> convective_flux_directional_jacobian_via_dfad (
        std::array<real,nstate> &conservative_soln,
        const dealii::Tensor<1,dim,real> &normal) const;

    real compute_scaled_viscosity_coefficient_derivative_wrt_temperature_via_dfad (
        std::array<real,nstate> &conservative_soln) const;

    template<typename real2>
    std::array<dealii::Tensor<1,dim,real2>,nstate> 
    dissipative_flux_given_velocities_viscous_stress_tensor_and_heat_flux (
        const dealii::Tensor<1,dim,real2> &vel,
        const dealii::Tensor<2,dim,real2> &viscous_stress_tensor,
        const dealii::Tensor<1,dim,real2> &heat_flux) const;

protected:

    template <typename real2>
    std::array<dealii::Tensor<1,dim,real2>,nstate> 
    dissipative_flux_templated (
        const std::array<real2,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &solution_gradient) const;

    void boundary_wall (
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;

    void boundary_manufactured_solution (
        const dealii::Point<dim, real> &pos,
        const dealii::Tensor<1,dim,real> &normal_int,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;

private:
    real get_tensor_magnitude_sqr (const dealii::Tensor<2,dim,real> &tensor) const;

};

} // Physics namespace
} // PHiLiP namespace

#endif