ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
DynamicUBBDoublePrecision.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2022-2025 The ESPResSo project
3 * Copyright (C) 2020-2025 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+13.gdfd203a, lbmpy
22// v1.3.7+15.g5018a18, sympy v1.12.1, lbmpy_walberla/pystencils_walberla from
23// waLBerla commit 191cf58b16b96d1d2f050dcbd9e88443995b2222
24
25/*
26 * Boundary class.
27 * Adapted from the waLBerla source file
28 * https://i10git.cs.fau.de/walberla/walberla/-/blob/3e54d4f2336e47168ad87e3caaf7b3b082d86ca7/python/pystencils_walberla/templates/Boundary.tmpl.h
29 */
30
31#pragma once
32
33#include <core/DataTypes.h>
34
35#include <blockforest/StructuredBlockForest.h>
36#include <core/debug/Debug.h>
37#include <domain_decomposition/BlockDataID.h>
38#include <domain_decomposition/IBlock.h>
39#include <field/FlagField.h>
40#include <field/GhostLayerField.h>
41
42#include <array>
43#include <cassert>
44#include <functional>
45#include <memory>
46#include <vector>
47
48#if defined(__clang__)
49#pragma clang diagnostic push
50#pragma clang diagnostic ignored "-Wunused-variable"
51#pragma clang diagnostic ignored "-Wunused-parameter"
52#elif defined(__GNUC__) or defined(__GNUG__)
53#pragma GCC diagnostic push
54#pragma GCC diagnostic ignored "-Wunused-variable"
55#pragma GCC diagnostic ignored "-Wunused-parameter"
56#endif
57
58#ifdef __GNUC__
59#define RESTRICT __restrict__
60#elif _MSC_VER
61#define RESTRICT __restrict
62#else
63#define RESTRICT
64#endif
65
66#ifdef WALBERLA_BUILD_WITH_HALF_PRECISION_SUPPORT
67using walberla::half;
68#endif
69
70namespace walberla {
71namespace lbm {
72
73class DynamicUBBDoublePrecision {
74public:
75 struct IndexInfo {
76 int32_t x;
77 int32_t y;
78 int32_t z;
79 int32_t dir;
80 double vel_0;
81 double vel_1;
82 double vel_2;
83 IndexInfo(int32_t x_, int32_t y_, int32_t z_, int32_t dir_)
84 : x(x_), y(y_), z(z_), dir(dir_), vel_0(), vel_1(), vel_2() {}
85 bool operator==(const IndexInfo &o) const {
86 return x == o.x && y == o.y && z == o.z && dir == o.dir &&
87 floatIsEqual(vel_0, o.vel_0) && floatIsEqual(vel_1, o.vel_1) &&
88 floatIsEqual(vel_2, o.vel_2);
89 }
90 };
91
92 class IndexVectors {
93 public:
94 using CpuIndexVector = std::vector<IndexInfo>;
95
96 enum Type { ALL = 0, INNER = 1, OUTER = 2, NUM_TYPES = 3 };
97
98 IndexVectors() = default;
99 bool operator==(IndexVectors const &other) const {
100 return other.cpuVectors_ == cpuVectors_;
101 }
102
103 auto &indexVector(Type t) { return cpuVectors_[t]; }
104 auto const &indexVector(Type t) const { return cpuVectors_[t]; }
105 IndexInfo *pointerCpu(Type t) {
106 return cpuVectors_[t].empty() ? nullptr : cpuVectors_[t].data();
107 }
108
109 void syncGPU() {}
110
111 private:
112 std::vector<CpuIndexVector> cpuVectors_{NUM_TYPES};
113 };
114
115 struct ForceStruct {
116 double F_0;
117 double F_1;
118 double F_2;
119 ForceStruct()
120 : F_0(double_c(0.0)), F_1(double_c(0.0)), F_2(double_c(0.0)) {}
121 bool operator==(const ForceStruct &o) const {
122 return floatIsEqual(F_0, o.F_0) && floatIsEqual(F_1, o.F_1) &&
123 floatIsEqual(F_2, o.F_2);
124 }
125 };
126
127 class ForceVector {
128 public:
129 ForceVector() = default;
130 bool operator==(ForceVector const &other) const {
131 return other.cpuVector_ == cpuVector_;
132 }
133
134 auto &forceVector() { return cpuVector_; }
135 auto const &forceVector() const { return cpuVector_; }
136 ForceStruct *pointerCpu() {
137 return cpuVector_.empty() ? nullptr : cpuVector_.data();
138 }
139 bool empty() const { return cpuVector_.empty(); }
140
141 Vector3<double> getForce() {
142 syncCPU();
143 Vector3<double> result(double_c(0.0));
144 for (auto const &force : cpuVector_) {
145 result[0] += force.F_0;
146 result[1] += force.F_1;
147 result[2] += force.F_2;
148 }
149 return result;
150 }
151
152 void syncGPU() {}
153
154 void syncCPU() {}
155
156 private:
157 std::vector<ForceStruct> cpuVector_;
158 };
159
160 DynamicUBBDoublePrecision(
161 const std::shared_ptr<StructuredBlockForest> &blocks, BlockDataID pdfsID_,
162 std::function<Vector3<double>(
163 const Cell &, const shared_ptr<StructuredBlockForest> &, IBlock &)>
164 &velocityCallbackDynamicUBBDoublePrecision)
165 : elementInitialiser(velocityCallbackDynamicUBBDoublePrecision),
166 pdfsID(pdfsID_) {
167 auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) {
168 return new IndexVectors();
169 };
170 indexVectorID = blocks->addStructuredBlockData<IndexVectors>(
171 createIdxVector, "IndexField_DynamicUBBDoublePrecision");
172 auto createForceVector = [](IBlock *const, StructuredBlockStorage *const) {
173 return new ForceVector();
174 };
175 forceVectorID = blocks->addStructuredBlockData<ForceVector>(
176 createForceVector, "forceVector_DynamicUBBDoublePrecision");
177 }
178
179 void run(IBlock *block);
180
181 void operator()(IBlock *block) { run(block); }
182
183 void inner(IBlock *block);
184
185 void outer(IBlock *block);
186
187 Vector3<double> getForce(IBlock *block) {
188 auto *forceVector = block->getData<ForceVector>(forceVectorID);
189 if (forceVector->empty())
190 return Vector3<double>(double_c(0.0));
191 return forceVector->getForce();
192 }
193
194 std::function<void(IBlock *)> getSweep() {
195 return [this](IBlock *b) { this->run(b); };
196 }
197
198 std::function<void(IBlock *)> getInnerSweep() {
199 return [this](IBlock *b) { this->inner(b); };
200 }
201
202 std::function<void(IBlock *)> getOuterSweep() {
203 return [this](IBlock *b) { this->outer(b); };
204 }
205
206 template <typename FlagField_T>
207 void fillFromFlagField(const std::shared_ptr<StructuredBlockForest> &blocks,
208 ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID,
209 FlagUID domainFlagUID) {
210 for (auto &block : *blocks)
211 fillFromFlagField<FlagField_T>(blocks, &block, flagFieldID,
212 boundaryFlagUID, domainFlagUID);
213 }
214
215 template <typename FlagField_T>
216 void fillFromFlagField(const shared_ptr<StructuredBlockForest> &blocks,
217 IBlock *block, ConstBlockDataID flagFieldID,
218 FlagUID boundaryFlagUID, FlagUID domainFlagUID) {
219 auto *indexVectors = block->getData<IndexVectors>(indexVectorID);
220 auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
221 auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
222 auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
223 auto *forceVector = block->getData<ForceVector>(forceVectorID);
224
225 auto *flagField = block->getData<FlagField_T>(flagFieldID);
226
227 if (!(flagField->flagExists(boundaryFlagUID) and
228 flagField->flagExists(domainFlagUID)))
229 return;
230
231 auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
232 auto domainFlag = flagField->getFlag(domainFlagUID);
233
234 auto inner = flagField->xyzSize();
235 inner.expand(cell_idx_t(-1));
236
237 indexVectorAll.clear();
238 indexVectorInner.clear();
239 indexVectorOuter.clear();
240
241 for (auto it = flagField->beginWithGhostLayerXYZ(
242 cell_idx_c(flagField->nrOfGhostLayers() - 1));
243 it != flagField->end(); ++it) {
244 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
245 continue;
246
247 if (isFlagSet(it.neighbor(0, 0, 0, 0), boundaryFlag)) {
248 auto element = IndexInfo(it.x(), it.y(), it.z(), 0);
249 auto const InitialisationAdditionalData = elementInitialiser(
250 Cell(it.x() + 0, it.y() + 0, it.z() + 0), blocks, *block);
251 element.vel_0 = InitialisationAdditionalData[0];
252 element.vel_1 = InitialisationAdditionalData[1];
253 element.vel_2 = InitialisationAdditionalData[2];
254 indexVectorAll.emplace_back(element);
255 if (inner.contains(it.x(), it.y(), it.z()))
256 indexVectorInner.emplace_back(element);
257 else
258 indexVectorOuter.emplace_back(element);
259 }
260 }
261
262 for (auto it = flagField->beginWithGhostLayerXYZ(
263 cell_idx_c(flagField->nrOfGhostLayers() - 1));
264 it != flagField->end(); ++it) {
265 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
266 continue;
267
268 if (isFlagSet(it.neighbor(0, 1, 0, 0), boundaryFlag)) {
269 auto element = IndexInfo(it.x(), it.y(), it.z(), 1);
270 auto const InitialisationAdditionalData = elementInitialiser(
271 Cell(it.x() + 0, it.y() + 1, it.z() + 0), blocks, *block);
272 element.vel_0 = InitialisationAdditionalData[0];
273 element.vel_1 = InitialisationAdditionalData[1];
274 element.vel_2 = InitialisationAdditionalData[2];
275 indexVectorAll.emplace_back(element);
276 if (inner.contains(it.x(), it.y(), it.z()))
277 indexVectorInner.emplace_back(element);
278 else
279 indexVectorOuter.emplace_back(element);
280 }
281 }
282
283 for (auto it = flagField->beginWithGhostLayerXYZ(
284 cell_idx_c(flagField->nrOfGhostLayers() - 1));
285 it != flagField->end(); ++it) {
286 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
287 continue;
288
289 if (isFlagSet(it.neighbor(0, -1, 0, 0), boundaryFlag)) {
290 auto element = IndexInfo(it.x(), it.y(), it.z(), 2);
291 auto const InitialisationAdditionalData = elementInitialiser(
292 Cell(it.x() + 0, it.y() - 1, it.z() + 0), blocks, *block);
293 element.vel_0 = InitialisationAdditionalData[0];
294 element.vel_1 = InitialisationAdditionalData[1];
295 element.vel_2 = InitialisationAdditionalData[2];
296 indexVectorAll.emplace_back(element);
297 if (inner.contains(it.x(), it.y(), it.z()))
298 indexVectorInner.emplace_back(element);
299 else
300 indexVectorOuter.emplace_back(element);
301 }
302 }
303
304 for (auto it = flagField->beginWithGhostLayerXYZ(
305 cell_idx_c(flagField->nrOfGhostLayers() - 1));
306 it != flagField->end(); ++it) {
307 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
308 continue;
309
310 if (isFlagSet(it.neighbor(-1, 0, 0, 0), boundaryFlag)) {
311 auto element = IndexInfo(it.x(), it.y(), it.z(), 3);
312 auto const InitialisationAdditionalData = elementInitialiser(
313 Cell(it.x() - 1, it.y() + 0, it.z() + 0), blocks, *block);
314 element.vel_0 = InitialisationAdditionalData[0];
315 element.vel_1 = InitialisationAdditionalData[1];
316 element.vel_2 = InitialisationAdditionalData[2];
317 indexVectorAll.emplace_back(element);
318 if (inner.contains(it.x(), it.y(), it.z()))
319 indexVectorInner.emplace_back(element);
320 else
321 indexVectorOuter.emplace_back(element);
322 }
323 }
324
325 for (auto it = flagField->beginWithGhostLayerXYZ(
326 cell_idx_c(flagField->nrOfGhostLayers() - 1));
327 it != flagField->end(); ++it) {
328 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
329 continue;
330
331 if (isFlagSet(it.neighbor(1, 0, 0, 0), boundaryFlag)) {
332 auto element = IndexInfo(it.x(), it.y(), it.z(), 4);
333 auto const InitialisationAdditionalData = elementInitialiser(
334 Cell(it.x() + 1, it.y() + 0, it.z() + 0), blocks, *block);
335 element.vel_0 = InitialisationAdditionalData[0];
336 element.vel_1 = InitialisationAdditionalData[1];
337 element.vel_2 = InitialisationAdditionalData[2];
338 indexVectorAll.emplace_back(element);
339 if (inner.contains(it.x(), it.y(), it.z()))
340 indexVectorInner.emplace_back(element);
341 else
342 indexVectorOuter.emplace_back(element);
343 }
344 }
345
346 for (auto it = flagField->beginWithGhostLayerXYZ(
347 cell_idx_c(flagField->nrOfGhostLayers() - 1));
348 it != flagField->end(); ++it) {
349 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
350 continue;
351
352 if (isFlagSet(it.neighbor(0, 0, 1, 0), boundaryFlag)) {
353 auto element = IndexInfo(it.x(), it.y(), it.z(), 5);
354 auto const InitialisationAdditionalData = elementInitialiser(
355 Cell(it.x() + 0, it.y() + 0, it.z() + 1), blocks, *block);
356 element.vel_0 = InitialisationAdditionalData[0];
357 element.vel_1 = InitialisationAdditionalData[1];
358 element.vel_2 = InitialisationAdditionalData[2];
359 indexVectorAll.emplace_back(element);
360 if (inner.contains(it.x(), it.y(), it.z()))
361 indexVectorInner.emplace_back(element);
362 else
363 indexVectorOuter.emplace_back(element);
364 }
365 }
366
367 for (auto it = flagField->beginWithGhostLayerXYZ(
368 cell_idx_c(flagField->nrOfGhostLayers() - 1));
369 it != flagField->end(); ++it) {
370 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
371 continue;
372
373 if (isFlagSet(it.neighbor(0, 0, -1, 0), boundaryFlag)) {
374 auto element = IndexInfo(it.x(), it.y(), it.z(), 6);
375 auto const InitialisationAdditionalData = elementInitialiser(
376 Cell(it.x() + 0, it.y() + 0, it.z() - 1), blocks, *block);
377 element.vel_0 = InitialisationAdditionalData[0];
378 element.vel_1 = InitialisationAdditionalData[1];
379 element.vel_2 = InitialisationAdditionalData[2];
380 indexVectorAll.emplace_back(element);
381 if (inner.contains(it.x(), it.y(), it.z()))
382 indexVectorInner.emplace_back(element);
383 else
384 indexVectorOuter.emplace_back(element);
385 }
386 }
387
388 for (auto it = flagField->beginWithGhostLayerXYZ(
389 cell_idx_c(flagField->nrOfGhostLayers() - 1));
390 it != flagField->end(); ++it) {
391 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
392 continue;
393
394 if (isFlagSet(it.neighbor(-1, 1, 0, 0), boundaryFlag)) {
395 auto element = IndexInfo(it.x(), it.y(), it.z(), 7);
396 auto const InitialisationAdditionalData = elementInitialiser(
397 Cell(it.x() - 1, it.y() + 1, it.z() + 0), blocks, *block);
398 element.vel_0 = InitialisationAdditionalData[0];
399 element.vel_1 = InitialisationAdditionalData[1];
400 element.vel_2 = InitialisationAdditionalData[2];
401 indexVectorAll.emplace_back(element);
402 if (inner.contains(it.x(), it.y(), it.z()))
403 indexVectorInner.emplace_back(element);
404 else
405 indexVectorOuter.emplace_back(element);
406 }
407 }
408
409 for (auto it = flagField->beginWithGhostLayerXYZ(
410 cell_idx_c(flagField->nrOfGhostLayers() - 1));
411 it != flagField->end(); ++it) {
412 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
413 continue;
414
415 if (isFlagSet(it.neighbor(1, 1, 0, 0), boundaryFlag)) {
416 auto element = IndexInfo(it.x(), it.y(), it.z(), 8);
417 auto const InitialisationAdditionalData = elementInitialiser(
418 Cell(it.x() + 1, it.y() + 1, it.z() + 0), blocks, *block);
419 element.vel_0 = InitialisationAdditionalData[0];
420 element.vel_1 = InitialisationAdditionalData[1];
421 element.vel_2 = InitialisationAdditionalData[2];
422 indexVectorAll.emplace_back(element);
423 if (inner.contains(it.x(), it.y(), it.z()))
424 indexVectorInner.emplace_back(element);
425 else
426 indexVectorOuter.emplace_back(element);
427 }
428 }
429
430 for (auto it = flagField->beginWithGhostLayerXYZ(
431 cell_idx_c(flagField->nrOfGhostLayers() - 1));
432 it != flagField->end(); ++it) {
433 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
434 continue;
435
436 if (isFlagSet(it.neighbor(-1, -1, 0, 0), boundaryFlag)) {
437 auto element = IndexInfo(it.x(), it.y(), it.z(), 9);
438 auto const InitialisationAdditionalData = elementInitialiser(
439 Cell(it.x() - 1, it.y() - 1, it.z() + 0), blocks, *block);
440 element.vel_0 = InitialisationAdditionalData[0];
441 element.vel_1 = InitialisationAdditionalData[1];
442 element.vel_2 = InitialisationAdditionalData[2];
443 indexVectorAll.emplace_back(element);
444 if (inner.contains(it.x(), it.y(), it.z()))
445 indexVectorInner.emplace_back(element);
446 else
447 indexVectorOuter.emplace_back(element);
448 }
449 }
450
451 for (auto it = flagField->beginWithGhostLayerXYZ(
452 cell_idx_c(flagField->nrOfGhostLayers() - 1));
453 it != flagField->end(); ++it) {
454 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
455 continue;
456
457 if (isFlagSet(it.neighbor(1, -1, 0, 0), boundaryFlag)) {
458 auto element = IndexInfo(it.x(), it.y(), it.z(), 10);
459 auto const InitialisationAdditionalData = elementInitialiser(
460 Cell(it.x() + 1, it.y() - 1, it.z() + 0), blocks, *block);
461 element.vel_0 = InitialisationAdditionalData[0];
462 element.vel_1 = InitialisationAdditionalData[1];
463 element.vel_2 = InitialisationAdditionalData[2];
464 indexVectorAll.emplace_back(element);
465 if (inner.contains(it.x(), it.y(), it.z()))
466 indexVectorInner.emplace_back(element);
467 else
468 indexVectorOuter.emplace_back(element);
469 }
470 }
471
472 for (auto it = flagField->beginWithGhostLayerXYZ(
473 cell_idx_c(flagField->nrOfGhostLayers() - 1));
474 it != flagField->end(); ++it) {
475 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
476 continue;
477
478 if (isFlagSet(it.neighbor(0, 1, 1, 0), boundaryFlag)) {
479 auto element = IndexInfo(it.x(), it.y(), it.z(), 11);
480 auto const InitialisationAdditionalData = elementInitialiser(
481 Cell(it.x() + 0, it.y() + 1, it.z() + 1), blocks, *block);
482 element.vel_0 = InitialisationAdditionalData[0];
483 element.vel_1 = InitialisationAdditionalData[1];
484 element.vel_2 = InitialisationAdditionalData[2];
485 indexVectorAll.emplace_back(element);
486 if (inner.contains(it.x(), it.y(), it.z()))
487 indexVectorInner.emplace_back(element);
488 else
489 indexVectorOuter.emplace_back(element);
490 }
491 }
492
493 for (auto it = flagField->beginWithGhostLayerXYZ(
494 cell_idx_c(flagField->nrOfGhostLayers() - 1));
495 it != flagField->end(); ++it) {
496 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
497 continue;
498
499 if (isFlagSet(it.neighbor(0, -1, 1, 0), boundaryFlag)) {
500 auto element = IndexInfo(it.x(), it.y(), it.z(), 12);
501 auto const InitialisationAdditionalData = elementInitialiser(
502 Cell(it.x() + 0, it.y() - 1, it.z() + 1), blocks, *block);
503 element.vel_0 = InitialisationAdditionalData[0];
504 element.vel_1 = InitialisationAdditionalData[1];
505 element.vel_2 = InitialisationAdditionalData[2];
506 indexVectorAll.emplace_back(element);
507 if (inner.contains(it.x(), it.y(), it.z()))
508 indexVectorInner.emplace_back(element);
509 else
510 indexVectorOuter.emplace_back(element);
511 }
512 }
513
514 for (auto it = flagField->beginWithGhostLayerXYZ(
515 cell_idx_c(flagField->nrOfGhostLayers() - 1));
516 it != flagField->end(); ++it) {
517 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
518 continue;
519
520 if (isFlagSet(it.neighbor(-1, 0, 1, 0), boundaryFlag)) {
521 auto element = IndexInfo(it.x(), it.y(), it.z(), 13);
522 auto const InitialisationAdditionalData = elementInitialiser(
523 Cell(it.x() - 1, it.y() + 0, it.z() + 1), blocks, *block);
524 element.vel_0 = InitialisationAdditionalData[0];
525 element.vel_1 = InitialisationAdditionalData[1];
526 element.vel_2 = InitialisationAdditionalData[2];
527 indexVectorAll.emplace_back(element);
528 if (inner.contains(it.x(), it.y(), it.z()))
529 indexVectorInner.emplace_back(element);
530 else
531 indexVectorOuter.emplace_back(element);
532 }
533 }
534
535 for (auto it = flagField->beginWithGhostLayerXYZ(
536 cell_idx_c(flagField->nrOfGhostLayers() - 1));
537 it != flagField->end(); ++it) {
538 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
539 continue;
540
541 if (isFlagSet(it.neighbor(1, 0, 1, 0), boundaryFlag)) {
542 auto element = IndexInfo(it.x(), it.y(), it.z(), 14);
543 auto const InitialisationAdditionalData = elementInitialiser(
544 Cell(it.x() + 1, it.y() + 0, it.z() + 1), blocks, *block);
545 element.vel_0 = InitialisationAdditionalData[0];
546 element.vel_1 = InitialisationAdditionalData[1];
547 element.vel_2 = InitialisationAdditionalData[2];
548 indexVectorAll.emplace_back(element);
549 if (inner.contains(it.x(), it.y(), it.z()))
550 indexVectorInner.emplace_back(element);
551 else
552 indexVectorOuter.emplace_back(element);
553 }
554 }
555
556 for (auto it = flagField->beginWithGhostLayerXYZ(
557 cell_idx_c(flagField->nrOfGhostLayers() - 1));
558 it != flagField->end(); ++it) {
559 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
560 continue;
561
562 if (isFlagSet(it.neighbor(0, 1, -1, 0), boundaryFlag)) {
563 auto element = IndexInfo(it.x(), it.y(), it.z(), 15);
564 auto const InitialisationAdditionalData = elementInitialiser(
565 Cell(it.x() + 0, it.y() + 1, it.z() - 1), blocks, *block);
566 element.vel_0 = InitialisationAdditionalData[0];
567 element.vel_1 = InitialisationAdditionalData[1];
568 element.vel_2 = InitialisationAdditionalData[2];
569 indexVectorAll.emplace_back(element);
570 if (inner.contains(it.x(), it.y(), it.z()))
571 indexVectorInner.emplace_back(element);
572 else
573 indexVectorOuter.emplace_back(element);
574 }
575 }
576
577 for (auto it = flagField->beginWithGhostLayerXYZ(
578 cell_idx_c(flagField->nrOfGhostLayers() - 1));
579 it != flagField->end(); ++it) {
580 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
581 continue;
582
583 if (isFlagSet(it.neighbor(0, -1, -1, 0), boundaryFlag)) {
584 auto element = IndexInfo(it.x(), it.y(), it.z(), 16);
585 auto const InitialisationAdditionalData = elementInitialiser(
586 Cell(it.x() + 0, it.y() - 1, it.z() - 1), blocks, *block);
587 element.vel_0 = InitialisationAdditionalData[0];
588 element.vel_1 = InitialisationAdditionalData[1];
589 element.vel_2 = InitialisationAdditionalData[2];
590 indexVectorAll.emplace_back(element);
591 if (inner.contains(it.x(), it.y(), it.z()))
592 indexVectorInner.emplace_back(element);
593 else
594 indexVectorOuter.emplace_back(element);
595 }
596 }
597
598 for (auto it = flagField->beginWithGhostLayerXYZ(
599 cell_idx_c(flagField->nrOfGhostLayers() - 1));
600 it != flagField->end(); ++it) {
601 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
602 continue;
603
604 if (isFlagSet(it.neighbor(-1, 0, -1, 0), boundaryFlag)) {
605 auto element = IndexInfo(it.x(), it.y(), it.z(), 17);
606 auto const InitialisationAdditionalData = elementInitialiser(
607 Cell(it.x() - 1, it.y() + 0, it.z() - 1), blocks, *block);
608 element.vel_0 = InitialisationAdditionalData[0];
609 element.vel_1 = InitialisationAdditionalData[1];
610 element.vel_2 = InitialisationAdditionalData[2];
611 indexVectorAll.emplace_back(element);
612 if (inner.contains(it.x(), it.y(), it.z()))
613 indexVectorInner.emplace_back(element);
614 else
615 indexVectorOuter.emplace_back(element);
616 }
617 }
618
619 for (auto it = flagField->beginWithGhostLayerXYZ(
620 cell_idx_c(flagField->nrOfGhostLayers() - 1));
621 it != flagField->end(); ++it) {
622 if (!isFlagSet(it, domainFlag) || isFlagSet(it, boundaryFlag))
623 continue;
624
625 if (isFlagSet(it.neighbor(1, 0, -1, 0), boundaryFlag)) {
626 auto element = IndexInfo(it.x(), it.y(), it.z(), 18);
627 auto const InitialisationAdditionalData = elementInitialiser(
628 Cell(it.x() + 1, it.y() + 0, it.z() - 1), blocks, *block);
629 element.vel_0 = InitialisationAdditionalData[0];
630 element.vel_1 = InitialisationAdditionalData[1];
631 element.vel_2 = InitialisationAdditionalData[2];
632 indexVectorAll.emplace_back(element);
633 if (inner.contains(it.x(), it.y(), it.z()))
634 indexVectorInner.emplace_back(element);
635 else
636 indexVectorOuter.emplace_back(element);
637 }
638 }
639
640 indexVectors->syncGPU();
641 forceVector->forceVector().resize(indexVectorAll.size());
642 forceVector->syncGPU();
643 }
644
645private:
646 void run_impl(IBlock *block, IndexVectors::Type type);
647
648 BlockDataID indexVectorID;
649 BlockDataID forceVectorID;
650 std::function<Vector3<double>(
651 const Cell &, const shared_ptr<StructuredBlockForest> &, IBlock &)>
652 elementInitialiser;
653
654public:
655 static constexpr std::array<std::array<int, 19u>, 3u> neighborOffset = {{
656 {0, 0, 0, -1, 1, 0, 0, -1, 1, -1, 1, 0, 0, -1, 1, 0, 0, -1, 1},
657 {0, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, 1, -1, 0, 0, 1, -1, 0, 0},
658 {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, 1, 1, 1, -1, -1, -1, -1},
659 }};
660
661 auto const &getForceVector(IBlock const *block) const {
662 auto const *forceVector = block->getData<ForceVector>(forceVectorID);
663 return forceVector->forceVector();
664 }
665
666 auto const &getIndexVector(IBlock const *block) const {
667 auto const *indexVectors = block->getData<IndexVectors>(indexVectorID);
668 return indexVectors->indexVector(IndexVectors::ALL);
669 }
670
671 BlockDataID getIndexVectorID() const { return indexVectorID; }
672 BlockDataID getForceVectorID() const { return forceVectorID; }
673
674public:
675 BlockDataID pdfsID;
676};
677
678#if defined(__clang__)
679#pragma clang diagnostic pop
680#elif defined(__GNUC__) or defined(__GNUG__)
681#pragma GCC diagnostic pop
682#endif
683
684} // namespace lbm
685} // namespace walberla
Cell
Definition Cell.hpp:96
walberla::lbm::DynamicUBBDoublePrecision::fillFromFlagField
void fillFromFlagField(const shared_ptr< StructuredBlockForest > &blocks, IBlock *block, ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, FlagUID domainFlagUID)
Definition DynamicUBBDoublePrecision.h:216
walberla::lbm::DynamicUBBDoublePrecision::fillFromFlagField
void fillFromFlagField(const std::shared_ptr< StructuredBlockForest > &blocks, ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, FlagUID domainFlagUID)
Definition DynamicUBBDoublePrecision.h:207
walberla::lbm::DynamicUBBDoublePrecision::DynamicUBBDoublePrecision
DynamicUBBDoublePrecision(const std::shared_ptr< StructuredBlockForest > &blocks, BlockDataID pdfsID_, std::function< Vector3< double >(const Cell &, const shared_ptr< StructuredBlockForest > &, IBlock &)> &velocityCallbackDynamicUBBDoublePrecision)
Definition DynamicUBBDoublePrecision.h:160
walberla::lbm::DynamicUBBDoublePrecision::neighborOffset
static constexpr std::array< std::array< int, 19u >, 3u > neighborOffset
Definition DynamicUBBDoublePrecision.h:655
walberla::lbm::DynamicUBBDoublePrecision::getForceVector
auto const & getForceVector(IBlock const *block) const
Definition DynamicUBBDoublePrecision.h:661
walberla::lbm::DynamicUBBDoublePrecision::getIndexVector
auto const & getIndexVector(IBlock const *block) const
Definition DynamicUBBDoublePrecision.h:666
block
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:176
walberla
\file PackInfoPdfDoublePrecision.cpp \author pystencils
Definition EKWalberla.hpp:36
walberla::lbm::DynamicUBBDoublePrecision::IndexInfo::IndexInfo
IndexInfo(int32_t x_, int32_t y_, int32_t z_, int32_t dir_)
Definition DynamicUBBDoublePrecision.h:83