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