VertexCFD  0.0-dev
VertexCFD_Closure_IncompressibleKTauSource_impl.hpp
1 #ifndef VERTEXCFD_CLOSURE_INCOMPRESSIBLEKTAUSOURCE_IMPL_HPP
2 #define VERTEXCFD_CLOSURE_INCOMPRESSIBLEKTAUSOURCE_IMPL_HPP
3 
4 #include "utils/VertexCFD_Utils_SmoothMath.hpp"
5 #include "utils/VertexCFD_Utils_VectorField.hpp"
6 
7 #include <Panzer_HierarchicParallelism.hpp>
8 
9 #include <math.h>
10 
11 namespace VertexCFD
12 {
13 namespace ClosureModel
14 {
15 //---------------------------------------------------------------------------//
16 template<class EvalType, class Traits, int NumSpaceDim>
17 IncompressibleKTauSource<EvalType, Traits, NumSpaceDim>::IncompressibleKTauSource(
18  const panzer::IntegrationRule& ir,
19  const Teuchos::ParameterList& user_params)
20  : _nu("kinematic_viscosity", ir.dl_scalar)
21  , _nu_t("turbulent_eddy_viscosity", ir.dl_scalar)
22  , _turb_kinetic_energy("turb_kinetic_energy", ir.dl_scalar)
23  , _turb_specific_dissipation_rate("turb_specific_dissipation_rate",
24  ir.dl_scalar)
25  , _grad_turb_kinetic_energy("GRAD_turb_kinetic_energy", ir.dl_vector)
26  , _grad_turb_specific_dissipation_rate(
27  "GRAD_turb_specific_dissipation_rate", ir.dl_vector)
28  , _beta_star(0.09)
29  , _gamma(0.52)
30  , _beta_0(0.0708)
31  , _sigma_d(0.125)
32  , _sigma_t(0.5)
33  , _limit_production(true)
34  , _limit_destruction(true)
35  , _k_source("SOURCE_turb_kinetic_energy_equation", ir.dl_scalar)
36  , _k_prod("PRODUCTION_turb_kinetic_energy_equation", ir.dl_scalar)
37  , _k_dest("DESTRUCTION_turb_kinetic_energy_equation", ir.dl_scalar)
38  , _t_source("SOURCE_turb_specific_dissipation_rate_equation", ir.dl_scalar)
39  , _t_prod("PRODUCTION_turb_specific_dissipation_rate_equation",
40  ir.dl_scalar)
41  , _t_dest_nut(
42  "DESTRUCTION_eddy_viscosity_turb_specific_dissipation_rate_equation",
43  ir.dl_scalar)
44  , _t_dest_nu(
45  "DESTRUCTION_viscosity_turb_specific_dissipation_rate_equation",
46  ir.dl_scalar)
47  , _t_cross("CROSS_DIFFUSION_turb_specific_dissipation_rate_equation",
48  ir.dl_scalar)
49 {
50  // Check for user-defined coefficients or parameters
51  if (user_params.isSublist("Turbulence Parameters"))
52  {
53  Teuchos::ParameterList turb_list
54  = user_params.sublist("Turbulence Parameters");
55 
56  if (turb_list.isType<double>("beta_star"))
57  {
58  _beta_star = turb_list.get<double>("beta_star");
59  }
60 
61  if (turb_list.isType<double>("gamma"))
62  {
63  _gamma = turb_list.get<double>("gamma");
64  }
65 
66  if (turb_list.isType<double>("beta_0"))
67  {
68  _beta_0 = turb_list.get<double>("beta_0");
69  }
70 
71  if (turb_list.isType<double>("sigma_d"))
72  {
73  _sigma_d = turb_list.get<double>("sigma_d");
74  }
75 
76  if (turb_list.isType<double>("sigma_t"))
77  {
78  _sigma_t = turb_list.get<double>("sigma_t");
79  }
80 
81  if (turb_list.isType<bool>("Limit Production Term"))
82  {
83  _limit_production = turb_list.get<bool>("Limit Production Term");
84  }
85 
86  if (turb_list.isType<bool>("Limit Destruction Term"))
87  {
88  _limit_destruction = turb_list.get<bool>("Limit Destruction Term");
89  }
90  }
91 
92  // Add dependent fields
93  this->addDependentField(_nu);
94  this->addDependentField(_nu_t);
95  this->addDependentField(_turb_kinetic_energy);
96  this->addDependentField(_turb_specific_dissipation_rate);
97  this->addDependentField(_grad_turb_kinetic_energy);
98  this->addDependentField(_grad_turb_specific_dissipation_rate);
99 
100  Utils::addDependentVectorField(
101  *this, ir.dl_vector, _grad_velocity, "GRAD_velocity_");
102 
103  // Add evaluated fields
104  this->addEvaluatedField(_k_source);
105  this->addEvaluatedField(_k_prod);
106  this->addEvaluatedField(_k_dest);
107  this->addEvaluatedField(_t_source);
108  this->addEvaluatedField(_t_prod);
109  this->addEvaluatedField(_t_dest_nut);
110  this->addEvaluatedField(_t_dest_nu);
111  this->addEvaluatedField(_t_cross);
112 
113  // Closure model name
114  this->setName("K-Tau Incompressible Source "
115  + std::to_string(num_space_dim) + "D");
116 }
117 
118 //---------------------------------------------------------------------------//
119 template<class EvalType, class Traits, int NumSpaceDim>
120 void IncompressibleKTauSource<EvalType, Traits, NumSpaceDim>::evaluateFields(
121  typename Traits::EvalData workset)
122 {
123  auto policy = panzer::HP::inst().teamPolicy<scalar_type, PHX::Device>(
124  workset.num_cells);
125  Kokkos::parallel_for(this->getName(), policy, *this);
126 }
127 
128 //---------------------------------------------------------------------------//
129 template<class EvalType, class Traits, int NumSpaceDim>
130 KOKKOS_INLINE_FUNCTION void
131 IncompressibleKTauSource<EvalType, Traits, NumSpaceDim>::operator()(
132  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
133 {
134  const int cell = team.league_rank();
135  const int num_point = _nu_t.extent(1);
136  const double max_tol = 1.0e-10;
137  using std::pow;
138 
139  Kokkos::parallel_for(
140  Kokkos::TeamThreadRange(team, 0, num_point), [&](const int point) {
141  auto&& Sij_Sij = _t_prod(cell, point);
142  Sij_Sij = 0.0;
143  for (int i = 0; i < num_space_dim; ++i)
144  {
145  Sij_Sij += pow(_grad_velocity[i](cell, point, i), 2.0);
146  for (int j = i + 1; j < num_space_dim; ++j)
147  {
148  Sij_Sij += 0.5
149  * pow(_grad_velocity[i](cell, point, j)
150  + _grad_velocity[j](cell, point, i),
151  2.0);
152  }
153  }
154 
155  // Turbulent kinetic energy terms
156  _k_prod(cell, point) = 2.0 * _nu_t(cell, point) * Sij_Sij;
157 
158  _k_dest(cell, point)
159  = _beta_star * _turb_kinetic_energy(cell, point)
160  / SmoothMath::max(
161  _turb_specific_dissipation_rate(cell, point),
162  max_tol,
163  0.0);
164 
165  if (_limit_production)
166  {
167  // The limiter term is 10 times the k destruction term
168  _k_prod(cell, point) = SmoothMath::min(
169  _k_prod(cell, point), 10.0 * _k_dest(cell, point), 0.0);
170  }
171 
172  _k_source(cell, point)
173  = (_k_prod(cell, point) - _k_dest(cell, point));
174 
175  // Turbulent dissipation rate terms
176  _t_prod(cell, point)
177  = _beta_0
178  - _gamma
179  * pow(_turb_specific_dissipation_rate(cell, point), 2.0)
180  * 2.0 * Sij_Sij;
181 
182  // Compute the destruction and cross diffusion term
183  auto&& dkdxj_dwdxj = _t_cross(cell, point);
184  dkdxj_dwdxj = 0.0;
185  auto&& limit_source = _t_dest_nu(cell, point);
186  limit_source = 0.0;
187  for (int j = 0; j < num_space_dim; ++j)
188  {
189  limit_source += pow(
190  _grad_turb_specific_dissipation_rate(cell, point, j), 2.0);
191 
192  dkdxj_dwdxj
193  += _grad_turb_kinetic_energy(cell, point, j)
194  * _grad_turb_specific_dissipation_rate(cell, point, j);
195  }
196 
197  _t_cross(cell, point)
198  = _sigma_d * _turb_specific_dissipation_rate(cell, point)
199  * SmoothMath::min(dkdxj_dwdxj, 0.0, 0.0);
200 
201  _t_dest_nut(cell, point) = 2.0 * _turb_kinetic_energy(cell, point)
202  * _sigma_t * limit_source;
203  _t_dest_nu(cell, point)
204  = 2.0 * _nu(cell, point) * limit_source
205  * SmoothMath::max(
206  1
207  / (SmoothMath::max(
208  _turb_specific_dissipation_rate(cell, point),
209  max_tol,
210  0.0)),
211  max_tol,
212  0.0);
213  if (_limit_destruction
214  && _t_dest_nu(cell, point) >= _t_dest_nut(cell, point))
215  {
216  _t_source(cell, point)
217  = _t_prod(cell, point)
218  - (SmoothMath::min(
219  4.0 / 3.0 * _beta_0,
220  _t_dest_nut(cell, point) + _t_dest_nu(cell, point),
221  0.0))
222  + _t_cross(cell, point);
223  }
224  else if (_limit_destruction
225  && _t_dest_nut(cell, point) > _t_dest_nu(cell, point))
226  {
227  _t_source(cell, point)
228  = _t_prod(cell, point)
229  - (SmoothMath::min(
230  4.0 / 3.0 * _beta_0, _t_dest_nu(cell, point), 0.0)
231  + _t_dest_nut(cell, point))
232  + _t_cross(cell, point);
233  }
234  else
235  {
236  _t_source(cell, point)
237  = _t_prod(cell, point)
238  - (_t_dest_nut(cell, point) + _t_dest_nu(cell, point))
239  + _t_cross(cell, point);
240  }
241  });
242 }
243 
244 //---------------------------------------------------------------------------//
245 
246 } // end namespace ClosureModel
247 } // end namespace VertexCFD
248 
249 #endif // end
250  // VERTEXCFD_CLOSURE_INCOMPRESSIBLEKTAUSOURCE_IMPL_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20