VertexCFD  0.0-dev
VertexCFD_Closure_IncompressibleTauSUPG_impl.hpp
1 #ifndef VERTEXCFD_CLOSURE_INCOMPRESSIBLETAUSUPG_IMPL_HPP
2 #define VERTEXCFD_CLOSURE_INCOMPRESSIBLETAUSUPG_IMPL_HPP
3 
4 #include "utils/VertexCFD_Utils_VectorField.hpp"
5 
6 #include <Panzer_HierarchicParallelism.hpp>
7 #include <Panzer_Workset_Utilities.hpp>
8 
9 #include <Teuchos_StandardParameterEntryValidators.hpp>
10 
11 #include <string>
12 
13 namespace VertexCFD
14 {
15 namespace ClosureModel
16 {
17 //---------------------------------------------------------------------------//
18 template<class EvalType, class Traits, int NumSpaceDim>
19 IncompressibleTauSUPG<EvalType, Traits, NumSpaceDim>::IncompressibleTauSUPG(
20  const panzer::IntegrationRule& ir,
21  const Teuchos::ParameterList& fluid_params,
22  const Teuchos::ParameterList& closure_params)
23  : _tau_supg_cont("tau_supg_continuity", ir.dl_scalar)
24  , _tau_supg_mom("tau_supg_momentum", ir.dl_scalar)
25  , _tau_supg_energy("tau_supg_energy", ir.dl_scalar)
26  , _element_length("element_length", ir.dl_vector)
27  , _rho("density", ir.dl_scalar)
28  , _nu("kinematic_viscosity", ir.dl_scalar)
29  , _k("thermal_conductivity", ir.dl_scalar)
30  , _cp("specific_heat_capacity", ir.dl_scalar)
31  , _alpha_cont(0.5)
32  , _alpha_mom(0.5)
33  , _alpha_ene(0.5)
34  , _solve_temp(fluid_params.get<bool>("Build Temperature Equation"))
35  , _ir_degree(ir.cubature_degree)
36  , _tau_model_temp(TauModelTemp::TempOneDXNodal)
37 {
38  // Get tau model type for Navier-Stokes equations
39  {
40  const auto type_validator_mom = Teuchos::rcp(
41  new Teuchos::StringToIntegralParameterEntryValidator<TauModel>(
42  Teuchos::tuple<std::string>("Steady", "Transient", "NoSUPG"),
43  "Transient"));
44 
45  _tau_model = type_validator_mom->getIntegralValue(
46  closure_params.get<std::string>("Tau model for Navier-Stokes "
47  "equations"));
48 
49  // SUPG coefficient
50  if (closure_params.isType<double>("SUPG coefficient for continuity "
51  "equation"))
52  _alpha_cont = closure_params.get<double>(
53  "SUPG coefficient for continuity equation");
54  if (closure_params.isType<double>("SUPG coefficient for momentum "
55  "equations"))
56  _alpha_mom = closure_params.get<double>(
57  "SUPG coefficient for momentum equations");
58  }
59 
60  // Get tau model type for temperature equation
61  if (_solve_temp)
62  {
63  const auto type_validator_temp = Teuchos::rcp(
64  new Teuchos::StringToIntegralParameterEntryValidator<TauModelTemp>(
65  Teuchos::tuple<std::string>("TempMultiDSUPGSteady",
66  "TempMultiDSUPGTransient",
67  "TempNoSUPG",
68  "TempOneDXNodal",
69  "TempOneDXUpwind"),
70  "TempOneDXNodal"));
71  _tau_model_temp = type_validator_temp->getIntegralValue(
72  closure_params.get<std::string>("Tau model for temperature "
73  "equation"));
74 
75  // SUPG coefficient
76  if (closure_params.isType<double>("SUPG coefficient for temperature "
77  "equation"))
78  _alpha_ene = closure_params.get<double>(
79  "SUPG coefficient for temperature equation");
80  }
81 
82  // Evaluated fields
83  this->addEvaluatedField(_tau_supg_cont);
84  this->addEvaluatedField(_tau_supg_mom);
85  if (_solve_temp)
86  this->addEvaluatedField(_tau_supg_energy);
87 
88  // Dependent fields
89  Utils::addDependentVectorField(*this, ir.dl_scalar, _velocity, "velocity_");
90  this->addDependentField(_element_length);
91  this->addDependentField(_nu);
92  if (_solve_temp)
93  {
94  this->addDependentField(_rho);
95  this->addDependentField(_k);
96  this->addDependentField(_cp);
97  }
98 
99  this->setName("Incompressible Tau SUPG " + std::to_string(num_space_dim)
100  + "D");
101 }
102 
103 //---------------------------------------------------------------------------//
104 template<class EvalType, class Traits, int NumSpaceDim>
105 void IncompressibleTauSUPG<EvalType, Traits, NumSpaceDim>::postRegistrationSetup(
106  typename Traits::SetupData sd, PHX::FieldManager<Traits>&)
107 {
108  _ir_index = panzer::getIntegrationRuleIndex(_ir_degree, (*sd.worksets_)[0]);
109 }
110 
111 //---------------------------------------------------------------------------//
112 template<class EvalType, class Traits, int NumSpaceDim>
113 void IncompressibleTauSUPG<EvalType, Traits, NumSpaceDim>::evaluateFields(
114  typename Traits::EvalData workset)
115 {
116  // Get FEM basis order
117  _basis_order = workset.bases[_ir_index]->basis_layout->getBasis()->order();
118 
119  // Get time step
120  const double dt = workset.step_size;
121  _time_constant = 4.0 / (dt * dt);
122 
123  // Allocate space for thread-local temporaries in parallel for
124  size_t bytes;
125  if constexpr (Sacado::IsADType<scalar_type>::value)
126  {
127  const int fad_size = Kokkos::dimension_scalar(_tau_supg_mom.get_view());
128  bytes = scratch_view::shmem_size(
129  _tau_supg_mom.extent(1), NUM_TMPS, fad_size);
130  }
131  else
132  {
133  bytes = scratch_view::shmem_size(_tau_supg_mom.extent(1), NUM_TMPS);
134  }
135 
136  auto policy = panzer::HP::inst().teamPolicy<scalar_type, PHX::Device>(
137  workset.num_cells);
138  policy.set_scratch_size(0, Kokkos::PerTeam(bytes));
139  Kokkos::parallel_for(this->getName(), policy, *this);
140 }
141 
142 //---------------------------------------------------------------------------//
143 template<class EvalType, class Traits, int NumSpaceDim>
144 KOKKOS_INLINE_FUNCTION void
145 IncompressibleTauSUPG<EvalType, Traits, NumSpaceDim>::operator()(
146  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
147 {
148  const int cell = team.league_rank();
149  const int num_point = _tau_supg_mom.extent(1);
150  using Kokkos::cosh;
151  using Kokkos::sinh;
152  using Kokkos::sqrt;
153  const double tol = 1.0e-8;
154 
155  scratch_view tmp_field;
156  if constexpr (Sacado::IsADType<scalar_type>::value)
157  {
158  const int fad_size = Kokkos::dimension_scalar(_tau_supg_mom.get_view());
159  tmp_field
160  = scratch_view(team.team_shmem(), num_point, NUM_TMPS, fad_size);
161  }
162  else
163  {
164  tmp_field = scratch_view(team.team_shmem(), num_point, NUM_TMPS);
165  }
166 
167  Kokkos::parallel_for(
168  Kokkos::TeamThreadRange(team, 0, num_point), [&](const int point) {
169  // Initialize local variables
170  auto&& u2_over_h2 = tmp_field(point, U2_OVER_H2);
171  auto&& pe = tmp_field(point, PE);
172  auto&& d = tmp_field(point, D);
173  u2_over_h2 = 0.0;
174  double h2 = 0.0;
175 
176  // Compute terms `h^2` and '\frac{u^2}{h^2}`
177  for (int dim = 0; dim < num_space_dim; ++dim)
178  {
179  h2 += _element_length(cell, point, dim)
180  * _element_length(cell, point, dim);
181  u2_over_h2 += _velocity[dim](cell, point)
182  * _velocity[dim](cell, point)
183  / (_element_length(cell, point, dim)
184  * _element_length(cell, point, dim));
185  }
186 
187  // SUPG coefficient for momentum and continuity equations
188  if (_tau_model == TauModel::NoSUPG)
189  {
190  _tau_supg_cont(cell, point) = 0.0;
191  _tau_supg_mom(cell, point) = 0.0;
192  }
193  else if (_tau_model == TauModel::Steady)
194  {
195  _tau_supg_mom(cell, point)
196  = 1.0
197  / sqrt(4.0 * u2_over_h2
198  + 9.0 * 16.0 * _nu(cell, point) * _nu(cell, point)
199  / (h2 * h2));
200  }
201  else // Transient
202  {
203  _tau_supg_mom(cell, point)
204  = 1.0
205  / sqrt(_time_constant + 4.0 * u2_over_h2
206  + 9.0 * 16.0 * _nu(cell, point) * _nu(cell, point)
207  / (h2 * h2));
208  }
209  _tau_supg_cont(cell, point) = _tau_supg_mom(cell, point);
210  _tau_supg_cont(cell, point) *= _alpha_cont / _basis_order;
211  _tau_supg_mom(cell, point) *= _alpha_mom / _basis_order;
212 
213  // SUPG coefficient for temperature equation
214  if (_solve_temp)
215  {
216  if (_tau_model_temp == TauModelTemp::TempNoSUPG)
217  {
218  _tau_supg_energy(cell, point) = 0.0;
219  }
220  else if (_tau_model_temp == TauModelTemp::TempOneDXUpwind)
221  {
222  _tau_supg_energy(cell, point)
223  = _alpha_ene / _basis_order
224  * _element_length(cell, point, 0)
225  / (_velocity[0](cell, point) + tol);
226  }
227  else if (_tau_model_temp == TauModelTemp::TempOneDXNodal)
228  {
229  // Nodal value: exact solution in 1D - second-order
230  // accuracy
231  d = _k(cell, point)
232  / (_rho(cell, point) * _cp(cell, point));
233  pe = 0.5 * _velocity[0](cell, point)
234  * _element_length(cell, point, 0) / d;
235  _tau_supg_energy(cell, point) = cosh(pe) / sinh(pe)
236  - 1.0 / pe;
237  _tau_supg_energy(cell, point)
238  *= _alpha_ene / _basis_order
239  * _element_length(cell, point, 0)
240  / (_velocity[0](cell, point) + tol);
241  }
242  else if (_tau_model_temp == TauModelTemp::TempMultiDSUPGSteady)
243  {
244  d = _k(cell, point)
245  / (_rho(cell, point) * _cp(cell, point));
246  _tau_supg_energy(cell, point)
247  = _alpha_ene / _basis_order
248  / sqrt(4.0 * u2_over_h2
249  + 9.0 * 16.0 * d * d / (h2 * h2));
250  }
251  else if (_tau_model_temp
252  == TauModelTemp::TempMultiDSUPGTransient)
253  {
254  d = _k(cell, point)
255  / (_rho(cell, point) * _cp(cell, point));
256  _tau_supg_energy(cell, point)
257  = _alpha_ene / _basis_order
258  / sqrt(_time_constant + 4.0 * u2_over_h2
259  + 9.0 * 16.0 * d * d / (h2 * h2));
260  }
261  }
262  });
263 }
264 
265 //---------------------------------------------------------------------------//
266 
267 } // end namespace ClosureModel
268 } // end namespace VertexCFD
269 
270 #endif // end VERTEXCFD_CLOSURE_INCOMPRESSIBLETAUSUPG_IMPL_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20