physics_factory.cpp

#include <boost/preprocessor/seq/for_each.hpp>

#include "parameters/all_parameters.h"
#include "parameters/parameters_manufactured_solution.h"

#include <deal.II/base/tensor.h>

#include "ADTypes.hpp"

#include "physics_factory.h"
#include "manufactured_solution.h"
#include "physics.h"
#include "convection_diffusion.h"
#include "burgers.h"
#include "burgers_rewienski.h"
#include "euler.h"
#include "mhd.h"
#include "navier_stokes.h"
#include "physics_model.h"
#include "real_gas.h"

namespace PHiLiP {
namespace Physics {

template <int dim, int nspecies, int nstate, typename real>
std::shared_ptr < PhysicsBase<dim,nspecies,nstate,real> >
PhysicsFactory<dim,nspecies,nstate,real>
::create_Physics(const Parameters::AllParameters               *const parameters_input,
                 std::shared_ptr< ModelBase<dim,nspecies,nstate,real> > model_input)
{
    using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation;
    PDE_enum pde_type = parameters_input->pde_type;

    return create_Physics(parameters_input, pde_type, model_input);
}

template <int dim, int nspecies, int nstate, typename real>
std::shared_ptr < PhysicsBase<dim,nspecies,nstate,real> >
PhysicsFactory<dim,nspecies,nstate,real>
::create_Physics(const Parameters::AllParameters                              *const parameters_input,
                 const Parameters::AllParameters::PartialDifferentialEquation pde_type,
                 std::shared_ptr< ModelBase<dim,nspecies,nstate,real> >                model_input)
{
    using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation;

    if constexpr(nspecies==1) {
        // generating the manufactured solution from the manufactured solution factory
        std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> >  manufactured_solution_function 
            = ManufacturedSolutionFactory<dim,nspecies,real>::create_ManufacturedSolution(parameters_input, nstate);

        // setting the diffusion tensor and advection vectors from parameters (if needed)
        const dealii::Tensor<2,3,double> diffusion_tensor      = parameters_input->manufactured_convergence_study_param.manufactured_solution_param.diffusion_tensor;
        const dealii::Tensor<1,3,double> advection_vector      = parameters_input->manufactured_convergence_study_param.manufactured_solution_param.advection_vector;
        const double                     diffusion_coefficient = parameters_input->manufactured_convergence_study_param.manufactured_solution_param.diffusion_coefficient;

        if (pde_type == PDE_enum::advection || pde_type == PDE_enum::advection_vector) {
            if constexpr (nstate<=2) 
                return std::make_shared < ConvectionDiffusion<dim,nspecies,nstate,real> >(
                    parameters_input,
                    true, false,
                    diffusion_tensor, advection_vector, diffusion_coefficient,
                    manufactured_solution_function);
        } else if (pde_type == PDE_enum::diffusion) {
            if constexpr (nstate==1) 
                return std::make_shared < ConvectionDiffusion<dim,nspecies,nstate,real> >(
                    parameters_input,
                    false, true,
                    diffusion_tensor, advection_vector, diffusion_coefficient,
                    manufactured_solution_function,
                    parameters_input->test_type);
        } else if (pde_type == PDE_enum::convection_diffusion) {
            if constexpr (nstate==1) 
                return std::make_shared < ConvectionDiffusion<dim,nspecies,nstate,real> >(
                    parameters_input,
                    true, true,
                    diffusion_tensor, advection_vector, diffusion_coefficient,
                    manufactured_solution_function,
                    parameters_input->test_type);
        } else if (pde_type == PDE_enum::burgers_inviscid) {
            if constexpr (nstate==dim) 
                return std::make_shared < Burgers<dim,nspecies,nstate,real> >(
                    parameters_input,
                    parameters_input->burgers_param.diffusion_coefficient,
                    true, false,
                    diffusion_tensor, 
                    manufactured_solution_function,
                    parameters_input->test_type);
        } else if (pde_type == PDE_enum::burgers_viscous) {
            if constexpr (nstate==dim)
                return std::make_shared < Burgers<dim,nspecies,nstate,real> >(
                    parameters_input,
                    parameters_input->burgers_param.diffusion_coefficient,
                    true, true,
                    diffusion_tensor,
                    manufactured_solution_function);
        } else if (pde_type == PDE_enum::burgers_rewienski) {
            if constexpr (nstate==dim)
                return std::make_shared < BurgersRewienski<dim,nspecies,nstate,real> >(
                    parameters_input,
                    parameters_input->burgers_param.rewienski_a,
                    parameters_input->burgers_param.rewienski_b,
                    parameters_input->burgers_param.rewienski_manufactured_solution,
                    true,
                    false,
                    diffusion_tensor,
                    manufactured_solution_function);
        } else if (pde_type == PDE_enum::euler) {
            if constexpr (nstate==dim+2) {
                return std::make_shared < Euler<dim,nspecies,nstate,real> > (
                    parameters_input,
                    parameters_input->euler_param.ref_length,
                    parameters_input->euler_param.gamma_gas,
                    parameters_input->euler_param.mach_inf,
                    parameters_input->euler_param.angle_of_attack,
                    parameters_input->euler_param.side_slip_angle,
                    manufactured_solution_function,
                    parameters_input->two_point_num_flux_type);
            }
        } else if (pde_type == PDE_enum::mhd) {
            if constexpr (nstate == 8) 
                return std::make_shared < MHD<dim,nspecies,nstate,real> > (
                    parameters_input,
                    parameters_input->euler_param.gamma_gas,
                    diffusion_tensor, 
                    manufactured_solution_function);
        } else if (pde_type == PDE_enum::navier_stokes) {
            if constexpr (nstate==dim+2) {
                return std::make_shared < NavierStokes<dim,nspecies,nstate,real> > (
                    parameters_input,
                    parameters_input->euler_param.ref_length,
                    parameters_input->euler_param.gamma_gas,
                    parameters_input->euler_param.mach_inf,
                    parameters_input->euler_param.angle_of_attack,
                    parameters_input->euler_param.side_slip_angle,
                    parameters_input->navier_stokes_param.prandtl_number,
                    parameters_input->navier_stokes_param.reynolds_number_inf,
                    parameters_input->navier_stokes_param.use_constant_viscosity,
                    parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
                    parameters_input->navier_stokes_param.temperature_inf,
                    parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
                    parameters_input->navier_stokes_param.thermal_boundary_condition_type,
                    manufactured_solution_function,
                    parameters_input->two_point_num_flux_type);
            }
        } else if (pde_type == PDE_enum::navier_stokes_channel_flow_constant_source_term) {
            if constexpr (nstate==dim+2) {
                const double domain_length_y_direction = parameters_input->flow_solver_param.turbulent_channel_domain_length_y_direction;
                const double half_channel_height = domain_length_y_direction/2.0;
                return std::make_shared < NavierStokes_ChannelFlowConstantSourceTerm<dim,nspecies,nstate,real> > (
                    parameters_input,
                    parameters_input->euler_param.ref_length,
                    parameters_input->euler_param.gamma_gas,
                    parameters_input->euler_param.mach_inf,
                    parameters_input->euler_param.angle_of_attack,
                    parameters_input->euler_param.side_slip_angle,
                    parameters_input->navier_stokes_param.prandtl_number,
                    parameters_input->navier_stokes_param.reynolds_number_inf,
                    parameters_input->navier_stokes_param.use_constant_viscosity,
                    parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
                    parameters_input->flow_solver_param.turbulent_channel_friction_velocity_reynolds_number,
                    half_channel_height,
                    parameters_input->navier_stokes_param.temperature_inf,
                    parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
                    parameters_input->navier_stokes_param.thermal_boundary_condition_type,
                    manufactured_solution_function,
                    parameters_input->two_point_num_flux_type);
            }
        } else if (pde_type == PDE_enum::navier_stokes_channel_flow_constant_source_term_wall_model) {
            if constexpr (nstate==dim+2) {
                if(parameters_input->flow_solver_param.flow_case_type != Parameters::FlowSolverParam::FlowCaseType::channel_flow) {
                    std::cout << "Invalid flow case, can only create NavierStokes_ChannelFlowConstantSourceTerm_WallModel for ChannelFlow." << std::endl;
                    return nullptr;
                }
                const double domain_length_y_direction = parameters_input->flow_solver_param.turbulent_channel_domain_length_y_direction;
                const double half_channel_height = domain_length_y_direction/2.0;
                const int number_of_cells_y_direction = parameters_input->flow_solver_param.turbulent_channel_number_of_cells_y_direction;
                const double uniform_spacing_y_direction = domain_length_y_direction/double(number_of_cells_y_direction);
                return std::make_shared < NavierStokes_ChannelFlowConstantSourceTerm_WallModel<dim,nspecies,nstate,real> > (
                    parameters_input,
                    parameters_input->euler_param.ref_length,
                    parameters_input->euler_param.gamma_gas,
                    parameters_input->euler_param.mach_inf,
                    parameters_input->euler_param.angle_of_attack,
                    parameters_input->euler_param.side_slip_angle,
                    parameters_input->navier_stokes_param.prandtl_number,
                    parameters_input->navier_stokes_param.reynolds_number_inf,
                    parameters_input->navier_stokes_param.use_constant_viscosity,
                    parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
                    parameters_input->flow_solver_param.turbulent_channel_friction_velocity_reynolds_number,
                    half_channel_height,
                    uniform_spacing_y_direction,
                    parameters_input->navier_stokes_param.temperature_inf,
                    parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
                    parameters_input->navier_stokes_param.thermal_boundary_condition_type,
                    manufactured_solution_function,
                    parameters_input->two_point_num_flux_type);
            }
        } else if (pde_type == PDE_enum::physics_model || pde_type == PDE_enum::physics_model_filtered) {
            if constexpr (nstate>=dim+2) {
                return create_Physics_Model(parameters_input,
                                            manufactured_solution_function,
                                            model_input);
            }
        }
        // prevent warnings for dim=3,nstate=4, etc.
        (void) diffusion_tensor;
        (void) advection_vector;
        (void) diffusion_coefficient;
    } else if (pde_type == PDE_enum::real_gas) {
        if constexpr (nstate==dim+nspecies+1) {
            return std::make_shared < RealGas<dim,nspecies,nstate,real> > (parameters_input);
        }
    }
    std::cout << "Can't create PhysicsBase, invalid PDE type: " << pde_type << std::endl;
    assert(0==1 && "Can't create PhysicsBase, invalid PDE type");
    return nullptr;
}

template <int dim, int nspecies, int nstate, typename real>
std::shared_ptr < PhysicsBase<dim,nspecies,nstate,real> >
PhysicsFactory<dim,nspecies,nstate,real>
::create_Physics_Model(const Parameters::AllParameters                           *const parameters_input,
                       std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function,
                       std::shared_ptr< ModelBase<dim,nspecies,nstate,real> >             model_input)
{
    using PDE_enum = Parameters::AllParameters::PartialDifferentialEquation;
    const PDE_enum pde_type = parameters_input->pde_type;

    using Model_enum = Parameters::AllParameters::ModelType;
    const Model_enum model_type = parameters_input->model_type;

    using RANSModel_enum = Parameters::PhysicsModelParam::ReynoldsAveragedNavierStokesModel;
    const RANSModel_enum rans_model_type = parameters_input->physics_model_param.RANS_model_type;

    using FlowCase_enum = Parameters::FlowSolverParam::FlowCaseType;
    const FlowCase_enum flow_case_type = parameters_input->flow_solver_param.flow_case_type;

    // ===============================================================================
    // Physics Model
    // ===============================================================================
    
    // Create baseline physics object
    PDE_enum baseline_physics_type;

    // Flag to signal non-zero diffusion
    bool has_nonzero_diffusion;

    // Flag to signal non-zero physical source
    bool has_nonzero_physical_source;

    // -------------------------------------------------------------------------------
    // Large Eddy Simulation (LES)
    // -------------------------------------------------------------------------------
    if (model_type==Model_enum::large_eddy_simulation || model_type==Model_enum::navier_stokes_model) {
        has_nonzero_diffusion = true; // because of SGS model term (initialized as true)
        has_nonzero_physical_source = false; // no physical source terms by default
        if(flow_case_type==FlowCase_enum::channel_flow) {
            has_nonzero_physical_source = true; // forcing function
        }
        if constexpr ((nstate==dim+2) && (dim==3) && nspecies == 1) {
            // Assign baseline physics type (and corresponding nstates) based on the physics model type
            // -- Assign nstates for the baseline physics (constexpr because template parameter)
            constexpr int nstate_baseline_physics = dim+2;
            // -- Assign baseline physics type
            if(parameters_input->physics_model_param.euler_turbulence) {
                baseline_physics_type = PDE_enum::euler;
                if(model_type==Model_enum::navier_stokes_model) has_nonzero_diffusion = false; // no additional diffusion terms
            }
            else {
                baseline_physics_type = PDE_enum::navier_stokes;
            }

            // Create the physics model object in physics
            if (pde_type == PDE_enum::physics_model) {
                return std::make_shared < PhysicsModel<dim,nspecies,nstate,real,nstate_baseline_physics> > (
                    parameters_input,
                    baseline_physics_type,
                    model_input,
                    manufactured_solution_function,
                    has_nonzero_diffusion,
                    has_nonzero_physical_source);
            } else if(pde_type == PDE_enum::physics_model_filtered) {
                return std::make_shared < PhysicsModelFiltered<dim,nspecies,nstate,real,nstate_baseline_physics> > (
                    parameters_input,
                    baseline_physics_type,
                    model_input,
                    manufactured_solution_function,
                    has_nonzero_diffusion,
                    has_nonzero_physical_source);
            }
        }
        else {
            // LES does not exist for nstate!=(dim+2) || dim!=3
            (void) baseline_physics_type;
            (void) has_nonzero_diffusion;
            (void) has_nonzero_physical_source;
            return nullptr;
        }
    }
    // -------------------------------------------------------------------------------
    // Reynolds-Averaged Navier-Stokes (RANS) + RANS model
    // -------------------------------------------------------------------------------
    else if (model_type == Model_enum::reynolds_averaged_navier_stokes) {
        has_nonzero_diffusion = true; // RANS (baseline part) has diffusion terms
        has_nonzero_physical_source = true; // RANS (baseline part) has physical source terms
        if (rans_model_type == RANSModel_enum::SA_negative)
        {
            if constexpr (nstate==dim+3 && nspecies==1) {
                // Assign baseline physics type (and corresponding nstates) based on the physics model type
                // -- Assign nstates for the baseline physics (constexpr because template parameter)
                constexpr int nstate_baseline_physics = dim+2;
                // -- Assign baseline physics type
                if(parameters_input->physics_model_param.euler_turbulence) {
                    baseline_physics_type = PDE_enum::euler;
                }
                else {
                    baseline_physics_type = PDE_enum::navier_stokes;
                }

                // Create the physics model object in physics
                return std::make_shared < PhysicsModel<dim,nspecies,nstate,real,nstate_baseline_physics> > (
                    parameters_input,
                    baseline_physics_type,
                    model_input,
                    manufactured_solution_function,
                    has_nonzero_diffusion,
                    has_nonzero_physical_source);
            }
            else {
                // RANS+one-equation model does not exist for nstate!=(dim+3)
                (void) baseline_physics_type;
                (void) has_nonzero_physical_source;
                std::cout << "Can't create RANS + negative SA model for nstate!=(dim+3). " << std::endl;
                return nullptr;
            }
        }
    }
    else {
        // prevent warnings for dim=3,nstate=4, etc.
        (void) baseline_physics_type;
        (void) has_nonzero_diffusion;
        (void) has_nonzero_physical_source;
    }    
    std::cout << "Can't create PhysicsModel, invalid ModelType type: " << model_type << std::endl;
    assert(0==1 && "Can't create PhysicsModel, invalid ModelType type");
    return nullptr;
}

#if PHILIP_SPECIES==1
    // Define a sequence of indices representing the range of nstate
    #define POSSIBLE_NSTATE (1)(2)(3)(4)(5)(6)(8)

    // Define a macro to instantiate functions for a specific nstate
    #define INSTANTIATE_FOR_NSTATE(r, data, nstate) \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_FOR_NSTATE, _, POSSIBLE_NSTATE)
#else
    #define POSSIBLE_TYPE (double)(FadType)(RadType)(FadFadType)(RadFadType)
    #define INSTANTIATE_TYPES(r, data, type) \
        template class PhysicsFactory<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>;
    BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TYPES, _, POSSIBLE_TYPE)
#endif

} // Physics namespace
} // PHiLiP namespace