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