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 nspecies, int nstate, typename real>
class NavierStokes : public Euler <dim, nspecies, nstate, real>
{
protected:
    // For overloading the virtual functions defined in PhysicsBase
    using PhysicsBase<dim,nspecies,nstate,real>::dissipative_flux;
    using PhysicsBase<dim,nspecies,nstate,real>::source_term;
    using PhysicsBase<dim,nspecies,nstate,real>::boundary_face_values;
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,nspecies,real> > manufactured_solution_function = nullptr,
        const two_point_num_flux_enum                             two_point_num_flux_type = two_point_num_flux_enum::KG,
        const bool                                                has_nonzero_physical_source = false);

    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>
    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>
    dealii::Tensor<1,dim,real2> compute_velocities_parallel_to_wall(
        const std::array<real2,nstate> &conservative_soln,
        const dealii::Tensor<1,dim,real2> &normal_vector) const;

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

    template<typename real2>
    dealii::Tensor<1,dim,real2> compute_wall_tangent_vector_from_velocities_parallel_to_wall(
        const dealii::Tensor<1,dim,real2> &velocities_parallel_to_wall) const;

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

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

    template<typename real2>
    real2 compute_viscosity_coefficient_from_temperature (const real2 temperature) const;

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

    template<typename real2>
    real2 compute_viscosity_coefficient_sutherlands_law_from_temperature (const real2 temperature) 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_second_invariant (
        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_incompressible_enstrophy (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &conservative_soln_gradient) const;

    real compute_incompressible_palinstrophy (
        const std::array<real,nstate> &conservative_soln,
        const std::array<dealii::Tensor<1,dim,real>,3> &vorticity_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;

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

    real compute_density_gradient_magnitude (
        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;

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

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

    real compute_viscosity_times_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_viscosity_times_deviatoric_strain_rate_tensor_magnitude_sqr (
        const real integrated_viscosity_times_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_viscosity_times_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_viscosity_times_strain_rate_tensor_magnitude_sqr (
        const real integrated_viscosity_times_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;

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

    template<typename real2>
    dealii::Tensor<2,dim,real2>
    compute_viscous_stress_tensor_from_conservative_templated (
        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<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;

    real get_tensor_product_magnitude_sqr (
        const dealii::Tensor<2,dim,real> &tensor1,
        const dealii::Tensor<2,dim,real> &tensor2) const;

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

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

    std::array<real,nstate> dissipative_flux_dot_normal (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const bool on_boundary,
        const dealii::types::global_dof_index cell_index,
        const dealii::Tensor<1,dim,real> &normal,
        const int boundary_type) override;

    virtual std::array<real,nstate> dissipative_flux_dot_normal_on_adiabatic_boundary (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const dealii::types::global_dof_index cell_index,
        const dealii::Tensor<1,dim,real> &normal);

    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;

    virtual void boundary_face_values_viscous_flux (
        const int boundary_type,
        const dealii::Point<dim, real> &pos,
        const dealii::Tensor<1,dim,real> &normal,
        const std::array<real,nstate> &soln_int,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
        const std::array<real,nstate> &/*filtered_soln_int*/,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*filtered_soln_grad_int*/,
        std::array<real,nstate> &soln_bc,
        std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const override;

    void boundary_wall_viscous_flux (
        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;

    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;

public:
    dealii::Vector<double> post_compute_derived_quantities_vector (
        const dealii::Vector<double>              &uh,
        const std::vector<dealii::Tensor<1,dim> > &duh,
        const std::vector<dealii::Tensor<2,dim> > &dduh,
        const dealii::Tensor<1,dim>               &normals,
        const dealii::Point<dim>                  &evaluation_points) const override;
    
    std::vector<std::string> post_get_names () const override;
    
    std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const override;
    
    dealii::UpdateFlags post_get_needed_update_flags () const override;

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

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

};

template <int dim, int nspecies, int nstate, typename real>
class NavierStokes_ChannelFlowConstantSourceTerm : public NavierStokes <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    NavierStokes_ChannelFlowConstantSourceTerm( 
        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                                              reynolds_number_based_on_friction_velocity,
        const double                                              half_channel_height,
        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,nspecies,real> > manufactured_solution_function = nullptr,
        const two_point_num_flux_enum                             two_point_num_flux_type = two_point_num_flux_enum::KG);

    const double x_momentum_constant_source_term;

    ~NavierStokes_ChannelFlowConstantSourceTerm() {};

    std::array<real,nstate> physical_source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const dealii::types::global_dof_index cell_index) const override;
};

template <typename real>
class WallModelLookUpTable
{
    static const int NUMBER_OF_SAMPLE_POINTS = 38;
    static constexpr std::array<double,NUMBER_OF_SAMPLE_POINTS> yData = 
            {{0.0, 3.0, 5.0, 8.0, 10.0, 20.0, 35.0, 50.0, 75.0, 100.0, 125.0, 150.0,
              200.0, 250.0, 300.0, 350.0, 400.0, 500.0, 575.0, 650.0, 725.0, 800.0, 
              900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0, 1800.0, 
              2000.0, 2500.0, 3000.0, 3500.0, 4000.0, 4500.0, 5000.0}};
    static constexpr std::array<double,NUMBER_OF_SAMPLE_POINTS> xData = 
            {{0.0000000000000e+00, 8.7652737082527e+00, 2.3294562921324e+01, 5.4184774233564e+01,
              7.8809463502726e+01, 2.2429199943857e+02, 4.6798107411288e+02, 7.2296127461536e+02, 
              1.1658971129082e+03, 1.6286476514445e+03, 2.1076047703808e+03, 2.5997463219921e+03, 
              3.6154681137081e+03, 4.6644290903876e+03, 5.7398811747011e+03, 6.8373640504590e+03, 
              7.9537073133201e+03, 1.0234022966195e+04, 1.1979574602097e+04, 1.3750942900204e+04, 
              1.5545134905888e+04, 1.7359780637619e+04, 1.9807850068644e+04, 2.2285219477338e+04, 
              2.4788948528427e+04, 2.7316633948913e+04, 2.9866274540736e+04, 3.2436177975259e+04, 
              3.5024894352307e+04, 3.7631167553565e+04, 4.2892118831955e+04, 4.8211670442755e+04, 
              6.1730908478994e+04, 7.5515087179925e+04, 8.9519590674471e+04, 1.0371270861601e+05, 
              1.1807073118603e+05, 1.3257525485198e+05}};
public:
    WallModelLookUpTable(); 

    ~WallModelLookUpTable(){}; 

private:
    real interpolate(const real x, const bool extrapolate ) const; 

public:
    real get_wall_shear_stress_magnitude(
        const real wall_parallel_velocity, 
        const real distance, 
        const real viscosity_coefficient,
        const real density,
        const double reynolds_number_inf) const;
};

template <int dim, int nspecies, int nstate, typename real>
class NavierStokes_ChannelFlowConstantSourceTerm_WallModel : public NavierStokes_ChannelFlowConstantSourceTerm <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    NavierStokes_ChannelFlowConstantSourceTerm_WallModel( 
        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                                              reynolds_number_based_on_friction_velocity,
        const double                                              half_channel_height,
        const double                                              distance_from_wall_for_wall_model_input_velocity,
        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,nspecies,real> > manufactured_solution_function = nullptr,
        const two_point_num_flux_enum                             two_point_num_flux_type = two_point_num_flux_enum::KG);

    ~NavierStokes_ChannelFlowConstantSourceTerm_WallModel() {};

    const double distance_from_wall_for_wall_model_input_velocity;

    std::unique_ptr < WallModelLookUpTable<real> > wall_model_look_up_table;

    std::array<real,nstate> dissipative_flux_dot_normal_on_adiabatic_boundary (
        const std::array<real,nstate> &solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient,
        const std::array<real,nstate> &filtered_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &filtered_solution_gradient,
        const dealii::types::global_dof_index cell_index,
        const dealii::Tensor<1,dim,real> &normal) override;

    real get_velocity_component_parallel_to_wall_from_solution_and_normal_vector (
        const std::array<real,nstate> &conservative_soln,
        const dealii::Tensor<1,dim,real> &normal_vector) const;
};

} // Physics namespace
} // PHiLiP namespace

#endif