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.3.3, lbmpy v1.3.3,
22// lbmpy_walberla/pystencils_walberla from waLBerla commit
23// b0842e1a493ce19ef1bbb8d2cf382fc343970a7f
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:97
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