negative_spalart_allmaras_rans_model.h

#ifndef __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__
#define __NEGATIVE_SPALART_ALLMARAS_RANS_MODEL__

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

namespace PHiLiP {
namespace Physics {

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

    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;

    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;

    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;

    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;

    real compute_eddy_viscosity(
        const std::array<real,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<real,nstate_turbulence_model> &primitive_soln_turbulence_model) const;

    FadType compute_eddy_viscosity_fad(
        const std::array<FadType,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<FadType,nstate_turbulence_model> &primitive_soln_turbulence_model) const;

    std::array<real,nstate-(dim+2)> 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;

    std::array<FadType,nstate-(dim+2)> 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;

    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;

    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<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const override;

    std::vector<std::string> post_get_names () const override;

protected:
    template<typename real2> dealii::Tensor<2,dim,real2> compute_Reynolds_stress_tensor_templated (
        const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model) const;

    template<typename real2> dealii::Tensor<1,dim,real2> compute_Reynolds_heat_flux_templated (
        const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<dealii::Tensor<1,dim,real2>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model) const;

    template<typename real2> real2 scale_coefficient(
        const real2 coefficient) const;


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

    void boundary_outflow (
        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_inflow (
        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_farfield (
        std::array<real,nstate> &soln_bc) const override;

private:

    template<typename real2> real2 compute_eddy_viscosity_templated(
        const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model) const;


    template<typename real2> std::array<real2,nstate-(dim+2)> compute_effective_viscosity_turbulence_model_templated (
        const std::array<real2,nstate_navier_stokes> &primitive_soln_rans,
        const std::array<real2,nstate_turbulence_model> &primitive_soln_turbulence_model) const;


    template<typename real2> real2 compute_coefficient_Chi (
        const real2 nu_tilde,
        const real2 laminar_kinematic_viscosity) const;


    template<typename real2> real2 compute_coefficient_f_v1 (
        const real2 coefficient_Chi) const;


    real compute_coefficient_f_v2 (
        const real coefficient_Chi) const;


    template<typename real2> real2 compute_coefficient_f_n (
        const real2 nu_tilde,
        const real2 laminar_kinematic_viscosity) const;


    real compute_coefficient_f_t2 (
        const real coefficient_Chi) const;


    real compute_coefficient_f_w (
        const real coefficient_g) const;


    real compute_coefficient_r (
        const real nu_tilde,
        const real d_wall,
        const real s_tilde) const;


    real compute_coefficient_g (
        const real coefficient_r) const;

    real compute_s (
        const std::array<real,nstate_navier_stokes> &conservative_soln_rans,
        const std::array<dealii::Tensor<1,dim,real>,nstate_navier_stokes> &conservative_soln_gradient_rans) const;


    real compute_s_bar (
        const real coefficient_Chi,
        const real nu_tilde,
        const real d_wall) const;


    real compute_s_tilde (
        const real coefficient_Chi,
        const real nu_tilde,
        const real d_wall,
        const real s) const;


    std::array<real,nstate> compute_production_source (
        const real coefficient_f_t2,
        const real density,
        const real nu_tilde,
        const real s,
        const real s_tilde) const;


    std::array<real,nstate> compute_dissipation_source (
        const real coefficient_f_t2,
        const real density,
        const real nu_tilde,
        const real d_wall,
        const real s_tilde) const;


    std::array<real,nstate> compute_cross_source (
        const real density,
        const real nu_tilde,
        const real laminar_kinematic_viscosity,
        const std::array<dealii::Tensor<1,dim,real>,nstate_navier_stokes> &primitive_soln_gradient_rans,
        const std::array<dealii::Tensor<1,dim,real>,nstate_turbulence_model> &primitive_solution_gradient_turbulence_model) const;


    const real c_b1  = 0.1355;
    const real sigma = 2.0/3.0;
    const real c_b2  = 0.622;
    const real kappa = 0.41;
    const real c_w1  = c_b1/(kappa*kappa)+(1+c_b2)/sigma;
    const real c_w2  = 0.3;
    const real c_w3  = 2.0;
    const real c_v1  = 7.1;
    const real c_v2  = 0.7;
    const real c_v3  = 0.9;
    const real c_t3  = 1.2;
    const real c_t4  = 0.5;
    const real c_n1  = 16.0;
    const real r_lim = 10.0;

    const FadType sigma_fad = 2.0/3.0;
    const FadType c_v1_fad  = 7.1;
    const FadType c_n1_fad  = 16.0;
};


} // Physics namespace
} // PHiLiP namespace

#endif