VertexCFD  0.0-dev
VertexCFD_EquationSet_IncompressibleNavierStokes_impl.hpp
1 #ifndef VERTEXCFD_EQUATIONSET_INCOMPRESSIBLE_NAVIERSTOKES_IMPL_HPP
2 #define VERTEXCFD_EQUATIONSET_INCOMPRESSIBLE_NAVIERSTOKES_IMPL_HPP
3 
4 #include <Panzer_BasisIRLayout.hpp>
5 #include <Panzer_IntegrationRule.hpp>
6 
7 #include <Panzer_Integrator_BasisTimesScalar.hpp>
8 #include <Panzer_Integrator_GradBasisDotVector.hpp>
9 #include <Panzer_String_Utilities.hpp>
10 
11 #include <Phalanx_FieldManager.hpp>
12 
13 #include <Teuchos_ParameterList.hpp>
14 #include <Teuchos_RCP.hpp>
15 #include <Teuchos_StandardParameterEntryValidators.hpp>
16 
17 namespace VertexCFD
18 {
19 namespace EquationSet
20 {
21 //---------------------------------------------------------------------------//
22 template<class EvalType>
23 IncompressibleNavierStokes<EvalType>::IncompressibleNavierStokes(
24  const Teuchos::RCP<Teuchos::ParameterList>& params,
25  const int& default_integration_order,
26  const panzer::CellData& cell_data,
27  const Teuchos::RCP<panzer::GlobalData>& global_data,
28  const bool build_transient_support)
29  : panzer::EquationSet_DefaultImpl<EvalType>(params,
30  default_integration_order,
31  cell_data,
32  global_data,
33  build_transient_support)
34 {
35  // This equation set is always transient.
36  if (!this->buildTransientSupport())
37  {
38  throw std::logic_error(
39  "Transient support is required for solving all equation sets.");
40  }
41 
42  // Get the number of space dimensions.
43  _num_space_dim = cell_data.baseCellDimension();
44  if (!(_num_space_dim == 2 || _num_space_dim == 3))
45  {
46  throw std::runtime_error("Number of space dimensions not supported");
47  }
48 
49  // Set default parameter values and validate the inputs.
50  Teuchos::ParameterList valid_parameters;
51  this->setDefaultValidParameters(valid_parameters);
52  valid_parameters.set(
53  "Model ID", "", "Closure model id associated with this equation set");
54  valid_parameters.set("Basis Order", 1, "Order of the basis");
55  valid_parameters.set("Integration Order", 2, "Order of the integration");
56  valid_parameters.set("Continuity Model", "AC", "Continuity model string");
57  valid_parameters.set("Build Viscous Flux", true, "Viscous flux boolean");
58  valid_parameters.set(
59  "Build Temperature Equation", false, "Solve temperature equation");
60  valid_parameters.set("Build Inductionless MHD Equation",
61  false,
62  "Solve electric potential equation");
63  valid_parameters.set(
64  "Build Constant Source", false, "Constant source boolean");
65 
66  valid_parameters.set(
67  "Build Buoyancy Source", false, "Buoyancy source boolean");
68 
69  valid_parameters.set("Build Viscous Heat", false, "Viscous heat boolean");
70 
71  valid_parameters.set(
72  "Build Resistive Flux", false, "Resistive flux boolean");
73 
74  valid_parameters.set(
75  "Build Godunov-Powell Source", false, "Godunov-Powell source boolean");
76 
77  valid_parameters.set(
78  "Build Joule Heating Source", false, "Joule heating source boolean");
79 
80  const auto turbulence_validator = Teuchos::rcp(new Teuchos::StringValidator(
81  Teuchos::tuple<std::string>("No Turbulence Model",
82  "Spalart-Allmaras",
83  "Chien K-Epsilon",
84  "Standard K-Epsilon",
85  "Realizable K-Epsilon",
86  "K-Omega",
87  "SST K-Omega",
88  "K-Tau",
89  "WALE")));
90 
91  valid_parameters.set("Turbulence Model",
92  "No Turbulence Model",
93  "Turbulence model choice",
94  turbulence_validator);
95 
96  const auto stabilization_validator
97  = Teuchos::rcp(new Teuchos::StringValidator(
98  Teuchos::tuple<std::string>("No Stabilization Method",
99  "Taylor-Hood",
100  "Grad-Div",
101  "Taylor-Hood and Grad-Div")));
102 
103  valid_parameters.set("Stabilization Method",
104  "No Stabilization Method",
105  "Stabilization method choice",
106  stabilization_validator);
107 
108  params->validateParametersAndSetDefaults(valid_parameters);
109 
110  // Extract parameters.
111  const int basis_order = params->get<int>("Basis Order", 1);
112  const int integration_order
113  = params->get<int>("Integration Order", basis_order + 1);
114  const std::string model_id = params->get<std::string>("Model ID");
115  const std::string continuity_model
116  = params->get<std::string>("Continuity Model");
117  _build_viscous_flux = params->get<bool>("Build Viscous Flux", true);
118  _build_temp_equ = params->get<bool>("Build Temperature Equation", false);
119  _build_ind_less_equ
120  = params->get<bool>("Build Inductionless MHD Equation", false);
121  _build_constant_source = params->get<bool>("Build Constant Source", false);
122  _build_buoyancy_source = params->get<bool>("Build Buoyancy Source", false);
123  _build_viscous_heat = params->get<bool>("Build Viscous Heat", false);
124  _turbulence_model = params->get<std::string>("Turbulence Model");
125  _stabilization_method = params->get<std::string>("Stabilization Method");
126  _build_godunov_powell_source
127  = params->get<bool>("Build Godunov-Powell Source", false);
128 #ifndef VERTEXCFD_ENABLE_FULL_INDUCTION_MHD
129  if (_build_godunov_powell_source)
130  {
131  throw std::runtime_error(
132  "\"Build Godunov-Powell Source\" set to \"true\", "
133  "but sovler was built without full induction model. Set\n"
134  " cmake -DVertexCFD_ENABLE_FULL_INDUCTION_MHD=ON.");
135  }
136 #endif
137  _build_joule_heating_source
138  = params->get<bool>("Build Joule Heating Source", false);
139 
140  if (continuity_model == "AC")
141  {
142  _continuity_model = ConModel::AC;
143  }
144  else if (continuity_model == "EDAC")
145  {
146  _continuity_model = ConModel::EDAC;
147  }
148  else if (continuity_model == "EDACTempNC")
149  {
150  _continuity_model = ConModel::EDACTempNC;
151  }
152  else
153  {
154  const std::string msg
155  = "Continuity Model' must be AC, EDAC, or EDACTempNC.";
156  throw std::runtime_error(msg);
157  }
158 
159  if (_build_buoyancy_source && !_build_temp_equ)
160  {
161  std::string msg = "Build Buoyancy Source set to TRUE,\n";
162  msg += "but Build Temperature Equation set to FALSE.\n";
163  msg += "Please enable temperature equation to solve with buoyancy source.";
164  throw std::runtime_error(msg);
165  }
166 
167  if (_build_viscous_heat && !_build_temp_equ)
168  {
169  std::string msg = "Build Viscous Heat set to TRUE.\n";
170  msg += "but Build Temperature Equation set to FALSE.\n";
171  msg += "Please enable temperature equation to solve with buoyancy source.";
172  throw std::runtime_error(msg);
173  }
174 
175  // Initialize equation names and variable names for NS equations
176  _equ_dof_ns_pair.insert({"continuity", "lagrange_pressure"});
177  for (int d = 0; d < _num_space_dim; ++d)
178  {
179  const std::string ds = std::to_string(d);
180  _equ_dof_ns_pair.insert({"momentum_" + ds, "velocity_" + ds});
181  }
182  if (_build_temp_equ)
183  _equ_dof_ns_pair.insert({"energy", "temperature"});
184 
185  // Initialize equation name and variable name for EP equation
186  if (_build_ind_less_equ)
187  _equ_dof_ep_pair.insert(
188  {"electric_potential_equation", "electric_potential"});
189 
190  // Initialize equation names and variable names for turbulence model (TM)
191  // equations
192  if (std::string::npos != _turbulence_model.find("Spalart-Allmaras"))
193  {
194  _equ_dof_tm_pair.insert(
195  {"spalart_allmaras_equation", "spalart_allmaras_variable"});
196  }
197  else if (std::string::npos != _turbulence_model.find("K-Epsilon"))
198  {
199  _equ_dof_tm_pair.insert(
200  {"turb_kinetic_energy_equation", "turb_kinetic_energy"});
201  _equ_dof_tm_pair.insert(
202  {"turb_dissipation_rate_equation", "turb_dissipation_rate"});
203  }
204  else if (std::string::npos != _turbulence_model.find("K-Omega"))
205  {
206  _equ_dof_tm_pair.insert(
207  {"turb_kinetic_energy_equation", "turb_kinetic_energy"});
208  _equ_dof_tm_pair.insert({"turb_specific_dissipation_rate_equation",
209  "turb_specific_dissipation_rate"});
210  }
211  else if (std::string::npos != _turbulence_model.find("SST"))
212  {
213  _equ_dof_tm_pair.insert(
214  {"turb_kinetic_energy_equation", "turb_kinetic_energy"});
215  _equ_dof_tm_pair.insert({"turb_specific_dissipation_rate_equation",
216  "turb_specific_dissipation_rate"});
217  }
218  else if (std::string::npos != _turbulence_model.find("K-Tau"))
219  {
220  _equ_dof_tm_pair.insert(
221  {"turb_kinetic_energy_equation", "turb_kinetic_energy"});
222  _equ_dof_tm_pair.insert({"turb_specific_dissipation_rate_equation",
223  "turb_specific_dissipation_rate"});
224  }
225 
226  // Functions to set dofs
227  auto set_dofs = [&, this](const std::string& equ_name,
228  const std::string& dof_name) {
229  const std::string residual_name = "RESIDUAL_" + equ_name;
230  const std::string scatter_name = "SCATTER_" + equ_name;
231  const std::string basis_type = "HGrad";
232  int basis_order_dof = basis_order;
233  int integration_order_dof = integration_order;
234  if (std::string::npos != _stabilization_method.find("Taylor-Hood"))
235  {
236  if (std::string::npos != dof_name.find("velocity_"))
237  {
238  basis_order_dof++;
239  integration_order_dof++;
240  }
241  }
242  this->addDOF(dof_name,
243  basis_type,
244  basis_order_dof,
245  integration_order_dof,
246  residual_name,
247  scatter_name);
248  this->addDOFGrad(dof_name);
249  this->addDOFTimeDerivative(dof_name);
250  };
251 
252  // Setup degrees of freedom for NS equations.
253  for (auto it : _equ_dof_ns_pair)
254  {
255  const auto equ_name = it.first;
256  const auto dof_name = it.second;
257  set_dofs(equ_name, dof_name);
258  };
259 
260  // Setup degrees of freedom for EP equations.
261  for (auto it : _equ_dof_ep_pair)
262  {
263  const auto equ_name = it.first;
264  const auto dof_name = it.second;
265  set_dofs(equ_name, dof_name);
266  }
267 
268  // Setup degrees of freedom for TM equations.
269  for (auto it : _equ_dof_tm_pair)
270  {
271  const auto equ_name = it.first;
272  const auto dof_name = it.second;
273  set_dofs(equ_name, dof_name);
274  }
275 
276  // Add closure models and set-up DOFs
277  this->addClosureModel(model_id);
278  this->setupDOFs();
279 }
280 
281 //---------------------------------------------------------------------------//
282 template<class EvalType>
283 void IncompressibleNavierStokes<EvalType>::buildAndRegisterEquationSetEvaluators(
284  PHX::FieldManager<panzer::Traits>& fm,
285  const panzer::FieldLibrary&,
286  const Teuchos::ParameterList&) const
287 {
288  // Add basis time residuals
289  auto add_basis_time_scalar_residual = [&, this](
290  const std::string& equ_name,
291  const std::string& residual_name,
292  const std::string& dof_name,
293  const double& multiplier,
294  std::vector<std::string>& op_names) {
295  const std::string closure_model_residual = residual_name + "_"
296  + equ_name;
297  const std::string full_residual = "RESIDUAL_" + closure_model_residual;
298  const auto ir = this->getIntRuleForDOF(dof_name);
299  const auto basis = this->getBasisIRLayoutForDOF(dof_name);
300  auto op = Teuchos::rcp(
301  new panzer::Integrator_BasisTimesScalar<EvalType, panzer::Traits>(
302  panzer::EvaluatorStyle::EVALUATES,
303  full_residual,
304  closure_model_residual,
305  *basis,
306  *ir,
307  multiplier));
308  this->template registerEvaluator<EvalType>(fm, op);
309  op_names.push_back(full_residual);
310  };
311 
312  // Add grad-basis time residuals
313  auto add_grad_basis_time_residual = [&, this](
314  const std::string& equ_name,
315  const std::string& residual_name,
316  const std::string& dof_name,
317  const double& multiplier,
318  std::vector<std::string>& op_names) {
319  const std::string closure_model_residual = residual_name + "_"
320  + equ_name;
321  const std::string full_residual = "RESIDUAL_" + closure_model_residual;
322  const auto ir = this->getIntRuleForDOF(dof_name);
323  const auto basis = this->getBasisIRLayoutForDOF(dof_name);
324  auto op = Teuchos::rcp(
325  new panzer::Integrator_GradBasisDotVector<EvalType, panzer::Traits>(
326  panzer::EvaluatorStyle::EVALUATES,
327  full_residual,
328  closure_model_residual,
329  *basis,
330  *ir,
331  multiplier));
332  this->template registerEvaluator<EvalType>(fm, op);
333  op_names.push_back(full_residual);
334  };
335 
336  // Build total residuals for NS equations
337  for (auto it : _equ_dof_ns_pair)
338  {
339  // Define local variables
340  const auto equ_name = it.first;
341  const auto dof_name = it.second;
342  std::vector<std::string> residual_operator_names;
343 
344  // Time derivative residual
345  add_basis_time_scalar_residual(
346  equ_name, "DQDT", dof_name, 1.0, residual_operator_names);
347 
348  // Non-conservative source term residual
349  if (_continuity_model == ConModel::EDACTempNC)
350  {
351  if (std::string::npos != equ_name.find("energy"))
352  {
353  add_basis_time_scalar_residual(equ_name,
354  "NON_CONSERVATIVE_SOURCE",
355  dof_name,
356  1.0,
357  residual_operator_names);
358  }
359  else
360  {
361  // Convective flux residuals for momentum and continuity
362  add_grad_basis_time_residual(equ_name,
363  "CONVECTIVE_FLUX",
364  dof_name,
365  -1.0,
366  residual_operator_names);
367  }
368  }
369  else
370  {
371  // Convective flux residual
372  add_grad_basis_time_residual(equ_name,
373  "CONVECTIVE_FLUX",
374  dof_name,
375  -1.0,
376  residual_operator_names);
377  }
378 
379  // Viscous flux residual
380  if (_build_viscous_flux)
381  {
382  add_grad_basis_time_residual(equ_name,
383  "VISCOUS_FLUX",
384  dof_name,
385  1.0,
386  residual_operator_names);
387  }
388 
389  // Constant source residual
390  if (_build_constant_source)
391  {
392  // Constant source for momentum equations
393  if (std::string::npos != equ_name.find("momentum"))
394  {
395  add_basis_time_scalar_residual(equ_name,
396  "CONSTANT_SOURCE",
397  dof_name,
398  -1.0,
399  residual_operator_names);
400  }
401 
402  // Constant source for energy equation
403  if (std::string::npos != equ_name.find("energy"))
404  {
405  add_basis_time_scalar_residual(equ_name,
406  "CONSTANT_SOURCE",
407  dof_name,
408  -1.0,
409  residual_operator_names);
410  }
411  }
412 
413  // Buoyancy source residual
414  if (_build_buoyancy_source)
415  {
416  add_basis_time_scalar_residual(equ_name,
417  "BUOYANCY_SOURCE",
418  dof_name,
419  -1.0,
420  residual_operator_names);
421  }
422 
423  // Viscous heating source residual
424  if (_build_viscous_heat)
425  {
426  add_basis_time_scalar_residual(equ_name,
427  "VISCOUS_HEAT_SOURCE",
428  dof_name,
429  -1.0,
430  residual_operator_names);
431  }
432 
433  // Inductionless MHD source residuals
434  if (_build_ind_less_equ)
435  {
436  // Lorentz force for momentum equations
437  if (std::string::npos != equ_name.find("momentum"))
438  {
439  add_basis_time_scalar_residual(equ_name,
440  "VOLUMETRIC_SOURCE",
441  dof_name,
442  -1.0,
443  residual_operator_names);
444  }
445 
446  // Joule heating for energy equation
447  if (_build_joule_heating_source
448  && std::string::npos != equ_name.find("energy"))
449  {
450  add_basis_time_scalar_residual(equ_name,
451  "VOLUMETRIC_SOURCE",
452  dof_name,
453  -1.0,
454  residual_operator_names);
455  }
456  }
457 
458  // Godunov-Powell source for momentum equations in full induction MHD
459  if (_build_godunov_powell_source
460  && std::string::npos != equ_name.find("momentum"))
461  {
462  add_basis_time_scalar_residual(equ_name,
463  "GODUNOV_POWELL_SOURCE",
464  dof_name,
465  -1.0,
466  residual_operator_names);
467  }
468 
469  // Build and register residuals
470  this->buildAndRegisterResidualSummationEvaluator(
471  fm, dof_name, residual_operator_names, "RESIDUAL_" + equ_name);
472  }
473 
474  // Build total residual for EP equation
475  for (auto it : _equ_dof_ep_pair)
476  {
477  // Define local variables
478  const auto equ_name = it.first;
479  const auto dof_name = it.second;
480  std::vector<std::string> residual_operator_names;
481 
482  // Cross-product flux and diffusion flux residuals
483  add_grad_basis_time_residual(equ_name,
484  "ELECTRIC_POTENTIAL_FLUX",
485  dof_name,
486  1.0,
487  residual_operator_names);
488 
489  // Build and register residuals
490  this->buildAndRegisterResidualSummationEvaluator(
491  fm, dof_name, residual_operator_names, "RESIDUAL_" + equ_name);
492  }
493 
494  // Build total residuals for TM equations
495  for (auto it : _equ_dof_tm_pair)
496  {
497  // Define local variables
498  const auto eq_name = it.first;
499  const auto dof_name = it.second;
500  std::vector<std::string> residual_operator_names;
501 
502  // Add time, convective, viscous and source residuals
503  add_basis_time_scalar_residual(
504  eq_name, "DQDT", dof_name, 1.0, residual_operator_names);
505  add_grad_basis_time_residual(
506  eq_name, "CONVECTIVE_FLUX", dof_name, -1.0, residual_operator_names);
507  add_grad_basis_time_residual(
508  eq_name, "DIFFUSION_FLUX", dof_name, 1.0, residual_operator_names);
509  add_basis_time_scalar_residual(
510  eq_name, "SOURCE", dof_name, -1.0, residual_operator_names);
511 
512  // Build and register residuals
513  this->buildAndRegisterResidualSummationEvaluator(
514  fm, dof_name, residual_operator_names, "RESIDUAL_" + eq_name);
515  }
516 }
517 
518 //---------------------------------------------------------------------------//
519 
520 } // end namespace EquationSet
521 } // end namespace VertexCFD
522 
523 #endif // end VERTEXCFD_EQUATIONSET_INCOMPRESSIBLE_NAVIERSTOKES_IMPL_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20