VertexCFD  0.0-dev
VertexCFD_Closure_IncompressibleChienKEpsilonSource_impl.hpp
1 #ifndef VERTEXCFD_CLOSURE_INCOMPRESSIBLECHIENKEPSILONSOURCE_IMPL_HPP
2 #define VERTEXCFD_CLOSURE_INCOMPRESSIBLECHIENKEPSILONSOURCE_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 IncompressibleChienKEpsilonSource<EvalType, Traits, NumSpaceDim>::
18  IncompressibleChienKEpsilonSource(
19  const panzer::IntegrationRule& ir,
20  const Teuchos::RCP<panzer::GlobalData>& global_data,
21  const Teuchos::ParameterList& user_params)
22  : _nu("kinematic_viscosity", ir.dl_scalar)
23  , _nu_t("turbulent_eddy_viscosity", ir.dl_scalar)
24  , _turb_kinetic_energy("turb_kinetic_energy", ir.dl_scalar)
25  , _turb_dissipation_rate("turb_dissipation_rate", ir.dl_scalar)
26  , _distance("distance", ir.dl_scalar)
27  , _global_data(global_data)
28  , _C_nu(0.09)
29  , _C_1(1.35)
30  , _C_2(1.8)
31  , _f_one(1.0)
32  , _k_source("SOURCE_turb_kinetic_energy_equation", ir.dl_scalar)
33  , _k_prod("PRODUCTION_turb_kinetic_energy_equation", ir.dl_scalar)
34  , _k_dest("DESTRUCTION_turb_kinetic_energy_equation", ir.dl_scalar)
35  , _e_source("SOURCE_turb_dissipation_rate_equation", ir.dl_scalar)
36  , _e_prod("PRODUCTION_turb_dissipation_rate_equation", ir.dl_scalar)
37  , _e_dest("DESTRUCTION_turb_dissipation_rate_equation", ir.dl_scalar)
38 {
39  const Teuchos::ParameterList turb_list
40  = user_params.sublist("Turbulence Parameters");
41 
42  // Get boundary surface area for wall shear stress
43  _area = turb_list.get<double>("Boundary Surface Area");
44 
45  // Add dependent fields
46  this->addDependentField(_nu);
47  this->addDependentField(_nu_t);
48  this->addDependentField(_turb_kinetic_energy);
49  this->addDependentField(_turb_dissipation_rate);
50  this->addDependentField(_distance);
51 
52  Utils::addDependentVectorField(
53  *this, ir.dl_vector, _grad_velocity, "GRAD_velocity_");
54 
55  // Add evaluated fields
56  this->addEvaluatedField(_k_source);
57  this->addEvaluatedField(_k_prod);
58  this->addEvaluatedField(_k_dest);
59  this->addEvaluatedField(_e_source);
60  this->addEvaluatedField(_e_prod);
61  this->addEvaluatedField(_e_dest);
62 
63  // Closure model name
64  this->setName("Chien K-Epsilon Incompressible Source "
65  + std::to_string(num_space_dim) + "D");
66 }
67 
68 //---------------------------------------------------------------------------//
69 template<class EvalType, class Traits, int NumSpaceDim>
70 void IncompressibleChienKEpsilonSource<EvalType, Traits, NumSpaceDim>::evaluateFields(
71  typename Traits::EvalData workset)
72 {
73  const auto& pl = *_global_data->pl;
74  const std::string wall_shear_stress_name
75  = "Friction Velocity - friction_velocity";
76  _wall_shear_stress = pl.getValue<EvalType>(wall_shear_stress_name) / _area;
77 
78  auto policy = panzer::HP::inst().teamPolicy<scalar_type, PHX::Device>(
79  workset.num_cells);
80  Kokkos::parallel_for(this->getName(), policy, *this);
81 }
82 
83 //---------------------------------------------------------------------------//
84 template<class EvalType, class Traits, int NumSpaceDim>
85 KOKKOS_INLINE_FUNCTION void
86 IncompressibleChienKEpsilonSource<EvalType, Traits, NumSpaceDim>::operator()(
87  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
88 {
89  const int cell = team.league_rank();
90  const int num_point = _nu_t.extent(1);
91  const double max_tol = 1.0e-14;
92  using Kokkos::exp;
93  using Kokkos::pow;
94 
95  Kokkos::parallel_for(
96  Kokkos::TeamThreadRange(team, 0, num_point), [&](const int point) {
97  // Production term
98  _k_prod(cell, point) = 0.0;
99 
100  for (int i = 0; i < num_space_dim; ++i)
101  {
102  for (int j = 0; j < num_space_dim; ++j)
103  {
104  _k_prod(cell, point)
105  += pow(0.5
106  * (_grad_velocity[i](cell, point, j)
107  + _grad_velocity[j](cell, point, i)),
108  2.0);
109  }
110  }
111 
112  // Turbulent kinetic energy terms
113  _k_prod(cell, point) *= (2.0 * _nu_t(cell, point));
114 
115  _k_dest(cell, point)
116  = _turb_dissipation_rate(cell, point)
117  + 2 * _nu(cell, point) * _turb_kinetic_energy(cell, point)
118  / SmoothMath::max(
119  pow(_distance(cell, point), 2.0), max_tol, 0.0);
120 
121  _k_source(cell, point) = _k_prod(cell, point)
122  - _k_dest(cell, point);
123 
124  // Turbulent dissipation rate terms
125  _e_prod(cell, point)
126  = _k_prod(cell, point) * _f_one * _C_1
127  * _turb_dissipation_rate(cell, point)
128  / SmoothMath::max(
129  _turb_kinetic_energy(cell, point), max_tol, 0.0);
130 
131  if (_turb_dissipation_rate(cell, point) > 0.0)
132  {
133  _e_dest(cell, point)
134  = _C_2
135  * (1.0
136  - 0.22
137  * exp(-pow(
138  pow(_turb_kinetic_energy(cell, point),
139  2.0)
140  / SmoothMath::max(
141  6.0 * _nu(cell, point)
142  * _turb_dissipation_rate(
143  cell, point),
144  max_tol,
145  0.0),
146  2.0)))
147  * pow(_turb_dissipation_rate(cell, point), 2.0)
148  / SmoothMath::max(
149  _turb_kinetic_energy(cell, point), max_tol, 0.0)
150  + 2.0 * _nu(cell, point)
151  * _turb_dissipation_rate(cell, point)
152  / SmoothMath::max(
153  pow(_distance(cell, point), 2.0), max_tol, 0.0)
154  * exp(-0.5 * _wall_shear_stress
155  * _distance(cell, point) / _nu(cell, point));
156  }
157  else
158  {
159  _e_dest(cell, point) = 0.0;
160  }
161 
162  _e_source(cell, point) = _e_prod(cell, point)
163  - _e_dest(cell, point);
164  });
165 }
166 
167 //---------------------------------------------------------------------------//
168 
169 } // end namespace ClosureModel
170 } // end namespace VertexCFD
171 
172 #endif // end
173  // VERTEXCFD_CLOSURE_INCOMPRESSIBLECHIENKEPSILONSOURCE_IMPL_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20