reynolds_averaged_navier_stokes.h

#ifndef __REYNOLDS_AVERAGED_NAVIER_STOKES__
#define __REYNOLDS_AVERAGED_NAVIER_STOKES__

#include "model.h"
#include "navier_stokes.h"
#include "euler.h"

namespace PHiLiP {
namespace Physics {

template <int dim, int nspecies, int nstate, typename real>
class ReynoldsAveragedNavierStokesBase : public ModelBase <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    ReynoldsAveragedNavierStokesBase(
        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                                              turbulent_prandtl_number,
        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);

    static const int nstate_navier_stokes = dim+2;

    static const int nstate_turbulence_model = nstate-(dim+2);

    const double turbulent_prandtl_number;

    std::unique_ptr< NavierStokes<dim,nspecies,nstate_navier_stokes,real> > navier_stokes_physics;

    std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
        const std::array<real,nstate> &conservative_soln) 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 dealii::types::global_dof_index cell_index) 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) const override;


    std::array<real,nstate> convective_eigenvalues (
        const std::array<real,nstate> &/*conservative_soln*/,
        const dealii::Tensor<1,dim,real> &/*normal*/) const;


    real max_convective_eigenvalue (const std::array<real,nstate> &soln) const;


    real max_convective_normal_eigenvalue (
        const std::array<real,nstate> &soln,
        const dealii::Tensor<1,dim,real> &normal) const;

    std::array<real,nstate> source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_solution,
        const real current_time,
        const dealii::types::global_dof_index cell_index) const;

    std::array<real,nstate> convective_dissipative_source_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_solution,
        const dealii::types::global_dof_index cell_index) const;

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

    virtual dealii::Tensor<2,dim,real> compute_Reynolds_stress_tensor (
        const std::array<real,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,real>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<real,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual dealii::Tensor<1,dim,real> compute_Reynolds_heat_flux (
        const std::array<real,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,real>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<real,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual dealii::Tensor<2,dim,FadType> compute_Reynolds_stress_tensor_fad (
        const std::array<FadType,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<FadType,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual dealii::Tensor<1,dim,FadType> compute_Reynolds_heat_flux_fad (
        const std::array<FadType,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<FadType,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual std::array<real,nstate_turbulence_model> compute_effective_viscosity_turbulence_model (
        const std::array<real,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<real,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual std::array<FadType,nstate_turbulence_model> compute_effective_viscosity_turbulence_model_fad (
        const std::array<FadType,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<FadType,nstate_turbulence_model> &primitive_soln_turbulence_model) const = 0;

    virtual std::array<real,nstate> compute_production_dissipation_cross_term (
        const dealii::Point<dim,real> &pos,
        const std::array<real,nstate> &conservative_solution,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const = 0;

protected:
    template<typename real2> 
    real2 get_vector_magnitude_sqr (const dealii::Tensor<1,3,real2> &vector) const;

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

    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 dealii::types::global_dof_index cell_index) const;

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

    template <typename real2>
    std::array<real2,dim+2> extract_rans_conservative_solution (
        const std::array<real2,nstate> &conservative_soln) const;

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

    template <typename real2>
    std::array<dealii::Tensor<1,dim,real2>,nstate-(dim+2)> dissipative_flux_turbulence_model (
        const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model,
        const std::array<dealii::Tensor<1,dim,real2>,nstate_turbulence_model> &primitive_solution_gradient_turbulence_model) const;

    template <typename real2>
    std::array<real2,nstate-(dim+2)> convert_conservative_to_primitive_turbulence_model (
        const std::array<real2,nstate> &conservative_soln) const;

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


    std::array<real,nstate-(dim+2)> compute_mean_turbulence_property (
        const std::array<real,nstate> &conservative_soln1,
        const std::array<real,nstate> &conservative_soln2) const;

    dealii::Tensor<2,nstate,real> convective_flux_directional_jacobian (
        const std::array<real,nstate> &conservative_soln,
        const dealii::Tensor<1,dim,real> &normal) 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::types::global_dof_index cell_index) 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 dealii::types::global_dof_index cell_index) const;

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

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

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

    std::array<real,nstate> dissipative_source_term_computed_from_manufactured_solution (
        const dealii::Point<dim,real> &pos,
        const dealii::types::global_dof_index cell_index) const;

    std::array<real,nstate> physical_source_term_computed_from_manufactured_solution (
        const dealii::Point<dim,real> &pos,
        const dealii::types::global_dof_index cell_index) 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;
};

} // Physics namespace
} // PHiLiP namespace

#endif