ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
FieldAccessorsDoublePrecision.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2021-2023 The ESPResSo project
3 * Copyright (C) 2020 The waLBerla project
4 *
5 * This file is part of ESPResSo.
6 *
7 * ESPResSo is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * ESPResSo is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21// kernel generated with pystencils v1.2, lbmpy v1.2,
22// lbmpy_walberla/pystencils_walberla from waLBerla commit
23// 4d10e7f2358fc4a4f7e99195d0f67f0b759ecb6f
24
25/**
26 * @file
27 * Lattice field accessors.
28 * Adapted from the waLBerla source file
29 * https://i10git.cs.fau.de/walberla/walberla/-/blob/a16141524c58ab88386e2a0f8fdd7c63c5edd704/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
30 */
31
32#pragma once
33
34#include <core/DataTypes.h>
35#include <core/cell/Cell.h>
36#include <core/cell/CellInterval.h>
37#include <core/math/Matrix3.h>
38#include <core/math/Vector3.h>
39
40#include <field/GhostLayerField.h>
41#include <stencil/D3Q19.h>
42
43#include <array>
44#include <cassert>
45#include <tuple>
46#include <vector>
47
48#ifdef WALBERLA_CXX_COMPILER_IS_GNU
49#pragma GCC diagnostic push
50#pragma GCC diagnostic ignored "-Wunused-variable"
51#pragma GCC diagnostic ignored "-Wunused-parameter"
52#endif
53
54#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
55#pragma clang diagnostic push
56#pragma clang diagnostic ignored "-Wunused-variable"
57#pragma clang diagnostic ignored "-Wunused-parameter"
58#endif
59
60namespace walberla {
61namespace lbm {
62namespace accessor {
63
64namespace Population {
65inline std::array<double, 19u>
66get(GhostLayerField<double, uint_t{19u}> const *pdf_field, Cell const &cell) {
67 double const &xyz0 = pdf_field->get(cell, uint_t{0u});
68 std::array<double, 19u> pop;
69 pop[0u] = pdf_field->getF(&xyz0, uint_t{0u});
70 pop[1u] = pdf_field->getF(&xyz0, uint_t{1u});
71 pop[2u] = pdf_field->getF(&xyz0, uint_t{2u});
72 pop[3u] = pdf_field->getF(&xyz0, uint_t{3u});
73 pop[4u] = pdf_field->getF(&xyz0, uint_t{4u});
74 pop[5u] = pdf_field->getF(&xyz0, uint_t{5u});
75 pop[6u] = pdf_field->getF(&xyz0, uint_t{6u});
76 pop[7u] = pdf_field->getF(&xyz0, uint_t{7u});
77 pop[8u] = pdf_field->getF(&xyz0, uint_t{8u});
78 pop[9u] = pdf_field->getF(&xyz0, uint_t{9u});
79 pop[10u] = pdf_field->getF(&xyz0, uint_t{10u});
80 pop[11u] = pdf_field->getF(&xyz0, uint_t{11u});
81 pop[12u] = pdf_field->getF(&xyz0, uint_t{12u});
82 pop[13u] = pdf_field->getF(&xyz0, uint_t{13u});
83 pop[14u] = pdf_field->getF(&xyz0, uint_t{14u});
84 pop[15u] = pdf_field->getF(&xyz0, uint_t{15u});
85 pop[16u] = pdf_field->getF(&xyz0, uint_t{16u});
86 pop[17u] = pdf_field->getF(&xyz0, uint_t{17u});
87 pop[18u] = pdf_field->getF(&xyz0, uint_t{18u});
88 return pop;
89}
90
91inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
92 std::array<double, 19u> const &pop, Cell const &cell) {
93 double &xyz0 = pdf_field->get(cell, uint_t{0u});
94 pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
95 pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
96 pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
97 pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
98 pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
99 pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
100 pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
101 pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
102 pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
103 pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
104 pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
105 pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
106 pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
107 pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
108 pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
109 pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
110 pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
111 pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
112 pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
113}
114
115inline void initialize(GhostLayerField<double, uint_t{19u}> *pdf_field,
116 std::array<double, 19u> const &pop) {
117 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdf_field, {
118 double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
119 pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
120 pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
121 pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
122 pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
123 pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
124 pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
125 pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
126 pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
127 pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
128 pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
129 pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
130 pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
131 pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
132 pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
133 pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
134 pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
135 pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
136 pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
137 pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
138 });
139}
140
141inline std::vector<double>
142get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
143 CellInterval const &ci) {
144 std::vector<double> out;
145 out.reserve(ci.numCells() * uint_t(19u));
146 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
147 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
148 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
149 double const &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
150 out.emplace_back(pdf_field->getF(&xyz0, uint_t{0u}));
151 out.emplace_back(pdf_field->getF(&xyz0, uint_t{1u}));
152 out.emplace_back(pdf_field->getF(&xyz0, uint_t{2u}));
153 out.emplace_back(pdf_field->getF(&xyz0, uint_t{3u}));
154 out.emplace_back(pdf_field->getF(&xyz0, uint_t{4u}));
155 out.emplace_back(pdf_field->getF(&xyz0, uint_t{5u}));
156 out.emplace_back(pdf_field->getF(&xyz0, uint_t{6u}));
157 out.emplace_back(pdf_field->getF(&xyz0, uint_t{7u}));
158 out.emplace_back(pdf_field->getF(&xyz0, uint_t{8u}));
159 out.emplace_back(pdf_field->getF(&xyz0, uint_t{9u}));
160 out.emplace_back(pdf_field->getF(&xyz0, uint_t{10u}));
161 out.emplace_back(pdf_field->getF(&xyz0, uint_t{11u}));
162 out.emplace_back(pdf_field->getF(&xyz0, uint_t{12u}));
163 out.emplace_back(pdf_field->getF(&xyz0, uint_t{13u}));
164 out.emplace_back(pdf_field->getF(&xyz0, uint_t{14u}));
165 out.emplace_back(pdf_field->getF(&xyz0, uint_t{15u}));
166 out.emplace_back(pdf_field->getF(&xyz0, uint_t{16u}));
167 out.emplace_back(pdf_field->getF(&xyz0, uint_t{17u}));
168 out.emplace_back(pdf_field->getF(&xyz0, uint_t{18u}));
169 }
170 }
171 }
172 return out;
173}
174
175inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
176 std::vector<double> const &values, CellInterval const &ci) {
177 assert(uint_c(values.size()) == ci.numCells() * uint_t(19u));
178 auto values_ptr = values.data();
179 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
180 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
181 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
182 double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
183 pdf_field->getF(&xyz0, uint_t{0u}) = values_ptr[0u];
184 pdf_field->getF(&xyz0, uint_t{1u}) = values_ptr[1u];
185 pdf_field->getF(&xyz0, uint_t{2u}) = values_ptr[2u];
186 pdf_field->getF(&xyz0, uint_t{3u}) = values_ptr[3u];
187 pdf_field->getF(&xyz0, uint_t{4u}) = values_ptr[4u];
188 pdf_field->getF(&xyz0, uint_t{5u}) = values_ptr[5u];
189 pdf_field->getF(&xyz0, uint_t{6u}) = values_ptr[6u];
190 pdf_field->getF(&xyz0, uint_t{7u}) = values_ptr[7u];
191 pdf_field->getF(&xyz0, uint_t{8u}) = values_ptr[8u];
192 pdf_field->getF(&xyz0, uint_t{9u}) = values_ptr[9u];
193 pdf_field->getF(&xyz0, uint_t{10u}) = values_ptr[10u];
194 pdf_field->getF(&xyz0, uint_t{11u}) = values_ptr[11u];
195 pdf_field->getF(&xyz0, uint_t{12u}) = values_ptr[12u];
196 pdf_field->getF(&xyz0, uint_t{13u}) = values_ptr[13u];
197 pdf_field->getF(&xyz0, uint_t{14u}) = values_ptr[14u];
198 pdf_field->getF(&xyz0, uint_t{15u}) = values_ptr[15u];
199 pdf_field->getF(&xyz0, uint_t{16u}) = values_ptr[16u];
200 pdf_field->getF(&xyz0, uint_t{17u}) = values_ptr[17u];
201 pdf_field->getF(&xyz0, uint_t{18u}) = values_ptr[18u];
202 values_ptr += 19u;
203 }
204 }
205 }
206}
207} // namespace Population
208
209namespace Vector {
210inline Vector3<double> get(GhostLayerField<double, uint_t{3u}> const *vec_field,
211 Cell const &cell) {
212 const double &xyz0 = vec_field->get(cell, uint_t{0u});
213 Vector3<double> vec;
214 vec[0] = vec_field->getF(&xyz0, uint_t{0u});
215 vec[1] = vec_field->getF(&xyz0, uint_t{1u});
216 vec[2] = vec_field->getF(&xyz0, uint_t{2u});
217 return vec;
218}
219
220inline void set(GhostLayerField<double, uint_t{3u}> *vec_field,
221 Vector3<double> const &vec, Cell const &cell) {
222 double &xyz0 = vec_field->get(cell, uint_t{0u});
223 vec_field->getF(&xyz0, uint_t{0u}) = vec[0u];
224 vec_field->getF(&xyz0, uint_t{1u}) = vec[1u];
225 vec_field->getF(&xyz0, uint_t{2u}) = vec[2u];
226}
227
228inline void add(GhostLayerField<double, uint_t{3u}> *vec_field,
229 Vector3<double> const &vec, Cell const &cell) {
230 double &xyz0 = vec_field->get(cell, uint_t{0u});
231 vec_field->getF(&xyz0, uint_t{0u}) += vec[0u];
232 vec_field->getF(&xyz0, uint_t{1u}) += vec[1u];
233 vec_field->getF(&xyz0, uint_t{2u}) += vec[2u];
234}
235
236inline void initialize(GhostLayerField<double, uint_t{3u}> *vec_field,
237 Vector3<double> const &vec) {
238 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
239 double &xyz0 = vec_field->get(x, y, z, uint_t{0u});
240 vec_field->getF(&xyz0, uint_t{0u}) = vec[0u];
241 vec_field->getF(&xyz0, uint_t{1u}) = vec[1u];
242 vec_field->getF(&xyz0, uint_t{2u}) = vec[2u];
243 });
244}
245
246inline void add_to_all(GhostLayerField<double, uint_t{3u}> *vec_field,
247 Vector3<double> const &vec) {
248 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
249 double &xyz0 = vec_field->get(x, y, z, uint_t{0u});
250 vec_field->getF(&xyz0, uint_t{0u}) += vec[0u];
251 vec_field->getF(&xyz0, uint_t{1u}) += vec[1u];
252 vec_field->getF(&xyz0, uint_t{2u}) += vec[2u];
253 });
254}
255
256inline std::vector<double>
257get(GhostLayerField<double, uint_t{3u}> const *vec_field,
258 CellInterval const &ci) {
259 std::vector<double> out;
260 out.reserve(ci.numCells() * uint_t(3u));
261 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
262 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
263 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
264 const double &xyz0 = vec_field->get(x, y, z, uint_t{0u});
265 out.emplace_back(vec_field->getF(&xyz0, uint_t{0u}));
266 out.emplace_back(vec_field->getF(&xyz0, uint_t{1u}));
267 out.emplace_back(vec_field->getF(&xyz0, uint_t{2u}));
268 }
269 }
270 }
271 return out;
272}
273
274inline void set(GhostLayerField<double, uint_t{3u}> *vec_field,
275 std::vector<double> const &values, CellInterval const &ci) {
276 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
277 auto values_ptr = values.data();
278 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
279 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
280 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
281 double &xyz0 = vec_field->get(x, y, z, uint_t{0u});
282 vec_field->getF(&xyz0, uint_t{0u}) = values_ptr[0u];
283 vec_field->getF(&xyz0, uint_t{1u}) = values_ptr[1u];
284 vec_field->getF(&xyz0, uint_t{2u}) = values_ptr[2u];
285 values_ptr += 3u;
286 }
287 }
288 }
289}
290} // namespace Vector
291
292namespace EquilibriumDistribution {
293inline double get(stencil::Direction const direction,
294 Vector3<double> const &u = Vector3<double>(double(0.0)),
295 double rho = double(1.0)) {
296
297 using namespace stencil;
298 switch (direction) {
299 case C:
300 return rho * -0.33333333333333331 * (u[0] * u[0]) +
301 rho * -0.33333333333333331 * (u[1] * u[1]) +
302 rho * -0.33333333333333331 * (u[2] * u[2]) +
303 rho * 0.33333333333333331;
304 case N:
305 return rho * -0.16666666666666666 * (u[0] * u[0]) +
306 rho * -0.16666666666666666 * (u[2] * u[2]) +
307 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[1] +
308 rho * 0.16666666666666666 * (u[1] * u[1]);
309 case S:
310 return rho * -0.16666666666666666 * u[1] +
311 rho * -0.16666666666666666 * (u[0] * u[0]) +
312 rho * -0.16666666666666666 * (u[2] * u[2]) +
313 rho * 0.055555555555555552 +
314 rho * 0.16666666666666666 * (u[1] * u[1]);
315 case W:
316 return rho * -0.16666666666666666 * u[0] +
317 rho * -0.16666666666666666 * (u[1] * u[1]) +
318 rho * -0.16666666666666666 * (u[2] * u[2]) +
319 rho * 0.055555555555555552 +
320 rho * 0.16666666666666666 * (u[0] * u[0]);
321 case E:
322 return rho * -0.16666666666666666 * (u[1] * u[1]) +
323 rho * -0.16666666666666666 * (u[2] * u[2]) +
324 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[0] +
325 rho * 0.16666666666666666 * (u[0] * u[0]);
326 case T:
327 return rho * -0.16666666666666666 * (u[0] * u[0]) +
328 rho * -0.16666666666666666 * (u[1] * u[1]) +
329 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[2] +
330 rho * 0.16666666666666666 * (u[2] * u[2]);
331 case B:
332 return rho * -0.16666666666666666 * u[2] +
333 rho * -0.16666666666666666 * (u[0] * u[0]) +
334 rho * -0.16666666666666666 * (u[1] * u[1]) +
335 rho * 0.055555555555555552 +
336 rho * 0.16666666666666666 * (u[2] * u[2]);
337 case NW:
338 return rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[1] +
339 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
340 rho * 0.083333333333333329 * (u[0] * u[0]) +
341 rho * 0.083333333333333329 * (u[1] * u[1]);
342 case NE:
343 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
344 rho * 0.083333333333333329 * u[1] +
345 rho * 0.083333333333333329 * (u[0] * u[0]) +
346 rho * 0.083333333333333329 * (u[1] * u[1]) +
347 rho * 0.25 * u[0] * u[1];
348 case SW:
349 return rho * -0.083333333333333329 * u[0] +
350 rho * -0.083333333333333329 * u[1] + rho * 0.027777777777777776 +
351 rho * 0.083333333333333329 * (u[0] * u[0]) +
352 rho * 0.083333333333333329 * (u[1] * u[1]) +
353 rho * 0.25 * u[0] * u[1];
354 case SE:
355 return rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[0] * u[1] +
356 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
357 rho * 0.083333333333333329 * (u[0] * u[0]) +
358 rho * 0.083333333333333329 * (u[1] * u[1]);
359 case TN:
360 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
361 rho * 0.083333333333333329 * u[2] +
362 rho * 0.083333333333333329 * (u[1] * u[1]) +
363 rho * 0.083333333333333329 * (u[2] * u[2]) +
364 rho * 0.25 * u[1] * u[2];
365 case TS:
366 return rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[1] * u[2] +
367 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
368 rho * 0.083333333333333329 * (u[1] * u[1]) +
369 rho * 0.083333333333333329 * (u[2] * u[2]);
370 case TW:
371 return rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[2] +
372 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
373 rho * 0.083333333333333329 * (u[0] * u[0]) +
374 rho * 0.083333333333333329 * (u[2] * u[2]);
375 case TE:
376 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
377 rho * 0.083333333333333329 * u[2] +
378 rho * 0.083333333333333329 * (u[0] * u[0]) +
379 rho * 0.083333333333333329 * (u[2] * u[2]) +
380 rho * 0.25 * u[0] * u[2];
381 case BN:
382 return rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[1] * u[2] +
383 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
384 rho * 0.083333333333333329 * (u[1] * u[1]) +
385 rho * 0.083333333333333329 * (u[2] * u[2]);
386 case BS:
387 return rho * -0.083333333333333329 * u[1] +
388 rho * -0.083333333333333329 * u[2] + rho * 0.027777777777777776 +
389 rho * 0.083333333333333329 * (u[1] * u[1]) +
390 rho * 0.083333333333333329 * (u[2] * u[2]) +
391 rho * 0.25 * u[1] * u[2];
392 case BW:
393 return rho * -0.083333333333333329 * u[0] +
394 rho * -0.083333333333333329 * u[2] + rho * 0.027777777777777776 +
395 rho * 0.083333333333333329 * (u[0] * u[0]) +
396 rho * 0.083333333333333329 * (u[2] * u[2]) +
397 rho * 0.25 * u[0] * u[2];
398 case BE:
399 return rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[0] * u[2] +
400 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
401 rho * 0.083333333333333329 * (u[0] * u[0]) +
402 rho * 0.083333333333333329 * (u[2] * u[2]);
403 default:
404 WALBERLA_ABORT("Invalid Direction")
405 }
406}
407} // namespace EquilibriumDistribution
408
409namespace Equilibrium {
410inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
411 Vector3<double> const &u, double const rho, Cell const &cell) {
412
413 double &xyz0 = pdf_field->get(cell, uint_t{0u});
414 pdf_field->getF(&xyz0, uint_t{0u}) =
415 rho * -0.33333333333333331 * (u[0] * u[0]) +
416 rho * -0.33333333333333331 * (u[1] * u[1]) +
417 rho * -0.33333333333333331 * (u[2] * u[2]) + rho * 0.33333333333333331;
418 pdf_field->getF(&xyz0, uint_t{1u}) =
419 rho * -0.16666666666666666 * (u[0] * u[0]) +
420 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
421 rho * 0.16666666666666666 * u[1] +
422 rho * 0.16666666666666666 * (u[1] * u[1]);
423 pdf_field->getF(&xyz0, uint_t{2u}) =
424 rho * -0.16666666666666666 * u[1] +
425 rho * -0.16666666666666666 * (u[0] * u[0]) +
426 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
427 rho * 0.16666666666666666 * (u[1] * u[1]);
428 pdf_field->getF(&xyz0, uint_t{3u}) =
429 rho * -0.16666666666666666 * u[0] +
430 rho * -0.16666666666666666 * (u[1] * u[1]) +
431 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
432 rho * 0.16666666666666666 * (u[0] * u[0]);
433 pdf_field->getF(&xyz0, uint_t{4u}) =
434 rho * -0.16666666666666666 * (u[1] * u[1]) +
435 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
436 rho * 0.16666666666666666 * u[0] +
437 rho * 0.16666666666666666 * (u[0] * u[0]);
438 pdf_field->getF(&xyz0, uint_t{5u}) =
439 rho * -0.16666666666666666 * (u[0] * u[0]) +
440 rho * -0.16666666666666666 * (u[1] * u[1]) + rho * 0.055555555555555552 +
441 rho * 0.16666666666666666 * u[2] +
442 rho * 0.16666666666666666 * (u[2] * u[2]);
443 pdf_field->getF(&xyz0, uint_t{6u}) =
444 rho * -0.16666666666666666 * u[2] +
445 rho * -0.16666666666666666 * (u[0] * u[0]) +
446 rho * -0.16666666666666666 * (u[1] * u[1]) + rho * 0.055555555555555552 +
447 rho * 0.16666666666666666 * (u[2] * u[2]);
448 pdf_field->getF(&xyz0, uint_t{7u}) =
449 rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[1] +
450 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
451 rho * 0.083333333333333329 * (u[0] * u[0]) +
452 rho * 0.083333333333333329 * (u[1] * u[1]);
453 pdf_field->getF(&xyz0, uint_t{8u}) =
454 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
455 rho * 0.083333333333333329 * u[1] +
456 rho * 0.083333333333333329 * (u[0] * u[0]) +
457 rho * 0.083333333333333329 * (u[1] * u[1]) + rho * 0.25 * u[0] * u[1];
458 pdf_field->getF(&xyz0, uint_t{9u}) =
459 rho * -0.083333333333333329 * u[0] + rho * -0.083333333333333329 * u[1] +
460 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[0] * u[0]) +
461 rho * 0.083333333333333329 * (u[1] * u[1]) + rho * 0.25 * u[0] * u[1];
462 pdf_field->getF(&xyz0, uint_t{10u}) =
463 rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[0] * u[1] +
464 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
465 rho * 0.083333333333333329 * (u[0] * u[0]) +
466 rho * 0.083333333333333329 * (u[1] * u[1]);
467 pdf_field->getF(&xyz0, uint_t{11u}) =
468 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
469 rho * 0.083333333333333329 * u[2] +
470 rho * 0.083333333333333329 * (u[1] * u[1]) +
471 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[1] * u[2];
472 pdf_field->getF(&xyz0, uint_t{12u}) =
473 rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[1] * u[2] +
474 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
475 rho * 0.083333333333333329 * (u[1] * u[1]) +
476 rho * 0.083333333333333329 * (u[2] * u[2]);
477 pdf_field->getF(&xyz0, uint_t{13u}) =
478 rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[2] +
479 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
480 rho * 0.083333333333333329 * (u[0] * u[0]) +
481 rho * 0.083333333333333329 * (u[2] * u[2]);
482 pdf_field->getF(&xyz0, uint_t{14u}) =
483 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
484 rho * 0.083333333333333329 * u[2] +
485 rho * 0.083333333333333329 * (u[0] * u[0]) +
486 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[0] * u[2];
487 pdf_field->getF(&xyz0, uint_t{15u}) =
488 rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[1] * u[2] +
489 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
490 rho * 0.083333333333333329 * (u[1] * u[1]) +
491 rho * 0.083333333333333329 * (u[2] * u[2]);
492 pdf_field->getF(&xyz0, uint_t{16u}) =
493 rho * -0.083333333333333329 * u[1] + rho * -0.083333333333333329 * u[2] +
494 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[1] * u[1]) +
495 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[1] * u[2];
496 pdf_field->getF(&xyz0, uint_t{17u}) =
497 rho * -0.083333333333333329 * u[0] + rho * -0.083333333333333329 * u[2] +
498 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[0] * u[0]) +
499 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[0] * u[2];
500 pdf_field->getF(&xyz0, uint_t{18u}) =
501 rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[0] * u[2] +
502 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
503 rho * 0.083333333333333329 * (u[0] * u[0]) +
504 rho * 0.083333333333333329 * (u[2] * u[2]);
505}
506} // namespace Equilibrium
507
508namespace Density {
509inline double get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
510 Cell const &cell) {
511 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
512 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
513 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
514 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
515 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
516 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
517 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
518 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
519 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
520 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
521 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
522 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
523 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
524 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
525 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
526 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
527 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
528 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
529 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
530 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
531 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
532 const double vel1Term = f_1 + f_11 + f_15 + f_7;
533 const double vel2Term = f_12 + f_13 + f_5;
534 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
535 vel1Term + vel2Term;
536 return rho;
537}
538
539inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
540 double const rho_in, Cell const &cell) {
541 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
542 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
543 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
544 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
545 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
546 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
547 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
548 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
549 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
550 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
551 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
552 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
553 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
554 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
555 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
556 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
557 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
558 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
559 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
560 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
561 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
562 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
563 const double vel1Term = f_1 + f_11 + f_15 + f_7;
564 const double momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
565 const double vel2Term = f_12 + f_13 + f_5;
566 const double momdensity_2 =
567 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
568 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
569 vel1Term + vel2Term;
570
571 // calculate current velocity (before density change)
572 const double conversion = double(1) / rho;
573 Vector3<double> velocity;
574 velocity[0u] = momdensity_0 * conversion;
575 velocity[1u] = momdensity_1 * conversion;
576 velocity[2u] = momdensity_2 * conversion;
577
578 Equilibrium::set(pdf_field, velocity, rho_in, cell);
579}
580
581inline std::vector<double>
582get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
583 CellInterval const &ci) {
584 std::vector<double> out;
585 out.reserve(ci.numCells());
586 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
587 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
588 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
589 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
590 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
591 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
592 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
593 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
594 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
595 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
596 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
597 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
598 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
599 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
600 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
601 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
602 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
603 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
604 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
605 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
606 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
607 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
608 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
609 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
610 const double vel1Term = f_1 + f_11 + f_15 + f_7;
611 const double vel2Term = f_12 + f_13 + f_5;
612 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
613 vel0Term + vel1Term + vel2Term;
614 out.emplace_back(rho);
615 }
616 }
617 }
618 return out;
619}
620
621inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
622 std::vector<double> const &values, CellInterval const &ci) {
623 assert(uint_c(values.size()) == ci.numCells());
624 auto values_it = values.begin();
625 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
626 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
627 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
628 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
629 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
630 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
631 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
632 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
633 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
634 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
635 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
636 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
637 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
638 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
639 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
640 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
641 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
642 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
643 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
644 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
645 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
646 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
647 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
648 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
649 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
650 const double vel1Term = f_1 + f_11 + f_15 + f_7;
651 const double momdensity_1 =
652 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
653 const double vel2Term = f_12 + f_13 + f_5;
654 const double momdensity_2 =
655 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
656 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
657 vel0Term + vel1Term + vel2Term;
658
659 // calculate current velocity (before density change)
660 const double conversion = double(1) / rho;
661 Vector3<double> velocity;
662 velocity[0u] = momdensity_0 * conversion;
663 velocity[1u] = momdensity_1 * conversion;
664 velocity[2u] = momdensity_2 * conversion;
665
666 Equilibrium::set(pdf_field, velocity, *values_it, Cell{x, y, z});
667 ++values_it;
668 }
669 }
670 }
671}
672} // namespace Density
673
674namespace Velocity {
675inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
676 GhostLayerField<double, uint_t{3u}> const *force_field,
677 Vector3<double> const &u, Cell const &cell) {
678 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
679 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
680 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
681 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
682 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
683 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
684 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
685 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
686 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
687 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
688 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
689 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
690 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
691 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
692 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
693 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
694 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
695 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
696 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
697 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
698 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
699 const double vel1Term = f_1 + f_11 + f_15 + f_7;
700 const double vel2Term = f_12 + f_13 + f_5;
701 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
702 vel1Term + vel2Term;
703
704 const auto x = cell.x();
705 const auto y = cell.y();
706 const auto z = cell.z();
707 const double u_0 =
708 -force_field->get(x, y, z, 0) * 0.50000000000000000 / rho + u[0];
709 const double u_1 =
710 -force_field->get(x, y, z, 1) * 0.50000000000000000 / rho + u[1];
711 const double u_2 =
712 -force_field->get(x, y, z, 2) * 0.50000000000000000 / rho + u[2];
713
714 Equilibrium::set(pdf_field, Vector3<double>(u_0, u_1, u_2), rho, cell);
715}
716} // namespace Velocity
717
718namespace MomentumDensity {
719inline Vector3<double>
720reduce(GhostLayerField<double, uint_t{19u}> const *pdf_field,
721 GhostLayerField<double, uint_t{3u}> const *force_field) {
722 Vector3<double> momentumDensity(double{0});
723 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
724 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
725 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
726 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
727 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
728 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
729 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
730 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
731 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
732 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
733 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
734 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
735 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
736 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
737 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
738 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
739 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
740 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
741 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
742 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
743 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
744 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
745 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
746 const double vel1Term = f_1 + f_11 + f_15 + f_7;
747 const double momdensity_1 =
748 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
749 const double vel2Term = f_12 + f_13 + f_5;
750 const double momdensity_2 =
751 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
752 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
753 vel1Term + vel2Term;
754 const double md_0 =
755 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
756 const double md_1 =
757 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
758 const double md_2 =
759 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
760
761 momentumDensity[0u] += md_0;
762 momentumDensity[1u] += md_1;
763 momentumDensity[2u] += md_2;
764 });
765 return momentumDensity;
766}
767} // namespace MomentumDensity
768
769namespace PressureTensor {
770inline Matrix3<double>
771get(GhostLayerField<double, uint_t{19u}> const *pdf_field, Cell const &cell) {
772 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
773 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
774 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
775 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
776 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
777 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
778 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
779 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
780 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
781 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
782 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
783 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
784 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
785 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
786 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
787 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
788 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
789 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
790 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
791 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
792 const double p_0 =
793 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
794 const double p_1 = -f_10 - f_7 + f_8 + f_9;
795 const double p_2 = -f_13 + f_14 + f_17 - f_18;
796 const double p_3 = -f_10 - f_7 + f_8 + f_9;
797 const double p_4 =
798 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
799 const double p_5 = f_11 - f_12 - f_15 + f_16;
800 const double p_6 = -f_13 + f_14 + f_17 - f_18;
801 const double p_7 = f_11 - f_12 - f_15 + f_16;
802 const double p_8 =
803 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
804
805 Matrix3<double> pressureTensor;
806 pressureTensor[0u] = p_0;
807 pressureTensor[1u] = p_1;
808 pressureTensor[2u] = p_2;
809
810 pressureTensor[3u] = p_3;
811 pressureTensor[4u] = p_4;
812 pressureTensor[5u] = p_5;
813
814 pressureTensor[6u] = p_6;
815 pressureTensor[7u] = p_7;
816 pressureTensor[8u] = p_8;
817
818 return pressureTensor;
819}
820} // namespace PressureTensor
821
822} // namespace accessor
823} // namespace lbm
824} // namespace walberla
825
826#ifdef WALBERLA_CXX_COMPILER_IS_GNU
827#pragma GCC diagnostic pop
828#endif
829
830#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
831#pragma clang diagnostic pop
832#endif
N
float N[3]
Definition barnes_hut_gpu_cuda.cu:690
u
float u[3]
Definition barnes_hut_gpu_cuda.cu:690
Cell
Definition Cell.hpp:98
walberla::lbm::accessor::Density::set
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, double const rho_in, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:539
walberla::lbm::accessor::Density::get
double get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:509
walberla::lbm::accessor::EquilibriumDistribution::get
double get(stencil::Direction const direction, Vector3< double > const &u=Vector3< double >(double(0.0)), double rho=double(1.0))
Definition FieldAccessorsDoublePrecision.h:293
walberla::lbm::accessor::Equilibrium::set
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, Vector3< double > const &u, double const rho, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:410
walberla::lbm::accessor::MomentumDensity::reduce
Vector3< double > reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field)
Definition FieldAccessorsDoublePrecision.h:720
walberla::lbm::accessor::Population::get
std::array< double, 19u > get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:66
walberla::lbm::accessor::Population::set
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:91
walberla::lbm::accessor::Population::initialize
void initialize(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop)
Definition FieldAccessorsDoublePrecision.h:115
walberla::lbm::accessor::PressureTensor::get
Matrix3< double > get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:771
walberla::lbm::accessor::Vector::initialize
void initialize(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec)
Definition FieldAccessorsDoublePrecision.h:236
walberla::lbm::accessor::Vector::add
void add(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:228
walberla::lbm::accessor::Vector::set
void set(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:220
walberla::lbm::accessor::Vector::add_to_all
void add_to_all(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec)
Definition FieldAccessorsDoublePrecision.h:246
walberla::lbm::accessor::Vector::get
Vector3< double > get(GhostLayerField< double, uint_t{3u}> const *vec_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:210
walberla::lbm::accessor::Velocity::set
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field, Vector3< double > const &u, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:675
velocity
static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs)
Velocity of the virtual site.
Definition relative.cpp:64