ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
FieldAccessorsDoublePrecision.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2021-2023 The ESPResSo project
3 * Copyright (C) 2020 The waLBerla project
4 *
5 * This file is part of ESPResSo.
6 *
7 * ESPResSo is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * ESPResSo is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21// kernel generated with pystencils v1.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<double, uint_t{19u}> const *pdf_field,
64 Cell const &cell) {
65 double const &xyz0 = pdf_field->get(cell, uint_t{0u});
66 std::array<double, 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<double, uint_t{19u}> *pdf_field,
90 std::array<double, 19u> const &pop, Cell const &cell) {
91 double &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<double, uint_t{19u}> *pdf_field,
114 GhostLayerField<double, uint_t{3u}> *velocity_field,
115 GhostLayerField<double, uint_t{3u}> const *force_field,
116 std::array<double, 19u> const &pop, Cell const &cell) {
117 auto &xyz0 = pdf_field->get(cell, uint_t{0u});
118 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
119 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
120 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
121 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
122 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
123 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
124 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
125 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
126 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
127 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
128 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
129 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
130 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
131 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
132 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
133 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
134 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
135 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
136 const double 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 double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
141 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
142 const double vel1Term = f_1 + f_11 + f_15 + f_7;
143 const double momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
144 const double vel2Term = f_12 + f_13 + f_5;
145 const double momdensity_2 =
146 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
147 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
148 vel1Term + vel2Term;
149 const double md_0 =
150 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
151 const double md_1 =
152 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
153 const double md_2 =
154 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
155 const auto rho_inv = double{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<double, uint_t{19u}> *pdf_field,
162 std::array<double, 19u> const &pop) {
163 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdf_field, {
164 double &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<double, uint_t{19u}> const *pdf_field,
188 CellInterval const &ci) {
189 std::vector<double> 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 double 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<double, uint_t{19u}> *pdf_field,
221 std::vector<double> 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 double &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<double, uint_t{19u}> *pdf_field,
254 GhostLayerField<double, uint_t{3u}> *velocity_field,
255 GhostLayerField<double, uint_t{3u}> const *force_field,
256 std::vector<double> 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 double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
263 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u}) = pop[0u];
264 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u}) = pop[1u];
265 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u}) = pop[2u];
266 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u}) = pop[3u];
267 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u}) = pop[4u];
268 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u}) = pop[5u];
269 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u}) = pop[6u];
270 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u}) = pop[7u];
271 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u}) = pop[8u];
272 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u}) = pop[9u];
273 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u}) = pop[10u];
274 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u}) = pop[11u];
275 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u}) = pop[12u];
276 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u}) = pop[13u];
277 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u}) = pop[14u];
278 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u}) = pop[15u];
279 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u}) = pop[16u];
280 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u}) = pop[17u];
281 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u}) = pop[18u];
282 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
283 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
284 const double vel1Term = f_1 + f_11 + f_15 + f_7;
285 const double momdensity_1 =
286 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
287 const double vel2Term = f_12 + f_13 + f_5;
288 const double momdensity_2 =
289 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
290 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
291 vel0Term + vel1Term + vel2Term;
292 const double md_0 =
293 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
294 const double md_1 =
295 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
296 const double md_2 =
297 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
298 const auto rho_inv = double{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<double, uint_t{3u}> const *vec_field,
311 Cell const &cell) {
312 const double &xyz0 = vec_field->get(cell, uint_t{0u});
313 Vector3<double> 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<double, uint_t{3u}> *vec_field,
321 Vector3<double> const &vec, Cell const &cell) {
322 double &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<double, uint_t{3u}> *vec_field,
329 Vector3<double> const &vec, Cell const &cell) {
330 double &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<double, uint_t{3u}> *vec_field,
337 Vector3<double> const &vec) {
338 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
339 double &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<double, uint_t{3u}> *vec_field,
347 Vector3<double> const &vec) {
348 WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(vec_field, {
349 double &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<double, uint_t{3u}> const *vec_field,
357 CellInterval const &ci) {
358 std::vector<double> 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 double &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<double, uint_t{3u}> *vec_field,
374 std::vector<double> 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 double &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 double get(stencil::Direction const direction,
393 Vector3<double> const &u = Vector3<double>(double{0}),
394 double rho = double{1}) {
395
396 using namespace stencil;
397 switch (direction) {
398 case C:
399 return rho * -0.33333333333333331 * (u[0] * u[0]) +
400 rho * -0.33333333333333331 * (u[1] * u[1]) +
401 rho * -0.33333333333333331 * (u[2] * u[2]) +
402 rho * 0.33333333333333331;
403 case N:
404 return rho * -0.16666666666666666 * (u[0] * u[0]) +
405 rho * -0.16666666666666666 * (u[2] * u[2]) +
406 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[1] +
407 rho * 0.16666666666666666 * (u[1] * u[1]);
408 case S:
409 return rho * -0.16666666666666666 * u[1] +
410 rho * -0.16666666666666666 * (u[0] * u[0]) +
411 rho * -0.16666666666666666 * (u[2] * u[2]) +
412 rho * 0.055555555555555552 +
413 rho * 0.16666666666666666 * (u[1] * u[1]);
414 case W:
415 return rho * -0.16666666666666666 * u[0] +
416 rho * -0.16666666666666666 * (u[1] * u[1]) +
417 rho * -0.16666666666666666 * (u[2] * u[2]) +
418 rho * 0.055555555555555552 +
419 rho * 0.16666666666666666 * (u[0] * u[0]);
420 case E:
421 return rho * -0.16666666666666666 * (u[1] * u[1]) +
422 rho * -0.16666666666666666 * (u[2] * u[2]) +
423 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[0] +
424 rho * 0.16666666666666666 * (u[0] * u[0]);
425 case T:
426 return rho * -0.16666666666666666 * (u[0] * u[0]) +
427 rho * -0.16666666666666666 * (u[1] * u[1]) +
428 rho * 0.055555555555555552 + rho * 0.16666666666666666 * u[2] +
429 rho * 0.16666666666666666 * (u[2] * u[2]);
430 case B:
431 return rho * -0.16666666666666666 * u[2] +
432 rho * -0.16666666666666666 * (u[0] * u[0]) +
433 rho * -0.16666666666666666 * (u[1] * u[1]) +
434 rho * 0.055555555555555552 +
435 rho * 0.16666666666666666 * (u[2] * u[2]);
436 case NW:
437 return rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[1] +
438 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
439 rho * 0.083333333333333329 * (u[0] * u[0]) +
440 rho * 0.083333333333333329 * (u[1] * u[1]);
441 case NE:
442 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
443 rho * 0.083333333333333329 * u[1] +
444 rho * 0.083333333333333329 * (u[0] * u[0]) +
445 rho * 0.083333333333333329 * (u[1] * u[1]) +
446 rho * 0.25 * u[0] * u[1];
447 case SW:
448 return rho * -0.083333333333333329 * u[0] +
449 rho * -0.083333333333333329 * u[1] + rho * 0.027777777777777776 +
450 rho * 0.083333333333333329 * (u[0] * u[0]) +
451 rho * 0.083333333333333329 * (u[1] * u[1]) +
452 rho * 0.25 * u[0] * u[1];
453 case SE:
454 return rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[0] * u[1] +
455 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
456 rho * 0.083333333333333329 * (u[0] * u[0]) +
457 rho * 0.083333333333333329 * (u[1] * u[1]);
458 case TN:
459 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
460 rho * 0.083333333333333329 * u[2] +
461 rho * 0.083333333333333329 * (u[1] * u[1]) +
462 rho * 0.083333333333333329 * (u[2] * u[2]) +
463 rho * 0.25 * u[1] * u[2];
464 case TS:
465 return rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[1] * u[2] +
466 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
467 rho * 0.083333333333333329 * (u[1] * u[1]) +
468 rho * 0.083333333333333329 * (u[2] * u[2]);
469 case TW:
470 return rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[2] +
471 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
472 rho * 0.083333333333333329 * (u[0] * u[0]) +
473 rho * 0.083333333333333329 * (u[2] * u[2]);
474 case TE:
475 return rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
476 rho * 0.083333333333333329 * u[2] +
477 rho * 0.083333333333333329 * (u[0] * u[0]) +
478 rho * 0.083333333333333329 * (u[2] * u[2]) +
479 rho * 0.25 * u[0] * u[2];
480 case BN:
481 return rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[1] * u[2] +
482 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
483 rho * 0.083333333333333329 * (u[1] * u[1]) +
484 rho * 0.083333333333333329 * (u[2] * u[2]);
485 case BS:
486 return rho * -0.083333333333333329 * u[1] +
487 rho * -0.083333333333333329 * u[2] + rho * 0.027777777777777776 +
488 rho * 0.083333333333333329 * (u[1] * u[1]) +
489 rho * 0.083333333333333329 * (u[2] * u[2]) +
490 rho * 0.25 * u[1] * u[2];
491 case BW:
492 return rho * -0.083333333333333329 * u[0] +
493 rho * -0.083333333333333329 * u[2] + rho * 0.027777777777777776 +
494 rho * 0.083333333333333329 * (u[0] * u[0]) +
495 rho * 0.083333333333333329 * (u[2] * u[2]) +
496 rho * 0.25 * u[0] * u[2];
497 case BE:
498 return rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[0] * u[2] +
499 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
500 rho * 0.083333333333333329 * (u[0] * u[0]) +
501 rho * 0.083333333333333329 * (u[2] * u[2]);
502 default:
503 WALBERLA_ABORT("Invalid Direction")
504 }
505}
506} // namespace EquilibriumDistribution
507
508namespace Equilibrium {
509inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
510 Vector3<double> const &u, double const rho, Cell const &cell) {
511
512 double &xyz0 = pdf_field->get(cell, uint_t{0u});
513 pdf_field->getF(&xyz0, uint_t{0u}) =
514 rho * -0.33333333333333331 * (u[0] * u[0]) +
515 rho * -0.33333333333333331 * (u[1] * u[1]) +
516 rho * -0.33333333333333331 * (u[2] * u[2]) + rho * 0.33333333333333331;
517 pdf_field->getF(&xyz0, uint_t{1u}) =
518 rho * -0.16666666666666666 * (u[0] * u[0]) +
519 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
520 rho * 0.16666666666666666 * u[1] +
521 rho * 0.16666666666666666 * (u[1] * u[1]);
522 pdf_field->getF(&xyz0, uint_t{2u}) =
523 rho * -0.16666666666666666 * u[1] +
524 rho * -0.16666666666666666 * (u[0] * u[0]) +
525 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
526 rho * 0.16666666666666666 * (u[1] * u[1]);
527 pdf_field->getF(&xyz0, uint_t{3u}) =
528 rho * -0.16666666666666666 * u[0] +
529 rho * -0.16666666666666666 * (u[1] * u[1]) +
530 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
531 rho * 0.16666666666666666 * (u[0] * u[0]);
532 pdf_field->getF(&xyz0, uint_t{4u}) =
533 rho * -0.16666666666666666 * (u[1] * u[1]) +
534 rho * -0.16666666666666666 * (u[2] * u[2]) + rho * 0.055555555555555552 +
535 rho * 0.16666666666666666 * u[0] +
536 rho * 0.16666666666666666 * (u[0] * u[0]);
537 pdf_field->getF(&xyz0, uint_t{5u}) =
538 rho * -0.16666666666666666 * (u[0] * u[0]) +
539 rho * -0.16666666666666666 * (u[1] * u[1]) + rho * 0.055555555555555552 +
540 rho * 0.16666666666666666 * u[2] +
541 rho * 0.16666666666666666 * (u[2] * u[2]);
542 pdf_field->getF(&xyz0, uint_t{6u}) =
543 rho * -0.16666666666666666 * u[2] +
544 rho * -0.16666666666666666 * (u[0] * u[0]) +
545 rho * -0.16666666666666666 * (u[1] * u[1]) + rho * 0.055555555555555552 +
546 rho * 0.16666666666666666 * (u[2] * u[2]);
547 pdf_field->getF(&xyz0, uint_t{7u}) =
548 rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[1] +
549 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
550 rho * 0.083333333333333329 * (u[0] * u[0]) +
551 rho * 0.083333333333333329 * (u[1] * u[1]);
552 pdf_field->getF(&xyz0, uint_t{8u}) =
553 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
554 rho * 0.083333333333333329 * u[1] +
555 rho * 0.083333333333333329 * (u[0] * u[0]) +
556 rho * 0.083333333333333329 * (u[1] * u[1]) + rho * 0.25 * u[0] * u[1];
557 pdf_field->getF(&xyz0, uint_t{9u}) =
558 rho * -0.083333333333333329 * u[0] + rho * -0.083333333333333329 * u[1] +
559 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[0] * u[0]) +
560 rho * 0.083333333333333329 * (u[1] * u[1]) + rho * 0.25 * u[0] * u[1];
561 pdf_field->getF(&xyz0, uint_t{10u}) =
562 rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[0] * u[1] +
563 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
564 rho * 0.083333333333333329 * (u[0] * u[0]) +
565 rho * 0.083333333333333329 * (u[1] * u[1]);
566 pdf_field->getF(&xyz0, uint_t{11u}) =
567 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
568 rho * 0.083333333333333329 * u[2] +
569 rho * 0.083333333333333329 * (u[1] * u[1]) +
570 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[1] * u[2];
571 pdf_field->getF(&xyz0, uint_t{12u}) =
572 rho * -0.083333333333333329 * u[1] + rho * -0.25 * u[1] * u[2] +
573 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
574 rho * 0.083333333333333329 * (u[1] * u[1]) +
575 rho * 0.083333333333333329 * (u[2] * u[2]);
576 pdf_field->getF(&xyz0, uint_t{13u}) =
577 rho * -0.083333333333333329 * u[0] + rho * -0.25 * u[0] * u[2] +
578 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[2] +
579 rho * 0.083333333333333329 * (u[0] * u[0]) +
580 rho * 0.083333333333333329 * (u[2] * u[2]);
581 pdf_field->getF(&xyz0, uint_t{14u}) =
582 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
583 rho * 0.083333333333333329 * u[2] +
584 rho * 0.083333333333333329 * (u[0] * u[0]) +
585 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[0] * u[2];
586 pdf_field->getF(&xyz0, uint_t{15u}) =
587 rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[1] * u[2] +
588 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[1] +
589 rho * 0.083333333333333329 * (u[1] * u[1]) +
590 rho * 0.083333333333333329 * (u[2] * u[2]);
591 pdf_field->getF(&xyz0, uint_t{16u}) =
592 rho * -0.083333333333333329 * u[1] + rho * -0.083333333333333329 * u[2] +
593 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[1] * u[1]) +
594 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[1] * u[2];
595 pdf_field->getF(&xyz0, uint_t{17u}) =
596 rho * -0.083333333333333329 * u[0] + rho * -0.083333333333333329 * u[2] +
597 rho * 0.027777777777777776 + rho * 0.083333333333333329 * (u[0] * u[0]) +
598 rho * 0.083333333333333329 * (u[2] * u[2]) + rho * 0.25 * u[0] * u[2];
599 pdf_field->getF(&xyz0, uint_t{18u}) =
600 rho * -0.083333333333333329 * u[2] + rho * -0.25 * u[0] * u[2] +
601 rho * 0.027777777777777776 + rho * 0.083333333333333329 * u[0] +
602 rho * 0.083333333333333329 * (u[0] * u[0]) +
603 rho * 0.083333333333333329 * (u[2] * u[2]);
604}
605} // namespace Equilibrium
606
607namespace Density {
608inline double get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
609 Cell const &cell) {
610 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
611 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
612 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
613 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
614 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
615 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
616 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
617 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
618 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
619 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
620 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
621 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
622 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
623 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
624 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
625 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
626 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
627 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
628 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
629 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
630 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
631 const double vel1Term = f_1 + f_11 + f_15 + f_7;
632 const double vel2Term = f_12 + f_13 + f_5;
633 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
634 vel1Term + vel2Term;
635 return rho;
636}
637
638inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
639 double const rho_in, Cell const &cell) {
640 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
641 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
642 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
643 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
644 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
645 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
646 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
647 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
648 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
649 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
650 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
651 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
652 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
653 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
654 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
655 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
656 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
657 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
658 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
659 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
660 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
661 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
662 const double vel1Term = f_1 + f_11 + f_15 + f_7;
663 const double momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
664 const double vel2Term = f_12 + f_13 + f_5;
665 const double momdensity_2 =
666 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
667 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
668 vel1Term + vel2Term;
669
670 // calculate current velocity (before density change)
671 const double conversion = double{1} / rho;
672 Vector3<double> velocity;
673 velocity[0u] = momdensity_0 * conversion;
674 velocity[1u] = momdensity_1 * conversion;
675 velocity[2u] = momdensity_2 * conversion;
676
677 Equilibrium::set(pdf_field, velocity, rho_in, cell);
678}
679
680inline std::vector<double>
681get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
682 CellInterval const &ci) {
683 std::vector<double> out;
684 out.reserve(ci.numCells());
685 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
686 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
687 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
688 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
689 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
690 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
691 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
692 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
693 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
694 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
695 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
696 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
697 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
698 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
699 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
700 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
701 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
702 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
703 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
704 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
705 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
706 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
707 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
708 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
709 const double vel1Term = f_1 + f_11 + f_15 + f_7;
710 const double vel2Term = f_12 + f_13 + f_5;
711 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
712 vel0Term + vel1Term + vel2Term;
713 out.emplace_back(rho);
714 }
715 }
716 }
717 return out;
718}
719
720inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
721 std::vector<double> const &values, CellInterval const &ci) {
722 assert(uint_c(values.size()) == ci.numCells());
723 auto values_it = values.begin();
724 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
725 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
726 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
727 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
728 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
729 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
730 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
731 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
732 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
733 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
734 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
735 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
736 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
737 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
738 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
739 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
740 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
741 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
742 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
743 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
744 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
745 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
746 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
747 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
748 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
749 const double vel1Term = f_1 + f_11 + f_15 + f_7;
750 const double momdensity_1 =
751 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
752 const double vel2Term = f_12 + f_13 + f_5;
753 const double momdensity_2 =
754 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
755 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
756 vel0Term + vel1Term + vel2Term;
757
758 // calculate current velocity (before density change)
759 const double conversion = double{1} / rho;
760 Vector3<double> velocity;
761 velocity[0u] = momdensity_0 * conversion;
762 velocity[1u] = momdensity_1 * conversion;
763 velocity[2u] = momdensity_2 * conversion;
764
765 Equilibrium::set(pdf_field, velocity, *values_it, Cell{x, y, z});
766 ++values_it;
767 }
768 }
769 }
770}
771} // namespace Density
772
773namespace Velocity {
774inline auto get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
775 GhostLayerField<double, uint_t{3u}> const *force_field,
776 Cell const &cell) {
777 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
778 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
779 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
780 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
781 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
782 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
783 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
784 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
785 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
786 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
787 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
788 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
789 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
790 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
791 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
792 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
793 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
794 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
795 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
796 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
797 const auto x = cell.x();
798 const auto y = cell.y();
799 const auto z = cell.z();
800 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
801 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
802 const double vel1Term = f_1 + f_11 + f_15 + f_7;
803 const double momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
804 const double vel2Term = f_12 + f_13 + f_5;
805 const double momdensity_2 =
806 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
807 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
808 vel1Term + vel2Term;
809 const double md_0 =
810 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
811 const double md_1 =
812 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
813 const double md_2 =
814 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
815 const double rho_inv = double{1} / rho;
816
817 return Vector3<double>(md_0 * rho_inv, md_1 * rho_inv, md_2 * rho_inv);
818}
819
820inline auto get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
821 GhostLayerField<double, uint_t{3u}> const *force_field,
822 CellInterval const &ci) {
823 std::vector<double> out;
824 out.reserve(ci.numCells() * uint_t(3u));
825 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
826 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
827 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
828 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
829 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
830 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
831 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
832 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
833 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
834 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
835 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
836 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
837 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
838 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
839 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
840 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
841 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
842 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
843 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
844 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
845 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
846 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
847 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
848 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
849 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
850 const double vel1Term = f_1 + f_11 + f_15 + f_7;
851 const double momdensity_1 =
852 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
853 const double vel2Term = f_12 + f_13 + f_5;
854 const double momdensity_2 =
855 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
856 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
857 vel0Term + vel1Term + vel2Term;
858 const double md_0 =
859 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
860 const double md_1 =
861 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
862 const double md_2 =
863 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
864 const double rho_inv = double{1} / rho;
865 out.emplace_back(md_0 * rho_inv);
866 out.emplace_back(md_1 * rho_inv);
867 out.emplace_back(md_2 * rho_inv);
868 }
869 }
870 }
871 return out;
872}
873
874inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
875 GhostLayerField<double, uint_t{3u}> *velocity_field,
876 GhostLayerField<double, uint_t{3u}> const *force_field,
877 Vector3<double> const &u, Cell const &cell) {
878 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
879 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
880 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
881 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
882 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
883 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
884 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
885 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
886 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
887 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
888 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
889 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
890 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
891 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
892 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
893 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
894 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
895 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
896 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
897 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
898 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
899 const double vel1Term = f_1 + f_11 + f_15 + f_7;
900 const double vel2Term = f_12 + f_13 + f_5;
901 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
902 vel1Term + vel2Term;
903
904 const auto x = cell.x();
905 const auto y = cell.y();
906 const auto z = cell.z();
907 const double u_0 =
908 -force_field->get(x, y, z, 0) * 0.50000000000000000 / rho + u[0];
909 const double u_1 =
910 -force_field->get(x, y, z, 1) * 0.50000000000000000 / rho + u[1];
911 const double u_2 =
912 -force_field->get(x, y, z, 2) * 0.50000000000000000 / rho + u[2];
913 velocity_field->get(x, y, z, uint_t{0u}) = u[0u];
914 velocity_field->get(x, y, z, uint_t{1u}) = u[1u];
915 velocity_field->get(x, y, z, uint_t{2u}) = u[2u];
916
917 Equilibrium::set(pdf_field, Vector3<double>(u_0, u_1, u_2), rho, cell);
918}
919
920inline void set(GhostLayerField<double, uint_t{19u}> *pdf_field,
921 GhostLayerField<double, uint_t{3u}> *velocity_field,
922 GhostLayerField<double, uint_t{3u}> const *force_field,
923 std::vector<double> const &values, CellInterval const &ci) {
924 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
925 auto u = values.data();
926 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
927 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
928 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
929 double &pdf_xyz0 = pdf_field->get(x, y, z, uint_t{0u});
930 double &vel_xyz0 = velocity_field->get(x, y, z, uint_t{0u});
931 const double f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
932 const double f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
933 const double f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
934 const double f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
935 const double f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
936 const double f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
937 const double f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
938 const double f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
939 const double f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
940 const double f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
941 const double f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
942 const double f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
943 const double f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
944 const double f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
945 const double f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
946 const double f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
947 const double f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
948 const double f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
949 const double f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
950 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
951 const double vel1Term = f_1 + f_11 + f_15 + f_7;
952 const double vel2Term = f_12 + f_13 + f_5;
953 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
954 vel0Term + vel1Term + vel2Term;
955
956 const double u_0 =
957 -force_field->get(x, y, z, 0) * 0.50000000000000000 / rho + u[0];
958 const double u_1 =
959 -force_field->get(x, y, z, 1) * 0.50000000000000000 / rho + u[1];
960 const double u_2 =
961 -force_field->get(x, y, z, 2) * 0.50000000000000000 / rho + u[2];
962 velocity_field->getF(&vel_xyz0, uint_t{0u}) = u[0u];
963 velocity_field->getF(&vel_xyz0, uint_t{1u}) = u[1u];
964 velocity_field->getF(&vel_xyz0, uint_t{2u}) = u[2u];
965
966 std::advance(u, 3);
967
968 Equilibrium::set(pdf_field, Vector3<double>(u_0, u_1, u_2), rho,
969 Cell{x, y, z});
970 }
971 }
972 }
973}
974} // namespace Velocity
975
976namespace Force {
977inline void set(GhostLayerField<double, uint_t{19u}> const *pdf_field,
978 GhostLayerField<double, uint_t{3u}> *velocity_field,
979 GhostLayerField<double, uint_t{3u}> *force_field,
980 Vector3<double> const &force, Cell const &cell) {
981 double const &pdf_xyz0 = pdf_field->get(cell, uint_t{0u});
982 double &vel_xyz0 = velocity_field->get(cell, uint_t{0u});
983 double &laf_xyz0 = force_field->get(cell, uint_t{0u});
984 const double f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
985 const double f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
986 const double f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
987 const double f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
988 const double f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
989 const double f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
990 const double f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
991 const double f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
992 const double f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
993 const double f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
994 const double f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
995 const double f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
996 const double f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
997 const double f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
998 const double f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
999 const double f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
1000 const double f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
1001 const double f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
1002 const double f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
1003 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1004 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1005 const double vel1Term = f_1 + f_11 + f_15 + f_7;
1006 const double momdensity_1 = -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1007 const double vel2Term = f_12 + f_13 + f_5;
1008 const double momdensity_2 =
1009 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1010 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
1011 vel1Term + vel2Term;
1012 const double md_0 = force[0u] * 0.50000000000000000 + momdensity_0;
1013 const double md_1 = force[1u] * 0.50000000000000000 + momdensity_1;
1014 const double md_2 = force[2u] * 0.50000000000000000 + momdensity_2;
1015 auto const rho_inv = double{1} / rho;
1016
1017 force_field->getF(&laf_xyz0, uint_t{0u}) = force[0u];
1018 force_field->getF(&laf_xyz0, uint_t{1u}) = force[1u];
1019 force_field->getF(&laf_xyz0, uint_t{2u}) = force[2u];
1020
1021 velocity_field->getF(&vel_xyz0, uint_t{0u}) = md_0 * rho_inv;
1022 velocity_field->getF(&vel_xyz0, uint_t{1u}) = md_1 * rho_inv;
1023 velocity_field->getF(&vel_xyz0, uint_t{2u}) = md_2 * rho_inv;
1024}
1025
1026inline void set(GhostLayerField<double, uint_t{19u}> const *pdf_field,
1027 GhostLayerField<double, uint_t{3u}> *velocity_field,
1028 GhostLayerField<double, uint_t{3u}> *force_field,
1029 std::vector<double> const &values, CellInterval const &ci) {
1030 assert(uint_c(values.size()) == ci.numCells() * uint_t(3u));
1031 auto force = values.data();
1032 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
1033 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
1034 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
1035 double const &pdf_xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1036 double &vel_xyz0 = velocity_field->get(x, y, z, uint_t{0u});
1037 double &laf_xyz0 = force_field->get(x, y, z, uint_t{0u});
1038 const double f_0 = pdf_field->getF(&pdf_xyz0, uint_t{0u});
1039 const double f_1 = pdf_field->getF(&pdf_xyz0, uint_t{1u});
1040 const double f_2 = pdf_field->getF(&pdf_xyz0, uint_t{2u});
1041 const double f_3 = pdf_field->getF(&pdf_xyz0, uint_t{3u});
1042 const double f_4 = pdf_field->getF(&pdf_xyz0, uint_t{4u});
1043 const double f_5 = pdf_field->getF(&pdf_xyz0, uint_t{5u});
1044 const double f_6 = pdf_field->getF(&pdf_xyz0, uint_t{6u});
1045 const double f_7 = pdf_field->getF(&pdf_xyz0, uint_t{7u});
1046 const double f_8 = pdf_field->getF(&pdf_xyz0, uint_t{8u});
1047 const double f_9 = pdf_field->getF(&pdf_xyz0, uint_t{9u});
1048 const double f_10 = pdf_field->getF(&pdf_xyz0, uint_t{10u});
1049 const double f_11 = pdf_field->getF(&pdf_xyz0, uint_t{11u});
1050 const double f_12 = pdf_field->getF(&pdf_xyz0, uint_t{12u});
1051 const double f_13 = pdf_field->getF(&pdf_xyz0, uint_t{13u});
1052 const double f_14 = pdf_field->getF(&pdf_xyz0, uint_t{14u});
1053 const double f_15 = pdf_field->getF(&pdf_xyz0, uint_t{15u});
1054 const double f_16 = pdf_field->getF(&pdf_xyz0, uint_t{16u});
1055 const double f_17 = pdf_field->getF(&pdf_xyz0, uint_t{17u});
1056 const double f_18 = pdf_field->getF(&pdf_xyz0, uint_t{18u});
1057 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1058 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1059 const double vel1Term = f_1 + f_11 + f_15 + f_7;
1060 const double momdensity_1 =
1061 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1062 const double vel2Term = f_12 + f_13 + f_5;
1063 const double momdensity_2 =
1064 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1065 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 +
1066 vel0Term + vel1Term + vel2Term;
1067 const double md_0 = force[0u] * 0.50000000000000000 + momdensity_0;
1068 const double md_1 = force[1u] * 0.50000000000000000 + momdensity_1;
1069 const double md_2 = force[2u] * 0.50000000000000000 + momdensity_2;
1070 auto const rho_inv = double{1} / rho;
1071
1072 force_field->getF(&laf_xyz0, uint_t{0u}) = force[0u];
1073 force_field->getF(&laf_xyz0, uint_t{1u}) = force[1u];
1074 force_field->getF(&laf_xyz0, uint_t{2u}) = force[2u];
1075
1076 velocity_field->getF(&vel_xyz0, uint_t{0u}) = md_0 * rho_inv;
1077 velocity_field->getF(&vel_xyz0, uint_t{1u}) = md_1 * rho_inv;
1078 velocity_field->getF(&vel_xyz0, uint_t{2u}) = md_2 * rho_inv;
1079
1080 std::advance(force, 3);
1081 }
1082 }
1083 }
1084}
1085} // namespace Force
1086
1087namespace MomentumDensity {
1088inline auto reduce(GhostLayerField<double, uint_t{19u}> const *pdf_field,
1089 GhostLayerField<double, uint_t{3u}> const *force_field) {
1090 Vector3<double> momentumDensity(double{0});
1091 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
1092 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1093 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1094 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1095 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1096 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1097 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1098 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1099 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1100 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1101 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1102 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1103 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1104 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1105 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1106 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1107 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1108 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1109 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1110 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1111 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1112 const double vel0Term = f_10 + f_14 + f_18 + f_4 + f_8;
1113 const double momdensity_0 = -f_13 - f_17 - f_3 - f_7 - f_9 + vel0Term;
1114 const double vel1Term = f_1 + f_11 + f_15 + f_7;
1115 const double momdensity_1 =
1116 -f_10 - f_12 - f_16 - f_2 + f_8 - f_9 + vel1Term;
1117 const double vel2Term = f_12 + f_13 + f_5;
1118 const double momdensity_2 =
1119 f_11 + f_14 - f_15 - f_16 - f_17 - f_18 - f_6 + vel2Term;
1120 const double rho = f_0 + f_16 + f_17 + f_2 + f_3 + f_6 + f_9 + vel0Term +
1121 vel1Term + vel2Term;
1122 const double md_0 =
1123 force_field->get(x, y, z, 0) * 0.50000000000000000 + momdensity_0;
1124 const double md_1 =
1125 force_field->get(x, y, z, 1) * 0.50000000000000000 + momdensity_1;
1126 const double md_2 =
1127 force_field->get(x, y, z, 2) * 0.50000000000000000 + momdensity_2;
1128
1129 momentumDensity[0u] += md_0;
1130 momentumDensity[1u] += md_1;
1131 momentumDensity[2u] += md_2;
1132 });
1133 return momentumDensity;
1134}
1135} // namespace MomentumDensity
1136
1137namespace PressureTensor {
1138inline auto get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
1139 Cell const &cell) {
1140 const double &xyz0 = pdf_field->get(cell, uint_t{0u});
1141 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1142 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1143 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1144 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1145 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1146 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1147 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1148 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1149 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1150 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1151 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1152 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1153 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1154 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1155 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1156 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1157 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1158 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1159 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1160 const double p_0 =
1161 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1162 const double p_1 = -f_10 - f_7 + f_8 + f_9;
1163 const double p_2 = -f_13 + f_14 + f_17 - f_18;
1164 const double p_3 = -f_10 - f_7 + f_8 + f_9;
1165 const double p_4 =
1166 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1167 const double p_5 = f_11 - f_12 - f_15 + f_16;
1168 const double p_6 = -f_13 + f_14 + f_17 - f_18;
1169 const double p_7 = f_11 - f_12 - f_15 + f_16;
1170 const double p_8 =
1171 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1172
1173 Matrix3<double> pressureTensor;
1174 pressureTensor[0u] = p_0;
1175 pressureTensor[1u] = p_1;
1176 pressureTensor[2u] = p_2;
1177
1178 pressureTensor[3u] = p_3;
1179 pressureTensor[4u] = p_4;
1180 pressureTensor[5u] = p_5;
1181
1182 pressureTensor[6u] = p_6;
1183 pressureTensor[7u] = p_7;
1184 pressureTensor[8u] = p_8;
1185
1186 return pressureTensor;
1187}
1188
1189inline auto get(GhostLayerField<double, uint_t{19u}> const *pdf_field,
1190 CellInterval const &ci) {
1191 std::vector<double> out;
1192 out.reserve(ci.numCells() * uint_t(9u));
1193 for (auto x = ci.xMin(); x <= ci.xMax(); ++x) {
1194 for (auto y = ci.yMin(); y <= ci.yMax(); ++y) {
1195 for (auto z = ci.zMin(); z <= ci.zMax(); ++z) {
1196 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1197 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1198 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1199 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1200 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1201 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1202 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1203 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1204 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1205 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1206 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1207 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1208 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1209 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1210 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1211 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1212 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1213 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1214 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1215 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1216 const double p_0 =
1217 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1218 const double p_1 = -f_10 - f_7 + f_8 + f_9;
1219 const double p_2 = -f_13 + f_14 + f_17 - f_18;
1220 const double p_3 = -f_10 - f_7 + f_8 + f_9;
1221 const double p_4 =
1222 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1223 const double p_5 = f_11 - f_12 - f_15 + f_16;
1224 const double p_6 = -f_13 + f_14 + f_17 - f_18;
1225 const double p_7 = f_11 - f_12 - f_15 + f_16;
1226 const double p_8 =
1227 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1228
1229 out.emplace_back(p_0);
1230 out.emplace_back(p_1);
1231 out.emplace_back(p_2);
1232
1233 out.emplace_back(p_3);
1234 out.emplace_back(p_4);
1235 out.emplace_back(p_5);
1236
1237 out.emplace_back(p_6);
1238 out.emplace_back(p_7);
1239 out.emplace_back(p_8);
1240 }
1241 }
1242 }
1243 return out;
1244}
1245
1246inline auto reduce(GhostLayerField<double, uint_t{19u}> const *pdf_field) {
1247 Matrix3<double> pressureTensor(double{0});
1248 WALBERLA_FOR_ALL_CELLS_XYZ(pdf_field, {
1249 const double &xyz0 = pdf_field->get(x, y, z, uint_t{0u});
1250 const double f_0 = pdf_field->getF(&xyz0, uint_t{0u});
1251 const double f_1 = pdf_field->getF(&xyz0, uint_t{1u});
1252 const double f_2 = pdf_field->getF(&xyz0, uint_t{2u});
1253 const double f_3 = pdf_field->getF(&xyz0, uint_t{3u});
1254 const double f_4 = pdf_field->getF(&xyz0, uint_t{4u});
1255 const double f_5 = pdf_field->getF(&xyz0, uint_t{5u});
1256 const double f_6 = pdf_field->getF(&xyz0, uint_t{6u});
1257 const double f_7 = pdf_field->getF(&xyz0, uint_t{7u});
1258 const double f_8 = pdf_field->getF(&xyz0, uint_t{8u});
1259 const double f_9 = pdf_field->getF(&xyz0, uint_t{9u});
1260 const double f_10 = pdf_field->getF(&xyz0, uint_t{10u});
1261 const double f_11 = pdf_field->getF(&xyz0, uint_t{11u});
1262 const double f_12 = pdf_field->getF(&xyz0, uint_t{12u});
1263 const double f_13 = pdf_field->getF(&xyz0, uint_t{13u});
1264 const double f_14 = pdf_field->getF(&xyz0, uint_t{14u});
1265 const double f_15 = pdf_field->getF(&xyz0, uint_t{15u});
1266 const double f_16 = pdf_field->getF(&xyz0, uint_t{16u});
1267 const double f_17 = pdf_field->getF(&xyz0, uint_t{17u});
1268 const double f_18 = pdf_field->getF(&xyz0, uint_t{18u});
1269 const double p_0 =
1270 f_10 + f_13 + f_14 + f_17 + f_18 + f_3 + f_4 + f_7 + f_8 + f_9;
1271 const double p_1 = -f_10 - f_7 + f_8 + f_9;
1272 const double p_2 = -f_13 + f_14 + f_17 - f_18;
1273 const double p_3 = -f_10 - f_7 + f_8 + f_9;
1274 const double p_4 =
1275 f_1 + f_10 + f_11 + f_12 + f_15 + f_16 + f_2 + f_7 + f_8 + f_9;
1276 const double p_5 = f_11 - f_12 - f_15 + f_16;
1277 const double p_6 = -f_13 + f_14 + f_17 - f_18;
1278 const double p_7 = f_11 - f_12 - f_15 + f_16;
1279 const double p_8 =
1280 f_11 + f_12 + f_13 + f_14 + f_15 + f_16 + f_17 + f_18 + f_5 + f_6;
1281
1282 pressureTensor[0u] += p_0;
1283 pressureTensor[1u] += p_1;
1284 pressureTensor[2u] += p_2;
1285
1286 pressureTensor[3u] += p_3;
1287 pressureTensor[4u] += p_4;
1288 pressureTensor[5u] += p_5;
1289
1290 pressureTensor[6u] += p_6;
1291 pressureTensor[7u] += p_7;
1292 pressureTensor[8u] += p_8;
1293 });
1294 return pressureTensor;
1295}
1296} // namespace PressureTensor
1297
1298} // namespace accessor
1299} // namespace lbm
1300} // namespace walberla
1301
1302#ifdef WALBERLA_CXX_COMPILER_IS_GNU
1303#pragma GCC diagnostic pop
1304#endif
1305
1306#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
1307#pragma clang diagnostic pop
1308#endif
Cell
Definition Cell.hpp:97
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