Loading [MathJax]/extensions/TeX/AMSmath.js
ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages Concepts
DynamicUBBDoublePrecisionCUDA.h
Go to the documentation of this file.
1//======================================================================================================================
2//
3// This file is part of waLBerla. waLBerla is free software: you can
4// redistribute it and/or modify it under the terms of the GNU General Public
5// License as published by the Free Software Foundation, either version 3 of
6// the License, or (at your option) any later version.
7//
8// waLBerla is distributed in the hope that it will be useful, but WITHOUT
9// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
11// for more details.
12//
13// You should have received a copy of the GNU General Public License along
14// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
15//
16//! \\file DynamicUBBDoublePrecisionCUDA.h
17//! \\author pystencils
18//======================================================================================================================
19
20// kernel generated with pystencils v1.3.7, lbmpy v1.3.7, sympy v1.12.1,
21// lbmpy_walberla/pystencils_walberla from waLBerla commit
22// f36fa0a68bae59f0b516f6587ea8fa7c24a41141
23
24#pragma once
25#include "core/DataTypes.h"
26#include "core/logging/Logging.h"
27
28#include "blockforest/StructuredBlockForest.h"
29#include "core/debug/Debug.h"
30#include "domain_decomposition/BlockDataID.h"
31#include "domain_decomposition/IBlock.h"
32#include "field/FlagField.h"
33#include "gpu/FieldCopy.h"
34#include "gpu/GPUField.h"
35#include "gpu/GPUWrapper.h"
36
37#include <set>
38#include <vector>
39
40#ifdef __GNUC__
41#define RESTRICT __restrict__
42#elif _MSC_VER
43#define RESTRICT __restrict
44#else
45#define RESTRICT
46#endif
47
48#ifdef WALBERLA_BUILD_WITH_HALF_PRECISION_SUPPORT
49using walberla::half;
50#endif
51
52namespace walberla {
53namespace lbm {
54
55class DynamicUBBDoublePrecisionCUDA {
56public:
57 struct IndexInfo {
58 int32_t x;
59 int32_t y;
60 int32_t z;
61 int32_t dir;
62 double vel_0;
63 double vel_1;
64 double vel_2;
65 IndexInfo(int32_t x_, int32_t y_, int32_t z_, int32_t dir_)
66 : x(x_), y(y_), z(z_), dir(dir_), vel_0(), vel_1(), vel_2() {}
67 bool operator==(const IndexInfo &o) const {
68 return x == o.x && y == o.y && z == o.z && dir == o.dir &&
69 floatIsEqual(vel_0, o.vel_0) && floatIsEqual(vel_1, o.vel_1) &&
70 floatIsEqual(vel_2, o.vel_2);
71 }
72 };
73
74 class IndexVectors {
75 public:
76 using CpuIndexVector = std::vector<IndexInfo>;
77
78 enum Type { ALL = 0, INNER = 1, OUTER = 2, NUM_TYPES = 3 };
79
80 IndexVectors() = default;
81 bool operator==(IndexVectors const &other) const {
82 return other.cpuVectors_ == cpuVectors_;
83 }
84
85 ~IndexVectors() {
86 for (auto &gpuVec : gpuVectors_)
87 WALBERLA_GPU_CHECK(gpuFree(gpuVec));
88 }
89 CpuIndexVector &indexVector(Type t) { return cpuVectors_[t]; }
90 IndexInfo *pointerCpu(Type t) { return cpuVectors_[t].data(); }
91
92 IndexInfo *pointerGpu(Type t) { return gpuVectors_[t]; }
93 void syncGPU() {
94 for (auto &gpuVec : gpuVectors_)
95 WALBERLA_GPU_CHECK(gpuFree(gpuVec));
96 gpuVectors_.resize(cpuVectors_.size());
97
98 WALBERLA_ASSERT_EQUAL(cpuVectors_.size(), NUM_TYPES);
99 for (size_t i = 0; i < cpuVectors_.size(); ++i) {
100 auto &gpuVec = gpuVectors_[i];
101 auto &cpuVec = cpuVectors_[i];
102 WALBERLA_GPU_CHECK(
103 gpuMalloc(&gpuVec, sizeof(IndexInfo) * cpuVec.size()));
104 WALBERLA_GPU_CHECK(gpuMemcpy(gpuVec, &cpuVec[0],
105 sizeof(IndexInfo) * cpuVec.size(),
106 gpuMemcpyHostToDevice));
107 }
108 }
109
110 private:
111 std::vector<CpuIndexVector> cpuVectors_{NUM_TYPES};
112
113 using GpuIndexVector = IndexInfo *;
114 std::vector<GpuIndexVector> gpuVectors_;
115 };
116
117 DynamicUBBDoublePrecisionCUDA(
118 const shared_ptr<StructuredBlockForest> &blocks, BlockDataID pdfsID_,
119 std::function<Vector3<float64>(const Cell &,
120 const shared_ptr<StructuredBlockForest> &,
121 IBlock &)> &velocityCallback)
122 : elementInitialiser(velocityCallback), pdfsID(pdfsID_) {
123 auto createIdxVector = [](IBlock *const, StructuredBlockStorage *const) {
124 return new IndexVectors();
125 };
126 indexVectorID = blocks->addStructuredBlockData<IndexVectors>(
127 createIdxVector, "IndexField_DynamicUBBDoublePrecisionCUDA");
128 }
129
130 void run(IBlock *block, gpuStream_t stream = nullptr);
131
132 void operator()(IBlock *block, gpuStream_t stream = nullptr) {
133 run(block, stream);
134 }
135
136 void inner(IBlock *block, gpuStream_t stream = nullptr);
137
138 void outer(IBlock *block, gpuStream_t stream = nullptr);
139
140 Vector3<double> getForce(IBlock * /*block*/) {
141
142 WALBERLA_ABORT(
143 "Boundary condition was not generated including force calculation.")
144 return Vector3<double>(double_c(0.0));
145 }
146
147 std::function<void(IBlock *)> getSweep(gpuStream_t stream = nullptr) {
148 return [this, stream](IBlock *b) { this->run(b, stream); };
149 }
150
151 std::function<void(IBlock *)> getInnerSweep(gpuStream_t stream = nullptr) {
152 return [this, stream](IBlock *b) { this->inner(b, stream); };
153 }
154
155 std::function<void(IBlock *)> getOuterSweep(gpuStream_t stream = nullptr) {
156 return [this, stream](IBlock *b) { this->outer(b, stream); };
157 }
158
159 template <typename FlagField_T>
160 void fillFromFlagField(const shared_ptr<StructuredBlockForest> &blocks,
161 ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID,
162 FlagUID domainFlagUID) {
163 for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
164 fillFromFlagField<FlagField_T>(blocks, &*blockIt, flagFieldID,
165 boundaryFlagUID, domainFlagUID);
166 }
167
168 template <typename FlagField_T>
169 void fillFromFlagField(const shared_ptr<StructuredBlockForest> &blocks,
170 IBlock *block, ConstBlockDataID flagFieldID,
171 FlagUID boundaryFlagUID, FlagUID domainFlagUID) {
172 auto *indexVectors = block->getData<IndexVectors>(indexVectorID);
173 auto &indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
174 auto &indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
175 auto &indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
176
177 auto *flagField = block->getData<FlagField_T>(flagFieldID);
178
179 if (!(flagField->flagExists(boundaryFlagUID) &&
180 flagField->flagExists(domainFlagUID)))
181 return;
182
183 auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
184 auto domainFlag = flagField->getFlag(domainFlagUID);
185
186 auto inner = flagField->xyzSize();
187 inner.expand(cell_idx_t(-1));
188
189 indexVectorAll.clear();
190 indexVectorInner.clear();
191 indexVectorOuter.clear();
192
193 for (auto it = flagField->beginWithGhostLayerXYZ(
194 cell_idx_c(flagField->nrOfGhostLayers() - 1));
195 it != flagField->end(); ++it) {
196 if (!isFlagSet(it, domainFlag))
197 continue;
198
199 if (isFlagSet(it.neighbor(0, 0, 0, 0), boundaryFlag)) {
200 auto element = IndexInfo(it.x(), it.y(), it.z(), 0);
201 auto const InitialisationAdditionalData = elementInitialiser(
202 Cell(it.x() + 0, it.y() + 0, it.z() + 0), blocks, *block);
203 element.vel_0 = InitialisationAdditionalData[0];
204 element.vel_1 = InitialisationAdditionalData[1];
205 element.vel_2 = InitialisationAdditionalData[2];
206 indexVectorAll.push_back(element);
207 if (inner.contains(it.x(), it.y(), it.z()))
208 indexVectorInner.push_back(element);
209 else
210 indexVectorOuter.push_back(element);
211 }
212 }
213
214 for (auto it = flagField->beginWithGhostLayerXYZ(
215 cell_idx_c(flagField->nrOfGhostLayers() - 1));
216 it != flagField->end(); ++it) {
217 if (!isFlagSet(it, domainFlag))
218 continue;
219
220 if (isFlagSet(it.neighbor(0, 1, 0, 0), boundaryFlag)) {
221 auto element = IndexInfo(it.x(), it.y(), it.z(), 1);
222 auto const InitialisationAdditionalData = elementInitialiser(
223 Cell(it.x() + 0, it.y() + 1, it.z() + 0), blocks, *block);
224 element.vel_0 = InitialisationAdditionalData[0];
225 element.vel_1 = InitialisationAdditionalData[1];
226 element.vel_2 = InitialisationAdditionalData[2];
227 indexVectorAll.push_back(element);
228 if (inner.contains(it.x(), it.y(), it.z()))
229 indexVectorInner.push_back(element);
230 else
231 indexVectorOuter.push_back(element);
232 }
233 }
234
235 for (auto it = flagField->beginWithGhostLayerXYZ(
236 cell_idx_c(flagField->nrOfGhostLayers() - 1));
237 it != flagField->end(); ++it) {
238 if (!isFlagSet(it, domainFlag))
239 continue;
240
241 if (isFlagSet(it.neighbor(0, -1, 0, 0), boundaryFlag)) {
242 auto element = IndexInfo(it.x(), it.y(), it.z(), 2);
243 auto const InitialisationAdditionalData = elementInitialiser(
244 Cell(it.x() + 0, it.y() - 1, it.z() + 0), blocks, *block);
245 element.vel_0 = InitialisationAdditionalData[0];
246 element.vel_1 = InitialisationAdditionalData[1];
247 element.vel_2 = InitialisationAdditionalData[2];
248 indexVectorAll.push_back(element);
249 if (inner.contains(it.x(), it.y(), it.z()))
250 indexVectorInner.push_back(element);
251 else
252 indexVectorOuter.push_back(element);
253 }
254 }
255
256 for (auto it = flagField->beginWithGhostLayerXYZ(
257 cell_idx_c(flagField->nrOfGhostLayers() - 1));
258 it != flagField->end(); ++it) {
259 if (!isFlagSet(it, domainFlag))
260 continue;
261
262 if (isFlagSet(it.neighbor(-1, 0, 0, 0), boundaryFlag)) {
263 auto element = IndexInfo(it.x(), it.y(), it.z(), 3);
264 auto const InitialisationAdditionalData = elementInitialiser(
265 Cell(it.x() - 1, it.y() + 0, it.z() + 0), blocks, *block);
266 element.vel_0 = InitialisationAdditionalData[0];
267 element.vel_1 = InitialisationAdditionalData[1];
268 element.vel_2 = InitialisationAdditionalData[2];
269 indexVectorAll.push_back(element);
270 if (inner.contains(it.x(), it.y(), it.z()))
271 indexVectorInner.push_back(element);
272 else
273 indexVectorOuter.push_back(element);
274 }
275 }
276
277 for (auto it = flagField->beginWithGhostLayerXYZ(
278 cell_idx_c(flagField->nrOfGhostLayers() - 1));
279 it != flagField->end(); ++it) {
280 if (!isFlagSet(it, domainFlag))
281 continue;
282
283 if (isFlagSet(it.neighbor(1, 0, 0, 0), boundaryFlag)) {
284 auto element = IndexInfo(it.x(), it.y(), it.z(), 4);
285 auto const InitialisationAdditionalData = elementInitialiser(
286 Cell(it.x() + 1, it.y() + 0, it.z() + 0), blocks, *block);
287 element.vel_0 = InitialisationAdditionalData[0];
288 element.vel_1 = InitialisationAdditionalData[1];
289 element.vel_2 = InitialisationAdditionalData[2];
290 indexVectorAll.push_back(element);
291 if (inner.contains(it.x(), it.y(), it.z()))
292 indexVectorInner.push_back(element);
293 else
294 indexVectorOuter.push_back(element);
295 }
296 }
297
298 for (auto it = flagField->beginWithGhostLayerXYZ(
299 cell_idx_c(flagField->nrOfGhostLayers() - 1));
300 it != flagField->end(); ++it) {
301 if (!isFlagSet(it, domainFlag))
302 continue;
303
304 if (isFlagSet(it.neighbor(0, 0, 1, 0), boundaryFlag)) {
305 auto element = IndexInfo(it.x(), it.y(), it.z(), 5);
306 auto const InitialisationAdditionalData = elementInitialiser(
307 Cell(it.x() + 0, it.y() + 0, it.z() + 1), blocks, *block);
308 element.vel_0 = InitialisationAdditionalData[0];
309 element.vel_1 = InitialisationAdditionalData[1];
310 element.vel_2 = InitialisationAdditionalData[2];
311 indexVectorAll.push_back(element);
312 if (inner.contains(it.x(), it.y(), it.z()))
313 indexVectorInner.push_back(element);
314 else
315 indexVectorOuter.push_back(element);
316 }
317 }
318
319 for (auto it = flagField->beginWithGhostLayerXYZ(
320 cell_idx_c(flagField->nrOfGhostLayers() - 1));
321 it != flagField->end(); ++it) {
322 if (!isFlagSet(it, domainFlag))
323 continue;
324
325 if (isFlagSet(it.neighbor(0, 0, -1, 0), boundaryFlag)) {
326 auto element = IndexInfo(it.x(), it.y(), it.z(), 6);
327 auto const InitialisationAdditionalData = elementInitialiser(
328 Cell(it.x() + 0, it.y() + 0, it.z() - 1), blocks, *block);
329 element.vel_0 = InitialisationAdditionalData[0];
330 element.vel_1 = InitialisationAdditionalData[1];
331 element.vel_2 = InitialisationAdditionalData[2];
332 indexVectorAll.push_back(element);
333 if (inner.contains(it.x(), it.y(), it.z()))
334 indexVectorInner.push_back(element);
335 else
336 indexVectorOuter.push_back(element);
337 }
338 }
339
340 for (auto it = flagField->beginWithGhostLayerXYZ(
341 cell_idx_c(flagField->nrOfGhostLayers() - 1));
342 it != flagField->end(); ++it) {
343 if (!isFlagSet(it, domainFlag))
344 continue;
345
346 if (isFlagSet(it.neighbor(-1, 1, 0, 0), boundaryFlag)) {
347 auto element = IndexInfo(it.x(), it.y(), it.z(), 7);
348 auto const InitialisationAdditionalData = elementInitialiser(
349 Cell(it.x() - 1, it.y() + 1, it.z() + 0), blocks, *block);
350 element.vel_0 = InitialisationAdditionalData[0];
351 element.vel_1 = InitialisationAdditionalData[1];
352 element.vel_2 = InitialisationAdditionalData[2];
353 indexVectorAll.push_back(element);
354 if (inner.contains(it.x(), it.y(), it.z()))
355 indexVectorInner.push_back(element);
356 else
357 indexVectorOuter.push_back(element);
358 }
359 }
360
361 for (auto it = flagField->beginWithGhostLayerXYZ(
362 cell_idx_c(flagField->nrOfGhostLayers() - 1));
363 it != flagField->end(); ++it) {
364 if (!isFlagSet(it, domainFlag))
365 continue;
366
367 if (isFlagSet(it.neighbor(1, 1, 0, 0), boundaryFlag)) {
368 auto element = IndexInfo(it.x(), it.y(), it.z(), 8);
369 auto const InitialisationAdditionalData = elementInitialiser(
370 Cell(it.x() + 1, it.y() + 1, it.z() + 0), blocks, *block);
371 element.vel_0 = InitialisationAdditionalData[0];
372 element.vel_1 = InitialisationAdditionalData[1];
373 element.vel_2 = InitialisationAdditionalData[2];
374 indexVectorAll.push_back(element);
375 if (inner.contains(it.x(), it.y(), it.z()))
376 indexVectorInner.push_back(element);
377 else
378 indexVectorOuter.push_back(element);
379 }
380 }
381
382 for (auto it = flagField->beginWithGhostLayerXYZ(
383 cell_idx_c(flagField->nrOfGhostLayers() - 1));
384 it != flagField->end(); ++it) {
385 if (!isFlagSet(it, domainFlag))
386 continue;
387
388 if (isFlagSet(it.neighbor(-1, -1, 0, 0), boundaryFlag)) {
389 auto element = IndexInfo(it.x(), it.y(), it.z(), 9);
390 auto const InitialisationAdditionalData = elementInitialiser(
391 Cell(it.x() - 1, it.y() - 1, it.z() + 0), blocks, *block);
392 element.vel_0 = InitialisationAdditionalData[0];
393 element.vel_1 = InitialisationAdditionalData[1];
394 element.vel_2 = InitialisationAdditionalData[2];
395 indexVectorAll.push_back(element);
396 if (inner.contains(it.x(), it.y(), it.z()))
397 indexVectorInner.push_back(element);
398 else
399 indexVectorOuter.push_back(element);
400 }
401 }
402
403 for (auto it = flagField->beginWithGhostLayerXYZ(
404 cell_idx_c(flagField->nrOfGhostLayers() - 1));
405 it != flagField->end(); ++it) {
406 if (!isFlagSet(it, domainFlag))
407 continue;
408
409 if (isFlagSet(it.neighbor(1, -1, 0, 0), boundaryFlag)) {
410 auto element = IndexInfo(it.x(), it.y(), it.z(), 10);
411 auto const InitialisationAdditionalData = elementInitialiser(
412 Cell(it.x() + 1, it.y() - 1, it.z() + 0), blocks, *block);
413 element.vel_0 = InitialisationAdditionalData[0];
414 element.vel_1 = InitialisationAdditionalData[1];
415 element.vel_2 = InitialisationAdditionalData[2];
416 indexVectorAll.push_back(element);
417 if (inner.contains(it.x(), it.y(), it.z()))
418 indexVectorInner.push_back(element);
419 else
420 indexVectorOuter.push_back(element);
421 }
422 }
423
424 for (auto it = flagField->beginWithGhostLayerXYZ(
425 cell_idx_c(flagField->nrOfGhostLayers() - 1));
426 it != flagField->end(); ++it) {
427 if (!isFlagSet(it, domainFlag))
428 continue;
429
430 if (isFlagSet(it.neighbor(0, 1, 1, 0), boundaryFlag)) {
431 auto element = IndexInfo(it.x(), it.y(), it.z(), 11);
432 auto const InitialisationAdditionalData = elementInitialiser(
433 Cell(it.x() + 0, it.y() + 1, it.z() + 1), blocks, *block);
434 element.vel_0 = InitialisationAdditionalData[0];
435 element.vel_1 = InitialisationAdditionalData[1];
436 element.vel_2 = InitialisationAdditionalData[2];
437 indexVectorAll.push_back(element);
438 if (inner.contains(it.x(), it.y(), it.z()))
439 indexVectorInner.push_back(element);
440 else
441 indexVectorOuter.push_back(element);
442 }
443 }
444
445 for (auto it = flagField->beginWithGhostLayerXYZ(
446 cell_idx_c(flagField->nrOfGhostLayers() - 1));
447 it != flagField->end(); ++it) {
448 if (!isFlagSet(it, domainFlag))
449 continue;
450
451 if (isFlagSet(it.neighbor(0, -1, 1, 0), boundaryFlag)) {
452 auto element = IndexInfo(it.x(), it.y(), it.z(), 12);
453 auto const InitialisationAdditionalData = elementInitialiser(
454 Cell(it.x() + 0, it.y() - 1, it.z() + 1), blocks, *block);
455 element.vel_0 = InitialisationAdditionalData[0];
456 element.vel_1 = InitialisationAdditionalData[1];
457 element.vel_2 = InitialisationAdditionalData[2];
458 indexVectorAll.push_back(element);
459 if (inner.contains(it.x(), it.y(), it.z()))
460 indexVectorInner.push_back(element);
461 else
462 indexVectorOuter.push_back(element);
463 }
464 }
465
466 for (auto it = flagField->beginWithGhostLayerXYZ(
467 cell_idx_c(flagField->nrOfGhostLayers() - 1));
468 it != flagField->end(); ++it) {
469 if (!isFlagSet(it, domainFlag))
470 continue;
471
472 if (isFlagSet(it.neighbor(-1, 0, 1, 0), boundaryFlag)) {
473 auto element = IndexInfo(it.x(), it.y(), it.z(), 13);
474 auto const InitialisationAdditionalData = elementInitialiser(
475 Cell(it.x() - 1, it.y() + 0, it.z() + 1), blocks, *block);
476 element.vel_0 = InitialisationAdditionalData[0];
477 element.vel_1 = InitialisationAdditionalData[1];
478 element.vel_2 = InitialisationAdditionalData[2];
479 indexVectorAll.push_back(element);
480 if (inner.contains(it.x(), it.y(), it.z()))
481 indexVectorInner.push_back(element);
482 else
483 indexVectorOuter.push_back(element);
484 }
485 }
486
487 for (auto it = flagField->beginWithGhostLayerXYZ(
488 cell_idx_c(flagField->nrOfGhostLayers() - 1));
489 it != flagField->end(); ++it) {
490 if (!isFlagSet(it, domainFlag))
491 continue;
492
493 if (isFlagSet(it.neighbor(1, 0, 1, 0), boundaryFlag)) {
494 auto element = IndexInfo(it.x(), it.y(), it.z(), 14);
495 auto const InitialisationAdditionalData = elementInitialiser(
496 Cell(it.x() + 1, it.y() + 0, it.z() + 1), blocks, *block);
497 element.vel_0 = InitialisationAdditionalData[0];
498 element.vel_1 = InitialisationAdditionalData[1];
499 element.vel_2 = InitialisationAdditionalData[2];
500 indexVectorAll.push_back(element);
501 if (inner.contains(it.x(), it.y(), it.z()))
502 indexVectorInner.push_back(element);
503 else
504 indexVectorOuter.push_back(element);
505 }
506 }
507
508 for (auto it = flagField->beginWithGhostLayerXYZ(
509 cell_idx_c(flagField->nrOfGhostLayers() - 1));
510 it != flagField->end(); ++it) {
511 if (!isFlagSet(it, domainFlag))
512 continue;
513
514 if (isFlagSet(it.neighbor(0, 1, -1, 0), boundaryFlag)) {
515 auto element = IndexInfo(it.x(), it.y(), it.z(), 15);
516 auto const InitialisationAdditionalData = elementInitialiser(
517 Cell(it.x() + 0, it.y() + 1, it.z() - 1), blocks, *block);
518 element.vel_0 = InitialisationAdditionalData[0];
519 element.vel_1 = InitialisationAdditionalData[1];
520 element.vel_2 = InitialisationAdditionalData[2];
521 indexVectorAll.push_back(element);
522 if (inner.contains(it.x(), it.y(), it.z()))
523 indexVectorInner.push_back(element);
524 else
525 indexVectorOuter.push_back(element);
526 }
527 }
528
529 for (auto it = flagField->beginWithGhostLayerXYZ(
530 cell_idx_c(flagField->nrOfGhostLayers() - 1));
531 it != flagField->end(); ++it) {
532 if (!isFlagSet(it, domainFlag))
533 continue;
534
535 if (isFlagSet(it.neighbor(0, -1, -1, 0), boundaryFlag)) {
536 auto element = IndexInfo(it.x(), it.y(), it.z(), 16);
537 auto const InitialisationAdditionalData = elementInitialiser(
538 Cell(it.x() + 0, it.y() - 1, it.z() - 1), blocks, *block);
539 element.vel_0 = InitialisationAdditionalData[0];
540 element.vel_1 = InitialisationAdditionalData[1];
541 element.vel_2 = InitialisationAdditionalData[2];
542 indexVectorAll.push_back(element);
543 if (inner.contains(it.x(), it.y(), it.z()))
544 indexVectorInner.push_back(element);
545 else
546 indexVectorOuter.push_back(element);
547 }
548 }
549
550 for (auto it = flagField->beginWithGhostLayerXYZ(
551 cell_idx_c(flagField->nrOfGhostLayers() - 1));
552 it != flagField->end(); ++it) {
553 if (!isFlagSet(it, domainFlag))
554 continue;
555
556 if (isFlagSet(it.neighbor(-1, 0, -1, 0), boundaryFlag)) {
557 auto element = IndexInfo(it.x(), it.y(), it.z(), 17);
558 auto const InitialisationAdditionalData = elementInitialiser(
559 Cell(it.x() - 1, it.y() + 0, it.z() - 1), blocks, *block);
560 element.vel_0 = InitialisationAdditionalData[0];
561 element.vel_1 = InitialisationAdditionalData[1];
562 element.vel_2 = InitialisationAdditionalData[2];
563 indexVectorAll.push_back(element);
564 if (inner.contains(it.x(), it.y(), it.z()))
565 indexVectorInner.push_back(element);
566 else
567 indexVectorOuter.push_back(element);
568 }
569 }
570
571 for (auto it = flagField->beginWithGhostLayerXYZ(
572 cell_idx_c(flagField->nrOfGhostLayers() - 1));
573 it != flagField->end(); ++it) {
574 if (!isFlagSet(it, domainFlag))
575 continue;
576
577 if (isFlagSet(it.neighbor(1, 0, -1, 0), boundaryFlag)) {
578 auto element = IndexInfo(it.x(), it.y(), it.z(), 18);
579 auto const InitialisationAdditionalData = elementInitialiser(
580 Cell(it.x() + 1, it.y() + 0, it.z() - 1), blocks, *block);
581 element.vel_0 = InitialisationAdditionalData[0];
582 element.vel_1 = InitialisationAdditionalData[1];
583 element.vel_2 = InitialisationAdditionalData[2];
584 indexVectorAll.push_back(element);
585 if (inner.contains(it.x(), it.y(), it.z()))
586 indexVectorInner.push_back(element);
587 else
588 indexVectorOuter.push_back(element);
589 }
590 }
591
592 indexVectors->syncGPU();
593 }
594
595private:
596 void run_impl(IBlock *block, IndexVectors::Type type,
597 gpuStream_t stream = nullptr);
598
599 BlockDataID indexVectorID;
600
601 std::function<Vector3<float64>(
602 const Cell &, const shared_ptr<StructuredBlockForest> &, IBlock &)>
603 elementInitialiser;
604
605public:
606 BlockDataID pdfsID;
607};
608
609} // namespace lbm
610} // namespace walberla
Cell
Definition Cell.hpp:96
walberla::lbm::DynamicUBBDoublePrecisionCUDA::fillFromFlagField
void fillFromFlagField(const shared_ptr< StructuredBlockForest > &blocks, IBlock *block, ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, FlagUID domainFlagUID)
Definition DynamicUBBDoublePrecisionCUDA.h:169
walberla::lbm::DynamicUBBDoublePrecisionCUDA::getSweep
std::function< void(IBlock *)> getSweep(gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.h:147
walberla::lbm::DynamicUBBDoublePrecisionCUDA::getInnerSweep
std::function< void(IBlock *)> getInnerSweep(gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.h:151
walberla::lbm::DynamicUBBDoublePrecisionCUDA::run
void run(IBlock *block, gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.cu:177
walberla::lbm::DynamicUBBDoublePrecisionCUDA::fillFromFlagField
void fillFromFlagField(const shared_ptr< StructuredBlockForest > &blocks, ConstBlockDataID flagFieldID, FlagUID boundaryFlagUID, FlagUID domainFlagUID)
Definition DynamicUBBDoublePrecisionCUDA.h:160
walberla::lbm::DynamicUBBDoublePrecisionCUDA::getOuterSweep
std::function< void(IBlock *)> getOuterSweep(gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.h:155
walberla::lbm::DynamicUBBDoublePrecisionCUDA::DynamicUBBDoublePrecisionCUDA
DynamicUBBDoublePrecisionCUDA(const shared_ptr< StructuredBlockForest > &blocks, BlockDataID pdfsID_, std::function< Vector3< float64 >(const Cell &, const shared_ptr< StructuredBlockForest > &, IBlock &)> &velocityCallback)
Definition DynamicUBBDoublePrecisionCUDA.h:117
walberla::lbm::DynamicUBBDoublePrecisionCUDA::operator()
void operator()(IBlock *block, gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.h:132
walberla::lbm::DynamicUBBDoublePrecisionCUDA::inner
void inner(IBlock *block, gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.cu:181
walberla::lbm::DynamicUBBDoublePrecisionCUDA::outer
void outer(IBlock *block, gpuStream_t stream=nullptr)
Definition DynamicUBBDoublePrecisionCUDA.cu:185
stream
cudaStream_t stream[1]
CUDA streams for parallel computing on CPU and GPU.
Definition common_cuda.cu:34
block
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
walberla
\file PackInfoPdfDoublePrecision.cpp \author pystencils
Definition EKWalberla.hpp:36
walberla::lbm::DynamicUBBDoublePrecisionCUDA::IndexInfo::IndexInfo
IndexInfo(int32_t x_, int32_t y_, int32_t z_, int32_t dir_)
Definition DynamicUBBDoublePrecisionCUDA.h:65