ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
FieldAccessorsSinglePrecision.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<float, 19u>
66get(GhostLayerField<float, uint_t{19u}> const *pdf_field, Cell const &cell) {
67 float const &xyz0 = pdf_field->get(cell, uint_t{0u});
68 std::array<float, 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<float, uint_t{19u}> *pdf_field,
92 std::array<float, 19u> const &pop, Cell const &cell) {
93 float &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<float, uint_t{19u}> *pdf_field,
116 std::array<float, 19u> const &pop) {
117 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdf_field, {
118 float &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<float>
142get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
143 CellInterval const &ci) {
144 std::vector<float> 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 float 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<float, uint_t{19u}> *pdf_field,
176 std::vector<float> 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 float &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<float> get(GhostLayerField<float, uint_t{3u}> const *vec_field,
211 Cell const &cell) {
212 const float &xyz0 = vec_field->get(cell, uint_t{0u});
213 Vector3<float> 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<float, uint_t{3u}> *vec_field,
221 Vector3<float> const &vec, Cell const &cell) {
222 float &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<float, uint_t{3u}> *vec_field,
229 Vector3<float> const &vec, Cell const &cell) {
230 float &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<float, uint_t{3u}> *vec_field,
237 Vector3<float> const &vec) {
238 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
239 float &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<float, uint_t{3u}> *vec_field,
247 Vector3<float> const &vec) {
248 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
249 float &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<float>
257get(GhostLayerField<float, uint_t{3u}> const *vec_field,
258 CellInterval const &ci) {
259 std::vector<float> 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 float &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<float, uint_t{3u}> *vec_field,
275 std::vector<float> 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 float &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 float get(stencil::Direction const direction,
294 Vector3<float> const &u = Vector3<float>(float(0.0)),
295 float rho = float(1.0)) {
296
297 using namespace stencil;
298 switch (direction) {
299 case C:
300 return rho * -0.33333333333333331f * (u[0] * u[0]) +
301 rho * -0.33333333333333331f * (u[1] * u[1]) +
302 rho * -0.33333333333333331f * (u[2] * u[2]) +
303 rho * 0.33333333333333331f;
304 case N:
305 return rho * -0.16666666666666666f * (u[0] * u[0]) +
306 rho * -0.16666666666666666f * (u[2] * u[2]) +
307 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[1] +
308 rho * 0.16666666666666666f * (u[1] * u[1]);
309 case S:
310 return rho * -0.16666666666666666f * u[1] +
311 rho * -0.16666666666666666f * (u[0] * u[0]) +
312 rho * -0.16666666666666666f * (u[2] * u[2]) +
313 rho * 0.055555555555555552f +
314 rho * 0.16666666666666666f * (u[1] * u[1]);
315 case W:
316 return rho * -0.16666666666666666f * u[0] +
317 rho * -0.16666666666666666f * (u[1] * u[1]) +
318 rho * -0.16666666666666666f * (u[2] * u[2]) +
319 rho * 0.055555555555555552f +
320 rho * 0.16666666666666666f * (u[0] * u[0]);
321 case E:
322 return rho * -0.16666666666666666f * (u[1] * u[1]) +
323 rho * -0.16666666666666666f * (u[2] * u[2]) +
324 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[0] +
325 rho * 0.16666666666666666f * (u[0] * u[0]);
326 case T:
327 return rho * -0.16666666666666666f * (u[0] * u[0]) +
328 rho * -0.16666666666666666f * (u[1] * u[1]) +
329 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[2] +
330 rho * 0.16666666666666666f * (u[2] * u[2]);
331 case B:
332 return rho * -0.16666666666666666f * u[2] +
333 rho * -0.16666666666666666f * (u[0] * u[0]) +
334 rho * -0.16666666666666666f * (u[1] * u[1]) +
335 rho * 0.055555555555555552f +
336 rho * 0.16666666666666666f * (u[2] * u[2]);
337 case NW:
338 return rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[1] +
339 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
340 rho * 0.083333333333333329f * (u[0] * u[0]) +
341 rho * 0.083333333333333329f * (u[1] * u[1]);
342 case NE:
343 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
344 rho * 0.083333333333333329f * u[1] +
345 rho * 0.083333333333333329f * (u[0] * u[0]) +
346 rho * 0.083333333333333329f * (u[1] * u[1]) +
347 rho * 0.25f * u[0] * u[1];
348 case SW:
349 return rho * -0.083333333333333329f * u[0] +
350 rho * -0.083333333333333329f * u[1] + rho * 0.027777777777777776f +
351 rho * 0.083333333333333329f * (u[0] * u[0]) +
352 rho * 0.083333333333333329f * (u[1] * u[1]) +
353 rho * 0.25f * u[0] * u[1];
354 case SE:
355 return rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[0] * u[1] +
356 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
357 rho * 0.083333333333333329f * (u[0] * u[0]) +
358 rho * 0.083333333333333329f * (u[1] * u[1]);
359 case TN:
360 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
361 rho * 0.083333333333333329f * u[2] +
362 rho * 0.083333333333333329f * (u[1] * u[1]) +
363 rho * 0.083333333333333329f * (u[2] * u[2]) +
364 rho * 0.25f * u[1] * u[2];
365 case TS:
366 return rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[1] * u[2] +
367 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
368 rho * 0.083333333333333329f * (u[1] * u[1]) +
369 rho * 0.083333333333333329f * (u[2] * u[2]);
370 case TW:
371 return rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[2] +
372 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
373 rho * 0.083333333333333329f * (u[0] * u[0]) +
374 rho * 0.083333333333333329f * (u[2] * u[2]);
375 case TE:
376 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
377 rho * 0.083333333333333329f * u[2] +
378 rho * 0.083333333333333329f * (u[0] * u[0]) +
379 rho * 0.083333333333333329f * (u[2] * u[2]) +
380 rho * 0.25f * u[0] * u[2];
381 case BN:
382 return rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[1] * u[2] +
383 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
384 rho * 0.083333333333333329f * (u[1] * u[1]) +
385 rho * 0.083333333333333329f * (u[2] * u[2]);
386 case BS:
387 return rho * -0.083333333333333329f * u[1] +
388 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
389 rho * 0.083333333333333329f * (u[1] * u[1]) +
390 rho * 0.083333333333333329f * (u[2] * u[2]) +
391 rho * 0.25f * u[1] * u[2];
392 case BW:
393 return rho * -0.083333333333333329f * u[0] +
394 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
395 rho * 0.083333333333333329f * (u[0] * u[0]) +
396 rho * 0.083333333333333329f * (u[2] * u[2]) +
397 rho * 0.25f * u[0] * u[2];
398 case BE:
399 return rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[0] * u[2] +
400 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
401 rho * 0.083333333333333329f * (u[0] * u[0]) +
402 rho * 0.083333333333333329f * (u[2] * u[2]);
403 default:
404 WALBERLA_ABORT("Invalid Direction")
405 }
406}
407} // namespace EquilibriumDistribution
408
409namespace Equilibrium {
410inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
411 Vector3<float> const &u, float const rho, Cell const &cell) {
412
413 float &xyz0 = pdf_field->get(cell, uint_t{0u});
414 pdf_field->getF(&xyz0, uint_t{0u}) =
415 rho * -0.33333333333333331f * (u[0] * u[0]) +
416 rho * -0.33333333333333331f * (u[1] * u[1]) +
417 rho * -0.33333333333333331f * (u[2] * u[2]) + rho * 0.33333333333333331f;
418 pdf_field->getF(&xyz0, uint_t{1u}) =
419 rho * -0.16666666666666666f * (u[0] * u[0]) +
420 rho * -0.16666666666666666f * (u[2] * u[2]) +
421 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[1] +
422 rho * 0.16666666666666666f * (u[1] * u[1]);
423 pdf_field->getF(&xyz0, uint_t{2u}) =
424 rho * -0.16666666666666666f * u[1] +
425 rho * -0.16666666666666666f * (u[0] * u[0]) +
426 rho * -0.16666666666666666f * (u[2] * u[2]) +
427 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[1] * u[1]);
428 pdf_field->getF(&xyz0, uint_t{3u}) =
429 rho * -0.16666666666666666f * u[0] +
430 rho * -0.16666666666666666f * (u[1] * u[1]) +
431 rho * -0.16666666666666666f * (u[2] * u[2]) +
432 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[0] * u[0]);
433 pdf_field->getF(&xyz0, uint_t{4u}) =
434 rho * -0.16666666666666666f * (u[1] * u[1]) +
435 rho * -0.16666666666666666f * (u[2] * u[2]) +
436 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[0] +
437 rho * 0.16666666666666666f * (u[0] * u[0]);
438 pdf_field->getF(&xyz0, uint_t{5u}) =
439 rho * -0.16666666666666666f * (u[0] * u[0]) +
440 rho * -0.16666666666666666f * (u[1] * u[1]) +
441 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[2] +
442 rho * 0.16666666666666666f * (u[2] * u[2]);
443 pdf_field->getF(&xyz0, uint_t{6u}) =
444 rho * -0.16666666666666666f * u[2] +
445 rho * -0.16666666666666666f * (u[0] * u[0]) +
446 rho * -0.16666666666666666f * (u[1] * u[1]) +
447 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[2] * u[2]);
448 pdf_field->getF(&xyz0, uint_t{7u}) =
449 rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[1] +
450 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
451 rho * 0.083333333333333329f * (u[0] * u[0]) +
452 rho * 0.083333333333333329f * (u[1] * u[1]);
453 pdf_field->getF(&xyz0, uint_t{8u}) =
454 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
455 rho * 0.083333333333333329f * u[1] +
456 rho * 0.083333333333333329f * (u[0] * u[0]) +
457 rho * 0.083333333333333329f * (u[1] * u[1]) + rho * 0.25f * u[0] * u[1];
458 pdf_field->getF(&xyz0, uint_t{9u}) =
459 rho * -0.083333333333333329f * u[0] +
460 rho * -0.083333333333333329f * u[1] + rho * 0.027777777777777776f +
461 rho * 0.083333333333333329f * (u[0] * u[0]) +
462 rho * 0.083333333333333329f * (u[1] * u[1]) + rho * 0.25f * u[0] * u[1];
463 pdf_field->getF(&xyz0, uint_t{10u}) =
464 rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[0] * u[1] +
465 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
466 rho * 0.083333333333333329f * (u[0] * u[0]) +
467 rho * 0.083333333333333329f * (u[1] * u[1]);
468 pdf_field->getF(&xyz0, uint_t{11u}) =
469 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
470 rho * 0.083333333333333329f * u[2] +
471 rho * 0.083333333333333329f * (u[1] * u[1]) +
472 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[1] * u[2];
473 pdf_field->getF(&xyz0, uint_t{12u}) =
474 rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[1] * u[2] +
475 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
476 rho * 0.083333333333333329f * (u[1] * u[1]) +
477 rho * 0.083333333333333329f * (u[2] * u[2]);
478 pdf_field->getF(&xyz0, uint_t{13u}) =
479 rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[2] +
480 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
481 rho * 0.083333333333333329f * (u[0] * u[0]) +
482 rho * 0.083333333333333329f * (u[2] * u[2]);
483 pdf_field->getF(&xyz0, uint_t{14u}) =
484 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
485 rho * 0.083333333333333329f * u[2] +
486 rho * 0.083333333333333329f * (u[0] * u[0]) +
487 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[0] * u[2];
488 pdf_field->getF(&xyz0, uint_t{15u}) =
489 rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[1] * u[2] +
490 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
491 rho * 0.083333333333333329f * (u[1] * u[1]) +
492 rho * 0.083333333333333329f * (u[2] * u[2]);
493 pdf_field->getF(&xyz0, uint_t{16u}) =
494 rho * -0.083333333333333329f * u[1] +
495 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
496 rho * 0.083333333333333329f * (u[1] * u[1]) +
497 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[1] * u[2];
498 pdf_field->getF(&xyz0, uint_t{17u}) =
499 rho * -0.083333333333333329f * u[0] +
500 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
501 rho * 0.083333333333333329f * (u[0] * u[0]) +
502 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[0] * u[2];
503 pdf_field->getF(&xyz0, uint_t{18u}) =
504 rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[0] * u[2] +
505 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
506 rho * 0.083333333333333329f * (u[0] * u[0]) +
507 rho * 0.083333333333333329f * (u[2] * u[2]);
508}
509} // namespace Equilibrium
510
511namespace Density {
512inline float get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
513 Cell const &cell) {
514 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
515 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
516 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
517 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
518 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
519 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
520 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
521 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
522 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
523 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
524 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
525 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
526 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
527 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
528 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
529 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
530 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
531 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
532 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
533 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
534 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
535 const float vel1Term = f_1 + f_11 + f_15 + f_7;
536 const float vel2Term = f_12 + f_13 + f_5;
537 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
538 vel1Term + vel2Term;
539 return rho;
540}
541
542inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
543 float const rho_in, Cell const &cell) {
544 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
545 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
546 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
547 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
548 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
549 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
550 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
551 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
552 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
553 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
554 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
555 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
556 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
557 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
558 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
559 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
560 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
561 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
562 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
563 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
564 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
565 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
566 const float vel1Term = f_1 + f_11 + f_15 + f_7;
567 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
568 const float vel2Term = f_12 + f_13 + f_5;
569 const float momdensity_2 =
570 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
571 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
572 vel1Term + vel2Term;
573
574 // calculate current velocity (before density change)
575 const float conversion = float(1) / rho;
576 Vector3<float> velocity;
577 velocity[0u] = momdensity_0 * conversion;
578 velocity[1u] = momdensity_1 * conversion;
579 velocity[2u] = momdensity_2 * conversion;
580
581 Equilibrium::set(pdf_field, velocity, rho_in, cell);
582}
583
584inline std::vector<float>
585get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
586 CellInterval const &ci) {
587 std::vector<float> out;
588 out.reserve(ci.numCells());
589 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
590 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
591 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
592 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
593 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
594 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
595 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
596 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
597 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
598 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
599 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
600 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
601 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
602 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
603 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
604 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
605 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
606 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
607 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
608 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
609 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
610 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
611 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
612 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
613 const float vel1Term = f_1 + f_11 + f_15 + f_7;
614 const float vel2Term = f_12 + f_13 + f_5;
615 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
616 vel1Term + vel2Term;
617 out.emplace_back(rho);
618 }
619 }
620 }
621 return out;
622}
623
624inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
625 std::vector<float> const &values, CellInterval const &ci) {
626 assert(uint_c(values.size()) == ci.numCells());
627 auto values_it = values.begin();
628 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
629 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
630 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
631 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
632 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
633 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
634 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
635 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
636 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
637 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
638 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
639 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
640 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
641 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
642 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
643 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
644 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
645 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
646 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
647 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
648 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
649 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
650 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
651 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
652 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
653 const float vel1Term = f_1 + f_11 + f_15 + f_7;
654 const float momdensity_1 =
655 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
656 const float vel2Term = f_12 + f_13 + f_5;
657 const float momdensity_2 =
658 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
659 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
660 vel1Term + vel2Term;
661
662 // calculate current velocity (before density change)
663 const float conversion = float(1) / rho;
664 Vector3<float> velocity;
665 velocity[0u] = momdensity_0 * conversion;
666 velocity[1u] = momdensity_1 * conversion;
667 velocity[2u] = momdensity_2 * conversion;
668
669 Equilibrium::set(pdf_field, velocity, *values_it, Cell{x, y, z});
670 ++values_it;
671 }
672 }
673 }
674}
675} // namespace Density
676
677namespace Velocity {
678inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
679 GhostLayerField<float, uint_t{3u}> const *force_field,
680 Vector3<float> const &u, Cell const &cell) {
681 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
682 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
683 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
684 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
685 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
686 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
687 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
688 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
689 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
690 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
691 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
692 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
693 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
694 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
695 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
696 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
697 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
698 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
699 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
700 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
701 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
702 const float vel1Term = f_1 + f_11 + f_15 + f_7;
703 const float vel2Term = f_12 + f_13 + f_5;
704 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
705 vel1Term + vel2Term;
706
707 const auto x = cell.x();
708 const auto y = cell.y();
709 const auto z = cell.z();
710 const float u_0 =
711 -force_field->get(x, y, z, 0) * 0.50000000000000000f / rho + u[0];
712 const float u_1 =
713 -force_field->get(x, y, z, 1) * 0.50000000000000000f / rho + u[1];
714 const float u_2 =
715 -force_field->get(x, y, z, 2) * 0.50000000000000000f / rho + u[2];
716
717 Equilibrium::set(pdf_field, Vector3<float>(u_0, u_1, u_2), rho, cell);
718}
719} // namespace Velocity
720
721namespace MomentumDensity {
722inline Vector3<float>
723reduce(GhostLayerField<float, uint_t{19u}> const *pdf_field,
724 GhostLayerField<float, uint_t{3u}> const *force_field) {
725 Vector3<float> momentumDensity(float{0});
726 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
727 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
728 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
729 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
730 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
731 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
732 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
733 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
734 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
735 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
736 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
737 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
738 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
739 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
740 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
741 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
742 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
743 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
744 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
745 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
746 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
747 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
748 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
749 const float vel1Term = f_1 + f_11 + f_15 + f_7;
750 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
751 const float vel2Term = f_12 + f_13 + f_5;
752 const float momdensity_2 =
753 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
754 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
755 vel1Term + vel2Term;
756 const float md_0 =
757 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
758 const float md_1 =
759 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
760 const float md_2 =
761 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
762
763 momentumDensity[0u] += md_0;
764 momentumDensity[1u] += md_1;
765 momentumDensity[2u] += md_2;
766 });
767 return momentumDensity;
768}
769} // namespace MomentumDensity
770
771namespace PressureTensor {
772inline Matrix3<float> get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
773 Cell const &cell) {
774 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
775 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
776 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
777 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
778 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
779 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
780 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
781 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
782 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
783 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
784 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
785 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
786 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
787 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
788 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
789 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
790 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
791 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
792 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
793 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
794 const float p_0 =
795 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
796 const float p_1 = -f_10 - f_7 + f_8 + f_9;
797 const float p_2 = -f_13 + f_14 + f_17 - f_18;
798 const float p_3 = -f_10 - f_7 + f_8 + f_9;
799 const float p_4 =
800 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
801 const float p_5 = f_11 - f_12 - f_15 + f_16;
802 const float p_6 = -f_13 + f_14 + f_17 - f_18;
803 const float p_7 = f_11 - f_12 - f_15 + f_16;
804 const float p_8 =
805 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
806
807 Matrix3<float> pressureTensor;
808 pressureTensor[0u] = p_0;
809 pressureTensor[1u] = p_1;
810 pressureTensor[2u] = p_2;
811
812 pressureTensor[3u] = p_3;
813 pressureTensor[4u] = p_4;
814 pressureTensor[5u] = p_5;
815
816 pressureTensor[6u] = p_6;
817 pressureTensor[7u] = p_7;
818 pressureTensor[8u] = p_8;
819
820 return pressureTensor;
821}
822} // namespace PressureTensor
823
824} // namespace accessor
825} // namespace lbm
826} // namespace walberla
827
828#ifdef WALBERLA_CXX_COMPILER_IS_GNU
829#pragma GCC diagnostic pop
830#endif
831
832#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
833#pragma clang diagnostic pop
834#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
ScriptInterface::Observables::PressureTensor
ParamlessObservableInterface<::Observables::PressureTensor > PressureTensor
Definition ParamlessObservable.hpp:67
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