MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
lbsolver.h
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#ifndef LBSOLVER_H
8#define LBSOLVER_H
9
10#include "COUPLER/surfacecoupling.h"
11#include "GEOM/geometryintersection.h"
12#include "GRID/cartesiangrid.h"
13#include "GRID/cartesiangridpoint.h"
14#include "GRID/cartesianpointbasedcell.h"
15#include "POST/postprocessing.h"
16#include "UTIL/functions.h"
17#include "UTIL/maiafftw.h"
18#include "cartesiansolver.h"
19#include "lbbndcnd.h"
20#include "lbcellcollector.h"
21#include "variables.h"
22
23class LbInterfaceCell;
24class LbParentCell;
25
26template <MInt nDim, class ppType>
27class PostProcessing;
28
29// TODO labels:LB miro: This is only introduced temporarily. For now, macroscopic
30// vairbales are calculated within the collision step (INCOLLISION). Some require a
31// calculation before the collision (PRECOLLISION) and other want an update
32// right after the propagation (POSTPROPAGATION). Until all is unified, it is
33// indicated by this enum which mode is used.
34enum LbUpdateMacroscopic { INCOLLISION, PRECOLLISION, POSTPROPAGATION };
35
39template <MInt nDim>
40class LbSolver : public maia::CartesianSolver<nDim, LbSolver<nDim>> {
41 template <MInt nDim_>
42 friend class LbBndCnd;
43
44 template <MInt nDim_>
45 friend class LbInterface;
46
47 friend class maia::CartesianSolver<nDim, LbSolver<nDim>>;
48
49 template <MInt nDim_, MInt nDist, class SysEqn>
50 friend class CouplingLB;
51
52 template <MInt nDim_, MInt nDist, class SysEqn>
53 friend class LbLpt;
54
55 template <MInt nDim_, class ppType>
56 friend class PostProcessing;
57 template <MInt nDim_>
58 friend class PostProcessingLb;
59
60 public:
61 // Type for cell properties
62 using Base = LbSolver<nDim>;
63 using Cell = typename maia::grid::tree::Tree<nDim>::Cell;
64 using SolverCell = LbCell;
65 using CartesianSolver = typename maia::CartesianSolver<nDim, LbSolver>;
66 using Grid = typename CartesianSolver::Grid;
67 using GridProxy = typename CartesianSolver::GridProxy;
68
69 // used CartesianSolver
70 using CartesianSolver::c_globalGridId;
71 using CartesianSolver::domainId;
72 using CartesianSolver::domainOffset;
73 using CartesianSolver::exchangeData;
74 using CartesianSolver::getIdentifier;
75 using CartesianSolver::grid;
76 using CartesianSolver::haloCell;
77 using CartesianSolver::haloCellId;
78 using CartesianSolver::isActive;
79 using CartesianSolver::m_bandWidth;
80 using CartesianSolver::m_freeIndices;
81 using CartesianSolver::m_initFromRestartFile;
82 using CartesianSolver::m_innerBandWidth;
83 using CartesianSolver::m_Ma;
84 using CartesianSolver::m_outerBandWidth;
85 using CartesianSolver::m_Re;
86 using CartesianSolver::m_residualInterval;
87 using CartesianSolver::m_restartFile;
88 using CartesianSolver::m_restartInterval;
89 using CartesianSolver::m_restartTimeStep;
90 using CartesianSolver::m_solutionInterval;
91 using CartesianSolver::m_solverId;
92 using CartesianSolver::maxLevel;
93 using CartesianSolver::maxNoGridCells;
94 using CartesianSolver::maxRefinementLevel;
95 using CartesianSolver::maxUniformRefinementLevel;
96 using CartesianSolver::minLevel;
97 using CartesianSolver::mpiComm;
98 using CartesianSolver::neighborDomain;
99 using CartesianSolver::noDomains;
100 using CartesianSolver::noHaloCells;
101 using CartesianSolver::noNeighborDomains;
102 using CartesianSolver::noWindowCells;
103 using CartesianSolver::outputDir;
104 using CartesianSolver::readSolverSamplingVarNames;
105 using CartesianSolver::reductionFactor;
106 using CartesianSolver::restartDir;
107 using CartesianSolver::solverId;
108 using CartesianSolver::solverMethod;
109 using CartesianSolver::updateDomainInfo;
110 using CartesianSolver::windowCell;
111 using CartesianSolver::windowCellId;
112
113
114 MInt noInternalCells() const override { return grid().noInternalCells(); }
115
116 MBool isCompressible() { return m_isCompressible; }
117
119 MBool a_isBndryCell(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::IsBndryCell); }
120
122 maia::lb::cell::BitsetType::reference a_isBndryCell(const MInt cellId) {
123 return a_hasProperty(cellId, SolverCell::IsBndryCell);
124 }
125
127 MBool a_isBndryGhostCell(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::IsGhost); }
128
130 MBool a_isInterface(const MInt cellId) const { return m_cells.hasProperty(cellId, SolverCell::IsInterface); }
131
133 maia::lb::cell::BitsetType::reference a_isInterface(const MInt cellId) {
134 return m_cells.hasProperty(cellId, SolverCell::IsInterface);
135 }
136
138 maia::lb::cell::BitsetType::reference a_onlyBoundary(const MInt cellId) {
139 return a_hasProperty(cellId, SolverCell::OnlyBoundary);
140 }
141
143 MBool a_onlyBoundary(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::OnlyBoundary); }
144
146 MFloat& a_kappa(const MInt cellId) { return m_cells.kappa(cellId); }
147
149 MFloat a_kappa(const MInt cellId) const { return m_cells.kappa(cellId); }
150
152 MFloat& a_diffusivity(const MInt cellId) { return m_cells.diffusivity(cellId); }
153
155 MFloat a_diffusivity(const MInt cellId) const { return m_cells.diffusivity(cellId); }
156
158 MFloat& a_nu(const MInt cellId) { return m_cells.nu(cellId); }
159
161 MFloat& a_oldNu(const MInt cellId) { return m_cells.oldNu(cellId); }
162
164 const MFloat& a_coordinate(const MInt cellId, const MInt dir) const { return grid().tree().coordinate(cellId, dir); }
165
167 MFloat c_coordinate(const MInt cellId, const MInt dir) const { return grid().tree().coordinate(cellId, dir); }
168
170 MFloat& a_distribution(const MInt cellId, const MInt dir) { return m_cells.distributions(cellId, dir); }
171
173 MFloat a_distribution(const MInt cellId, const MInt dir) const { return m_cells.distributions(cellId, dir); }
174
176 MFloat& a_oldDistribution(const MInt cellId, const MInt dir) { return m_cells.oldDistributions(cellId, dir); }
177
179 MFloat a_oldDistribution(const MInt cellId, const MInt dir) const { return m_cells.oldDistributions(cellId, dir); }
180
182 MFloat& a_externalForces(const MInt cellId, const MInt dim) { return m_cells.externalForces(cellId, dim); }
183
185 MFloat a_externalForces(const MInt cellId, const MInt dim) const { return m_cells.externalForces(cellId, dim); }
186
187 MFloat& a_previousDistribution(const MInt cellId, const MInt distr) {
188 return m_cells.previousDistribution(cellId, distr);
189 }
190
191 MFloat a_previousDistribution(const MInt cellId, const MInt distr) const {
192 return m_cells.previousDistribution(cellId, distr);
193 }
194
195 MFloat& a_previousVariable(const MInt cellId, const MInt varId) { return m_cells.previousVariable(cellId, varId); }
196
197 MFloat a_previousVariable(const MInt cellId, const MInt varId) const {
198 return m_cells.previousVariable(cellId, varId);
199 }
200
202 MFloat& a_nuT(const MInt cellId) { return m_cells.nuT(cellId); }
203
205 MFloat a_nuT(const MInt cellId) const { return m_cells.nuT(cellId); }
206
208 MFloat& a_oldNuT(const MInt cellId) { return m_cells.oldNuT(cellId); }
209
211 MFloat a_oldNuT(const MInt cellId) const { return m_cells.oldNuT(cellId); }
212
214 MFloat a_nu(const MInt cellId) const { return m_cells.nu(cellId); }
215
217 MFloat a_oldNu(const MInt cellId) const { return m_cells.oldNu(cellId); }
218
220 maia::lb::cell::BitsetType::reference a_isInterfaceChild(const MInt cellId) {
221 return a_hasProperty(cellId, SolverCell::IsInterfaceChild);
222 }
223
225 MBool a_isInterfaceChild(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::IsInterfaceChild); }
226
228 maia::lb::cell::BitsetType::reference a_isInterfaceParent(const MInt cellId) {
229 return a_hasProperty(cellId, SolverCell::IsInterfaceParent);
230 }
231
233 MBool a_isInterfaceParent(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::IsInterfaceParent); }
234
236 MFloat& a_distributionThermal(const MInt cellId, const MInt dir) { return m_cells.distributionsThermal(cellId, dir); }
237
239 MFloat a_distributionThermal(const MInt cellId, const MInt dir) const {
240 return m_cells.distributionsThermal(cellId, dir);
241 }
242
244 MFloat& a_oldDistributionThermal(const MInt cellId, const MInt dir) {
245 return m_cells.oldDistributionsThermal(cellId, dir);
246 }
247
249 MFloat a_oldDistributionThermal(const MInt cellId, const MInt dir) const {
250 return m_cells.oldDistributionsThermal(cellId, dir);
251 }
252
254 MFloat& a_distributionTransport(const MInt cellId, const MInt dir) {
255 return m_cells.distributionsTransport(cellId, dir);
256 }
257
259 MFloat a_distributionTransport(const MInt cellId, const MInt dir) const {
260 return m_cells.distributionsTransport(cellId, dir);
261 }
262
264 MFloat& a_oldDistributionTransport(const MInt cellId, const MInt dir) {
265 return m_cells.oldDistributionsTransport(cellId, dir);
266 }
267
269 MFloat a_oldDistributionTransport(const MInt cellId, const MInt dir) const {
270 return m_cells.oldDistributionsTransport(cellId, dir);
271 }
272
274 MFloat& a_variable(const MInt cellId, const MInt varId) { return m_cells.variables(cellId, varId); }
275
277 MFloat a_variable(const MInt cellId, const MInt varId) const { return m_cells.variables(cellId, varId); }
278
280 MFloat& a_oldVariable(const MInt cellId, const MInt varId) { return m_cells.oldVariables(cellId, varId); }
281
283 MFloat a_oldVariable(const MInt cellId, const MInt varId) const { return m_cells.oldVariables(cellId, varId); }
284
289 // template <MBool interpolate/* = true */>
290 template <MBool interpolate = true>
291 MFloat a_interpolatedVariable(const MInt cellId, const MInt varId, const MFloat dt = 0.0) const {
292 IF_CONSTEXPR(interpolate) { return (1.0 + dt) * a_variable(cellId, varId) - dt * a_oldVariable(cellId, varId); }
293 else { // this is identical to a_variable
294 return a_variable(cellId, varId);
295 }
296 }
297
299 MBool a_hasProperty(const MInt cellId, const Cell p) const { return grid().tree().hasProperty(cellId, p); }
300
302 maia::lb::cell::BitsetType::reference a_hasProperty(const MInt cellId, const SolverCell p) {
303 return m_cells.hasProperty(cellId, p);
304 }
305
307 MBool a_hasProperty(const MInt cellId, const SolverCell p) const { return m_cells.hasProperty(cellId, p); }
308
310 MInt& a_bndId(const MInt cellId) { return m_cells.bndId(cellId); }
311
313 MInt a_bndId(const MInt cellId) const { return m_cells.bndId(cellId); }
314
315 // sampling data methods
316 MInt getCellIdByIndex(const MInt index) { return index; }
317 virtual void initInterpolationForCell(const MInt NotUsed(cellId)){};
318
319 MFloat c_cellLengthAtLevel(const MInt level) const { return grid().cellLengthAtLevel(level); }
320
321 MFloat c_cellLengthAtCell(const MInt cellId) const { return grid().cellLengthAtLevel(a_level(cellId)); }
322
324 MInt& c_neighborId(const MInt cellId, const MInt dir) {
325 // Note: use neighbor list instead of grid().m_tree, also includes diagonal neighbors
326 return grid().neighborList(cellId, dir);
327 }
328
330 MInt c_neighborId(const MInt cellId, const MInt dir) const {
331 // Note: use neighbor list instead of grid().m_tree, also includes diagonal neighbors
332 return grid().neighborList(cellId, dir);
333 }
334
336 MInt a_hasNeighbor(const MInt cellId, const MInt dir) const {
337 // Note: use neighbor list instead of grid().m_tree, also includes diagonal neighbors
338 return static_cast<MInt>(grid().neighborList(cellId, dir) > -1);
339 }
340
341 MBool a_isHalo(const MInt cellId) const { return m_cells.hasProperty(cellId, SolverCell::IsHalo); }
342
343 maia::lb::cell::BitsetType::reference a_isHalo(const MInt cellId) {
344 // ENSURE_VALID_GCELLID(gCellId);
345 return m_cells.hasProperty(cellId, SolverCell::IsHalo);
346 }
347
348 MBool a_isWindow(const MInt cellId) const { return m_cells.hasProperty(cellId, SolverCell::IsWindow); }
349
350 maia::lb::cell::BitsetType::reference a_isWindow(const MInt cellId) {
351 // ENSURE_VALID_GCELLID(gCellId);
352 return m_cells.hasProperty(cellId, SolverCell::IsWindow);
353 }
354
355 void a_resetPropertiesSolver(const MInt cellId) { m_cells.resetProperties(cellId); }
356
357 MFloat* a_variables_ptr(MInt pCellId) { return &a_variable(pCellId, 0); }
358
359 MFloat* a_oldVariables_ptr(MInt pCellId) { return &a_oldVariable(pCellId, 0); }
360
361 MInt c_globalId(const MInt cellId) const { return grid().tree().globalId(cellId); }
362
363 MInt a_level(const MInt cellId) const { return m_cells.level(cellId); }
364
365 MInt& a_level(const MInt cellId) { return m_cells.level(cellId); }
366
367 MFloat a_alphaGas(const MInt cellId) const { return m_cells.invVolumeFraction(cellId); }
368
369 MFloat a_alphaGasLim(const MInt cellId) const { return mMax(mMin(m_cells.invVolumeFraction(cellId), 1.0), 0.0); }
370
371 MFloat& a_alphaGas(const MInt cellId) { return m_cells.invVolumeFraction(cellId); }
372
373 MFloat a_uOtherPhase(const MInt cellId, const MInt dir) const { return m_cells.uOtherPhase(cellId, dir); }
374
375 MFloat& a_uOtherPhase(const MInt cellId, const MInt dir) { return m_cells.uOtherPhase(cellId, dir); }
376
377 MInt getIdAtPoint(const MFloat* point, [[maybe_unused]] MBool globalUnique = false) {
378 return grid().findContainingLeafCell(point);
379 }
380
382 MInt c_noCells() const { return grid().tree().size(); }
383
384 MInt c_level(const MInt cellId) const { return grid().tree().level(cellId); }
385
386 MInt c_noChildren(const MInt cellId) const { return grid().tree().noChildren(cellId); }
387
388 MBool c_isLeafCell(const MInt cellId) const { return grid().tree().isLeafCell(cellId); }
389
390 MInt c_childId(const MInt cellId, const MInt childNumber) const { return grid().tree().child(cellId, childNumber); }
391
392 MInt c_parentId(const MInt cellId) const { return grid().tree().parent(cellId); }
393
394 // Coupling Level-Set Accessors
395 //-------------------------------------------------------------------------
397 MFloat& a_levelSetFunctionMB(const MInt cellId, const MInt set) {
398#ifdef WAR_NVHPC_PSTL
399 const MInt pos = (cellId * m_noLevelSetsUsedForMb) + set;
400 return m_levelSetValues[pos];
401#else
402 return m_levelSetValues[IDX_LSSETMB(cellId, set)];
403#endif
404 }
405
407 MFloat a_levelSetFunctionMB(const MInt cellId, const MInt set) const {
408#ifdef WAR_NVHPC_PSTL
409 const MInt pos = (cellId * m_noLevelSetsUsedForMb) + set;
410 return m_levelSetValues[pos];
411#else
412 return m_levelSetValues[IDX_LSSETMB(cellId, set)];
413#endif
414 }
415
417 MInt& a_associatedBodyIds(const MInt cellId, const MInt set) {
418#ifdef WAR_NVHPC_PSTL
419 const MInt pos = (cellId * m_noLevelSetsUsedForMb) + set;
420 return m_associatedBodyIds[pos];
421#else
422 return m_associatedBodyIds[IDX_LSSETMB(cellId, set)];
423#endif
424 }
425
427 MInt a_associatedBodyIds(const MInt cellId, const MInt set) const {
428#ifdef WAR_NVHPC_PSTL
429 const MInt pos = (cellId * m_noLevelSetsUsedForMb) + set;
430 return m_associatedBodyIds[pos];
431#else
432 return m_associatedBodyIds[IDX_LSSETMB(cellId, set)];
433#endif
434 }
435
437 MBool& a_isG0CandidateOfSet(const MInt cellId, const MInt set) { return m_isG0CandidateOfSet[cellId][set]; }
438
440 MBool a_isG0CandidateOfSet(const MInt cellId, const MInt set) const { return m_isG0CandidateOfSet[cellId][set]; }
441
443 MBool a_isActive(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::IsActive); }
444
446 maia::lb::cell::BitsetType::reference a_isActive(const MInt cellId) {
447 return a_hasProperty(cellId, SolverCell::IsActive);
448 }
449
451 MBool a_wasActive(const MInt cellId) const { return a_hasProperty(cellId, SolverCell::WasActive); }
452
454 maia::lb::cell::BitsetType::reference a_wasActive(const MInt cellId) {
455 return a_hasProperty(cellId, SolverCell::WasActive);
456 }
457 // end Coupling Level-Set Accessors
458 //-------------------------------------------------------------------------
459
460
461 // Swap fields (TODO labels:LB only exchange pointers)
462 void swap_variables(MInt pCellId) {
463#ifdef WAR_NVHPC_PSTL
464 for(MInt varId = 0; varId < m_noVariables; ++varId) {
465#else
466 for(MInt varId = 0; varId < noVariables(); ++varId) {
467#endif
468 std::swap(a_variable(pCellId, varId), a_oldVariable(pCellId, varId));
469 }
470 }
471
472 LbSolver(MInt id, MInt dist, GridProxy& gridProxy_, Geometry<nDim>& geometry_, const MPI_Comm comm);
473 ~LbSolver() override;
474 MBool isActive() const override { return grid().isActive(); }
475 void updateCellCollectorFromGrid();
476 MInt getCurrentTimeStep() const override { return globalTimeStep; }
477 virtual void activateAllCells(){};
478 virtual void activateInnerCells();
479 virtual void activateWindowCells();
480 virtual void activateAllButHaloCells();
481 virtual void propagation_step() = 0;
482 virtual void propagation_step_vol() = 0;
483 virtual void propagation_step_thermal() = 0;
484 virtual void propagation_step_thermal_vol() = 0;
485 virtual void propagation_step_transport() = 0;
486 virtual void propagation_step_transport_vol() = 0;
487 virtual void propagation_step_thermaltransport() = 0;
488 virtual void propagation_step_thermaltransport_vol() = 0;
489
490 virtual void volumeForces() = 0;
491
492 virtual void exchange(MInt mode = 1);
493 virtual void exchangeLb(MInt mode);
494 virtual void exchangeLbNB(MInt mode);
495 virtual void exchangeOldDistributions();
496
497 void (LbSolver::*m_sendMethod)();
498 void (LbSolver::*m_receiveMethod)();
499 virtual void sendNormal();
500 virtual void receiveNormal();
501 virtual void sendReduced();
502 virtual void receiveReduced();
503
504 // Member functions needed for adaptation
505 virtual void removeChildsLb(const MInt parentId) = 0;
506 virtual void refineCellLb(const MInt parentId, const MInt* childIds) = 0;
507 virtual void restartBndCnd() = 0;
508 virtual void prolongation() = 0;
509 void writeInfo();
510 void initializeRefinedCellsPerLevel();
511 void initializeNewInterfaceParents();
512 void saveOutput();
513 void computeFFTStatistics();
514 void saveAdaptedGridFile(MInt* const p_recalcCellIds);
515
516 void (LbSolver::*m_gatherMethod)();
517 void (LbSolver::*m_scatterMethod)();
518 inline void gatherNormal();
519 inline void scatterNormal();
520 inline void gatherReduced();
521 inline void scatterReduced();
522
523 void printCommunicationMethod();
524 // Coupling Level-Set Functions
525 //-------------------------------------------------------------------------
526 void exchangeG0Candidates();
527 void preCoupleLs(std::vector<MInt>& maxGCellLevels);
528 void createBndryToBodyMapping(maia::coupling::Mapping& bndryToBodyMapping,
529 maia::coupling::Mapping& bodyToBndryMapping);
530 void findG0Candidates(std::vector<MInt>& maxGCellLevels);
531 MInt setUpLbInterpolationStencil(const MInt cellId, MInt* interpolationCells, MFloat* point);
532 MFloat interpolateFieldDataLb(MInt*, MFloat*, MInt varId, std::function<MFloat(MInt, MInt)> scalarField,
533 std::function<MFloat(MInt, MInt)> coordinate);
534 void initializeMovingBoundaries();
535 // End Coupling Level-Set Functions
536 //-------------------------------------------------------------------------
537
538 void resetActiveCellList(MInt mode = 0);
539 void resetExternalSources();
540 void exchangeExternalSources();
541
542 void prepareCommunication();
543 void prepareCommunicationNormal();
544 void prepareCommunicationReduced();
545 void markCellsForAdditionalComm();
546 void resetComm();
547 void resetCellLists(MBool resize = true);
548 void loadRestartFile() override;
549 virtual void saveRestartFile();
550 void copyGridProperties();
551 void returnCellInfo(MInt cellId);
552
553 void initSolver() override;
554 void finalizeInitSolver() override;
555 void preTimeStep() override;
556 MBool solutionStep() override;
557 void saveSolverSolution(MBool forceOutput, const MBool finalTimeStep) override;
558 virtual MBool maxResidual() {
559 m_log << "Called LbSolver::maxResidual without Implementation" << std::endl;
560 return true;
561 }
562 void outputInitSummary();
563
564 virtual void cleanUp(){};
565
566 void postTimeStep() {} // until now it is left unused
567
568 void storeOldDistributions();
569 void storeOldVariables();
570 void initNu(const MInt cellId, const MFloat nu);
571
572 protected:
573 virtual void clb_collision_step() = 0;
574 virtual void clb_smagorinsky_collision_step() = 0;
575 virtual void cumulant_collision_step() = 0;
576 virtual void bgki_collision_step(){};
577 virtual void bgki_collision_step_mb(){};
578 virtual void bgki_collision_step_mb_thermal(){};
579 virtual void bgki_collision_step_Guo_forcing(){};
580 virtual void bgki_init_collision_step(){};
581 virtual void bgkc_collision_step(){};
582 virtual void bgki_smagorinsky_collision_step(){};
583 virtual void bgki_smagorinsky_collision_step2(){};
584 virtual void bgki_smago_wall_collision_step(){};
585 virtual void bgki_dynamic_smago_collision_step(){};
586 virtual void bgki_euler_collision_step(){};
587 virtual void bgki_thermal_collision_step(){};
588 virtual void bgki_innerEnergy_collision_step(){};
589 virtual void bgki_totalEnergy_collision_step(){};
590 virtual void bgkc_transport_collision_step(){};
591 virtual void bgkc_thermal_transport_collision_step(){};
592 virtual void bgkc_innerenergy_transport_collision_step(){};
593 virtual void bgkc_totalenergy_transport_collision_step(){};
594 virtual void mrt_collision_step() = 0;
595 virtual void mrt2_collision_step() = 0;
596 virtual void mrt_smagorinsky_collision_step() = 0;
597 virtual void rbgk_collision_step() = 0;
598 virtual void rbgk_dynamic_smago_collision_step() = 0;
599 virtual void rbgk_smagorinsky_collision_step() = 0;
600
601 virtual void initializeLatticeBgk(){};
602 virtual void initializeLatticeBgkThermal(){};
603 virtual void initializeLatticeBgkTransport(){};
604 virtual void initializeLatticeBgkThermalTransport(){};
605 virtual void restartInitLb(){};
606
607 virtual void initSrcTermController() = 0;
608 virtual void initSrcTerms() = 0;
609 virtual void preCollisionSrcTerm() = 0;
610 virtual void postCollisionSrcTerm() = 0;
611 virtual void postPropagationSrcTerm() = 0;
612 virtual void postCollisionBc() = 0;
613 virtual void postPropagationBc() = 0;
614 virtual void updateVariablesFromOldDist() = 0;
615 virtual void updateVariablesFromOldDist_preCollision() = 0;
616 virtual void updateViscosity() = 0;
617
618 virtual void calcNodalLsValues(){};
619
620 inline void setIsActiveDefaultStates();
621 inline void setInActiveBndryCells();
622 inline void setInActiveMBCells();
623 inline void fillActiveCellList();
624 void getReLambdaAndUrmsInit();
625
626 // sampling variables
627 MBool m_isInitSamplingVars = false;
628 std::vector<MFloat**> m_samplingVariables{};
629
630 private:
631 void (LbSolver::*m_propagationStepMethod)();
632 void (LbSolver::*m_solutionStepMethod)();
633 void (LbSolver::*m_initializeMethod)();
634 void (LbSolver::*m_restartInitLbMethod)();
635
636 // Function pointer for the exchange of data
637 void (LbSolver::*m_exchangeMethod)(MInt mode);
638
639 void initTimer();
640 void averageTimer();
641
642 // TODO labels:LB move this function to something like src/UTIL/numerics.h ?
650 void getDerivativeStencilAndCoefficient(std::array<MInt, 5>& cellIds, std::array<MFloat, 5>& coef) {
651 // ?-?-x-?-?
652 if(cellIds[0] > -1) {
653 // x-x-x-?-?
654 if(cellIds[4] > -1) {
655 // x-x-0-x-x -> central 4th order
656 coef[0] = F1B12;
657 coef[1] = -F2B3;
658 coef[2] = 0.0;
659 coef[3] = F2B3;
660 coef[4] = -F1B12;
661 } else if(cellIds[3] > -1) {
662 // x-x-x-x-0 -> backward + 1 forward difference 3rd order
663 coef[0] = F1B6;
664 coef[1] = -1.0;
665 coef[2] = F1B2;
666 coef[3] = F1B3;
667 coef[4] = 0.0;
668 } else {
669 // x-x-x-0-0 -> backward difference 2nd order
670 coef[0] = F1B2;
671 coef[1] = -2.0;
672 coef[2] = F3B2;
673 coef[3] = 0.0;
674 coef[4] = 0.0;
675 }
676 } else if(cellIds[4] > -1) {
677 // 0-?-x-x-x
678 if(cellIds[1] > -1) {
679 // 0-x-x-x-x -> forward + 1 backward difference 3rd order
680 coef[0] = 0.0;
681 coef[1] = -F1B3;
682 coef[2] = -F1B2;
683 coef[3] = 1.0;
684 coef[4] = -F1B6;
685 } else {
686 // 0-0-x-x-x -> forward difference 2nd order
687 coef[0] = 0.0;
688 coef[1] = 0.0;
689 coef[2] = -F3B2;
690 coef[3] = 2.0;
691 coef[4] = -F1B2;
692 }
693 } else {
694 // 0-?-x-?-0
695 if(cellIds[1] > -1) {
696 // 0-x-x-?-0
697 if(cellIds[3] > -1) {
698 // 0-x-0-x-0 -> central difference 2nd order
699 coef[0] = 0.0;
700 coef[1] = -F1B2;
701 coef[2] = 0.0;
702 coef[3] = F1B2;
703 coef[4] = 0.0;
704 } else {
705 // 0-x-x-0-0 -> backward difference 1st order
706 coef[0] = 0.0;
707 coef[1] = -1.0;
708 coef[2] = 1.0;
709 coef[3] = 0.0;
710 coef[4] = 0.0;
711 }
712 } else if(cellIds[3] > -1) {
713 // 0-0-x-x-0 -> forward difference 1st order
714 coef[0] = 0.0;
715 coef[1] = 0.0;
716 coef[2] = -1.0;
717 coef[3] = 1.0;
718 coef[4] = 0.0;
719 } else {
720 // 0-0-x-0-0
721 // Really not possible -> all coef == 0.0 -> no gradient here
722 coef[0] = 0.0;
723 coef[1] = 0.0;
724 coef[2] = 0.0;
725 coef[3] = 0.0;
726 coef[4] = 0.0;
727 }
728 }
729 }
730
731 void printScalingVariables();
732
733 public:
734 MBool m_calculateDissipation;
735
736 MString m_initMethod;
737
738 // restart
739 MBool m_initFromCoarse;
740
741 //(non-blocking) communication
742 MFloat** m_sendBuffers = nullptr;
743 MFloat** m_receiveBuffers = nullptr;
744 MInt** m_nghbrOffsetsWindow = nullptr;
745 MInt** m_nghbrOffsetsHalo = nullptr;
746 MFloat** m_baseAddresses = nullptr;
747 MInt m_noVarsTransfer{};
748 MInt m_noElementsTransfer{};
749 MInt m_noDistsTransfer{};
750 MInt m_dataBlockSizeTotal{};
751 MInt* m_dataBlockSizes = nullptr;
752
753 // reduced communication
754 MBool m_reducedComm{};
755 MInt* m_noWindowDistDataPerDomain{};
756 MInt* m_noHaloDistDataPerDomain{};
757 MFloat*** m_commPtWindow{};
758 MFloat*** m_commPtHalo{};
759 MInt* m_orderedNeighbors{};
760 MInt* m_needsFurtherExchange{};
761 MInt** m_windowDistsForExchange{};
762 MInt** m_haloDistsForExchange{};
763
764 MBool m_isRefined = false;
765 MFloat m_interfaceCellSize{};
766 MBool nonBlockingComm() const { return m_nonBlockingComm; }
767
768 // Load balancing
769 MBool m_setCellDataFinished = false;
770
771 // Coupling Level-Set
772 MInt m_currentNoG0Cells;
773 MFloat** m_nodalGValues{};
774 MInt* m_G0CellList{};
775 MInt* m_G0CellMapping{};
776 MFloat** m_G0WallDist{};
777 std::vector<CutCandidate<nDim>> m_G0Candidates;
778 MInt* m_sendBufferMB{};
779 MInt* m_receiveBufferMB{};
780 MInt m_maxNoSets = -1;
781 MInt m_levelSetId = 0;
782 MFloat* m_levelSetValues{};
783 MInt* m_associatedBodyIds{};
784 MBool** m_isG0CandidateOfSet{};
785 MInt m_noLevelSetsUsedForMb;
786 MFloat m_maxBodyRadius;
787 MInt m_noEmbeddedBodies;
788 MBool m_constructGField;
789 MBool m_useOnlyCollectedLS = false;
790 MBool m_allowBndryAsG0 = false;
791
792 MFloat* m_bodyVelocity = nullptr;
793 MBool m_trackMovingBndry;
794 MInt m_trackMbEnd;
795 MInt m_trackMbStart;
796 MInt m_noG0CandidatesTotal;
797 MInt* m_initialActiveCells = nullptr;
798
799 // Residual
800 MInt m_resTimestepPos;
801 MInt m_resRePos;
802 std::ofstream mRes; // output stream for the residual file
803
804 MInt m_solutionOffset;
805 MFloat m_referenceLength;
806 MFloat m_referenceLengthSTL;
807
808 MFloat m_noOuterBndryCells;
809
810 // Non-Newtonian
811 MFloat m_n = F0;
812 MBool m_nonNewtonian = false;
813
815 MInt noVariables() const override { return m_noVariables; }
816 MInt noDistributions() const { return m_noDistributions; }
817 MFloat time() const override { return m_time; }
818
819 using Geom = Geometry<nDim>;
820
822 constexpr const Geom& geometry() const { return *m_geometry; }
823
824 GeometryIntersection<nDim>* m_geometryIntersection;
825
826
827 protected:
828 MBool m_nonBlockingComm = false;
830 Geom& geometry() { return *m_geometry; }
831
832 private:
833 // Number of variables
834 MInt m_noVariables;
835
836
837 public:
845 Collector<PointBasedCell<nDim>>* m_extractedCells = nullptr;
846 Collector<CartesianGridPoint<nDim>>* m_gridPoints = nullptr;
847
849
850 MInt a_noCells() const { return m_cells.size(); }
851
852 MFloat a_Nu() const { return m_nu; }
853 MInt a_FFTInterval() const { return m_fftInterval; }
854
855 private:
856 MPI_Request* mpi_request{};
857 MPI_Request* mpi_requestR{};
858 MPI_Request* mpi_requestS{};
859
860 protected:
861 Geometry<nDim>* m_geometry;
863 maia::lb::collector::LbCellCollector<nDim> m_cells;
864
865 // Pointer to the class MConservativeVariables
866 MConservativeVariables<nDim>* CV;
867 // Pointer to the class MPrimitiveVariables
868 MPrimitiveVariables<nDim>* PV;
869
870 // domain dimensions
871 MInt m_arraySize[nDim]; // in cell units
872 MFloat* m_domainBoundaries{}; // in macroscopic units
873
874 LbBndCnd<nDim>* m_bndCnd;
875
876 // subgrid properties
877 MFloat m_Cs;
878 MFloat m_deltaX;
879
880 // time (non-dimensional by using L=1m and u=u_infinity)
881 MFloat m_time = F0;
882 MFloat m_dt = -1.0;
883
884 MFloat m_densityGradient{};
885
886 MFloat m_domainLength;
887 MInt maxTwoPower{};
888 MInt m_currentMaxNoCells{};
889 MInt m_initialNoCells[3]{};
890
891 MBool m_densityFluctuations;
892 MInt m_noPeakModes;
893 MBool m_FftInit;
894
895 // FFT output and number of timesteps before wrinting next FFT output
896 MInt m_fftInterval;
897
898 // Defines the discretization Model
899 MInt m_noDistributions;
900 MInt m_methodId;
901
902 // Viscosity
903 MFloat m_nu;
904
905 // for BCs
906 MFloat m_rho1;
907 MFloat m_rho2;
908
909 // pulsatile
910 MFloat m_alpha;
911 MFloat m_pulsatileFrequency;
912
913 MBool m_isInitRun = false;
914
915 MBool m_tanhInit;
916 MFloat m_initRe;
917 MInt m_initTime;
918 MInt m_initStartTime;
919 MFloat m_finalRe;
920 MFloat m_tanhScaleFactor;
921
922 // for isotropic turbulence
923 MFloat m_UrmsInit;
924 MFloat m_ReLambda;
925
926 // Relaxation factor
927 MFloat m_omega;
928
929 // external Forcing terms
930 MFloat* m_Fext{};
931
932 MFloat m_Ga{};
933
934 // dimension independent data
936 MInt* m_startLevel{};
937
938 // Use an external pressure force
939 MBool m_externalForcing;
940
941 MBool m_particleMomentumCoupling;
942
943 MBool m_saveExternalForces;
944
945 MBool m_updateAfterPropagation;
946 LbUpdateMacroscopic m_updateMacroscopicLocation = INCOLLISION;
947 MBool m_initDensityGradient;
948 MFloat m_volumeAccel[3] = {};
949 MFloat m_volumeAccelBase[3] = {};
950 MBool m_controlVelocity = false;
951 struct {
952 MInt dir;
953 MInt interval;
954 MFloat KT;
955 MFloat KI;
956 MFloat KD;
957 MFloat lastGlobalAvgV = -1.0;
958 MFloat integratedError = 0.0;
959 MFloat derivedError = 0.0;
960 MFloat previousError = 0.0;
961 MBool restart = false;
962 } m_velocityControl;
963
964 MBool m_isEELiquid;
965 struct {
966 MBool restartWithoutAlpha;
967 MFloat alphaInf;
968 MFloat initialAlpha;
969 MBool gravity;
970 MFloat gravityAccelM[3];
971 MFloat* Fg{};
972 } m_EELiquid;
973
974 MFloat* m_residual = nullptr;
975 MFloat* m_tmpResidual = nullptr;
976 MInt* m_tmpResidualLvl = nullptr;
977 MInt* m_maxResId = nullptr;
978 MFloat** m_rescoordinates = nullptr;
979 MString m_resFileName = "";
980
981 MFloat** m_momentumFlux = nullptr;
982 MFloat* m_MijMij = nullptr;
983 MFloat* m_MijLij = nullptr;
984
985 MFloat m_ReTau;
986
987 // the id of the segment used to get m_referenceLength
988 MInt m_referenceLengthSegId;
989 MBool m_refineDiagonals = true;
990
991 // Thermal variables
992 MFloat m_Pr;
993 MFloat m_omegaT;
994 MFloat m_kappa;
995 MFloat m_initTemperatureKelvin;
996 MFloat m_blasiusPos;
997
998 MBool m_isCompressible = false;
999 MInt m_isThermal = 0;
1000 MInt m_innerEnergy = 0;
1001 MInt m_totalEnergy = 0;
1002 MFloat m_CouettePoiseuilleRatio = F0;
1003 MInt m_calcTotalPressureGradient = 0;
1004
1005 // Transport variables
1006 MFloat m_Pe;
1007 MFloat m_omegaD;
1008 MFloat m_diffusivity;
1009 MFloat m_initCon;
1010 MInt m_isTransport = 0;
1011
1012 MFloat* m_cellLength = nullptr;
1013
1014 virtual void treatInterfaceCells();
1015 void correctInterfaceBcCells();
1016 virtual void setActiveCellList();
1017
1019 virtual void initializeInterfaceCells();
1020 virtual void buildInterfaceCells();
1021 virtual void resetInterfaceCells();
1022
1023 LbInterpolationType m_interpolationType;
1024 virtual void setInterpolationNeighbors();
1025 virtual void setInterpolationCoefficients();
1026 void setInterpolationNeighborsBC();
1027 void setInterpolationCoefficientsBC();
1028
1029 std::vector<Collector<LbInterfaceCell>*> m_interfaceChildren;
1030 std::vector<Collector<LbParentCell>*> m_interfaceParents;
1031
1032 MInt* m_noInterfaceChildren = nullptr;
1033 MInt* m_noInterfaceParents = nullptr;
1034
1035 MBool m_correctInterfaceBcCells = false;
1036 struct InterfaceLink {
1037 MInt levelId;
1038 MInt interfaceId;
1039 };
1040 std::vector<InterfaceLink> m_interfaceChildrenBc;
1041
1042 MBool m_saveDerivatives;
1043 MBool m_initRestart;
1044
1045 // The number of dimensions
1046 MInt* m_activeCellList = nullptr;
1047 MInt* m_activeCellListLvlOffset = nullptr;
1048
1049 protected:
1050 virtual void writeGridToTecFile(const MChar* fileName);
1051 virtual void writeGridToVtkFile(const MChar* fileName);
1052 virtual void writeGridToVtkFileAscii(const MChar* fileName);
1053 virtual void initPressureForce() = 0;
1054 virtual void initVolumeForces() = 0;
1055 virtual void initRunCorrection() = 0;
1056 void determineEqualDiagonals();
1057
1058 public: // parallel IO
1059 void saveUVWRhoTOnlyPar(const MChar* fileName, const MChar* gridInputFileName, MInt* recalcIdTree = nullptr);
1060 void saveRestartWithDistributionsPar(const MChar* fileName, const MChar* gridInputFileName,
1061 MInt* recalcIdTree = nullptr);
1062 void loadRestartWithDistributionsPar(const MChar* fileName);
1063 void loadRestartWithoutDistributionsPar(const MChar* fileName);
1064 void loadRestartWithoutDistributionsParFromCoarse(const MChar* fileName);
1065 // derivatives
1066 void calculateVelocityDerivative(const MInt cellId, MFloat (&gradient)[nDim][nDim]);
1067 template <MBool tI = false> // interpolation in time
1068 MFloat calculateDerivative(const MInt cellId, const MInt velComp, const MInt spaceDir, const MFloat dt = 0.0);
1069 MFloat calculateInterpolatedDerivative(const MInt cellId, const MInt velComp, const MInt spaceDir, const MFloat dt) {
1070 return calculateDerivative<true>(cellId, velComp, spaceDir, dt);
1071 }
1072 inline MFloat calculatePressureDerivative(MInt cellId, MInt spaceDir);
1073 inline void getStrainTensor(MFloatScratchSpace& derivatives, MInt cellId, MFloatScratchSpace& strain);
1074 inline void getVorticityTensor(MFloatScratchSpace& derivatives, MInt cellId, MFloatScratchSpace& strain);
1075
1076 // postprocessing
1077 void accessSampleVariables(MInt cellId, MFloat*& vars);
1078 void interpolateToParentCells(MInt parentlevel);
1079 void saveCoarseSolution();
1080
1081 void createMPIComm(const MInt* ranks, MInt noRanks, MPI_Comm* comm);
1082 inline void writeSegmentBoundaryVTK(MFloat** bndVs, MInt num);
1083 void updateTime();
1084
1085 protected:
1086 virtual MFloat calculateReferenceLength(MInt segmentId);
1087 inline MFloat calcCharLenAll(MInt segmentId);
1088 inline MFloat calcCharLenParallelSplit(MInt segmentId);
1089
1090 public:
1091 virtual void getSolverSamplingProperties(std::vector<MInt>& samplingVars, std::vector<MInt>& noSamplingVars,
1092 std::vector<std::vector<MString>>& samplingVarNames,
1093 const MString featureName = "");
1094 virtual void initSolverSamplingVariables(const std::vector<MInt>& varIds, const std::vector<MInt>& noSamplingVars);
1095 virtual void calcSamplingVariables(const std::vector<MInt>& varIds, const MBool exchange);
1096 void interpolateVariablesInCell(const MInt cellId, const MFloat* position, MFloat* interpolationResult);
1097
1098 void calcSamplingVarAtPoint(const MFloat* point, const MInt id, const MInt sampleVarId, MFloat* state,
1099 const MBool interpolate = false);
1100 void loadSampleVariables(MInt timeStep);
1101 void getSampleVariables(MInt cellId, const MFloat*& vars);
1102 void getSampleVariables(const MInt cellId, std::vector<MFloat>& vars);
1103 MBool getSampleVarsDerivatives(const MInt cellId, std::vector<MFloat>& vars);
1104 void calculateHeatRelease() { std::cerr << "calculateHeatRelease LbSolver " << std::endl; }
1105 void getHeatRelease(MFloat*& heatRelease) { std::cerr << "getHeatRelease LbSolver " << heatRelease << std::endl; }
1106
1107 // adaptation
1108 void prepareAdaptation(std::vector<std::vector<MFloat>>& NotUsed(sensors),
1109 std::vector<MFloat>& NotUsed(sensorWeight),
1110 std::vector<std::bitset<64>>& NotUsed(sensorCellFlag),
1111 std::vector<MInt>& NotUsed(sensorSolverId)) override{};
1112 void prepareAdaptation() override;
1113 void setSensors(std::vector<std::vector<MFloat>>& sensors, std::vector<MFloat>& sensorWeight,
1114 std::vector<std::bitset<64>>& sensorCellFlag, std::vector<MInt>& sensorSolverId) override;
1115 void reinitAfterAdaptation() override{};
1116 void postAdaptation() override;
1117 void finalizeAdaptation() override;
1118 void removeChilds(const MInt) override;
1119 void removeCell(const MInt) override;
1120 void refineCell(const MInt) override;
1121 void swapCells(const MInt, const MInt) override;
1122 void swapProxy(const MInt cellId0, const MInt cellId1) override;
1123 void resizeGridMap() override;
1124 MBool prepareRestart(MBool writeRestart, MBool& writeGridRestart) override;
1125 void reIntAfterRestart(MBool doneRestart);
1126 void writeRestartFile(const MBool writeRestart, const MBool writeBackup, const MString gridFileName,
1127 MInt* recalcIdTree) override;
1128
1129 // Dynamic load balancing
1130 void balance(const MInt* const /*noCellsToReceiveByDomain*/, const MInt* const /*noCellsToSendByDomain*/,
1131 const MInt* const /*targetDomainsByCell*/, const MInt /*oldNoCells*/) override {
1132 TERMM(1, "Use split balancing methods for LB solver instead of balance().");
1133 };
1134 void setCellWeights(MFloat* solverCellWeight) override;
1135 void resetSolver() override{};
1136 void finalizeBalance() override{};
1137
1138 // DLB loads
1139 MInt noLoadTypes() const override;
1140 void getLoadQuantities(MInt* const loadQuantities) const override;
1141 void getDefaultWeights(MFloat* weights, std::vector<MString>& names) const override;
1142 MFloat getCellLoad(const MInt cellId, const MFloat* const weights) const override;
1143
1144 // DLB (Split balance)
1145 MBool hasSplitBalancing() const override { return true; }
1146 void balancePre() override;
1147 void balancePost() override;
1148
1149 // DLB exchange information (Split balance)
1150
1151 // Number of fields to exchange during DLB
1152 MInt noCellDataDlb() const override {
1153 if(grid().isActive()) {
1154 return noSolverCellData();
1155 } else {
1156 return 0;
1157 }
1158 };
1159
1160 // Datatype (enum) for a given field
1161 MInt cellDataTypeDlb(const MInt dataId) const {
1162 TRACE();
1163
1164 if(!isActive()) {
1165 TERMM(1, "Error: cellDataTypeDlb() might give wrong results on inactive ranks.");
1166 return -1;
1167 }
1168
1169 MInt dataType = -1;
1170 if(dataId > -1 && dataId < noSolverCellData()) {
1171 // LB solver cell data
1172 dataType = s_cellDataTypeDlb[dataId];
1173 } else {
1174 TERMM(1, "The requested dataId is not valid: " + std::to_string(dataId) + " ("
1175 + std::to_string(noSolverCellData()) + ", " + std::to_string(noCellDataDlb()) + ")");
1176 }
1177
1178 return dataType;
1179 }
1180
1181 // Data size per cell
1182 MInt cellDataSizeDlb(const MInt dataId, const MInt gridCellId) {
1183 TRACE();
1184
1185 // Inactive ranks do not have any data to communicate
1186 if(!isActive()) {
1187 return 0;
1188 }
1189
1190 // Convert to solver cell id and check
1191 const MInt cellId = grid().tree().grid2solver(gridCellId);
1192 if(cellId < 0) {
1193 return 0;
1194 }
1195
1196 MInt dataSize = 0;
1197
1198 if(dataId > -1 && dataId < noSolverCellData()) {
1199 // LB solver cell data
1200 switch(dataId) {
1201 case CellData::VARIABLES:
1202 case CellData::OLD_VARIABLES:
1203 dataSize = noVariables();
1204 break;
1205 case CellData::DISTRIBUTIONS:
1206 case CellData::OLD_DISTRIBUTIONS:
1207 dataSize = noDistributions();
1208 break;
1209 case CellData::NU:
1210 case CellData::IS_ACTIVE:
1211 dataSize = 1;
1212 break;
1213 default:
1214 TERMM(1, "Unknown data id. (" + std::to_string(dataId) + ")");
1215 break;
1216 }
1217 } else {
1218 TERMM(1, "The requested dataId is not valid.");
1219 }
1220
1221 return dataSize;
1222 }
1223
1224 // Copy cell data field to buffer (int version)
1225 void getCellDataDlb(const MInt dataId, const MInt oldNoCells, const MInt* const bufferIdToCellId,
1226 MInt* const data) override {
1227 getCellDataDlb_(dataId, oldNoCells, bufferIdToCellId, data);
1228 };
1229
1230 // Copy cell data field to buffer (float version)
1231 void getCellDataDlb(const MInt dataId, const MInt oldNoCells, const MInt* const bufferIdToCellId,
1232 MFloat* const data) override {
1233 getCellDataDlb_(dataId, oldNoCells, bufferIdToCellId, data);
1234 };
1235
1236 // Copy buffer data to data field (int version)
1237 void setCellDataDlb(const MInt dataId, const MInt* const data) override { setCellDataDlb_(dataId, data); };
1238
1239 // Copy buffer data to data field (float version)
1240 void setCellDataDlb(const MInt dataId, const MFloat* const data) override { setCellDataDlb_(dataId, data); };
1241
1251 template <typename DataType>
1252 void setCellDataDlb_(const MInt dataId, const DataType* const data) {
1253 TRACE();
1254
1255#ifdef LBCOLLECTOR_SOA_MEMORY_LAYOUT
1256 TERMM(1, "Balance data accessor not defined for SOA cell collector");
1257#endif
1258
1259 // This function can be called multiple times, prevent override
1260 if(m_setCellDataFinished) {
1261 return;
1262 }
1263
1264 // Nothing to do if solver is not active
1265 if(!grid().isActive()) {
1266 return;
1267 }
1268
1269 const MInt noCells = m_cells.size();
1270
1271 if(dataId > -1 && dataId < noSolverCellData()) {
1272 // LB solver cell data
1273 switch(dataId) {
1274 case CellData::VARIABLES:
1275 std::copy_n(data, noCells * noVariables(), &m_cells.variables(0, 0));
1276 break;
1277 case CellData::OLD_VARIABLES:
1278 std::copy_n(data, noCells * noVariables(), &m_cells.oldVariables(0, 0));
1279 break;
1280 case CellData::DISTRIBUTIONS:
1281 std::copy_n(data, noCells * noDistributions(), &m_cells.distributions(0, 0));
1282 break;
1283 case CellData::OLD_DISTRIBUTIONS:
1284 std::copy_n(data, noCells * noDistributions(), &m_cells.oldDistributions(0, 0));
1285 break;
1286 case CellData::NU:
1287 std::copy_n(data, noCells, &m_cells.nu(0));
1288 break;
1289 case CellData::IS_ACTIVE:
1290 for(MInt i = 0; i < noCells; i++) {
1291 a_isActive(i) = (MInt)data[i];
1292 }
1293 break;
1294 default:
1295 TERMM(1, "Unknown data id (" + std::to_string(dataId) + ").");
1296 break;
1297 }
1298 } else {
1299 TERMM(1, "Invalid dataId: " + std::to_string(dataId));
1300 }
1301 }
1302
1314 template <typename DataType>
1315 void getCellDataDlb_(const MInt dataId, const MInt oldNoCells, const MInt* const bufferIdToCellId,
1316 DataType* const data) {
1317 TRACE();
1318
1319#ifdef LBCOLLECTOR_SOA_MEMORY_LAYOUT
1320 TERMM(1, "Balance data accessor not defined for SOA cell collector");
1321#endif
1322
1323 ASSERT((dataId > -1 && dataId < noSolverCellData()), "Invalid dataId in getCellData!");
1324
1325 MInt localBufferId = 0;
1326 for(MInt i = 0; i < oldNoCells; i++) {
1327 const MInt gridCellId = bufferIdToCellId[i];
1328 if(gridCellId < 0) {
1329 continue;
1330 }
1331
1332 const MInt cellId = grid().tree().grid2solver(gridCellId);
1333 if(cellId < 0 || cellId >= noInternalCells()) {
1334 continue;
1335 }
1336
1337 // LB solver cell data
1338 switch(dataId) {
1339 case CellData::VARIABLES: {
1340 const MInt dataSize = noVariables();
1341 std::copy_n(&m_cells.variables(cellId, 0), dataSize, &data[localBufferId * dataSize]);
1342 break;
1343 }
1344 case CellData::OLD_VARIABLES: {
1345 const MInt dataSize = noVariables();
1346 std::copy_n(&m_cells.oldVariables(cellId, 0), dataSize, &data[localBufferId * dataSize]);
1347 break;
1348 }
1349 case CellData::DISTRIBUTIONS: {
1350 const MInt dataSize = noDistributions();
1351 std::copy_n(&m_cells.distributions(cellId, 0), dataSize, &data[localBufferId * dataSize]);
1352 break;
1353 }
1354 case CellData::OLD_DISTRIBUTIONS: {
1355 const MInt dataSize = noDistributions();
1356 std::copy_n(&m_cells.oldDistributions(cellId, 0), dataSize, &data[localBufferId * dataSize]);
1357 break;
1358 }
1359 case CellData::NU: {
1360 const MInt dataSize = 1;
1361 std::copy_n(&m_cells.nu(cellId), dataSize, &data[localBufferId * dataSize]);
1362 break;
1363 }
1364 case CellData::IS_ACTIVE: {
1365 const MInt dataSize = 1;
1366 data[localBufferId * dataSize] = a_isActive(cellId);
1367 break;
1368 }
1369 default:
1370 TERMM(1, "Unknown data id (" + std::to_string(dataId) + ").");
1371 break;
1372 }
1373 localBufferId++;
1374 }
1375 }
1376
1377 // Helper struct to differentiate cell data which needs to be exchanged during DLB
1378 // TODO labels:LB Does not support thermal, change if proper sysEqn is introduced.
1379 struct CellData {
1380 static constexpr const MInt count = 6;
1381
1382 static constexpr const MInt VARIABLES = 0;
1383 static constexpr const MInt OLD_VARIABLES = 1;
1384 static constexpr const MInt DISTRIBUTIONS = 2;
1385 static constexpr const MInt OLD_DISTRIBUTIONS = 3;
1386 static constexpr const MInt NU = 4;
1387 static constexpr const MInt IS_ACTIVE = 5;
1388 };
1389
1390 // Data types of cell data
1391 const std::array<MInt, CellData::count> s_cellDataTypeDlb = {{MFLOAT, MFLOAT, MFLOAT, MFLOAT, MFLOAT, MINT}};
1392
1393 // Number of different data fields
1394 MInt noSolverCellData() const { return CellData::count; };
1395
1396 protected:
1397 // Timers
1398 struct {
1399 MInt solver;
1400 MInt initSolver;
1401 MInt solutionStep;
1402 MInt collision;
1403 MInt collisionBC;
1404 MInt propagation;
1405 MInt propagationBC;
1406 MInt srcTerms;
1407 MInt exchange;
1408 MInt residual;
1409 MInt packing;
1410 MInt unpacking;
1411 MInt communication;
1412
1413 MInt findG0Cells;
1414 MInt resetListsMb;
1415 MInt findG0Candidates;
1416 MInt geomNodal;
1417 MInt geomExchange;
1418 MInt calcNodalValues;
1419
1420 MInt prepComm;
1421
1422 MInt fft;
1423 } m_t;
1424
1425 // Status flags for adaptation
1426 public:
1427 MBool m_adaptationSinceLastRestart;
1428 MBool m_adaptationSinceLastSolution;
1429 MBool m_adaptation;
1430
1431#ifdef WAR_NVHPC_PSTL
1432 std::array<MFloat, 50> m_faculty;
1433 std::array<MInt, 27> m_nFld;
1434 std::array<MInt, 27> m_pFld;
1435 MInt m_idFld[POWX(3, nDim)][nDim];
1436 std::array<MInt, 24> m_mFld1;
1437 std::array<MInt, 24> m_mFld2;
1438 std::array<MInt, 27> m_oppositeDist;
1439 std::array<MFloat, 81> m_ppdfDir;
1440 std::array<MFloat, 4> m_tp;
1441 std::array<MInt, 3> m_distFld;
1442 std::array<MInt, 27> m_distType;
1443#endif
1444
1445 // Members for adaptation of grid file name
1446 MString m_reinitFileName; // Current grid file name.
1447 MString m_reinitFilePath; // Path to current grid file.
1448 protected:
1449 // Properties for adaptation.
1450 MString m_adaptationInitMethod;
1451 std::set<MInt> m_refinedParents;
1452 MBool m_solidLayerExtension = false;
1453 MBool m_writeLsData = false;
1454};
1455
1475template <MInt nDim>
1476template <MBool tI> // interpolation in time, added by Daniel Lauwers
1477MFloat LbSolver<nDim>::calculateDerivative(const MInt cellId, const MInt comp, const MInt spaceDir, const MFloat dt) {
1478 TRACE();
1479 const MInt lbdir1 = 2 * spaceDir;
1480 const MInt lbdir2 = lbdir1 + 1;
1481
1482 const MInt left = c_neighborId(cellId, lbdir1);
1483 const MInt right = c_neighborId(cellId, lbdir2);
1484 const MInt leftleft = (left > -1) ? c_neighborId(left, lbdir1) : -1;
1485 const MInt rightright = (right > -1) ? c_neighborId(right, lbdir2) : -1;
1486
1487 MFloat gradient = F0;
1488
1489 if((left > -1) && (right > -1)) {
1490 if(leftleft > -1 && rightright > -1) { // use central differences 4th order
1491 gradient = (-a_interpolatedVariable<tI>(rightright, comp, dt) + 8.0 * a_interpolatedVariable<tI>(right, comp, dt)
1492 - 8.0 * a_interpolatedVariable<tI>(left, comp, dt) + a_interpolatedVariable<tI>(leftleft, comp, dt))
1493 / 12.0;
1494 } else {
1495 if(leftleft > -1) { // backward differences 3nd order
1496 gradient =
1497 (2.0 * a_interpolatedVariable<tI>(right, comp, dt) + 3.0 * a_interpolatedVariable<tI>(cellId, comp, dt)
1498 - 6.0 * a_interpolatedVariable<tI>(left, comp, dt) + a_interpolatedVariable<tI>(leftleft, comp, dt))
1499 / 6.0;
1500 } else if(rightright > -1) { // forward differences 3nd order
1501 gradient =
1502 (-a_interpolatedVariable<tI>(rightright, comp, dt) + 6.0 * a_interpolatedVariable<tI>(right, comp, dt)
1503 - 3.0 * a_interpolatedVariable<tI>(cellId, comp, dt) - 2.0 * a_interpolatedVariable<tI>(left, comp, dt))
1504 / 6.0;
1505 } else { // use central differences 2nd order
1506 gradient = (a_interpolatedVariable<tI>(right, comp, dt) - a_interpolatedVariable<tI>(left, comp, dt)) / 2.0;
1507 }
1508 }
1509 } else { // use forward or backward differences 1st order
1510 // backward
1511 if(left > -1) {
1512 gradient = a_interpolatedVariable<tI>(cellId, comp, dt) - a_interpolatedVariable<tI>(left, comp, dt);
1513 }
1514 // forward
1515 else if(right > -1) {
1516 gradient = a_interpolatedVariable<tI>(right, comp, dt) - a_interpolatedVariable<tI>(cellId, comp, dt);
1517 }
1518 }
1519 return gradient;
1520}
1521
1528template <MInt nDim>
1529void LbSolver<nDim>::calculateVelocityDerivative(const MInt cellId, MFloat (&gradient)[nDim][nDim]) {
1530 TRACE();
1531 const MFloat F1bdx = FFPOW2(maxLevel() - a_level(cellId)); // LB units
1532 for(MInt dir = 0; dir < nDim; dir++) {
1533 const MInt dir0 = 2 * dir;
1534 const MInt dir1 = 2 * dir + 1;
1535
1536 auto isValid = [&](const MInt l_cellId, const MInt l_dir) -> MInt {
1537 if(l_cellId != -1) {
1538 const MInt nghbrId = c_neighborId(l_cellId, l_dir);
1539 if(nghbrId > -1 && a_isActive(nghbrId)) {
1540 return nghbrId;
1541 }
1542 }
1543 return -1;
1544 };
1545 // 1) set coefficients
1546 constexpr MInt noCellsInStencil = 5;
1547 std::array<MFloat, noCellsInStencil> coef;
1548 std::array<MInt, noCellsInStencil> cellIds;
1549 // 1.1) get cellIds and coefficients
1550 cellIds[1] = isValid(cellId, dir0); // left
1551 cellIds[0] = isValid(cellIds[1], dir0); // leftleft
1552 cellIds[2] = cellId; // middle
1553 cellIds[3] = isValid(cellId, dir1); // right
1554 cellIds[4] = isValid(cellIds[3], dir1); // rightright
1555 getDerivativeStencilAndCoefficient(cellIds, coef);
1556 // 1.2) scale coefficients with correct dx in case of refined
1557 for(MInt i = 0; i < noCellsInStencil; i++) {
1558 coef[i] *= F1bdx;
1559 }
1560 // 2) calculate gradient
1561 for(MInt veloId = 0; veloId < nDim; veloId++) {
1562 gradient[veloId][dir] = 0.0;
1563 for(MInt i = 0; i < noCellsInStencil; i++) {
1564 gradient[veloId][dir] += cellIds[i] < 0 ? 0.0 : coef[i] * a_variable(cellIds[i], veloId);
1565 }
1566 }
1567 }
1568}
1569
1570#endif
GridCell
GridCell
Grid cell Property Labels.
Definition: cartesiangridcellproperties.h:14
Collector::size
MInt size()
Definition: collector.h:29
LbLpt
Definition: lblpt.h:30
MConservativeVariables
Storage of the Position of the Conservative Variables (RHO, RHO_VV, RHO_E) in the value vectors of th...
Definition: variables.h:31
MPrimitiveVariables
Storage of the Position of the Primitive Variables (u, v, w, T, p) in the value vectors of the solver...
Definition: variables.h:110
PostProcessing
Definition: postprocessing.h:85
PostProcessingLb
Definition: postprocessinglb.h:21
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
Solver::noInternalCells
virtual MInt noInternalCells() const =0
Return the number of internal cells within this solver.
maia::container::Container::size
constexpr MInt size() const
Return size (i.e., currently used number of nodes)
Definition: container.h:89
maia::lb::collector::LbCellCollector
Class that represents LB cell collector.
Definition: lbcellcollector.h:114
maia::lb::collector::LbCellCollector::variables
MFloat & variables(const MInt id, const MInt eid)
Definition: lbcellcollector.h:426
maia::lb::collector::LbCellCollector::nu
MFloat & nu(const MInt id)
Accessor for nu.
Definition: lbcellcollector.h:334
maia::lb::collector::LbCellCollector::oldVariables
MFloat & oldVariables(const MInt id, const MInt eid)
Definition: lbcellcollector.h:453
maia::lb::collector::LbCellCollector::distributions
MFloat & distributions(const MInt id, const MInt eid)
Definition: lbcellcollector.h:480
maia::lb::collector::LbCellCollector::oldDistributions
MFloat & oldDistributions(const MInt id, const MInt eid)
Definition: lbcellcollector.h:502
LbInterpolationType
LbInterpolationType
Definition: enums.h:265
MINT
@ MINT
Definition: enums.h:269
MFLOAT
@ MFLOAT
Definition: enums.h:269
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
POWX
T constexpr POWX(T base, MUint exponent)
Compile time power calculation.
Definition: functions.h:146
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
LbCell
LbCell
LB cell Property Labels.
Definition: lbcellproperties.h:14
LbUpdateMacroscopic
LbUpdateMacroscopic
Definition: lbsolver.h:34
POSTPROPAGATION
@ POSTPROPAGATION
Definition: lbsolver.h:34
INCOLLISION
@ INCOLLISION
Definition: lbsolver.h:34
PRECOLLISION
@ PRECOLLISION
Definition: lbsolver.h:34
FFPOW2
constexpr MFloat FFPOW2(MInt x)
Definition: maiaconstants.h:162
MInt
int32_t MInt
Definition: maiatypes.h:62
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
MChar
char MChar
Definition: maiatypes.h:56
maia::collector_memory::cellId
void const MInt cellId
Definition: collector.h:239
maia::grid::hilbert::index
IdType index(const FloatType *const x, const IdType level)
Return Hilbert index for given location and level in 2D or 3D.
Definition: hilbert.h:165
maia::grid::tree::Cell
GridCell Cell
Underlying enum type for property access.
Definition: cartesiangridtree.h:101
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54
maia::coupling::Mapping
Multi-to-multi mapping class.
Definition: surfacecoupling.h:186