Loading [MathJax]/extensions/TeX/AMSmath.js
ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages Concepts
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.7, lbmpy v1.3.7, sympy v1.12.1,
22// lbmpy_walberla/pystencils_walberla from waLBerla commit
23// f36fa0a68bae59f0b516f6587ea8fa7c24a41141
24
25/*
26 * Lattice field accessors.
27 * Adapted from the waLBerla source file
28 * https://i10git.cs.fau.de/walberla/walberla/-/blob/a16141524c58ab88386e2a0f8fdd7c63c5edd704/python/lbmpy_walberla/templates/LatticeModel.tmpl.h
29 */
30
31#pragma once
32
33#include <core/DataTypes.h>
34#include <core/cell/Cell.h>
35#include <core/cell/CellInterval.h>
36#include <core/math/Matrix3.h>
37#include <core/math/Vector3.h>
38
39#include <field/GhostLayerField.h>
40#include <stencil/D3Q19.h>
41
42#include <array>
43#include <cassert>
44#include <iterator>
45#include <tuple>
46#include <vector>
47
48#ifdef WALBERLA_CXX_COMPILER_IS_GNU
49#pragma GCC diagnostic push
50#pragma GCC diagnostic ignored "-Wunused-variable"
51#endif
52
53#ifdef WALBERLA_CXX_COMPILER_IS_CLANG
54#pragma clang diagnostic push
55#pragma clang diagnostic ignored "-Wunused-variable"
56#endif
57
58namespace walberla {
59namespace lbm {
60namespace accessor {
61
62namespace Population {
63inline auto get(GhostLayerField<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:96
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