Loading [MathJax]/extensions/TeX/AMSmath.js
ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages Concepts
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.3.7, lbmpy v1.3.7, sympy v1.12.1,
22// lbmpy_walberla/pystencils_walberla from waLBerla commit
23// f36fa0a68bae59f0b516f6587ea8fa7c24a41141
24
25/*
26 * Lattice field accessors.
27 * Adapted from the waLBerla source file
28 * https://i10git.cs.fau.de/walberla/walberla/-/blob/a16141524c58ab88386e2a0f8fdd7c63c5edd704/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
29 */
30
31#pragma once
32
33#include <core/DataTypes.h>
34#include <core/cell/Cell.h>
35#include <core/cell/CellInterval.h>
36#include <core/math/Matrix3.h>
37#include <core/math/Vector3.h>
38
39#include <field/GhostLayerField.h>
40#include <stencil/D3Q19.h>
41
42#include <array>
43#include <cassert>
44#include <iterator>
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#endif
52
53#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
54#pragma clang diagnostic push
55#pragma clang diagnostic ignored "-Wunused-variable"
56#endif
57
58namespace walberla {
59namespace lbm {
60namespace accessor {
61
62namespace Population {
63inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
64 Cell const &cell) {
65 float const &xyz0 = pdf_field->get(cell, uint_t{0u});
66 std::array<float, 19u> pop;
67 pop[0u] = pdf_field->getF(&xyz0, uint_t{0u});
68 pop[1u] = pdf_field->getF(&xyz0, uint_t{1u});
69 pop[2u] = pdf_field->getF(&xyz0, uint_t{2u});
70 pop[3u] = pdf_field->getF(&xyz0, uint_t{3u});
71 pop[4u] = pdf_field->getF(&xyz0, uint_t{4u});
72 pop[5u] = pdf_field->getF(&xyz0, uint_t{5u});
73 pop[6u] = pdf_field->getF(&xyz0, uint_t{6u});
74 pop[7u] = pdf_field->getF(&xyz0, uint_t{7u});
75 pop[8u] = pdf_field->getF(&xyz0, uint_t{8u});
76 pop[9u] = pdf_field->getF(&xyz0, uint_t{9u});
77 pop[10u] = pdf_field->getF(&xyz0, uint_t{10u});
78 pop[11u] = pdf_field->getF(&xyz0, uint_t{11u});
79 pop[12u] = pdf_field->getF(&xyz0, uint_t{12u});
80 pop[13u] = pdf_field->getF(&xyz0, uint_t{13u});
81 pop[14u] = pdf_field->getF(&xyz0, uint_t{14u});
82 pop[15u] = pdf_field->getF(&xyz0, uint_t{15u});
83 pop[16u] = pdf_field->getF(&xyz0, uint_t{16u});
84 pop[17u] = pdf_field->getF(&xyz0, uint_t{17u});
85 pop[18u] = pdf_field->getF(&xyz0, uint_t{18u});
86 return pop;
87}
88
89inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
90 std::array<float, 19u> const &pop, Cell const &cell) {
91 float &xyz0 = pdf_field->get(cell, uint_t{0u});
92 pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
93 pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
94 pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
95 pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
96 pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
97 pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
98 pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
99 pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
100 pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
101 pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
102 pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
103 pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
104 pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
105 pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
106 pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
107 pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
108 pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
109 pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
110 pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
111}
112
113inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
114 GhostLayerField<float, uint_t{3u}> *velocity_field,
115 GhostLayerField<float, uint_t{3u}> const *force_field,
116 std::array<float, 19u> const &pop, Cell const &cell) {
117 auto &xyz0 = pdf_field->get(cell, uint_t{0u});
118 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
119 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
120 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
121 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
122 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
123 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
124 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
125 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
126 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
127 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
128 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
129 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
130 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
131 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
132 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
133 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
134 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
135 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
136 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
137 const auto x = cell.x();
138 const auto y = cell.y();
139 const auto z = cell.z();
140 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
141 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
142 const float vel1Term = f_1 + f_11 + f_15 + f_7;
143 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
144 const float vel2Term = f_12 + f_13 + f_5;
145 const float momdensity_2 =
146 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
147 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
148 vel1Term + vel2Term;
149 const float md_0 =
150 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
151 const float md_1 =
152 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
153 const float md_2 =
154 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
155 const auto rho_inv = float{1} / rho;
156 velocity_field->get(cell, uint_t{0u}) = md_0 * rho_inv;
157 velocity_field->get(cell, uint_t{1u}) = md_1 * rho_inv;
158 velocity_field->get(cell, uint_t{2u}) = md_2 * rho_inv;
159}
160
161inline void initialize(GhostLayerField<float, uint_t{19u}> *pdf_field,
162 std::array<float, 19u> const &pop) {
163 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdf_field, {
164 float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
165 pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
166 pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
167 pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
168 pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
169 pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
170 pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
171 pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
172 pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
173 pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
174 pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
175 pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
176 pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
177 pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
178 pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
179 pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
180 pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
181 pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
182 pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
183 pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
184 });
185}
186
187inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
188 CellInterval const &ci) {
189 std::vector<float> out;
190 out.reserve(ci.numCells() * uint_t(19u));
191 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
192 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
193 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
194 float const &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
195 out.emplace_back(pdf_field->getF(&xyz0, uint_t{0u}));
196 out.emplace_back(pdf_field->getF(&xyz0, uint_t{1u}));
197 out.emplace_back(pdf_field->getF(&xyz0, uint_t{2u}));
198 out.emplace_back(pdf_field->getF(&xyz0, uint_t{3u}));
199 out.emplace_back(pdf_field->getF(&xyz0, uint_t{4u}));
200 out.emplace_back(pdf_field->getF(&xyz0, uint_t{5u}));
201 out.emplace_back(pdf_field->getF(&xyz0, uint_t{6u}));
202 out.emplace_back(pdf_field->getF(&xyz0, uint_t{7u}));
203 out.emplace_back(pdf_field->getF(&xyz0, uint_t{8u}));
204 out.emplace_back(pdf_field->getF(&xyz0, uint_t{9u}));
205 out.emplace_back(pdf_field->getF(&xyz0, uint_t{10u}));
206 out.emplace_back(pdf_field->getF(&xyz0, uint_t{11u}));
207 out.emplace_back(pdf_field->getF(&xyz0, uint_t{12u}));
208 out.emplace_back(pdf_field->getF(&xyz0, uint_t{13u}));
209 out.emplace_back(pdf_field->getF(&xyz0, uint_t{14u}));
210 out.emplace_back(pdf_field->getF(&xyz0, uint_t{15u}));
211 out.emplace_back(pdf_field->getF(&xyz0, uint_t{16u}));
212 out.emplace_back(pdf_field->getF(&xyz0, uint_t{17u}));
213 out.emplace_back(pdf_field->getF(&xyz0, uint_t{18u}));
214 }
215 }
216 }
217 return out;
218}
219
220inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
221 std::vector<float> const &values, CellInterval const &ci) {
222 assert(uint_c(values.size()) == ci.numCells() * uint_t(19u));
223 auto pop = values.data();
224 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
225 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
226 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
227 float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
228 pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
229 pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
230 pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
231 pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
232 pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
233 pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
234 pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
235 pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
236 pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
237 pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
238 pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
239 pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
240 pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
241 pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
242 pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
243 pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
244 pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
245 pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
246 pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
247 std::advance(pop, 19);
248 }
249 }
250 }
251}
252
253inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
254 GhostLayerField<float, uint_t{3u}> *velocity_field,
255 GhostLayerField<float, uint_t{3u}> const *force_field,
256 std::vector<float> const &values, CellInterval const &ci) {
257 assert(uint_c(values.size()) == ci.numCells() * uint_t(19u));
258 auto pop = values.data();
259 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
260 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
261 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
262 float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
263 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
264 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
265 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
266 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
267 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
268 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
269 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
270 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
271 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
272 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
273 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
274 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
275 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
276 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
277 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
278 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
279 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
280 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
281 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
282 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
283 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
284 const float vel1Term = f_1 + f_11 + f_15 + f_7;
285 const float momdensity_1 =
286 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
287 const float vel2Term = f_12 + f_13 + f_5;
288 const float momdensity_2 =
289 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
290 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
291 vel1Term + vel2Term;
292 const float md_0 =
293 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
294 const float md_1 =
295 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
296 const float md_2 =
297 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
298 const auto rho_inv = float{1} / rho;
299 velocity_field->get(x, y, z, uint_t{0u}) = md_0 * rho_inv;
300 velocity_field->get(x, y, z, uint_t{1u}) = md_1 * rho_inv;
301 velocity_field->get(x, y, z, uint_t{2u}) = md_2 * rho_inv;
302 std::advance(pop, 19);
303 }
304 }
305 }
306}
307} // namespace Population
308
309namespace Vector {
310inline auto get(GhostLayerField<float, uint_t{3u}> const *vec_field,
311 Cell const &cell) {
312 const float &xyz0 = vec_field->get(cell, uint_t{0u});
313 Vector3<float> vec;
314 vec[0] = vec_field->getF(&xyz0, uint_t{0u});
315 vec[1] = vec_field->getF(&xyz0, uint_t{1u});
316 vec[2] = vec_field->getF(&xyz0, uint_t{2u});
317 return vec;
318}
319
320inline void set(GhostLayerField<float, uint_t{3u}> *vec_field,
321 Vector3<float> const &vec, Cell const &cell) {
322 float &xyz0 = vec_field->get(cell, uint_t{0u});
323 vec_field->getF(&xyz0, uint_t{0u}) = vec[0u];
324 vec_field->getF(&xyz0, uint_t{1u}) = vec[1u];
325 vec_field->getF(&xyz0, uint_t{2u}) = vec[2u];
326}
327
328inline void add(GhostLayerField<float, uint_t{3u}> *vec_field,
329 Vector3<float> const &vec, Cell const &cell) {
330 float &xyz0 = vec_field->get(cell, uint_t{0u});
331 vec_field->getF(&xyz0, uint_t{0u}) += vec[0u];
332 vec_field->getF(&xyz0, uint_t{1u}) += vec[1u];
333 vec_field->getF(&xyz0, uint_t{2u}) += vec[2u];
334}
335
336inline void initialize(GhostLayerField<float, uint_t{3u}> *vec_field,
337 Vector3<float> const &vec) {
338 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
339 float &xyz0 = vec_field->get(x, y, z, uint_t{0u});
340 vec_field->getF(&xyz0, uint_t{0u}) = vec[0u];
341 vec_field->getF(&xyz0, uint_t{1u}) = vec[1u];
342 vec_field->getF(&xyz0, uint_t{2u}) = vec[2u];
343 });
344}
345
346inline void add_to_all(GhostLayerField<float, uint_t{3u}> *vec_field,
347 Vector3<float> const &vec) {
348 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
349 float &xyz0 = vec_field->get(x, y, z, uint_t{0u});
350 vec_field->getF(&xyz0, uint_t{0u}) += vec[0u];
351 vec_field->getF(&xyz0, uint_t{1u}) += vec[1u];
352 vec_field->getF(&xyz0, uint_t{2u}) += vec[2u];
353 });
354}
355
356inline auto get(GhostLayerField<float, uint_t{3u}> const *vec_field,
357 CellInterval const &ci) {
358 std::vector<float> out;
359 out.reserve(ci.numCells() * uint_t(3u));
360 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
361 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
362 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
363 const float &xyz0 = vec_field->get(x, y, z, uint_t{0u});
364 out.emplace_back(vec_field->getF(&xyz0, uint_t{0u}));
365 out.emplace_back(vec_field->getF(&xyz0, uint_t{1u}));
366 out.emplace_back(vec_field->getF(&xyz0, uint_t{2u}));
367 }
368 }
369 }
370 return out;
371}
372
373inline void set(GhostLayerField<float, uint_t{3u}> *vec_field,
374 std::vector<float> const &values, CellInterval const &ci) {
375 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
376 auto values_ptr = values.data();
377 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
378 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
379 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
380 float &xyz0 = vec_field->get(x, y, z, uint_t{0u});
381 vec_field->getF(&xyz0, uint_t{0u}) = values_ptr[0u];
382 vec_field->getF(&xyz0, uint_t{1u}) = values_ptr[1u];
383 vec_field->getF(&xyz0, uint_t{2u}) = values_ptr[2u];
384 std::advance(values_ptr, 3);
385 }
386 }
387 }
388}
389} // namespace Vector
390
391namespace EquilibriumDistribution {
392inline float get(stencil::Direction const direction,
393 Vector3<float> const &u = Vector3<float>(float{0}),
394 float rho = float{1}) {
395
396 using namespace stencil;
397 switch (direction) {
398 case C:
399 return rho * -0.33333333333333331f * (u[0] * u[0]) +
400 rho * -0.33333333333333331f * (u[1] * u[1]) +
401 rho * -0.33333333333333331f * (u[2] * u[2]) +
402 rho * 0.33333333333333331f;
403 case N:
404 return rho * -0.16666666666666666f * (u[0] * u[0]) +
405 rho * -0.16666666666666666f * (u[2] * u[2]) +
406 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[1] +
407 rho * 0.16666666666666666f * (u[1] * u[1]);
408 case S:
409 return rho * -0.16666666666666666f * u[1] +
410 rho * -0.16666666666666666f * (u[0] * u[0]) +
411 rho * -0.16666666666666666f * (u[2] * u[2]) +
412 rho * 0.055555555555555552f +
413 rho * 0.16666666666666666f * (u[1] * u[1]);
414 case W:
415 return rho * -0.16666666666666666f * u[0] +
416 rho * -0.16666666666666666f * (u[1] * u[1]) +
417 rho * -0.16666666666666666f * (u[2] * u[2]) +
418 rho * 0.055555555555555552f +
419 rho * 0.16666666666666666f * (u[0] * u[0]);
420 case E:
421 return rho * -0.16666666666666666f * (u[1] * u[1]) +
422 rho * -0.16666666666666666f * (u[2] * u[2]) +
423 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[0] +
424 rho * 0.16666666666666666f * (u[0] * u[0]);
425 case T:
426 return rho * -0.16666666666666666f * (u[0] * u[0]) +
427 rho * -0.16666666666666666f * (u[1] * u[1]) +
428 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[2] +
429 rho * 0.16666666666666666f * (u[2] * u[2]);
430 case B:
431 return rho * -0.16666666666666666f * u[2] +
432 rho * -0.16666666666666666f * (u[0] * u[0]) +
433 rho * -0.16666666666666666f * (u[1] * u[1]) +
434 rho * 0.055555555555555552f +
435 rho * 0.16666666666666666f * (u[2] * u[2]);
436 case NW:
437 return rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[1] +
438 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
439 rho * 0.083333333333333329f * (u[0] * u[0]) +
440 rho * 0.083333333333333329f * (u[1] * u[1]);
441 case NE:
442 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
443 rho * 0.083333333333333329f * u[1] +
444 rho * 0.083333333333333329f * (u[0] * u[0]) +
445 rho * 0.083333333333333329f * (u[1] * u[1]) +
446 rho * 0.25f * u[0] * u[1];
447 case SW:
448 return rho * -0.083333333333333329f * u[0] +
449 rho * -0.083333333333333329f * u[1] + rho * 0.027777777777777776f +
450 rho * 0.083333333333333329f * (u[0] * u[0]) +
451 rho * 0.083333333333333329f * (u[1] * u[1]) +
452 rho * 0.25f * u[0] * u[1];
453 case SE:
454 return rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[0] * u[1] +
455 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
456 rho * 0.083333333333333329f * (u[0] * u[0]) +
457 rho * 0.083333333333333329f * (u[1] * u[1]);
458 case TN:
459 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
460 rho * 0.083333333333333329f * u[2] +
461 rho * 0.083333333333333329f * (u[1] * u[1]) +
462 rho * 0.083333333333333329f * (u[2] * u[2]) +
463 rho * 0.25f * u[1] * u[2];
464 case TS:
465 return rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[1] * u[2] +
466 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
467 rho * 0.083333333333333329f * (u[1] * u[1]) +
468 rho * 0.083333333333333329f * (u[2] * u[2]);
469 case TW:
470 return rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[2] +
471 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
472 rho * 0.083333333333333329f * (u[0] * u[0]) +
473 rho * 0.083333333333333329f * (u[2] * u[2]);
474 case TE:
475 return rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
476 rho * 0.083333333333333329f * u[2] +
477 rho * 0.083333333333333329f * (u[0] * u[0]) +
478 rho * 0.083333333333333329f * (u[2] * u[2]) +
479 rho * 0.25f * u[0] * u[2];
480 case BN:
481 return rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[1] * u[2] +
482 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
483 rho * 0.083333333333333329f * (u[1] * u[1]) +
484 rho * 0.083333333333333329f * (u[2] * u[2]);
485 case BS:
486 return rho * -0.083333333333333329f * u[1] +
487 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
488 rho * 0.083333333333333329f * (u[1] * u[1]) +
489 rho * 0.083333333333333329f * (u[2] * u[2]) +
490 rho * 0.25f * u[1] * u[2];
491 case BW:
492 return rho * -0.083333333333333329f * u[0] +
493 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
494 rho * 0.083333333333333329f * (u[0] * u[0]) +
495 rho * 0.083333333333333329f * (u[2] * u[2]) +
496 rho * 0.25f * u[0] * u[2];
497 case BE:
498 return rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[0] * u[2] +
499 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
500 rho * 0.083333333333333329f * (u[0] * u[0]) +
501 rho * 0.083333333333333329f * (u[2] * u[2]);
502 default:
503 WALBERLA_ABORT("Invalid Direction")
504 }
505}
506} // namespace EquilibriumDistribution
507
508namespace Equilibrium {
509inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
510 Vector3<float> const &u, float const rho, Cell const &cell) {
511
512 float &xyz0 = pdf_field->get(cell, uint_t{0u});
513 pdf_field->getF(&xyz0, uint_t{0u}) =
514 rho * -0.33333333333333331f * (u[0] * u[0]) +
515 rho * -0.33333333333333331f * (u[1] * u[1]) +
516 rho * -0.33333333333333331f * (u[2] * u[2]) + rho * 0.33333333333333331f;
517 pdf_field->getF(&xyz0, uint_t{1u}) =
518 rho * -0.16666666666666666f * (u[0] * u[0]) +
519 rho * -0.16666666666666666f * (u[2] * u[2]) +
520 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[1] +
521 rho * 0.16666666666666666f * (u[1] * u[1]);
522 pdf_field->getF(&xyz0, uint_t{2u}) =
523 rho * -0.16666666666666666f * u[1] +
524 rho * -0.16666666666666666f * (u[0] * u[0]) +
525 rho * -0.16666666666666666f * (u[2] * u[2]) +
526 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[1] * u[1]);
527 pdf_field->getF(&xyz0, uint_t{3u}) =
528 rho * -0.16666666666666666f * u[0] +
529 rho * -0.16666666666666666f * (u[1] * u[1]) +
530 rho * -0.16666666666666666f * (u[2] * u[2]) +
531 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[0] * u[0]);
532 pdf_field->getF(&xyz0, uint_t{4u}) =
533 rho * -0.16666666666666666f * (u[1] * u[1]) +
534 rho * -0.16666666666666666f * (u[2] * u[2]) +
535 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[0] +
536 rho * 0.16666666666666666f * (u[0] * u[0]);
537 pdf_field->getF(&xyz0, uint_t{5u}) =
538 rho * -0.16666666666666666f * (u[0] * u[0]) +
539 rho * -0.16666666666666666f * (u[1] * u[1]) +
540 rho * 0.055555555555555552f + rho * 0.16666666666666666f * u[2] +
541 rho * 0.16666666666666666f * (u[2] * u[2]);
542 pdf_field->getF(&xyz0, uint_t{6u}) =
543 rho * -0.16666666666666666f * u[2] +
544 rho * -0.16666666666666666f * (u[0] * u[0]) +
545 rho * -0.16666666666666666f * (u[1] * u[1]) +
546 rho * 0.055555555555555552f + rho * 0.16666666666666666f * (u[2] * u[2]);
547 pdf_field->getF(&xyz0, uint_t{7u}) =
548 rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[1] +
549 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
550 rho * 0.083333333333333329f * (u[0] * u[0]) +
551 rho * 0.083333333333333329f * (u[1] * u[1]);
552 pdf_field->getF(&xyz0, uint_t{8u}) =
553 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
554 rho * 0.083333333333333329f * u[1] +
555 rho * 0.083333333333333329f * (u[0] * u[0]) +
556 rho * 0.083333333333333329f * (u[1] * u[1]) + rho * 0.25f * u[0] * u[1];
557 pdf_field->getF(&xyz0, uint_t{9u}) =
558 rho * -0.083333333333333329f * u[0] +
559 rho * -0.083333333333333329f * u[1] + rho * 0.027777777777777776f +
560 rho * 0.083333333333333329f * (u[0] * u[0]) +
561 rho * 0.083333333333333329f * (u[1] * u[1]) + rho * 0.25f * u[0] * u[1];
562 pdf_field->getF(&xyz0, uint_t{10u}) =
563 rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[0] * u[1] +
564 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
565 rho * 0.083333333333333329f * (u[0] * u[0]) +
566 rho * 0.083333333333333329f * (u[1] * u[1]);
567 pdf_field->getF(&xyz0, uint_t{11u}) =
568 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
569 rho * 0.083333333333333329f * u[2] +
570 rho * 0.083333333333333329f * (u[1] * u[1]) +
571 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[1] * u[2];
572 pdf_field->getF(&xyz0, uint_t{12u}) =
573 rho * -0.083333333333333329f * u[1] + rho * -0.25f * u[1] * u[2] +
574 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
575 rho * 0.083333333333333329f * (u[1] * u[1]) +
576 rho * 0.083333333333333329f * (u[2] * u[2]);
577 pdf_field->getF(&xyz0, uint_t{13u}) =
578 rho * -0.083333333333333329f * u[0] + rho * -0.25f * u[0] * u[2] +
579 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[2] +
580 rho * 0.083333333333333329f * (u[0] * u[0]) +
581 rho * 0.083333333333333329f * (u[2] * u[2]);
582 pdf_field->getF(&xyz0, uint_t{14u}) =
583 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
584 rho * 0.083333333333333329f * u[2] +
585 rho * 0.083333333333333329f * (u[0] * u[0]) +
586 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[0] * u[2];
587 pdf_field->getF(&xyz0, uint_t{15u}) =
588 rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[1] * u[2] +
589 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[1] +
590 rho * 0.083333333333333329f * (u[1] * u[1]) +
591 rho * 0.083333333333333329f * (u[2] * u[2]);
592 pdf_field->getF(&xyz0, uint_t{16u}) =
593 rho * -0.083333333333333329f * u[1] +
594 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
595 rho * 0.083333333333333329f * (u[1] * u[1]) +
596 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[1] * u[2];
597 pdf_field->getF(&xyz0, uint_t{17u}) =
598 rho * -0.083333333333333329f * u[0] +
599 rho * -0.083333333333333329f * u[2] + rho * 0.027777777777777776f +
600 rho * 0.083333333333333329f * (u[0] * u[0]) +
601 rho * 0.083333333333333329f * (u[2] * u[2]) + rho * 0.25f * u[0] * u[2];
602 pdf_field->getF(&xyz0, uint_t{18u}) =
603 rho * -0.083333333333333329f * u[2] + rho * -0.25f * u[0] * u[2] +
604 rho * 0.027777777777777776f + rho * 0.083333333333333329f * u[0] +
605 rho * 0.083333333333333329f * (u[0] * u[0]) +
606 rho * 0.083333333333333329f * (u[2] * u[2]);
607}
608} // namespace Equilibrium
609
610namespace Density {
611inline float get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
612 Cell const &cell) {
613 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
614 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
615 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
616 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
617 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
618 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
619 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
620 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
621 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
622 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
623 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
624 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
625 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
626 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
627 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
628 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
629 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
630 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
631 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
632 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
633 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
634 const float vel1Term = f_1 + f_11 + f_15 + f_7;
635 const float vel2Term = f_12 + f_13 + f_5;
636 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
637 vel1Term + vel2Term;
638 return rho;
639}
640
641inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
642 float const rho_in, Cell const &cell) {
643 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
644 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
645 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
646 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
647 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
648 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
649 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
650 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
651 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
652 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
653 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
654 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
655 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
656 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
657 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
658 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
659 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
660 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
661 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
662 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
663 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
664 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
665 const float vel1Term = f_1 + f_11 + f_15 + f_7;
666 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
667 const float vel2Term = f_12 + f_13 + f_5;
668 const float momdensity_2 =
669 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
670 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
671 vel1Term + vel2Term;
672
673 // calculate current velocity (before density change)
674 const float conversion = float{1} / rho;
675 Vector3<float> velocity;
676 velocity[0u] = momdensity_0 * conversion;
677 velocity[1u] = momdensity_1 * conversion;
678 velocity[2u] = momdensity_2 * conversion;
679
680 Equilibrium::set(pdf_field, velocity, rho_in, cell);
681}
682
683inline std::vector<float>
684get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
685 CellInterval const &ci) {
686 std::vector<float> out;
687 out.reserve(ci.numCells());
688 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
689 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
690 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
691 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
692 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
693 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
694 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
695 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
696 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
697 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
698 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
699 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
700 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
701 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
702 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
703 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
704 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
705 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
706 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
707 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
708 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
709 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
710 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
711 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
712 const float vel1Term = f_1 + f_11 + f_15 + f_7;
713 const float vel2Term = f_12 + f_13 + f_5;
714 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
715 vel1Term + vel2Term;
716 out.emplace_back(rho);
717 }
718 }
719 }
720 return out;
721}
722
723inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
724 std::vector<float> const &values, CellInterval const &ci) {
725 assert(uint_c(values.size()) == ci.numCells());
726 auto values_it = values.begin();
727 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
728 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
729 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
730 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
731 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
732 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
733 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
734 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
735 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
736 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
737 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
738 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
739 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
740 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
741 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
742 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
743 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
744 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
745 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
746 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
747 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
748 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
749 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
750 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
751 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
752 const float vel1Term = f_1 + f_11 + f_15 + f_7;
753 const float momdensity_1 =
754 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
755 const float vel2Term = f_12 + f_13 + f_5;
756 const float momdensity_2 =
757 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
758 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
759 vel1Term + vel2Term;
760
761 // calculate current velocity (before density change)
762 const float conversion = float{1} / rho;
763 Vector3<float> velocity;
764 velocity[0u] = momdensity_0 * conversion;
765 velocity[1u] = momdensity_1 * conversion;
766 velocity[2u] = momdensity_2 * conversion;
767
768 Equilibrium::set(pdf_field, velocity, *values_it, Cell{x, y, z});
769 ++values_it;
770 }
771 }
772 }
773}
774} // namespace Density
775
776namespace Velocity {
777inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
778 GhostLayerField<float, uint_t{3u}> const *force_field,
779 Cell const &cell) {
780 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
781 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
782 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
783 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
784 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
785 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
786 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
787 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
788 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
789 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
790 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
791 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
792 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
793 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
794 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
795 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
796 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
797 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
798 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
799 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
800 const auto x = cell.x();
801 const auto y = cell.y();
802 const auto z = cell.z();
803 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
804 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
805 const float vel1Term = f_1 + f_11 + f_15 + f_7;
806 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
807 const float vel2Term = f_12 + f_13 + f_5;
808 const float momdensity_2 =
809 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
810 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
811 vel1Term + vel2Term;
812 const float md_0 =
813 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
814 const float md_1 =
815 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
816 const float md_2 =
817 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
818 const float rho_inv = float{1} / rho;
819
820 return Vector3<float>(md_0 * rho_inv, md_1 * rho_inv, md_2 * rho_inv);
821}
822
823inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
824 GhostLayerField<float, uint_t{3u}> const *force_field,
825 CellInterval const &ci) {
826 std::vector<float> out;
827 out.reserve(ci.numCells() * uint_t(3u));
828 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
829 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
830 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
831 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
832 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
833 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
834 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
835 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
836 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
837 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
838 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
839 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
840 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
841 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
842 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
843 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
844 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
845 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
846 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
847 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
848 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
849 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
850 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
851 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
852 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
853 const float vel1Term = f_1 + f_11 + f_15 + f_7;
854 const float momdensity_1 =
855 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
856 const float vel2Term = f_12 + f_13 + f_5;
857 const float momdensity_2 =
858 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
859 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
860 vel1Term + vel2Term;
861 const float md_0 =
862 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
863 const float md_1 =
864 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
865 const float md_2 =
866 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
867 const float rho_inv = float{1} / rho;
868 out.emplace_back(md_0 * rho_inv);
869 out.emplace_back(md_1 * rho_inv);
870 out.emplace_back(md_2 * rho_inv);
871 }
872 }
873 }
874 return out;
875}
876
877inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
878 GhostLayerField<float, uint_t{3u}> *velocity_field,
879 GhostLayerField<float, uint_t{3u}> const *force_field,
880 Vector3<float> const &u, Cell const &cell) {
881 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
882 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
883 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
884 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
885 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
886 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
887 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
888 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
889 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
890 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
891 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
892 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
893 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
894 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
895 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
896 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
897 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
898 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
899 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
900 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
901 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
902 const float vel1Term = f_1 + f_11 + f_15 + f_7;
903 const float vel2Term = f_12 + f_13 + f_5;
904 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
905 vel1Term + vel2Term;
906
907 const auto x = cell.x();
908 const auto y = cell.y();
909 const auto z = cell.z();
910 const float u_0 =
911 -force_field->get(x, y, z, 0) * 0.50000000000000000f / rho + u[0];
912 const float u_1 =
913 -force_field->get(x, y, z, 1) * 0.50000000000000000f / rho + u[1];
914 const float u_2 =
915 -force_field->get(x, y, z, 2) * 0.50000000000000000f / rho + u[2];
916 velocity_field->get(x, y, z, uint_t{0u}) = u[0u];
917 velocity_field->get(x, y, z, uint_t{1u}) = u[1u];
918 velocity_field->get(x, y, z, uint_t{2u}) = u[2u];
919
920 Equilibrium::set(pdf_field, Vector3<float>(u_0, u_1, u_2), rho, cell);
921}
922
923inline void set(GhostLayerField<float, uint_t{19u}> *pdf_field,
924 GhostLayerField<float, uint_t{3u}> *velocity_field,
925 GhostLayerField<float, uint_t{3u}> const *force_field,
926 std::vector<float> const &values, CellInterval const &ci) {
927 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
928 auto u = values.data();
929 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
930 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
931 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
932 float &pdf_xyz0 = pdf_field->get(x, y, z, uint_t{0u});
933 float &vel_xyz0 = velocity_field->get(x, y, z, uint_t{0u});
934 const float f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
935 const float f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
936 const float f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
937 const float f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
938 const float f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
939 const float f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
940 const float f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
941 const float f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
942 const float f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
943 const float f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
944 const float f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
945 const float f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
946 const float f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
947 const float f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
948 const float f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
949 const float f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
950 const float f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
951 const float f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
952 const float f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
953 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
954 const float vel1Term = f_1 + f_11 + f_15 + f_7;
955 const float vel2Term = f_12 + f_13 + f_5;
956 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
957 vel1Term + vel2Term;
958
959 const float u_0 =
960 -force_field->get(x, y, z, 0) * 0.50000000000000000f / rho + u[0];
961 const float u_1 =
962 -force_field->get(x, y, z, 1) * 0.50000000000000000f / rho + u[1];
963 const float u_2 =
964 -force_field->get(x, y, z, 2) * 0.50000000000000000f / rho + u[2];
965 velocity_field->getF(&vel_xyz0, uint_t{0u}) = u[0u];
966 velocity_field->getF(&vel_xyz0, uint_t{1u}) = u[1u];
967 velocity_field->getF(&vel_xyz0, uint_t{2u}) = u[2u];
968
969 std::advance(u, 3);
970
971 Equilibrium::set(pdf_field, Vector3<float>(u_0, u_1, u_2), rho,
972 Cell{x, y, z});
973 }
974 }
975 }
976}
977} // namespace Velocity
978
979namespace Force {
980inline void set(GhostLayerField<float, uint_t{19u}> const *pdf_field,
981 GhostLayerField<float, uint_t{3u}> *velocity_field,
982 GhostLayerField<float, uint_t{3u}> *force_field,
983 Vector3<float> const &force, Cell const &cell) {
984 float const &pdf_xyz0 = pdf_field->get(cell, uint_t{0u});
985 float &vel_xyz0 = velocity_field->get(cell, uint_t{0u});
986 float &laf_xyz0 = force_field->get(cell, uint_t{0u});
987 const float f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
988 const float f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
989 const float f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
990 const float f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
991 const float f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
992 const float f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
993 const float f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
994 const float f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
995 const float f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
996 const float f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
997 const float f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
998 const float f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
999 const float f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
1000 const float f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
1001 const float f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
1002 const float f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
1003 const float f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
1004 const float f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
1005 const float f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
1006 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1007 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1008 const float vel1Term = f_1 + f_11 + f_15 + f_7;
1009 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1010 const float vel2Term = f_12 + f_13 + f_5;
1011 const float momdensity_2 =
1012 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1013 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
1014 vel1Term + vel2Term;
1015 const float md_0 = force[0u] * 0.50000000000000000f + momdensity_0;
1016 const float md_1 = force[1u] * 0.50000000000000000f + momdensity_1;
1017 const float md_2 = force[2u] * 0.50000000000000000f + momdensity_2;
1018 auto const rho_inv = float{1} / rho;
1019
1020 force_field->getF(&laf_xyz0, uint_t{0u}) = force[0u];
1021 force_field->getF(&laf_xyz0, uint_t{1u}) = force[1u];
1022 force_field->getF(&laf_xyz0, uint_t{2u}) = force[2u];
1023
1024 velocity_field->getF(&vel_xyz0, uint_t{0u}) = md_0 * rho_inv;
1025 velocity_field->getF(&vel_xyz0, uint_t{1u}) = md_1 * rho_inv;
1026 velocity_field->getF(&vel_xyz0, uint_t{2u}) = md_2 * rho_inv;
1027}
1028
1029inline void set(GhostLayerField<float, uint_t{19u}> const *pdf_field,
1030 GhostLayerField<float, uint_t{3u}> *velocity_field,
1031 GhostLayerField<float, uint_t{3u}> *force_field,
1032 std::vector<float> const &values, CellInterval const &ci) {
1033 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
1034 auto force = values.data();
1035 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
1036 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
1037 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
1038 float const &pdf_xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1039 float &vel_xyz0 = velocity_field->get(x, y, z, uint_t{0u});
1040 float &laf_xyz0 = force_field->get(x, y, z, uint_t{0u});
1041 const float f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
1042 const float f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
1043 const float f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
1044 const float f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
1045 const float f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
1046 const float f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
1047 const float f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
1048 const float f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
1049 const float f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
1050 const float f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
1051 const float f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
1052 const float f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
1053 const float f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
1054 const float f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
1055 const float f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
1056 const float f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
1057 const float f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
1058 const float f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
1059 const float f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
1060 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1061 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1062 const float vel1Term = f_1 + f_11 + f_15 + f_7;
1063 const float momdensity_1 =
1064 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1065 const float vel2Term = f_12 + f_13 + f_5;
1066 const float momdensity_2 =
1067 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1068 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
1069 vel1Term + vel2Term;
1070 const float md_0 = force[0u] * 0.50000000000000000f + momdensity_0;
1071 const float md_1 = force[1u] * 0.50000000000000000f + momdensity_1;
1072 const float md_2 = force[2u] * 0.50000000000000000f + momdensity_2;
1073 auto const rho_inv = float{1} / rho;
1074
1075 force_field->getF(&laf_xyz0, uint_t{0u}) = force[0u];
1076 force_field->getF(&laf_xyz0, uint_t{1u}) = force[1u];
1077 force_field->getF(&laf_xyz0, uint_t{2u}) = force[2u];
1078
1079 velocity_field->getF(&vel_xyz0, uint_t{0u}) = md_0 * rho_inv;
1080 velocity_field->getF(&vel_xyz0, uint_t{1u}) = md_1 * rho_inv;
1081 velocity_field->getF(&vel_xyz0, uint_t{2u}) = md_2 * rho_inv;
1082
1083 std::advance(force, 3);
1084 }
1085 }
1086 }
1087}
1088} // namespace Force
1089
1090namespace MomentumDensity {
1091inline auto reduce(GhostLayerField<float, uint_t{19u}> const *pdf_field,
1092 GhostLayerField<float, uint_t{3u}> const *force_field) {
1093 Vector3<float> momentumDensity(float{0});
1094 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
1095 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1096 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1097 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1098 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1099 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1100 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1101 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1102 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1103 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1104 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1105 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1106 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1107 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1108 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1109 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1110 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1111 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1112 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1113 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1114 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1115 const float vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1116 const float momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1117 const float vel1Term = f_1 + f_11 + f_15 + f_7;
1118 const float momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1119 const float vel2Term = f_12 + f_13 + f_5;
1120 const float momdensity_2 =
1121 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1122 const float rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
1123 vel1Term + vel2Term;
1124 const float md_0 =
1125 force_field->get(x, y, z, 0) * 0.50000000000000000f + momdensity_0;
1126 const float md_1 =
1127 force_field->get(x, y, z, 1) * 0.50000000000000000f + momdensity_1;
1128 const float md_2 =
1129 force_field->get(x, y, z, 2) * 0.50000000000000000f + momdensity_2;
1130
1131 momentumDensity[0u] += md_0;
1132 momentumDensity[1u] += md_1;
1133 momentumDensity[2u] += md_2;
1134 });
1135 return momentumDensity;
1136}
1137} // namespace MomentumDensity
1138
1139namespace PressureTensor {
1140inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
1141 Cell const &cell) {
1142 const float &xyz0 = pdf_field->get(cell, uint_t{0u});
1143 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1144 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1145 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1146 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1147 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1148 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1149 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1150 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1151 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1152 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1153 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1154 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1155 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1156 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1157 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1158 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1159 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1160 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1161 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1162 const float p_0 =
1163 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1164 const float p_1 = -f_10 - f_7 + f_8 + f_9;
1165 const float p_2 = -f_13 + f_14 + f_17 - f_18;
1166 const float p_3 = -f_10 - f_7 + f_8 + f_9;
1167 const float p_4 =
1168 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1169 const float p_5 = f_11 - f_12 - f_15 + f_16;
1170 const float p_6 = -f_13 + f_14 + f_17 - f_18;
1171 const float p_7 = f_11 - f_12 - f_15 + f_16;
1172 const float p_8 =
1173 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1174
1175 Matrix3<float> pressureTensor;
1176 pressureTensor[0u] = p_0;
1177 pressureTensor[1u] = p_1;
1178 pressureTensor[2u] = p_2;
1179
1180 pressureTensor[3u] = p_3;
1181 pressureTensor[4u] = p_4;
1182 pressureTensor[5u] = p_5;
1183
1184 pressureTensor[6u] = p_6;
1185 pressureTensor[7u] = p_7;
1186 pressureTensor[8u] = p_8;
1187
1188 return pressureTensor;
1189}
1190
1191inline auto get(GhostLayerField<float, uint_t{19u}> const *pdf_field,
1192 CellInterval const &ci) {
1193 std::vector<float> out;
1194 out.reserve(ci.numCells() * uint_t(9u));
1195 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
1196 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
1197 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
1198 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1199 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1200 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1201 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1202 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1203 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1204 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1205 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1206 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1207 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1208 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1209 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1210 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1211 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1212 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1213 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1214 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1215 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1216 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1217 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1218 const float p_0 =
1219 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1220 const float p_1 = -f_10 - f_7 + f_8 + f_9;
1221 const float p_2 = -f_13 + f_14 + f_17 - f_18;
1222 const float p_3 = -f_10 - f_7 + f_8 + f_9;
1223 const float p_4 =
1224 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1225 const float p_5 = f_11 - f_12 - f_15 + f_16;
1226 const float p_6 = -f_13 + f_14 + f_17 - f_18;
1227 const float p_7 = f_11 - f_12 - f_15 + f_16;
1228 const float p_8 =
1229 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1230
1231 out.emplace_back(p_0);
1232 out.emplace_back(p_1);
1233 out.emplace_back(p_2);
1234
1235 out.emplace_back(p_3);
1236 out.emplace_back(p_4);
1237 out.emplace_back(p_5);
1238
1239 out.emplace_back(p_6);
1240 out.emplace_back(p_7);
1241 out.emplace_back(p_8);
1242 }
1243 }
1244 }
1245 return out;
1246}
1247
1248inline auto reduce(GhostLayerField<float, uint_t{19u}> const *pdf_field) {
1249 Matrix3<float> pressureTensor(float{0});
1250 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
1251 const float &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1252 const float f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1253 const float f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1254 const float f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1255 const float f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1256 const float f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1257 const float f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1258 const float f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1259 const float f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1260 const float f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1261 const float f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1262 const float f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1263 const float f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1264 const float f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1265 const float f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1266 const float f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1267 const float f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1268 const float f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1269 const float f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1270 const float f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1271 const float p_0 =
1272 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1273 const float p_1 = -f_10 - f_7 + f_8 + f_9;
1274 const float p_2 = -f_13 + f_14 + f_17 - f_18;
1275 const float p_3 = -f_10 - f_7 + f_8 + f_9;
1276 const float p_4 =
1277 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1278 const float p_5 = f_11 - f_12 - f_15 + f_16;
1279 const float p_6 = -f_13 + f_14 + f_17 - f_18;
1280 const float p_7 = f_11 - f_12 - f_15 + f_16;
1281 const float p_8 =
1282 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1283
1284 pressureTensor[0u] += p_0;
1285 pressureTensor[1u] += p_1;
1286 pressureTensor[2u] += p_2;
1287
1288 pressureTensor[3u] += p_3;
1289 pressureTensor[4u] += p_4;
1290 pressureTensor[5u] += p_5;
1291
1292 pressureTensor[6u] += p_6;
1293 pressureTensor[7u] += p_7;
1294 pressureTensor[8u] += p_8;
1295 });
1296 return pressureTensor;
1297}
1298} // namespace PressureTensor
1299
1300} // namespace accessor
1301} // namespace lbm
1302} // namespace walberla
1303
1304#ifdef WALBERLA_CXX_COMPILER_IS_GNU
1305#pragma GCC diagnostic pop
1306#endif
1307
1308#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
1309#pragma clang diagnostic pop
1310#endif
Cell
Definition Cell.hpp:96
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:638
walberla::lbm::accessor::Density::get
double get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:608
walberla::lbm::accessor::EquilibriumDistribution::get
double get(stencil::Direction const direction, Vector3< double > const &u=Vector3< double >(double{0}), double rho=double{1})
Definition FieldAccessorsDoublePrecision.h:392
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:509
walberla::lbm::accessor::Force::set
void set(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> *velocity_field, GhostLayerField< double, uint_t{3u}> *force_field, Vector3< double > const &force, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:977
walberla::lbm::accessor::MomentumDensity::reduce
auto reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field)
Definition FieldAccessorsDoublePrecision.h:1088
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:89
walberla::lbm::accessor::Population::get
auto get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:63
walberla::lbm::accessor::Population::initialize
void initialize(GhostLayerField< double, uint_t{19u}> *pdf_field, std::array< double, 19u > const &pop)
Definition FieldAccessorsDoublePrecision.h:161
walberla::lbm::accessor::PressureTensor::get
auto get(GhostLayerField< double, uint_t{19u}> const *pdf_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:1138
walberla::lbm::accessor::PressureTensor::reduce
auto reduce(GhostLayerField< double, uint_t{19u}> const *pdf_field)
Definition FieldAccessorsDoublePrecision.h:1246
walberla::lbm::accessor::Vector::initialize
void initialize(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec)
Definition FieldAccessorsDoublePrecision.h:336
walberla::lbm::accessor::Vector::add
void add(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:328
walberla::lbm::accessor::Vector::set
void set(GhostLayerField< double, uint_t{3u}> *vec_field, Vector3< double > const &vec, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:320
walberla::lbm::accessor::Vector::get
auto get(GhostLayerField< double, uint_t{3u}> const *vec_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:310
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:346
walberla::lbm::accessor::Velocity::set
void set(GhostLayerField< double, uint_t{19u}> *pdf_field, GhostLayerField< double, uint_t{3u}> *velocity_field, GhostLayerField< double, uint_t{3u}> const *force_field, Vector3< double > const &u, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:874
walberla::lbm::accessor::Velocity::get
auto get(GhostLayerField< double, uint_t{19u}> const *pdf_field, GhostLayerField< double, uint_t{3u}> const *force_field, Cell const &cell)
Definition FieldAccessorsDoublePrecision.h:774
walberla
\file PackInfoPdfDoublePrecision.cpp \author pystencils
Definition EKWalberla.hpp:36
velocity
static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs)
Velocity of the virtual site.
Definition relative.cpp:64