VertexCFD  0.0-dev
VertexCFD_BoundaryState_MethodManufacturedSolution_impl.hpp
1 #ifndef VERTEXCFD_BOUNDARYSTATE_METHODMANUFACTUREDSOLUTION_IMPL_HPP
2 #define VERTEXCFD_BOUNDARYSTATE_METHODMANUFACTUREDSOLUTION_IMPL_HPP
3 
4 #include "utils/VertexCFD_Utils_Constants.hpp"
5 #include <utils/VertexCFD_Utils_VectorField.hpp>
6 
7 #include <Panzer_HierarchicParallelism.hpp>
8 #include <Panzer_Workset_Utilities.hpp>
9 
10 namespace VertexCFD
11 {
12 namespace BoundaryCondition
13 {
14 //---------------------------------------------------------------------------//
15 template<class EvalType, class Traits, int NumSpaceDim>
16 MethodManufacturedSolution<EvalType, Traits, NumSpaceDim>::MethodManufacturedSolution(
17  const panzer::IntegrationRule& ir)
18  : _boundary_lagrange_pressure("BOUNDARY_lagrange_pressure", ir.dl_scalar)
19  , _boundary_temperature("BOUNDARY_temperature", ir.dl_scalar)
20  , _boundary_grad_lagrange_pressure("BOUNDARY_GRAD_lagrange_pressure",
21  ir.dl_vector)
22  , _boundary_grad_temperature("BOUNDARY_GRAD_temperature", ir.dl_vector)
23  , _ir_degree(ir.cubature_degree)
24 {
25  // Add evaludated fields
26  this->addEvaluatedField(_boundary_lagrange_pressure);
27  Utils::addEvaluatedVectorField(
28  *this, ir.dl_scalar, _boundary_velocity, "BOUNDARY_velocity_");
29  this->addEvaluatedField(_boundary_temperature);
30 
31  this->addEvaluatedField(_boundary_grad_lagrange_pressure);
32  Utils::addEvaluatedVectorField(*this,
33  ir.dl_vector,
34  _boundary_grad_velocity,
35  "BOUNDARY_GRAD_velocity_");
36  this->addEvaluatedField(_boundary_grad_temperature);
37 
38  this->setName("Boundary State Method Manufactured Solution "
39  + std::to_string(num_space_dim) + "D");
40 
41  // Initialize coefficients for MMS
42  _phi_coeff[0] = 0.0125;
43  _phi_coeff[1] = 1.0;
44  _phi_coeff[2] = 0.25;
45  _phi_coeff[3] = 0.5;
46  _phi_coeff[4] = 0.125;
47  _phi_coeff[5] = 0.0;
48 
49  _vel_coeff[0][0] = 0.0125;
50  _vel_coeff[0][1] = 0.08;
51  _vel_coeff[0][2] = 0.125;
52  _vel_coeff[0][3] = 0.0;
53  _vel_coeff[0][4] = 0.125;
54  _vel_coeff[0][5] = 0.25;
55 
56  _vel_coeff[1][0] = 0.0375;
57  _vel_coeff[1][1] = 1.125;
58  _vel_coeff[1][2] = 0.25;
59  _vel_coeff[1][3] = 0.0;
60  _vel_coeff[1][4] = 0.375;
61  _vel_coeff[1][5] = 0.5;
62 
63  _T_coeff[0] = 0.0625;
64  _T_coeff[1] = 1.0;
65  _T_coeff[2] = 0.375;
66  _T_coeff[3] = 0.25;
67  _T_coeff[4] = 0.25;
68  _T_coeff[5] = 0.5;
69 
70  if (num_space_dim == 3)
71  {
72  _phi_coeff[6] = 0.375;
73  _phi_coeff[7] = 1.0;
74 
75  _vel_coeff[0][6] = 0.25;
76  _vel_coeff[0][7] = 0.0;
77 
78  _vel_coeff[1][6] = 0.25;
79  _vel_coeff[1][7] = 0.5;
80 
81  _vel_coeff[2][0] = 0.025;
82  _vel_coeff[2][1] = 0.0;
83  _vel_coeff[2][2] = 0.125;
84  _vel_coeff[2][3] = 1.0;
85  _vel_coeff[2][4] = 0.25;
86  _vel_coeff[2][5] = 0.25;
87  _vel_coeff[2][6] = 0.25;
88  _vel_coeff[2][7] = 0.25;
89 
90  _T_coeff[6] = 0.125;
91  _T_coeff[7] = 0.5;
92  }
93 }
94 
95 //---------------------------------------------------------------------------//
96 template<class EvalType, class Traits, int NumSpaceDim>
97 void MethodManufacturedSolution<EvalType, Traits, NumSpaceDim>::postRegistrationSetup(
98  typename Traits::SetupData sd, PHX::FieldManager<Traits>&)
99 {
100  _ir_index = panzer::getIntegrationRuleIndex(_ir_degree, (*sd.worksets_)[0]);
101 }
102 
103 //---------------------------------------------------------------------------//
104 template<class EvalType, class Traits, int NumSpaceDim>
105 void MethodManufacturedSolution<EvalType, Traits, NumSpaceDim>::evaluateFields(
106  typename Traits::EvalData workset)
107 {
108  _ip_coords = workset.int_rules[_ir_index]->ip_coordinates;
109  auto policy = panzer::HP::inst().teamPolicy<scalar_type, PHX::Device>(
110  workset.num_cells);
111  Kokkos::parallel_for(this->getName(), policy, *this);
112 }
113 
114 //---------------------------------------------------------------------------//
115 template<class EvalType, class Traits, int NumSpaceDim>
116 KOKKOS_INLINE_FUNCTION void
117 MethodManufacturedSolution<EvalType, Traits, NumSpaceDim>::operator()(
118  const Kokkos::TeamPolicy<PHX::exec_space>::member_type& team) const
119 {
120  const int cell = team.league_rank();
121  const int num_point = _boundary_lagrange_pressure.extent(1);
122  using Constants::pi;
123  using std::cos;
124  using std::sin;
125 
126  Kokkos::parallel_for(
127  Kokkos::TeamThreadRange(team, 0, num_point), [&](const int point) {
128  // function = B + A * sin(2*pi*f_x*(x-phi_x)) *
129  // sin(2*pi*f_y*(x-phi_y))
130  // * sin(2*pi*f_z*(z-phi_z)) for 3D
131 
132  auto set_function
133  = [=](const Kokkos::Array<double, num_coeff> coeff,
134  const Kokkos::Array<double, num_space_dim> coords) {
135  double return_val = coeff[0];
136  for (int i = 0; i < num_space_dim; ++i)
137  return_val
138  *= sin(2.0 * pi * coeff[2 * (i + 1)]
139  * (coords[i] - coeff[2 * (i + 1) + 1]));
140  return_val += coeff[1];
141  return return_val;
142  };
143 
144  // function = A * 2*pi*f_x* cos(2*pi*f_x*(x-phi_x)) *
145  // sin(2*pi*f_y*(y-phi_y))
146  // * sin(2*pi*f_z*(z-phi_z)) for 3D
147  auto set_gradX_function =
148  [=](const Kokkos::Array<double, num_coeff> coeff,
149  const Kokkos::Array<double, num_space_dim> coords) {
150  double return_val
151  = coeff[0] * 2.0 * pi * coeff[2]
152  * cos(2.0 * pi * coeff[2] * (coords[0] - coeff[3]))
153  * sin(2.0 * pi * coeff[4] * (coords[1] - coeff[5]));
154  if (num_space_dim == 3)
155  return_val *= sin(2.0 * pi * coeff[6]
156  * (coords[2] - coeff[7]));
157  return return_val;
158  };
159 
160  // function = A * sin(2*pi*f_x*(x-phi_x)) * 2*pi*f_y*
161  // cos(2*pi*f_y*(y-phi_y))
162  // * sin(2*pi*f_z*(z-phi_z)) for 3D
163  auto set_gradY_function =
164  [=](const Kokkos::Array<double, num_coeff> coeff,
165  const Kokkos::Array<double, num_space_dim> coords) {
166  double return_val
167  = coeff[0] * 2.0 * pi * coeff[4]
168  * sin(2.0 * pi * coeff[2] * (coords[0] - coeff[3]))
169  * cos(2.0 * pi * coeff[4] * (coords[1] - coeff[5]));
170  if (num_space_dim == 3)
171  return_val *= sin(2.0 * pi * coeff[6]
172  * (coords[2] - coeff[7]));
173  return return_val;
174  };
175 
176  // function = A * sin(2*pi*f_x*(x-phi_x)) * sin(2*pi*f_y*(y-phi_y))
177  // * 2*pi*f_Z * cos(2*pi*f_z*(z-phi_z)) for 3D
178  auto set_gradZ_function =
179  [=](const Kokkos::Array<double, num_coeff> coeff,
180  const Kokkos::Array<double, num_space_dim> coords) {
181  return coeff[0] * 2.0 * pi * coeff[6]
182  * sin(2.0 * pi * coeff[2] * (coords[0] - coeff[3]))
183  * sin(2.0 * pi * coeff[4] * (coords[1] - coeff[5]))
184  * cos(2.0 * pi * coeff[6] * (coords[2] - coeff[7]));
185  };
186 
187  Kokkos::Array<double, num_space_dim> x;
188  for (int dim = 0; dim < num_space_dim; ++dim)
189  x[dim] = _ip_coords(cell, point, dim);
190 
191  const double phi = set_function(_phi_coeff, x);
192  const double T = set_function(_T_coeff, x);
193 
194  _boundary_lagrange_pressure(cell, point) = phi;
195  Kokkos::Array<scalar_type, num_space_dim> vel;
196  for (int i = 0; i < num_space_dim; ++i)
197  {
198  _boundary_velocity[i](cell, point)
199  = set_function(_vel_coeff[i], x);
200  vel[i] = _boundary_velocity[i](cell, point);
201  }
202  _boundary_temperature(cell, point) = T;
203 
204  _boundary_grad_lagrange_pressure(cell, point, 0)
205  = set_gradX_function(_phi_coeff, x);
206  _boundary_grad_lagrange_pressure(cell, point, 1)
207  = set_gradY_function(_phi_coeff, x);
208 
209  _boundary_grad_velocity[0](cell, point, 0)
210  = set_gradX_function(_vel_coeff[0], x);
211  _boundary_grad_velocity[0](cell, point, 1)
212  = set_gradY_function(_vel_coeff[0], x);
213  _boundary_grad_velocity[1](cell, point, 0)
214  = set_gradX_function(_vel_coeff[1], x);
215  _boundary_grad_velocity[1](cell, point, 1)
216  = set_gradY_function(_vel_coeff[1], x);
217 
218  _boundary_grad_temperature(cell, point, 0)
219  = set_gradX_function(_T_coeff, x);
220  _boundary_grad_temperature(cell, point, 1)
221  = set_gradY_function(_T_coeff, x);
222 
223  if (num_space_dim == 3)
224  {
225  _boundary_grad_lagrange_pressure(cell, point, 2)
226  = set_gradZ_function(_phi_coeff, x);
227 
228  _boundary_grad_velocity[0](cell, point, 2)
229  = set_gradZ_function(_vel_coeff[0], x);
230  _boundary_grad_velocity[1](cell, point, 2)
231  = set_gradZ_function(_vel_coeff[1], x);
232  _boundary_grad_velocity[2](cell, point, 0)
233  = set_gradX_function(_vel_coeff[2], x);
234  _boundary_grad_velocity[2](cell, point, 1)
235  = set_gradY_function(_vel_coeff[2], x);
236  _boundary_grad_velocity[2](cell, point, 2)
237  = set_gradZ_function(_vel_coeff[2], x);
238 
239  _boundary_grad_temperature(cell, point, 2)
240  = set_gradZ_function(_T_coeff, x);
241  }
242  });
243 }
244 
245 //---------------------------------------------------------------------------//
246 
247 } // end namespace BoundaryCondition
248 } // end namespace VertexCFD
249 
250 #endif // VERTEXCFD_BOUNDARYSTATE_METHODMANUFACTUREDSOLUTION_IMPL_HPP
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20