large_eddy_simulation.h

#ifndef __LARGE_EDDY_SIMULATION__
#define __LARGE_EDDY_SIMULATION__

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

namespace PHiLiP {
namespace Physics {

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

    const double ratio_of_filter_width_to_cell_size;

    std::unique_ptr< NavierStokes<dim,nspecies,nstate,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;


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


    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> &solution,
        const real current_time,
        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;

    double get_filter_width (const dealii::types::global_dof_index cell_index) const;

    double get_filter_width_from_poly_degree (
        const dealii::types::global_dof_index cell_index,
        const int cell_poly_degree) const;

    virtual dealii::Tensor<2,dim,real> compute_SGS_stress_tensor (
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const = 0;

    virtual dealii::Tensor<1,dim,real> compute_SGS_heat_flux (
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const = 0;

    virtual dealii::Tensor<2,dim,FadType> compute_SGS_stress_tensor_fad (
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const = 0;

    virtual dealii::Tensor<1,dim,FadType> compute_SGS_heat_flux_fad (
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const = 0;

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

protected:
    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;

    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> dissipative_source_term (
        const dealii::Point<dim,real> &pos,
        const dealii::types::global_dof_index cell_index) const;


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

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_Smagorinsky : public LargeEddySimulationBase <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_Smagorinsky(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    const double model_constant; 
    const bool apply_low_reynolds_number_eddy_viscosity_correction; 

    ~LargeEddySimulation_Smagorinsky() {};

    void set_unfiltered_conservative_solution(const std::array<real,nstate> &unfiltered_conservative_solution_) override;

    virtual double get_model_constant_times_filter_width (const dealii::types::global_dof_index cell_index) const;

    virtual double get_model_constant_times_filter_width_squared (const dealii::types::global_dof_index cell_index) const;

    real get_corrected_eddy_viscosity_low_reynolds_number(
        const real uncorrected_eddy_viscosity) const;

    FadType get_corrected_eddy_viscosity_low_reynolds_number_fad(
        const FadType uncorrected_eddy_viscosity) const;

    real get_scaled_fluid_kinematic_viscosity_from_unfiltered_solution() const;

    dealii::Tensor<2,dim,real> compute_SGS_stress_tensor (
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    dealii::Tensor<1,dim,real> compute_SGS_heat_flux (
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    dealii::Tensor<2,dim,FadType> compute_SGS_stress_tensor_fad (
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    dealii::Tensor<1,dim,FadType> compute_SGS_heat_flux_fad (
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    virtual real compute_eddy_viscosity(
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    virtual FadType compute_eddy_viscosity_fad(
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

protected:
    template<typename real2> dealii::Tensor<2,dim,real2> compute_SGS_stress_tensor_templated (
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

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

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

    double scaled_fluid_kinematic_viscosity_from_unfiltered_solution;

private:
    template<typename real2> real2 compute_eddy_viscosity_templated(
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;

    template<typename real2> real2 get_corrected_eddy_viscosity_low_reynolds_number_templated(
        const real2 uncorrected_eddy_viscosity) const;
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_WALE : public LargeEddySimulation_Smagorinsky <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_WALE(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    real compute_eddy_viscosity(
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

    FadType compute_eddy_viscosity_fad(
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

private:
    template<typename real2> real2 compute_eddy_viscosity_templated(
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_Vreman : public LargeEddySimulation_Smagorinsky <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_Vreman(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    real compute_eddy_viscosity(
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

    FadType compute_eddy_viscosity_fad(
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

private:
    template<typename real2> real2 compute_eddy_viscosity_templated(
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_ShearImprovedSmagorinsky : public LargeEddySimulation_Smagorinsky <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_ShearImprovedSmagorinsky(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    ~LargeEddySimulation_ShearImprovedSmagorinsky() {};

    real compute_eddy_viscosity(
        const std::array<real,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

    FadType compute_eddy_viscosity_fad(
        const std::array<FadType,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,FadType>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const override;

private:
    template<typename real2> real2 compute_eddy_viscosity_templated(
        const std::array<real2,nstate> &primitive_soln,
        const std::array<dealii::Tensor<1,dim,real2>,nstate> &primitive_soln_gradient,
        const dealii::types::global_dof_index cell_index) const;
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_VMS : public LargeEddySimulation_Smagorinsky <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_VMS(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        const unsigned int                                        poly_degree,
        const unsigned int                                        poly_degree_large_scales,
        const double                                              mesh_size,
        const double                                              curve_fit_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    ~LargeEddySimulation_VMS() {};

    const double poly_degree;
    const double poly_degree_large_scales;
    const double mesh_size;
    const double curve_fit_constant;

    double get_model_constant_times_filter_width (const dealii::types::global_dof_index cell_index) const override;
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_SmallSmallVMS : public LargeEddySimulation_VMS <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_SmallSmallVMS(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        const unsigned int                                        poly_degree,
        const unsigned int                                        poly_degree_large_scales,
        const double                                              mesh_size,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    ~LargeEddySimulation_SmallSmallVMS() {};
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_AllAllVMS : public LargeEddySimulation_VMS <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_AllAllVMS(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        const unsigned int                                        poly_degree,
        const unsigned int                                        poly_degree_large_scales,
        const double                                              mesh_size,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    ~LargeEddySimulation_AllAllVMS() {};
};

template <int dim, int nspecies, int nstate, typename real>
class LargeEddySimulation_DynamicSmagorinsky : public LargeEddySimulation_Smagorinsky <dim, nspecies, nstate, real>
{
public:
    using thermal_boundary_condition_enum = Parameters::NavierStokesParam::ThermalBoundaryCondition;
    using two_point_num_flux_enum = Parameters::AllParameters::TwoPointNumericalFlux;
    LargeEddySimulation_DynamicSmagorinsky(
        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,
        const double                                              turbulent_prandtl_number,
        const double                                              ratio_of_filter_width_to_cell_size,
        const double                                              model_constant,
        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                                                apply_low_reynolds_number_eddy_viscosity_correction = false);

    ~LargeEddySimulation_DynamicSmagorinsky() {};

    double get_model_constant_times_filter_width_squared (const dealii::types::global_dof_index cell_index) const override;
};

} // Physics namespace
} // PHiLiP namespace

#endif