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