[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
model_factory.cpp
1 #include <boost/preprocessor/seq/for_each.hpp>
2 
3 #include "parameters/all_parameters.h"
4 #include "parameters/parameters_manufactured_solution.h"
5 
6 #include <deal.II/base/tensor.h>
7 
8 #include "ADTypes.hpp"
9 
10 #include "model_factory.h"
11 #include "manufactured_solution.h"
12 #include "large_eddy_simulation.h"
13 #include "reynolds_averaged_navier_stokes.h"
14 #include "negative_spalart_allmaras_rans_model.h"
15 #include "navier_stokes_model.h"
16 
17 namespace PHiLiP {
18 namespace Physics {
19 
20 template <int dim, int nspecies, int nstate, typename real>
21 std::shared_ptr < ModelBase<dim,nspecies,nstate,real> >
23 ::create_Model(const Parameters::AllParameters *const parameters_input)
24 {
26  PDE_enum pde_type = parameters_input->pde_type;
27 
28  if((pde_type == PDE_enum::physics_model || pde_type == PDE_enum::physics_model_filtered) && nspecies==1) {
29  // generating the manufactured solution from the manufactured solution factory
30  std::shared_ptr< ManufacturedSolutionFunction<dim,nspecies,real> > manufactured_solution_function
32 
33  using Model_enum = Parameters::AllParameters::ModelType;
34  Model_enum model_type = parameters_input->model_type;
35 
36  // ===============================================================================
37  // Model
38  // ===============================================================================
39  // -------------------------------------------------------------------------------
40  // Large Eddy Simulation (LES)
41  // -------------------------------------------------------------------------------
42  if (model_type == Model_enum::large_eddy_simulation) {
43  if constexpr ((nstate==dim+2) && (dim==3)) {
44  // Create Large Eddy Simulation (LES) model based on the SGS model type
46  SGS_enum sgs_model_type = parameters_input->physics_model_param.SGS_model_type;
47  if (sgs_model_type == SGS_enum::smagorinsky) {
48  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49  // Smagorinsky model
50  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51  return std::make_shared < LargeEddySimulation_Smagorinsky<dim,nspecies,nstate,real> > (
52  parameters_input,
53  parameters_input->euler_param.ref_length,
54  parameters_input->euler_param.gamma_gas,
55  parameters_input->euler_param.mach_inf,
56  parameters_input->euler_param.angle_of_attack,
57  parameters_input->euler_param.side_slip_angle,
58  parameters_input->navier_stokes_param.prandtl_number,
59  parameters_input->navier_stokes_param.reynolds_number_inf,
60  parameters_input->navier_stokes_param.use_constant_viscosity,
61  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
62  parameters_input->navier_stokes_param.temperature_inf,
63  parameters_input->physics_model_param.turbulent_prandtl_number,
64  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
65  parameters_input->physics_model_param.smagorinsky_model_constant,
66  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
67  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
68  manufactured_solution_function,
69  parameters_input->two_point_num_flux_type,
70  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
71  } else if (sgs_model_type == SGS_enum::wall_adaptive_local_eddy_viscosity) {
72  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73  // WALE (Wall-Adapting Local Eddy-viscosity) eddy viscosity model
74  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75  return std::make_shared < LargeEddySimulation_WALE<dim,nspecies,nstate,real> > (
76  parameters_input,
77  parameters_input->euler_param.ref_length,
78  parameters_input->euler_param.gamma_gas,
79  parameters_input->euler_param.mach_inf,
80  parameters_input->euler_param.angle_of_attack,
81  parameters_input->euler_param.side_slip_angle,
82  parameters_input->navier_stokes_param.prandtl_number,
83  parameters_input->navier_stokes_param.reynolds_number_inf,
84  parameters_input->navier_stokes_param.use_constant_viscosity,
85  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
86  parameters_input->navier_stokes_param.temperature_inf,
87  parameters_input->physics_model_param.turbulent_prandtl_number,
88  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
89  parameters_input->physics_model_param.WALE_model_constant,
90  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
91  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
92  manufactured_solution_function,
93  parameters_input->two_point_num_flux_type,
94  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
95  } else if (sgs_model_type == SGS_enum::vreman) {
96  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97  // Vreman eddy viscosity model
98  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99  return std::make_shared < LargeEddySimulation_Vreman<dim,nspecies,nstate,real> > (
100  parameters_input,
101  parameters_input->euler_param.ref_length,
102  parameters_input->euler_param.gamma_gas,
103  parameters_input->euler_param.mach_inf,
104  parameters_input->euler_param.angle_of_attack,
105  parameters_input->euler_param.side_slip_angle,
106  parameters_input->navier_stokes_param.prandtl_number,
107  parameters_input->navier_stokes_param.reynolds_number_inf,
108  parameters_input->navier_stokes_param.use_constant_viscosity,
109  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
110  parameters_input->navier_stokes_param.temperature_inf,
111  parameters_input->physics_model_param.turbulent_prandtl_number,
112  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
113  parameters_input->physics_model_param.vreman_model_constant,
114  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
115  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
116  manufactured_solution_function,
117  parameters_input->two_point_num_flux_type,
118  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
119  } else if (sgs_model_type == SGS_enum::shear_improved_smagorinsky) {
120  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121  // Shear-improved Smagorinsky model
122  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123  return std::make_shared < LargeEddySimulation_ShearImprovedSmagorinsky<dim,nspecies,nstate,real> > (
124  parameters_input,
125  parameters_input->euler_param.ref_length,
126  parameters_input->euler_param.gamma_gas,
127  parameters_input->euler_param.mach_inf,
128  parameters_input->euler_param.angle_of_attack,
129  parameters_input->euler_param.side_slip_angle,
130  parameters_input->navier_stokes_param.prandtl_number,
131  parameters_input->navier_stokes_param.reynolds_number_inf,
132  parameters_input->navier_stokes_param.use_constant_viscosity,
133  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
134  parameters_input->navier_stokes_param.temperature_inf,
135  parameters_input->physics_model_param.turbulent_prandtl_number,
136  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
137  parameters_input->physics_model_param.smagorinsky_model_constant,
138  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
139  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
140  manufactured_solution_function,
141  parameters_input->two_point_num_flux_type,
142  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
143  } else if ((sgs_model_type == SGS_enum::small_small_variational_multiscale) ||
144  (sgs_model_type == SGS_enum::all_all_variational_multiscale)) {
145  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146  // Variational multiscale (VMS) eddy viscosity models
147  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148  const double domain_left = parameters_input->flow_solver_param.grid_left_bound;
149  const double domain_right = parameters_input->flow_solver_param.grid_right_bound;
150  const int number_of_cells_per_direction = parameters_input->flow_solver_param.number_of_grid_elements_per_dimension;
151  const double mesh_size = (domain_right - domain_left)/((double)number_of_cells_per_direction);
152 
153  if (sgs_model_type == SGS_enum::small_small_variational_multiscale) {
154  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155  // Small-Small Variational multiscale (SmallSmallVMS) eddy viscosity model
156  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157  return std::make_shared < LargeEddySimulation_SmallSmallVMS<dim,nspecies,nstate,real> > (
158  parameters_input,
159  parameters_input->euler_param.ref_length,
160  parameters_input->euler_param.gamma_gas,
161  parameters_input->euler_param.mach_inf,
162  parameters_input->euler_param.angle_of_attack,
163  parameters_input->euler_param.side_slip_angle,
164  parameters_input->navier_stokes_param.prandtl_number,
165  parameters_input->navier_stokes_param.reynolds_number_inf,
166  parameters_input->navier_stokes_param.use_constant_viscosity,
167  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
168  parameters_input->navier_stokes_param.temperature_inf,
169  parameters_input->physics_model_param.turbulent_prandtl_number,
170  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
171  parameters_input->physics_model_param.smagorinsky_model_constant,
172  parameters_input->flow_solver_param.poly_degree,
173  parameters_input->physics_model_param.poly_degree_max_large_scales,
174  mesh_size,
175  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
176  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
177  manufactured_solution_function,
178  parameters_input->two_point_num_flux_type,
179  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
180  } else if (sgs_model_type == SGS_enum::all_all_variational_multiscale) {
181  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182  // All-All Variational multiscale (AllAllVMS) eddy viscosity model
183  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184  return std::make_shared < LargeEddySimulation_AllAllVMS<dim,nspecies,nstate,real> > (
185  parameters_input,
186  parameters_input->euler_param.ref_length,
187  parameters_input->euler_param.gamma_gas,
188  parameters_input->euler_param.mach_inf,
189  parameters_input->euler_param.angle_of_attack,
190  parameters_input->euler_param.side_slip_angle,
191  parameters_input->navier_stokes_param.prandtl_number,
192  parameters_input->navier_stokes_param.reynolds_number_inf,
193  parameters_input->navier_stokes_param.use_constant_viscosity,
194  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
195  parameters_input->navier_stokes_param.temperature_inf,
196  parameters_input->physics_model_param.turbulent_prandtl_number,
197  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
198  parameters_input->physics_model_param.smagorinsky_model_constant,
199  parameters_input->flow_solver_param.poly_degree,
200  parameters_input->physics_model_param.poly_degree_max_large_scales,
201  mesh_size,
202  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
203  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
204  manufactured_solution_function,
205  parameters_input->two_point_num_flux_type,
206  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
207  } else {
208  std::cout << "Can't create LargeEddySimulationVMS, invalid SGSModelType type: " << sgs_model_type << std::endl;
209  assert(0==1 && "Can't create LargeEddySimulationVMS, invalid SGSModelType type");
210  return nullptr;
211  }
212  } else if (sgs_model_type == SGS_enum::dynamic_smagorinsky) {
213  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
214  // Dynamic Smagorinsky model
215  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216  return std::make_shared < LargeEddySimulation_DynamicSmagorinsky<dim,nspecies,nstate,real> > (
217  parameters_input,
218  parameters_input->euler_param.ref_length,
219  parameters_input->euler_param.gamma_gas,
220  parameters_input->euler_param.mach_inf,
221  parameters_input->euler_param.angle_of_attack,
222  parameters_input->euler_param.side_slip_angle,
223  parameters_input->navier_stokes_param.prandtl_number,
224  parameters_input->navier_stokes_param.reynolds_number_inf,
225  parameters_input->navier_stokes_param.use_constant_viscosity,
226  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
227  parameters_input->navier_stokes_param.temperature_inf,
228  parameters_input->physics_model_param.turbulent_prandtl_number,
229  parameters_input->physics_model_param.ratio_of_filter_width_to_cell_size,
230  parameters_input->physics_model_param.smagorinsky_model_constant,
231  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
232  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
233  manufactured_solution_function,
234  parameters_input->two_point_num_flux_type,
235  parameters_input->physics_model_param.apply_low_reynolds_number_eddy_viscosity_correction);
236  }
237  else {
238  std::cout << "Can't create LargeEddySimulationBase, invalid SGSModelType type: " << sgs_model_type << std::endl;
239  assert(0==1 && "Can't create LargeEddySimulationBase, invalid SGSModelType type");
240  return nullptr;
241  }
242  }
243  else {
244  // LES does not exist for nstate!=(dim+2) || dim!=3
245  std::cout << "Can't create LES for nstate!=(dim+2) or dim!=3" << std::endl;
246  assert(0==1 && "Can't create LES for nstate!=(dim+2) or dim!=3");
247  manufactured_solution_function = nullptr;
248  return nullptr;
249  }
250  }
251  // -------------------------------------------------------------------------------
252  // Navier-Stokes model
253  // -------------------------------------------------------------------------------
254  else if (model_type == Model_enum::navier_stokes_model) {
255  if constexpr ((nstate==dim+2) && (dim==3)) {
256  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257  // Navier-Stokes with model source terms (e.g. channel flow)
258  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259  return std::make_shared < NavierStokesWithModelSourceTerms<dim,nspecies,nstate,real> > (
260  parameters_input,
261  parameters_input->euler_param.ref_length,
262  parameters_input->euler_param.gamma_gas,
263  parameters_input->euler_param.mach_inf,
264  parameters_input->euler_param.angle_of_attack,
265  parameters_input->euler_param.side_slip_angle,
266  parameters_input->navier_stokes_param.prandtl_number,
267  parameters_input->navier_stokes_param.reynolds_number_inf,
268  parameters_input->navier_stokes_param.use_constant_viscosity,
269  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
270  parameters_input->navier_stokes_param.temperature_inf,
271  parameters_input->flow_solver_param.relaxation_coefficient_for_turbulent_channel_flow_source_term,
272  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
273  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
274  manufactured_solution_function,
275  parameters_input->two_point_num_flux_type);
276  } else {
277  // Navier-Stokes model does not exist for nstate!=(dim+2) || dim!=3
278  std::cout << "Can't create Navier-Stokes model for nstate!=(dim+2) or dim!=3" << std::endl;
279  assert(0==1 && "Can't create Navier-Stokes model for nstate!=(dim+2) or dim!=3");
280  manufactured_solution_function = nullptr;
281  return nullptr;
282  }
283  }
284  // -------------------------------------------------------------------------------
285  // Reynolds-Averaged Navier-Stokes (RANS) + RANS model
286  // -------------------------------------------------------------------------------
287  else if (model_type == Model_enum::reynolds_averaged_navier_stokes) {
289  RANSModel_enum rans_model_type = parameters_input->physics_model_param.RANS_model_type;
290  // Create Reynolds-Averaged Navier-Stokes (RANS) model with one-equation model
291  if(rans_model_type == RANSModel_enum::SA_negative){
292  if constexpr (nstate==dim+3) {
293  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
294  // SA negative model
295  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
296  return std::make_shared < ReynoldsAveragedNavierStokes_SAneg<dim,nspecies,nstate,real> > (
297  parameters_input,
298  parameters_input->euler_param.ref_length,
299  parameters_input->euler_param.gamma_gas,
300  parameters_input->euler_param.mach_inf,
301  parameters_input->euler_param.angle_of_attack,
302  parameters_input->euler_param.side_slip_angle,
303  parameters_input->navier_stokes_param.prandtl_number,
304  parameters_input->navier_stokes_param.reynolds_number_inf,
305  parameters_input->navier_stokes_param.use_constant_viscosity,
306  parameters_input->navier_stokes_param.nondimensionalized_constant_viscosity,
307  parameters_input->physics_model_param.turbulent_prandtl_number,
308  parameters_input->navier_stokes_param.temperature_inf,
309  parameters_input->navier_stokes_param.nondimensionalized_isothermal_wall_temperature,
310  parameters_input->navier_stokes_param.thermal_boundary_condition_type,
311  manufactured_solution_function,
312  parameters_input->two_point_num_flux_type);
313  }
314  else {
315  // SA negative does not exist for nstate!=(dim+3)
316  std::cout << "Can't create RANS for nstate!=(dim+2) or dim!=3" << std::endl;
317  assert(0==1 && "Can't create RANS for nstate!=(dim+2) or dim!=3");
318  manufactured_solution_function = nullptr;
319  return nullptr;
320  }
321  }
322  else {
323  std::cout << "Can't create ReynoldsAveragedNavierStokesBase, invalid RANSModelType type: " << rans_model_type << std::endl;
324  assert(0==1 && "Can't create ReynoldsAveragedNavierStokesBase, invalid RANSModelType type");
325  return nullptr;
326  }
327  }
328  else {
329  // prevent warnings for dim=3,nstate=4, etc.
330  // to avoid "unused variable" warnings
331  std::cout << "Can't create ModelBase, invalid ModelType type: " << model_type << std::endl;
332  assert(0==1 && "Can't create ModelBase, invalid ModelType type");
333  manufactured_solution_function = nullptr;
334  return nullptr;
335  }
336  }
337  else {
338  return nullptr;
339  }
340 }
341 
342 //----------------------------------------------------------------
343 //----------------------------------------------------------------
344 //----------------------------------------------------------------
345 // Instantiate explicitly
346 #if PHILIP_SPECIES==1
347  // Define a sequence of indices representing the range of nstate
348  #define POSSIBLE_NSTATE (1)(2)(3)(4)(5)(6)(8)
349 
350  // Define a macro to instantiate functions for a specific nstate
351  #define INSTANTIATE_FOR_NSTATE(r, data, nstate) \
352  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, double>; \
353  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, FadType>; \
354  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, RadType>; \
355  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, FadFadType>; \
356  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, nstate, RadFadType>;
357  BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_FOR_NSTATE, _, POSSIBLE_NSTATE)
358 #else
359  #define POSSIBLE_TYPE (double)(FadType)(RadType)(FadFadType)(RadFadType)
360  #define INSTANTIATE_TYPES(r, data, type) \
361  template class ModelFactory<PHILIP_DIM, PHILIP_SPECIES, PHILIP_DIM+PHILIP_SPECIES+1, type>;
362  BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_TYPES, _, POSSIBLE_TYPE)
363 #endif
364 } // Physics namespace
365 } // PHiLiP namespace
366 
PartialDifferentialEquation pde_type
Store the PDE type to be solved.
unsigned int number_of_grid_elements_per_dimension
Number of grid elements per dimension for hyper_cube mesh based cases.
FlowSolverParam flow_solver_param
Contains the parameters for simulation cases (flow solver test)
PartialDifferentialEquation
Possible Partial Differential Equations to solve.
Files for the baseline physics.
Definition: ADTypes.hpp:10
EulerParam euler_param
Contains parameters for the Euler equations non-dimensionalization.
ModelType
Types of models available.
Main parameter class that contains the various other sub-parameter classes.
double grid_left_bound
Left bound of domain for hyper_cube mesh based cases.
SubGridScaleModel SGS_model_type
Store the SubGridScale (SGS) model type.
double ref_length
Reference length.
ReynoldsAveragedNavierStokesModel
Types of Reynolds-averaged Navier-Stokes (RANS) models that can be used.
ReynoldsAveragedNavierStokesModel RANS_model_type
Store the Reynolds-averaged Navier-Stokes (RANS) model type.
static std::shared_ptr< ManufacturedSolutionFunction< dim, nspecies, real > > create_ManufacturedSolution(Parameters::AllParameters const *const param, int nstate)
Construct Manufactured solution object from global parameter file.
SubGridScaleModel
Types of sub-grid scale (SGS) models that can be used.
double grid_right_bound
Right bound of domain for hyper_cube mesh based cases.
ModelType model_type
Store the model type.
static std::shared_ptr< ModelBase< dim, nspecies, nstate, real > > create_Model(const Parameters::AllParameters *const parameters_input)
Factory to return the correct model given input parameters.
PhysicsModelParam physics_model_param
Contains parameters for Physics Model.