MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
cartesiansolver.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 CARTESIANSOLVER_H_
8#define CARTESIANSOLVER_H_
9
10#include <algorithm>
11
12#include "COMM/mpioverride.h"
13#include "FV/fvcartesiancellcollector.h"
14#include "GEOM/geometry.h"
15#include "GEOM/geometryanalytic.h"
16#include "GEOM/geometryelement.h"
17#include "GRID/cartesiangrid.h"
18#include "GRID/cartesiangridpoint.h"
19#include "GRID/cartesiangridproxy.h"
20#include "GRID/cartesianpointbasedcell.h"
21#include "INCLUDE/maiaconstants.h"
22#include "INCLUDE/maiamacro.h"
23#include "INCLUDE/maiatypes.h"
24#include "IO/context.h"
25#include "IO/parallelio.h"
26#include "MEMORY/scratch.h"
27#include "UTIL/debug.h"
28#include "property.h"
29#include "solver.h"
30
31// Forward declarations
32template <MInt nDim>
33class CartesianGrid;
34
35namespace maia {
36
37struct PatchRefinement {
38 // refinement patch
39 MInt noLocalPatchRfnLvls() { return (MInt)localRfnLevelMethods.size(); };
40 MInt noPatchesPerLevel(MInt addLevel) { return (MInt)localRfnLevelMethods[addLevel].size(); };
41
42 std::vector<MString> localRfnLevelMethods;
43
44 MInt* localRfnLevelPropertiesOffset = nullptr;
45 MInt* noLocalRfnPatchProperties = nullptr;
46 MFloat** localRfnPatchProperties = nullptr;
47};
48
49template <MInt nDim, class SolverType>
50class CartesianSolver : public Solver {
51 public:
52 // Types
53 using Grid = CartesianGrid<nDim>;
54 using GridProxy = typename maia::grid::Proxy<nDim>;
55 using Geom = Geometry<nDim>;
56 using TreeProxy = maia::grid::tree::TreeProxy<nDim>;
57 // using SolverCell = FvCell;
58 using Cell = maia::grid::tree::Cell;
59
60 // Methods
61 // Constructor takes solverId, reference to grid
62 CartesianSolver(const MInt solverId, GridProxy& gridProxy_, const MPI_Comm comm, const MBool checkActive = false);
63
65 MInt minLevel() const { return grid().minLevel(); }
66 MInt maxLevel() const { return grid().maxLevel(); }
67 MInt maxNoGridCells() const { return grid().maxNoCells(); }
68 MInt maxRefinementLevel() const { return grid().maxRefinementLevel(); }
69 MInt maxUniformRefinementLevel() const { return grid().maxUniformRefinementLevel(); }
70 MInt noNeighborDomains() const { return grid().noNeighborDomains(); }
71 MInt neighborDomain(const MInt id) const { return grid().neighborDomain(id); }
72 MLong domainOffset(const MInt id) const { return grid().domainOffset(id); }
73 MInt noHaloLayers() const { return grid().noHaloLayers(); }
74 MInt noHaloCells(const MInt domainId) const { return grid().noHaloCells(domainId); }
75 MInt haloCellId(const MInt domainId, const MInt cellId) const { return grid().haloCell(domainId, cellId); }
76 MInt noWindowCells(const MInt domainId) const { return grid().noWindowCells(domainId); }
77 MInt windowCellId(const MInt domainId, const MInt cellId) const { return grid().windowCell(domainId, cellId); }
78 MString gridInputFileName() const { return grid().gridInputFileName(); }
79 MFloat reductionFactor() const { return grid().reductionFactor(); }
80 MFloat centerOfGravity(const MInt dir) const { return grid().centerOfGravity(dir); }
81 MInt neighborList(const MInt cellId, const MInt dir) const { return grid().neighborList(cellId, dir); }
82 const MLong& localPartitionCellGlobalIds(const MInt cellId) const {
83 return grid().localPartitionCellGlobalIds(cellId);
84 }
85 MLong localPartitionCellOffsets(const MInt index) const { return grid().localPartitionCellOffsets(index); }
86 MInt noMinCells() const { return grid().noMinCells(); }
87 MInt minCell(const MInt id) const { return grid().minCell(id); }
88 const MInt& haloCell(const MInt domainId, const MInt cellId) const { return grid().haloCell(domainId, cellId); }
89 const MInt& windowCell(const MInt domainId, const MInt cellId) const { return grid().windowCell(domainId, cellId); }
90
91 MBool isActive() const override { return grid().isActive(); }
92
93 // Grid and tree proxy accessors
94 constexpr GridProxy& grid() const { return m_gridProxy; }
95 GridProxy& grid() { return m_gridProxy; }
96
97 virtual void sensorDerivative(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
98 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
99 TERMM(1, "Not implemented for this solver");
100 };
101 virtual void sensorDivergence(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
102 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
103 TERMM(1, "Not implemented for this solver");
104 };
105 virtual void sensorTotalPressure(std::vector<std::vector<MFloat>>& /*unused*/,
106 std::vector<std::bitset<64>>& /*unused*/, std::vector<MFloat>& /*unused*/,
107 MInt /*unused*/, MInt /*unused*/) {
108 TERMM(1, "Not implemented for this solver");
109 };
110 virtual void sensorEntropyGrad(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
111 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
112 TERMM(1, "Not implemented for this solver");
113 };
114 virtual void sensorEntropyQuot(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
115 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
116 TERMM(1, "Not implemented for this solver");
117 };
118 virtual void sensorVorticity(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
119 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
120 TERMM(1, "Not implemented for this solver");
121 };
122 virtual void sensorInterface(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
123 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
124 TERMM(1, "Not implemented for this solver");
125 };
126 void sensorLimit(std::vector<std::vector<MFloat>>& /*sensors*/, std::vector<std::bitset<64>>& /*sensorCellFlag*/,
127 std::vector<MFloat>& /*sensorWeight*/, MInt /*sensorOffset*/, MInt /*sen*/,
128 std::function<MFloat(MInt)> /*value*/, const MFloat /*limit*/, const MInt* /*range*/,
129 const MBool /*refineDiagonals*/, const MBool allowCoarsening = true);
130 void sensorSmooth(std::vector<std::vector<MFloat>>& /*sensors*/, std::vector<std::bitset<64>>& /*sensorCellFlag*/,
131 std::vector<MFloat>& /*sensorWeight*/, MInt /*sensorOffset*/, MInt /*sen*/);
132 void sensorBand(std::vector<std::vector<MFloat>>& /*sensors*/, std::vector<std::bitset<64>>& /*sensorCellFlag*/,
133 std::vector<MFloat>& /*sensorWeight*/, MInt /*sensorOffset*/, MInt /*sen*/);
134 virtual void sensorMeanStress(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
135 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
136 TERMM(1, "Not implemented for this solver");
137 };
138 virtual void sensorParticle(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
139 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
140 TERMM(1, "Not implemented for this solver");
141 };
142 virtual void sensorSpecies(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
143 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/) {
144 TERMM(1, "Not implemented for this solver");
145 };
146 virtual void sensorPatch(std::vector<std::vector<MFloat>>& sensor, std::vector<std::bitset<64>>& sensorCellFlag,
147 std::vector<MFloat>& sensorWeight, MInt sensorOffset, MInt sen) {
148 patchRefinement(sensor, sensorCellFlag, sensorWeight, sensorOffset, sen);
149 };
150 virtual void sensorCutOff(std::vector<std::vector<MFloat>>& /*unused*/, std::vector<std::bitset<64>>& /*unused*/,
151 std::vector<MFloat>& /*unused*/, MInt /*unused*/, MInt /*unused*/)
152 {
153 TERMM(1, "Not implemented for this solver");
154 };
155 void saveSensorData(const std::vector<std::vector<MFloat>>& sensors, const MInt& level, const MString& gridFileName,
156 const MInt* const recalcIds) override;
157
158 // Asserts that cellId belongs to a valid grid cell
159 void assertValidGridCellId(const MInt) const {}
161 MLong c_parentId(const MInt cellId) const { return solver().grid().tree().parent(cellId); }
163 MLong c_neighborId(const MInt cellId, const MInt dir) const { return solver().grid().tree().neighbor(cellId, dir); }
164 MInt c_noCells() const { return grid().tree().size(); }
165 MInt c_level(const MInt cellId) const { return grid().tree().level(cellId); }
166
167 // \brief Return global grid id
168 MLong c_globalGridId(const MInt cellId) { return grid().raw().a_globalId(grid().tree().solver2grid(cellId)); }
169
170 // General Exchange Functions
171 template <typename T>
172 void exchangeData(T* data, const MInt dataBlockSize = 1);
173 template <typename T>
174 void exchangeLeafData(std::function<T&(MInt, MInt)> data, const MInt noDat = 1);
175
176
177 template <class G, class S, class M>
178 void exchangeSparseLeafValues(G getData, S setData, const MInt dataSize, M cellMapping);
179
180 template <typename T>
181 void exchangeAzimuthalPer(T* data, MInt dataBlockSize = 1, MInt firstBlock = 0);
182
183 template <typename T>
184 void collectVariables(T* variablesIn, ScratchSpace<T>& variablesOut, const std::vector<MString>& variablesNameIn,
185 std::vector<MString>& variablesNameOut, const MInt noVars, const MInt noCells,
186 const MBool reverseOrder = false);
187 template <typename T>
188 void collectVariables(T** variablesIn, ScratchSpace<T>& variablesOut, const std::vector<MString>& variablesNameIn,
189 std::vector<MString>& variablesNameOut, const MInt noVars, const MInt noCells);
190
191 void saveGridFlowVars(const MChar* fileName, const MChar* gridFileName, const MInt noTotalCells,
192 const MInt noInternal, MFloatScratchSpace& dbVariables, std::vector<MString>& dbVariablesName,
193 MInt noDbVars, MIntScratchSpace& idVariables, std::vector<MString>& idVariablesName,
194 MInt noIdVars, MFloatScratchSpace& dbParameters, std::vector<MString>& dbParametersName,
195 MIntScratchSpace& idParameters, std::vector<MString>& idParametersName, MInt* recalcIds,
196 MFloat time);
197
198 template <typename T>
199 void collectParameters(T, ScratchSpace<T>&, const MChar*, std::vector<MString>&);
200
201 void calcRecalcCellIdsSolver(const MInt* const recalcIdsTree, MInt& noCells, MInt& noInternalCellIds,
202 std::vector<MInt>& recalcCellIdsSolver, std::vector<MInt>& reorderedCellIds);
203
204 protected:
205 // Grid/geometry-related methods
206 // stl-boundary related:
207 void identifyBoundaryCells(MBool* const isInterface, const std::vector<MInt>& bndCndIds = std::vector<MInt>());
208 void identifyBoundaryCells();
209 MBool cellIsOnGeometry(MInt cellId, Geometry<nDim>* geom);
210 void setBoundaryDistance(const MBool* const interfaceCell, const MFloat* const outerBandWidth,
211 MFloatScratchSpace& distance);
212 void markSurrndCells(MIntScratchSpace& inList, const MInt bandWidth, const MInt level,
213 const MBool refineDiagonals = true);
214 void receiveWindowTriangles();
215
216 // adaptation related
217 void compactCells();
218 MInt createCellId(const MInt gridCellId);
219 void removeCellId(const MInt cellId);
220
221 MInt inCell(const MInt cellId, MFloat* point, MFloat fac = F1);
222
223 MInt setUpInterpolationStencil(const MInt cellId, MInt*, const MFloat*, std::function<MBool(MInt, MInt)>,
224 MBool allowIncompleteStencil);
225
226 template <MBool cartesianInterpolation>
227 MFloat interpolateFieldData(MInt*, MFloat*, MInt varId, std::function<MFloat(MInt, MInt)> scalarField,
228 std::function<MFloat(MInt, MInt)> coordinate);
229
230 MFloat leastSquaresInterpolation(MInt*, MFloat*, MInt varId, std::function<MFloat(MInt, MInt)> scalarField,
231 std::function<MFloat(MInt, MInt)> coordinate);
232
233 void checkNoHaloLayers();
234
235 void mapInterpolationCells(std::map<MInt, MInt>& cellMap);
236 void setHaloCellsOnInactiveRanks();
237
238 void patchRefinement(std::vector<std::vector<MFloat>>& sensors, std::vector<std::bitset<64>>& sensorCellFlag,
239 std::vector<MFloat>& sensorWeight, MInt sensorOffset, MInt sen);
240
241 void reOrderCellIds(std::vector<MInt>& reOrderedCells);
242 void recomputeGlobalIds(std::vector<MInt>&, std::vector<MLong>&);
243
244 void extractPointIdsFromGrid(Collector<PointBasedCell<nDim>>*&,
245 Collector<CartesianGridPoint<nDim>>*&,
246 const MBool,
247 const std::map<MInt, MInt>&,
248 MInt levelThreshold = 999999,
249 MFloat* bBox = nullptr,
250 MBool levelSetMb = false) const;
251
252 private:
253 // Methods
254 // CRTP accessors
255 SolverType& solver() { return static_cast<SolverType&>(*this); }
256 constexpr const SolverType& solver() const { return static_cast<const SolverType&>(*this); }
257
258 protected:
259 MInt* m_rfnBandWidth;
260 MInt m_noSensors;
261 MInt m_adaptationInterval;
262 MInt m_adaptationStep = F0;
263 std::vector<MInt> m_maxSensorRefinementLevel;
264 std::vector<MFloat> m_sensorWeight;
265 std::vector<MFloat> m_sensorDerivativeVariables;
266 MBool m_adaptation;
267 MBool m_adapts;
268 MInt m_noInitialSensors;
269 MBool m_resTriggeredAdapt = false;
270 MInt m_noSmoothingLayers = -1;
271 MInt m_sensorBandAdditionalLayers;
272
273 // adaptation
274 MBool m_sensorInterface;
275 MBool m_sensorParticle;
276 std::vector<MString> m_sensorType;
277 MInt* m_recalcIds = nullptr;
278
279 using fun = void (CartesianSolver<nDim, SolverType>::*)(std::vector<std::vector<MFloat>>&,
280 std::vector<std::bitset<64>>&, std::vector<MFloat>&, MInt,
281 MInt);
282 std::vector<fun> m_sensorFnPtr;
283
284 // necessary for generalised interpolation function!
285 MInt m_maxNoSets = -1;
286
287 // Azimuthal communication
288 std::vector<MFloat> m_azimuthalCartRecCoord;
289
290 // Holds the reverse-direction-Information
291 const MInt m_revDir[6] = {1, 0, 3, 2, 5, 4};
292 const MInt m_noDirs = 2 * nDim;
293
294 private:
295 void readPatchProperties();
296 void initAdaptation();
297 void addChildren(std::vector<MInt>& reOrderedIds, const MInt parentId);
298
299
300 // Data
301 MString m_gridCutTest{};
302 GridProxy& m_gridProxy;
303
304 PatchRefinement* m_patchRefinement = nullptr;
305 MBool m_testPatch = false;
306 MInt m_patchStartTimeStep = -1;
307 MInt m_patchStopTimeStep = -1;
308};
309
310template <MInt nDim, class SolverType>
311CartesianSolver<nDim, SolverType>::CartesianSolver(const MInt solverId_, GridProxy& gridProxy_, const MPI_Comm comm,
312 const MBool checkActive)
313 : Solver(solverId_, comm, checkActive ? gridProxy_.isActive() : true), m_gridProxy(gridProxy_) {
314 // Get necessary properties
315 m_gridCutTest = "SAT";
316 m_gridCutTest = Context::getSolverProperty<MString>("gridCutTest", solverId_, AT_, &m_gridCutTest);
317
318 initAdaptation();
319}
320
325template <MInt nDim, class SolverType>
326void CartesianSolver<nDim, SolverType>::identifyBoundaryCells(MBool* const isInterface,
327 const std::vector<MInt>& bndCndIds) {
328 // TODO labels:GEOM,DLB add mode to tag only halos during DLB if isInterface property is communicated?
329 TRACE();
330
331 MFloat target[6] = {0, 0, 0, 0, 0, 0};
332
333 // Check all cells for intersection with geometry
334 const MInt noCells = solver().grid().tree().size();
335
336 for(MInt i = 0; i < noCells; i++) {
337 const MFloat cellHalfLength = solver().grid().cellLengthAtLevel(solver().grid().tree().level(i) + 1);
338 for(MInt j = 0; j < nDim; j++) {
339 target[j] = solver().a_coordinate(i, j) - cellHalfLength;
340 target[j + nDim] = solver().a_coordinate(i, j) + cellHalfLength;
341 }
342
343 std::vector<MInt> nodeList;
344 // TODO labels:GEOM why not read gridCutTest once in the geometry and decide there which getIntersectionElements
345 // function to use?
346 if(m_gridCutTest == "SAT") {
347 solver().geometry().getIntersectionElements(target, nodeList, cellHalfLength, &solver().a_coordinate(i, 0));
348 } else {
349 solver().geometry().getIntersectionElements(target, nodeList);
350 }
351
352 isInterface[i] = false;
353
354 const MInt noNodes = nodeList.size();
355 if(noNodes > 0 && bndCndIds.empty()) {
356 isInterface[i] = true;
357 } else if(noNodes > 0) {
358 for(MInt n = 0; n < noNodes; n++) {
359 for(const auto& bnd : bndCndIds) {
360 if(bnd == solver().geometry().elements[nodeList[n]].m_bndCndId) {
361 isInterface[i] = true;
362 }
363 }
364 }
365 }
366 }
367}
368
369
374template <MInt nDim, class SolverType>
375void CartesianSolver<nDim, SolverType>::identifyBoundaryCells() {
376 TRACE();
377
378 const MInt noCells = solver().grid().tree().size();
379 MBoolScratchSpace isInterface(noCells, AT_, "isInterface");
380 identifyBoundaryCells(&isInterface[0]);
381 for(MInt i = 0; i < noCells; i++) {
382 solver().a_isInterface(i) = isInterface[i];
383 }
384}
385
386
398template <MInt nDim, class SolverType>
399MBool CartesianSolver<nDim, SolverType>::cellIsOnGeometry(MInt cellId, Geometry<nDim>* geom) {
400 TRACE();
401
402 MFloat target[6] = {0, 0, 0, 0, 0, 0};
403 const MFloat cellHalfLength = solver().grid().cellLengthAtLevel(solver().a_level(cellId) + 1);
404 std::vector<MInt> nodeList;
405
406 for(MInt dim = 0; dim < nDim; dim++) {
407 target[dim] = solver().a_coordinate(cellId, dim) - cellHalfLength;
408 target[dim + nDim] = solver().a_coordinate(cellId, dim) + cellHalfLength;
409 }
410
411 if(m_gridCutTest == "SAT") {
412 geom->getIntersectionElements(target, nodeList, cellHalfLength, &solver().a_coordinate(cellId, 0));
413 } else {
414 geom->getIntersectionElements(target, nodeList);
415 }
416
417 return nodeList.size() > 0;
418}
419
420
446template <MInt nDim, class SolverType>
447void CartesianSolver<nDim, SolverType>::receiveWindowTriangles() {
448 TRACE();
449
450 using namespace std;
451
452 m_log << " - exchanging triangles " << std::endl;
453 vector<set<pair<MInt, MInt>>> triangleIdsPerDomain;
454
455 MPI_Request* mpi_request = nullptr;
456 MPI_Status status;
457 mAlloc(mpi_request, solver().grid().noNeighborDomains(), "mpi_request", AT_);
458
459 m_log << " * collecting window triangles" << std::endl;
460 // 1. collect unique triangles per neighbor domain
461 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
462 set<pair<MInt, MInt>> triangles;
463 for(MInt w = 0; w < (signed)solver().grid().noWindowCells(n); w++) {
464 MInt currentId = solver().grid().windowCell(n, w);
465
466 // on coarsest level
467 if(solver().grid().tree().parent(currentId) < 0) {
468 // a. Create target for check
469 MFloat target[6];
470 MFloat cellHalfLength = solver().grid().cellLengthAtLevel(solver().grid().tree().level(currentId) + 1);
471 cellHalfLength += 0.005 * cellHalfLength;
472
473 std::vector<MInt> nodeList;
474
475 for(MInt j = 0; j < nDim; j++) {
476 target[j] = solver().a_coordinate(currentId, j) - cellHalfLength;
477 target[j + nDim] = solver().a_coordinate(currentId, j) + cellHalfLength;
478 }
479
480 // b. Do the intersection test
481 solver().geometry().getIntersectionElements(target, nodeList, cellHalfLength,
482 &solver().a_coordinate(currentId, 0));
483
484 for(MInt t = 0; t < (signed)nodeList.size(); t++) {
485 triangles.insert(pair<MInt, MInt>(nodeList[t], solver().geometry().elements[nodeList[t]].m_originalId));
486 }
487 }
488 }
489 triangleIdsPerDomain.push_back(triangles);
490 }
491
492 // 2. let me know about the uniqueIds of my neighbors
493 MIntScratchSpace numUniques(solver().grid().noNeighborDomains(), AT_, "numUniques");
494 MInt myNumUniques = solver().geometry().m_uniqueOriginalTriId.size();
495 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
496 MPI_Issend(&myNumUniques, 1, MPI_INT, solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &mpi_request[n], AT_,
497 "myNumUniques");
498 }
499 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
500 MPI_Recv(&numUniques[n], 1, MPI_INT, solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &status, AT_,
501 "numUniques[n]");
502 }
503 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
504 MPI_Wait(&mpi_request[n], &status, AT_);
505 }
506
507 MInt sumuniques = 0;
508 MIntScratchSpace uniquesDomOff(solver().grid().noNeighborDomains(), AT_, "uniquesDomOff");
509 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
510 uniquesDomOff[n] = sumuniques;
511 sumuniques += numUniques[n];
512 }
513
514 // send and receive the uniques
515 MIntScratchSpace myUniques(myNumUniques, AT_, "myUniques");
516 MInt ua = 0;
517 for(auto u = solver().geometry().m_uniqueOriginalTriId.begin(); u != solver().geometry().m_uniqueOriginalTriId.end();
518 ++u, ua++) {
519 myUniques[ua] = *u;
520 }
521 MIntScratchSpace alluniques(sumuniques, AT_, "alluniques");
522
523 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
524 if(!myUniques.empty())
525 MPI_Issend(myUniques.getPointer(), myUniques.size(), MPI_INT, solver().grid().neighborDomain(n), 0,
526 MPI_COMM_WORLD, &mpi_request[n], AT_, "myUniques.getPointer()");
527 }
528 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
529 if(numUniques[n] > 0) {
530 MPI_Recv(&alluniques[uniquesDomOff[n]], numUniques[n], MPI_INT, solver().grid().neighborDomain(n), 0,
531 MPI_COMM_WORLD, &status, AT_, "alluniques[uniquesDomOff[n]]");
532 }
533 }
534 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
535 if(numUniques[n] > 0) {
536 MPI_Wait(&mpi_request[n], &status, AT_);
537 }
538 }
539
540 m_log << " * sending and receiving unique originalIds" << std::endl;
541 m_log << " + sending " << myUniques.size() << " originalIds to: ";
542 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
543 m_log << solver().grid().neighborDomain(n) << " ";
544 }
545 m_log << std::endl;
546 m_log << " * receiving unique originalIds" << std::endl;
547 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
548 m_log << " + receiving from " << solver().grid().neighborDomain(n) << ": " << numUniques[n] << std::endl;
549 }
550
551 m_log << " * removing doubles from sender list" << std::endl;
552
553 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
554 // for searching
555 set<MInt> un;
556 for(MInt off = uniquesDomOff[n]; off < uniquesDomOff[n] + numUniques[n]; off++) {
557 un.insert(alluniques[off]);
558 }
559
560 for(auto it = triangleIdsPerDomain[n].begin(); it != triangleIdsPerDomain[n].end();) {
561 auto iun = un.find(get<1>(*it));
562 if(iun != un.end()) {
563 triangleIdsPerDomain[n].erase(it);
564 it = triangleIdsPerDomain[n].begin();
565 } else {
566 ++it;
567 }
568 }
569 }
570
571
572 MIntScratchSpace toReceive(solver().grid().noNeighborDomains(), AT_, "toReceive");
573
574
575 m_log << " * sending numbers of triangles to send to other domains" << std::endl;
576 // 3. send information how many triangles will be send
577 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
578 MInt numtris = triangleIdsPerDomain[n].size();
579 MPI_Issend(&numtris, 1, MPI_INT, solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &mpi_request[n], AT_,
580 "numtris");
581 }
582
583 // receive information
584 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
585 MPI_Recv(&toReceive[n], 1, MPI_INT, solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &status, AT_,
586 "toReceive[n]");
587 }
588 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
589 MPI_Wait(&mpi_request[n], &status, AT_);
590 }
591
592 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
593 if(toReceive[n] > 0) {
594 m_log << " + receive from domain " << solver().grid().neighborDomain(n) << " : " << toReceive[n]
595 << std::endl;
596 }
597 if(!triangleIdsPerDomain[n].empty()) {
598 m_log << " + send to domain " << solver().grid().neighborDomain(n) << " : "
599 << triangleIdsPerDomain[n].size() << std::endl;
600 }
601 }
602
603
604 MInt sumofreceive = 0;
605 MInt sumofsend = 0;
606 MIntScratchSpace offsetreceive(solver().grid().noNeighborDomains(), AT_, "offsetreceive");
607 MIntScratchSpace offsetsend(solver().grid().noNeighborDomains(), AT_, "offsetsend");
608 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
609 offsetsend[n] = sumofsend;
610 sumofsend += triangleIdsPerDomain[n].size();
611 offsetreceive[n] = sumofreceive;
612 sumofreceive += toReceive[n];
613 }
614
615 m_log << " * receive offsets:" << std::endl;
616 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
617 if(toReceive[n] > 0) {
618 m_log << " + from domain " << solver().grid().neighborDomain(n) << ": " << offsetreceive[n] << std::endl;
619 }
620 }
621
622 m_log << " * send offsets:" << endl;
623 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
624 if(!triangleIdsPerDomain[n].empty()) {
625 m_log << " + from domain " << solver().grid().neighborDomain(n) << ": " << offsetsend[n] << std::endl;
626 }
627 }
628
629
630 // create array holding triangles to receive
631 MInt trisize = nDim // normal
632 + (nDim * nDim) // vertcies
633 + 3; // segmentId, bndCndId, orininalId
634
635 MFloatScratchSpace recTris(sumofreceive * trisize, AT_, "recTris");
636 MFloatScratchSpace sndTris(sumofsend * trisize, AT_, "sndTris");
637
638 m_log << " * sending triangles" << std::endl;
639 // send the triangles
640 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
641 if(triangleIdsPerDomain[n].empty()) {
642 continue;
643 }
644 // create buffer
645
646 MInt j = offsetsend[n] * trisize;
647 for(const auto& tri : triangleIdsPerDomain[n]) {
648 sndTris[j++] = (MFloat)solver().geometry().elements[get<0>(tri)].m_originalId;
649 sndTris[j++] = (MFloat)solver().geometry().elements[get<0>(tri)].m_segmentId;
650 sndTris[j++] = (MFloat)solver().geometry().elements[get<0>(tri)].m_bndCndId;
651
652 // normal
653 for(MInt d = 0; d < nDim; d++, j++) {
654 sndTris[j] = solver().geometry().elements[get<0>(tri)].m_normal[d];
655 }
656
657 // vertcies
658 for(MInt d1 = 0; d1 < nDim; d1++) {
659 for(MInt d2 = 0; d2 < nDim; d2++, j++) {
660 sndTris[j] = solver().geometry().elements[get<0>(tri)].m_vertices[d1][d2];
661 }
662 }
663 }
664 MPI_Issend(sndTris.getPointer() + (offsetsend[n] * trisize), trisize * triangleIdsPerDomain[n].size(), MPI_DOUBLE,
665 solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &mpi_request[n], AT_,
666 "sndTris.getPointer()+(offsetsend[n]*trisize)");
667 }
668
669 m_log << " * receiving triangles" << endl;
670 // receive triangles
671 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
672 if(toReceive[n] > 0) {
673 MPI_Recv(recTris.getPointer() + (offsetreceive[n] * trisize), toReceive[n] * trisize, MPI_DOUBLE,
674 solver().grid().neighborDomain(n), 0, MPI_COMM_WORLD, &status, AT_,
675 "recTris.getPointer()+(offsetreceive[n]*trisize)");
676 }
677 }
678 for(MInt n = 0; n < solver().grid().noNeighborDomains(); n++) {
679 if(toReceive[n] > 0) {
680 MPI_Wait(&mpi_request[n], &status, AT_);
681 }
682 }
683
684 m_log << " * inserting received triangles" << endl;
685
686 // if(solver().geometry().GetNoElements() == 0)
687 if(sumofreceive > 0) {
688 solver().geometry().resizeCollector(sumofreceive);
689 }
690
691
692 // 4. insert the triangles
693 for(MInt tri = 0; tri < sumofreceive * trisize; tri += trisize) {
694 auto originalId = (MInt)recTris[tri];
695 // check if triangle already in my domain
696 // not found, great, insert the triangle
697 if(solver().geometry().m_uniqueOriginalTriId.find(originalId) == solver().geometry().m_uniqueOriginalTriId.end()) {
698 solver().geometry().addElement(&recTris[tri]);
699 }
700 }
701
702 // 5. update tree
703 solver().geometry().calculateBoundingBox();
704
705 if(solver().geometry().m_debugParGeom && solver().geometry().GetNoElements() > 0) {
706 solver().geometry().writeParallelGeometryVTK("allpluswindow");
707 }
708
709 if(solver().geometry().GetNoElements() > 0 && sumofreceive > 0) {
710 solver().geometry().rebuildAdtTree();
711 }
712}
713
714
715// --------------------------------------------------------------------------------------
716
717
722template <MInt nDim, class SolverType>
723void CartesianSolver<nDim, SolverType>::compactCells() {
724 MIntScratchSpace oldCellId(m_gridProxy.maxNoCells(), AT_, "oldCellId");
725 MIntScratchSpace isToDelete(m_gridProxy.maxNoCells(), AT_, "isToDelete");
726 oldCellId.fill(-1);
727 isToDelete.fill(0);
728 for(auto& i : solver().m_freeIndices) {
729 isToDelete[i] = 1;
730 }
731 solver().m_freeIndices.clear();
732
733 if(grid().azimuthalPeriodicity()) {
734 m_gridProxy.correctAzimuthalHaloCells();
735 }
736
737 m_gridProxy.resizeGridMap(solver().m_cells.size());
738
739 // 0. some sanity checks
740 ASSERT(solver().m_cells.size() == m_gridProxy.tree().size(),
741 "m_cells: " + std::to_string(solver().m_cells.size()) + " tree: " + std::to_string(m_gridProxy.tree().size()));
742 for(MInt cellId = 0; cellId < solver().m_cells.size(); cellId++) {
743 ASSERT(isToDelete(cellId)
744 || solver().a_isHalo(cellId) == solver().grid().raw().a_isHalo(grid().tree().solver2grid(cellId)),
745 "");
746 }
747 for(MInt gridCellId = 0; gridCellId < solver().grid().raw().treeb().size(); gridCellId++) {
748 if(solver().grid().raw().treeb().solver(gridCellId, solver().solverId())) {
749 ASSERT(
750 grid().tree().grid2solver(gridCellId) > -1 && grid().tree().grid2solver(gridCellId) < solver().m_cells.size(),
751 std::to_string(gridCellId) + " " + std::to_string(grid().tree().grid2solver(gridCellId)) + " "
752 + std::to_string(solver().m_cells.size()) + " " + std::to_string(solver().grid().raw().treeb().size()));
753 } else {
754 ASSERT(grid().tree().grid2solver(gridCellId) < 0,
755 std::to_string(gridCellId) + " " + std::to_string(grid().tree().grid2solver(gridCellId)) + " "
756 + std::to_string(solver().m_cells.size()) + " "
757 + std::to_string(solver().grid().raw().treeb().size()));
758 }
759 }
760
761 // 1. determine number of cells and internal cells
762 MInt noCells = 0;
763 MInt noInternalCells = 0;
764 oldCellId.fill(-1);
765 for(MInt cellId = 0; cellId < solver().m_cells.size(); cellId++) {
766 if(isToDelete(cellId) != 0) {
767 continue;
768 }
769 oldCellId(cellId) = cellId;
770 if(!solver().a_isHalo(cellId)) {
771 noInternalCells++;
772 }
773 noCells++;
774 }
775
776 if(solver().grid().raw().treeb().noSolvers() == 1) {
777 ASSERT(noCells == solver().grid().raw().treeb().size(), "");
778 ASSERT(noInternalCells == solver().grid().raw().m_noInternalCells, "");
779 }
780
781 // 2. remove holes created by previously deleted cells and move halo cells to the back
782 MInt otherId = solver().m_cells.size() - 1;
783 for(MInt cellId = 0; cellId < noInternalCells; cellId++) {
784 if(isToDelete(cellId) || solver().a_isHalo(cellId)) {
785 while(isToDelete(otherId) || solver().a_isHalo(otherId)) {
786 otherId--;
787 }
788 ASSERT(cellId < otherId, "");
789
790 solver().swapCells(cellId, otherId);
791 grid().swapSolverIds(cellId, otherId);
792 std::swap(oldCellId(cellId), oldCellId(otherId));
793 std::swap(isToDelete(cellId), isToDelete(otherId));
794
795 ASSERT(grid().tree().solver2grid(cellId) > -1, "");
796 ASSERT(grid().tree().solver2grid(otherId) < 0 || solver().a_isHalo(otherId), "");
797 ASSERT(!solver().a_isHalo(cellId), "");
798 ASSERT(isToDelete(otherId) || solver().a_isHalo(otherId), "");
799 }
800 ASSERT(!solver().a_isHalo(cellId) && !isToDelete(cellId), "");
801 }
802
803 // 3. remove holes in the range of halo cells
804 otherId = solver().m_cells.size() - 1;
805 for(MInt cellId = noInternalCells; cellId < noCells; cellId++) {
806 if(isToDelete(cellId) != 0) {
807 while(isToDelete(otherId) != 0) {
808 otherId--;
809 }
810 ASSERT(cellId < otherId, "");
811 ASSERT(otherId >= noCells, "");
812 ASSERT(solver().a_isHalo(otherId) && !isToDelete(otherId), "");
813 solver().swapCells(cellId, otherId);
814 grid().swapSolverIds(cellId, otherId);
815 std::swap(oldCellId(cellId), oldCellId(otherId));
816 std::swap(isToDelete(cellId), isToDelete(otherId));
817
818 ASSERT(grid().tree().solver2grid(cellId) > -1, "");
819 ASSERT(grid().tree().solver2grid(otherId) < 0, "");
820 }
821 ASSERT(solver().a_isHalo(cellId) && !isToDelete(cellId), "");
822 }
823
824 for(MInt cellId = noInternalCells; cellId < noCells; cellId++) {
825 ASSERT(grid().tree().solver2grid(cellId) > -1
826 && grid().tree().solver2grid(cellId) < solver().grid().raw().treeb().size(),
827 "");
828 ASSERT(solver().a_isHalo(cellId) == solver().grid().raw().a_isHalo(grid().tree().solver2grid(cellId)), "");
829 }
830
831 solver().m_cells.size(noCells);
832 m_gridProxy.resizeGridMap(solver().m_cells.size());
833 ASSERT(solver().m_cells.size() == m_gridProxy.tree().size(), "");
834
835
836 if(solver().grid().raw().treeb().noSolvers() == 1) {
837 for(MInt gridCellId = 0; gridCellId < noCells; gridCellId++) {
838 ASSERT(grid().tree().solver2grid(gridCellId) == gridCellId, "");
839 ASSERT(grid().tree().grid2solver(gridCellId) == gridCellId, "");
840 }
841 }
842
843
844 /*
845 MIntScratchSpace newCellId(m_gridProxy.maxNoCells(), AT_, "newCellId");
846 newCellId.fill(-1);
847 for( MInt cellId = 0; cellId < noCells; cellId++ ) {
848 if ( oldCellId( cellId ) < 0 ) continue;
849 newCellId( oldCellId( cellId ) ) = cellId;
850 }
851 for(MInt i=0; i < noNeighborDomains(); i++) {
852 for ( MInt j = 0; j < (signed)m_gridProxy.m_haloCells[i].size(); j++ ) {
853 MInt cellId = m_gridProxy.m_haloCells[i][j];
854 if ( newCellId(cellId) > -1 ) {
855 m_gridProxy.m_haloCells[i][j] = newCellId(cellId);
856 }
857 }
858 for ( MInt j = 0; j < (signed)m_gridProxy.m_windowCells[i].size(); j++ ) {
859 MInt cellId = m_gridProxy.m_windowCells[i][j];
860 if ( newCellId(cellId) > -1 ) {
861 m_gridProxy.m_windowCells[i][j] = newCellId(cellId);
862 }
863 }
864 }
865 */
866
867 // m_gridProxy.tree().size(noCells);
868 // m_gridProxy.m_noInternalCells = noInternalCells;
869}
870
871
872// --------------------------------------------------------------------------------------
873
874
879template <MInt nDim, class SolverType>
880MInt CartesianSolver<nDim, SolverType>::createCellId(const MInt gridCellId) {
881 MInt solverCellId = -1;
882 if(solver().m_freeIndices.size() > 0) {
883 auto it = solver().m_freeIndices.begin();
884 solverCellId = *(it);
885 solver().m_freeIndices.erase(it);
886 m_gridProxy.resizeGridMap(solver().m_cells.size());
887 } else {
888 solverCellId = solver().m_cells.size();
889 solver().m_cells.append();
890 m_gridProxy.resizeGridMap(solver().m_cells.size());
891 }
892 ASSERT(solverCellId > -1 && solverCellId < solver().m_cells.size(), "");
893 if(!g_multiSolverGrid) {
894 ASSERT(solverCellId == gridCellId, std::to_string(solverCellId) + " " + std::to_string(gridCellId));
895 }
896
897 solver().a_resetPropertiesSolver(solverCellId);
898 solver().a_isHalo(solverCellId) = solver().grid().raw().a_isHalo(gridCellId);
899
900 grid().setSolver2grid(solverCellId, gridCellId);
901 grid().setGrid2solver(gridCellId, solverCellId);
902
903 return solverCellId;
904}
905
906
911template <MInt nDim, class SolverType>
912void CartesianSolver<nDim, SolverType>::removeCellId(const MInt cellId) {
913 const MInt gridCellId = grid().tree().solver2grid(cellId);
914 ASSERT(gridCellId > -1 && gridCellId < solver().grid().raw().treeb().size(), "");
915
916 solver().a_resetPropertiesSolver(cellId);
917
918 grid().setSolver2grid(cellId, std::numeric_limits<MInt>::min());
919 grid().setGrid2solver(gridCellId, std::numeric_limits<MInt>::min());
920 solver().grid().raw().treeb().solver(gridCellId, solver().solverId()) = false;
921
922 if(cellId == (solver().m_cells.size() - 1)) {
923 solver().m_cells.size(solver().m_cells.size() - 1);
924 m_gridProxy.resizeGridMap(solver().m_cells.size());
925 } else {
926 solver().m_freeIndices.insert(cellId);
927 }
928
929 ASSERT(g_multiSolverGrid || cellId == gridCellId, "");
930}
931
933// and sets the approximate distance to the interfaceCells!
938
939template <MInt nDim, class SolverType>
940void CartesianSolver<nDim, SolverType>::setBoundaryDistance(const MBool* const interfaceCell,
941 const MFloat* const outerBandWidth,
942 MFloatScratchSpace& distance) {
943 TRACE();
944
945 distance.fill(std::numeric_limits<MFloat>::max());
946 std::vector<MInt> currentLayer;
947 MIntScratchSpace onCurrentLayer(solver().a_noCells(), AT_, "onCurrentLayer");
948 MIntScratchSpace accessed(solver().a_noCells(), AT_, "accessed");
949 onCurrentLayer.fill(0);
950 accessed.fill(0);
951
952 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
953 if(interfaceCell[cellId]) {
954 currentLayer.push_back(cellId);
955 onCurrentLayer(cellId) = 1;
956 accessed(cellId) = 1;
957 distance(cellId) = F0;
958 }
959 }
960 MInt proceed = 1;
961 while(proceed != 0) {
962 proceed = 0;
963 solver().exchangeData(distance.data());
964 solver().exchangeData(onCurrentLayer.data());
965 for(MInt i = 0; i < solver().noNeighborDomains(); i++) {
966 for(MInt c = 0; c < solver().noHaloCells(i); c++) {
967 if(onCurrentLayer(solver().haloCellId(i, c))) {
968 currentLayer.push_back(solver().haloCellId(i, c));
969 accessed(solver().haloCellId(i, c)) = 1;
970 }
971 }
972 }
973 onCurrentLayer.fill(0);
974 std::vector<MInt> currentLayerBak(currentLayer);
975 for(auto& cellId : currentLayerBak) {
976 for(MInt dir = 0; dir < solver().m_noDirs; dir++) {
977 if(solver().a_hasNeighbor(cellId, dir)) {
978 MInt nghbrId = solver().c_neighborId(cellId, dir);
979 if(solver().c_isLeafCell(cellId)) {
980 MFloat dx = F1B2 * (solver().c_cellLengthAtCell(cellId) + solver().c_cellLengthAtCell(nghbrId));
981 distance(nghbrId) = mMin(distance(nghbrId), distance(cellId) + dx);
982 if(accessed(nghbrId) == 0) {
983 currentLayer.push_back(nghbrId);
984 accessed(nghbrId) = 1;
985 onCurrentLayer(nghbrId) = 1;
986 if(distance(nghbrId) < outerBandWidth[mMax(solver().minLevel(), solver().a_level(nghbrId) - 1)]
987 + F2 * solver().c_cellLengthAtCell(nghbrId)) {
988 proceed = 1;
989 }
990 }
991 } else {
992 for(MInt c = 0; c < solver().grid().m_maxNoChilds; c++) {
993 if(!childCode[dir][c]) {
994 continue;
995 }
996 MInt childId = solver().c_childId(nghbrId, c);
997 if(childId < 0) {
998 continue;
999 }
1000 MFloat dx = F1B2 * (solver().c_cellLengthAtCell(cellId) + solver().c_cellLengthAtCell(childId));
1001 distance(childId) = mMin(distance(childId), distance(cellId) + dx);
1002 if(accessed(childId) == 0) {
1003 currentLayer.push_back(childId);
1004 accessed(childId) = 1;
1005 onCurrentLayer(childId) = 1;
1006 if(distance(childId) < outerBandWidth[mMax(solver().minLevel(), solver().a_level(childId) - 1)]
1007 + F2 * solver().c_cellLengthAtCell(childId)) {
1008 proceed = 1;
1009 }
1010 }
1011 }
1012 }
1013 } else {
1014 MInt parentId = solver().c_parentId(cellId);
1015 while(parentId > -1 && !solver().a_hasNeighbor(parentId, dir)) {
1016 parentId = solver().c_parentId(parentId);
1017 }
1018 if(parentId > -1 && solver().a_hasNeighbor(parentId, dir)) {
1019 MInt nghbrId = solver().c_neighborId(parentId, dir);
1020 if(!solver().c_isLeafCell(nghbrId)) {
1021 continue;
1022 }
1023 MFloat dx = F1B2 * (solver().c_cellLengthAtCell(cellId) + solver().c_cellLengthAtCell(nghbrId));
1024 distance(nghbrId) = mMin(distance(nghbrId), distance(cellId) + dx);
1025 if(accessed(nghbrId) == 0) {
1026 currentLayer.push_back(nghbrId);
1027 accessed(nghbrId) = 1;
1028 onCurrentLayer(nghbrId) = 1;
1029 if(distance(nghbrId) < outerBandWidth[mMax(solver().minLevel(), solver().a_level(nghbrId) - 1)]
1030 + F2 * solver().c_cellLengthAtCell(nghbrId)) {
1031 proceed = 1;
1032 }
1033 }
1034 }
1035 }
1036 }
1037 }
1038 if(currentLayer.empty()) {
1039 proceed = 0;
1040 }
1041 MPI_Allreduce(MPI_IN_PLACE, &proceed, 1, MPI_INT, MPI_SUM, solver().mpiComm(), AT_, "MPI_IN_PLACE", "proceed");
1042 }
1043}
1044
1045
1055
1056template <MInt nDim, class SolverType>
1057void CartesianSolver<nDim, SolverType>::markSurrndCells(MIntScratchSpace& inList, const MInt bandWidth,
1058 const MInt level, const MBool refineDiagonals) {
1059 TRACE();
1060 static constexpr MInt noDirs = 2 * nDim;
1061 static constexpr std::array<MInt, 4> diagDirs{3, 2, 0, 1};
1062
1063
1064 // Loop over the number of refinement-grid-cells and add those to the list:
1065 for(MInt loopMarker = 1; loopMarker < bandWidth; loopMarker++) {
1066 for(MInt cellId = 0; cellId < solver().a_noCells(); cellId++) {
1067 if(solver().grid().tree().level(cellId) != level) {
1068 continue;
1069 }
1070 if(inList(cellId) != loopMarker) {
1071 continue;
1072 }
1073
1074 const MInt cellDist = loopMarker + 1;
1075
1076 // direct neighbors +x-x+y-y+z-z
1077 for(MInt mainDir = 0; mainDir < noDirs; mainDir++) {
1078 const MInt nghbrId = solver().c_neighborId(cellId, mainDir);
1079 if(nghbrId < 0) {
1080 continue;
1081 }
1082 if(inList(nghbrId) == 0) {
1083 inList(nghbrId) = cellDist;
1084 }
1085
1086 if(refineDiagonals) {
1087 if(mainDir < 4) {
1088 // add diagonal neighbors in x-y plane
1089 const MInt nghbrId2 = solver().c_neighborId(nghbrId, diagDirs.at(mainDir));
1090 if(nghbrId2 >= 0 && inList(nghbrId2) == 0) {
1091 inList(nghbrId2) = cellDist;
1092 }
1093 }
1094
1095 IF_CONSTEXPR(nDim == 3) {
1096 for(MInt dirZ = 4; dirZ < 6; dirZ++) {
1097 const MInt nghbrIdZ = solver().c_neighborId(cellId, dirZ);
1098 if(nghbrIdZ < 0) {
1099 continue;
1100 }
1101 if(nghbrIdZ >= 0 && inList(nghbrIdZ) == 0) {
1102 inList(nghbrIdZ) = cellDist;
1103 }
1104 for(MInt dir = 0; dir < 4; dir++) {
1105 const MInt nghbrId2 = solver().c_neighborId(nghbrIdZ, dir);
1106 if(nghbrId2 < 0) {
1107 continue;
1108 }
1109 if(inList(nghbrId2) == 0) {
1110 inList(nghbrId2) = cellDist;
1111 }
1112 // tridiagonal neighbors
1113 const MInt triDiagNghbrId = solver().c_neighborId(nghbrId2, diagDirs.at(dir));
1114 if(triDiagNghbrId >= 0 && inList(triDiagNghbrId) == 0) {
1115 inList(triDiagNghbrId) = cellDist;
1116 }
1117 }
1118 }
1119 }
1120 }
1121 }
1122 }
1123 exchangeData(inList.data());
1124 if(grid().azimuthalPeriodicity()) {
1125 exchangeAzimuthalPer(&inList[0]);
1126 }
1127 }
1128}
1129
1130
1136template <MInt nDim, class SolverType>
1137void CartesianSolver<nDim, SolverType>::checkNoHaloLayers() {
1138#ifndef NDEBUG
1139
1140 if(Context::propertyExists("periodicCells")) return;
1141 // NOTE: not working for Alexejs periodicExchange! As periodic cells there are only created
1142 // on solver level! The functionality there should be moved into the cartesiangrid!
1143 // and once this is done there, the additional exchange Cells are part of the proxy!
1144
1145
1146 // NOTE: searching only on the same level is not correct if the
1147 // domainBoundary is at a lower level and the first haloLayer is refined,
1148 // the correct number of leafCell-halo-layers needs to be ensured!
1149
1150 m_log << "check HaloLayer-Count for solver " << solver().solverId() << " ... ";
1151
1152 MBool uncorrectLayerCount = false;
1153 static constexpr MInt cornerIndices[8][3] = {{-1, -1, -1}, {1, -1, -1}, {-1, 1, -1}, {1, 1, -1},
1154 {-1, -1, 1}, {1, -1, 1}, {-1, 1, 1}, {1, 1, 1}};
1155
1156 /*
1157 std::multimap< MInt, MInt > firstHaloLayerOld;
1158 firstHaloLayerOld.clear();
1159
1160 //loop over all haloCells and find cells which have a window cell as neighbor
1161 for(MInt cellId = solver().noInternalCells(); cellId < solver().c_noCells(); cellId++) {
1162 ASSERT( solver().a_isHalo(cellId), "");
1163 //only check for solver leaf cells
1164 if(!grid().tree().isLeafCell(cellId)) continue;
1165 for(MInt dir = 0; dir < 2*nDim; dir++) {
1166 MInt nghbrId = -1;
1167 if(grid().tree().hasNeighbor(cellId, dir)) {
1168 nghbrId = grid().tree().neighbor(cellId,dir);
1169 } else if ( grid().tree().hasParent(cellId) && grid().tree().parent(cellId) < grid().tree().size() ) {
1170 nghbrId = grid().tree().neighbor(grid().tree().parent(cellId),dir);
1171 }
1172
1173 if ( nghbrId < 0 ) continue;
1174
1175
1176 //if(!grid().tree().isLeafCell(nghbrId) ){
1177 //get the adjacent Childs
1178 // for ( MInt child = 0; child < ipow(2, nDim); child++) {
1179 // if ( !childCode[ dir ][ child ] ) continue;
1180 // if ( grid().tree().child( nghbrId , child ) > -1 ) {
1181 // const MInt childId = grid().tree().child( nghbrId , child );
1182 // if(grid().tree().hasProperty(childId, Cell::IsWindow)){
1183 // firstHaloLayerOld.insert( std::make_pair( cellId, m_revDir[dir]));
1184 // }
1185 // }
1186 // }
1187 //} else
1188 if(grid().tree().hasProperty(nghbrId, Cell::IsWindow)){
1189 firstHaloLayerOld.insert(std::make_pair(cellId, m_revDir[dir]));
1190 ASSERT(cellId > -1, "");
1191 ASSERT(m_revDir[dir] > -1 && m_revDir[dir] < 2*nDim, std::to_string(m_revDir[dir]));
1192 }
1193 }
1194 }
1195
1196 for ( std::multimap<MInt,MInt>::iterator it = firstHaloLayerOld.begin();
1197 it != firstHaloLayerOld.end(); it++ ) {
1198 MInt cellId = it->first;
1199 const MInt dir = it->second;
1200 //const MInt level = grid().tree().level(cellId);
1201 ASSERT(dir > -1 && dir < 2*nDim, std::to_string(dir));
1202
1203 if(!grid().tree().isLeafCell(cellId) ) continue;
1204
1205 for(MInt layer = 1; layer < grid().noHaloLayers(); layer++){
1206 MInt nghbrId = -1;
1207
1208 if(grid().tree().hasNeighbor(cellId, dir)) {
1209 nghbrId = grid().tree().neighbor(cellId,dir);
1210 } else if (grid().tree().hasParent(cellId) && grid().tree().parent(cellId) < grid().tree().size()) {
1211 nghbrId = grid().tree().neighbor(grid().tree().parent(cellId),dir);
1212 }
1213
1214 if(nghbrId > -1 && !grid().tree().isLeafCell(nghbrId) ){
1215 //check that the two childs are present!
1216 for ( MInt child = 0; child < ipow(2, nDim); child++) {
1217 if ( !childCode[ dir ][ child ] ) continue;
1218 if ( grid().tree().child( nghbrId , child ) < 0 ) {
1219 std::cerr << "Incorrect halo-Layer count: 1 " << cellId << " " << solver().domainId() << " "
1220 << grid().tree().globalId(cellId) << " " << solver().solverId() << " " << dir << std::endl;
1221 uncorrectLayerCount = true;
1222 break;
1223 } else {
1224 //continue the leayer check from one of the childs onward!
1225 cellId = grid().tree().child( nghbrId , child );
1226 }
1227 }
1228 } else if(nghbrId < 0 || nghbrId > grid().tree().size()) {
1229 //check if the cell is at the domain-bndry
1230 //meaning that the potential neighbor would have to be outside the geometry!
1231
1232 //prepare to call pointIsInside in the geometry class:
1233 MBool outside = true;
1234 const MFloat cellHalfLength = F1B2*grid().cellLengthAtCell(cellId);
1235 MFloat corner[ 3 ] = { 0,0,0 };
1236
1237 //computing the coordinates the neighbor should be having!
1238 MFloat coords[nDim];
1239 for ( MInt k = 0; k < nDim; k++ ) {
1240 coords[k] = grid().tree().coordinate( cellId, k );
1241 }
1242 coords[dir/2] += ((dir%2==0)?-F1:F1) * grid().cellLengthAtCell(cellId);
1243
1244 for(MInt node = 0; node < ipow(2, nDim); node++){
1245 for( MInt dim = 0; dim < nDim; dim++ ){
1246 corner[ dim ] = coords[dim] + cornerIndices[ node ][dim] * cellHalfLength;
1247 }
1248 IF_CONSTEXPR(nDim == 2) {
1249 if(! solver().geometry().pointIsInside( corner )) outside = false;
1250 } else {
1251 if ( ! solver().geometry().pointIsInside2( corner )) outside = false;
1252 }
1253 }
1254
1255 if(outside) {
1256 break;
1257 } else {
1258
1259 //check if the cell is at the cutOff
1260 if(grid().hasCutOff() && grid().tree().cutOff(cellId) > -1) {
1261 break;
1262 }
1263
1264
1265 std::cerr << "Incorrect halo-Layer count: 2 " << cellId << " " << solver().domainId() << " "
1266 << grid().tree().globalId(cellId) << " " << solver().solverId() << " " << dir << " "
1267 << layer << std::endl;
1268 uncorrectLayerCount = true;
1269 std::cerr << "Coord: " << coords[0] << " " << coords[1] << " " << coords[nDim-1] << std::endl;
1270 break;
1271 }
1272 } else {
1273 //check passed, continue with next layer
1274 cellId = nghbrId;
1275 }
1276
1277
1278
1279
1280 }
1281 }
1282
1283 if(uncorrectLayerCount) {
1284 mTerm(1,AT_, "Incorrect number of halo-layers!");
1285 }
1286*/
1287
1288 std::set<MInt> haloCells;
1289
1290 std::multimap<MInt, MInt> firstHaloLayer;
1291 firstHaloLayer.clear();
1292
1293 // Create a list of all halo cells to be checked.
1294 // Add regular halo cells
1295 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
1296 for(MInt j = 0; j < grid().noHaloCells(i); j++) {
1297 const MInt cellId = grid().haloCell(i, j);
1298 haloCells.insert(cellId);
1299 }
1300 }
1301 // Add azimuthal periodic halo cells
1302 for(MInt i = 0; i < grid().noAzimuthalNeighborDomains(); i++) {
1303 for(MInt j = 0; j < grid().noAzimuthalHaloCells(i); j++) {
1304 const MInt cellId = grid().azimuthalHaloCell(i, j);
1305 haloCells.insert(cellId);
1306 }
1307 }
1308 for(MInt i = 0; i < grid().noAzimuthalUnmappedHaloCells(); i++) {
1309 const MInt cellId = grid().azimuthalUnmappedHaloCell(i);
1310 haloCells.insert(cellId);
1311 }
1312
1313 // loop over all haloCells and find cells which have a window cell as neighbor
1314 // these cells are the first layer of haloCells.
1315 for(auto it = haloCells.begin(); it != haloCells.end(); ++it) {
1316 MInt cellId = *it;
1317 if(!grid().tree().isLeafCell(cellId)) continue;
1318 ASSERT(solver().a_isHalo(cellId), "");
1319
1320 // skip partition levell ancestors!
1321 if(grid().tree().hasProperty(cellId, Cell::IsPartLvlAncestor)) continue;
1322
1323 for(MInt dir = 0; dir < 2 * nDim; dir++) {
1324 MInt nghbrId = -1; // adjacent neighbor
1325 // MInt nghbrId1 = -1; //adjacent child1
1326 // MInt nghbrId2 = -1; //adjacent child2
1327 if(grid().tree().hasNeighbor(cellId, dir)) {
1328 // direct neighbor
1329 nghbrId = grid().tree().neighbor(cellId, dir);
1330 // if the direct neighbor has children, get the two closest childs
1331 /*
1332 if(!grid().tree().isLeafCell(nghbrId) ){
1333 for ( MInt child = 0; child < ipow(2, nDim); child++) {
1334 if( !childCode[ dir ][ child ] ) continue;
1335 if( grid().tree().child( nghbrId , child ) < 0 ) continue;
1336 if(nghbrId1 < 0) {
1337 nghbrId1 = grid().tree().child( nghbrId , child );
1338 } else {
1339 ASSERT(nghbrId2 < 0, "");
1340 nghbrId2 = grid().tree().child( nghbrId , child );
1341 }
1342 }
1343 }
1344 */
1345 } else if(grid().tree().hasParent(cellId) && grid().tree().parent(cellId) < grid().tree().size()) {
1346 nghbrId = grid().tree().neighbor(grid().tree().parent(cellId), dir);
1347 }
1348 if(nghbrId < 0) continue;
1349
1350 // use childs
1351 // if(nghbrId1 > 0 && grid().tree().hasProperty(nghbrId1, Cell::IsWindow) &&
1352 // grid().tree().isLeafCell(cellId)){
1353 // firstHaloLayer.insert( std::make_pair( std::make_pair(cellId, i), m_revDir[dir]));
1354 //}
1355 // if(nghbrId2 > 0 && grid().tree().hasProperty(nghbrId2, Cell::IsWindow) &&
1356 // grid().tree().isLeafCell(cellId)){
1357 // firstHaloLayer.insert( std::make_pair( std::make_pair(cellId, i), m_revDir[dir]));
1358 //}
1359 // if(/*nghbrId1 < 0 && */solver().a_isWindow(nghbrId) /*&&
1360 // grid().tree().isLeafCell(cellId)*/){
1361 if(!grid().tree().hasProperty(nghbrId, Cell::IsHalo)) {
1362 firstHaloLayer.insert(std::make_pair(cellId, m_revDir[dir]));
1363 }
1364 }
1365 }
1366
1367
1368 for(std::multimap<MInt, MInt>::iterator it = firstHaloLayer.begin(); it != firstHaloLayer.end(); it++) {
1369 MInt cellId = it->first;
1370 // const MInt nghbrDomain = it->first.second;
1371 const MInt dir = it->second;
1372 ASSERT(dir > -1 && dir < 2 * nDim, std::to_string(dir));
1373 ASSERT(grid().tree().isLeafCell(cellId), "");
1374 ASSERT(solver().a_isHalo(cellId), "");
1375
1376 for(MInt layer = 1; layer < grid().noHaloLayers(); layer++) {
1377 MInt nghbrId = -1; // adjacent neighbor
1378 MInt nghbrId1 = -1; // adjacent child1
1379 // MInt nghbrId2 = -1; //adjacent child2
1380 if(grid().tree().hasNeighbor(cellId, dir)) {
1381 nghbrId = grid().tree().neighbor(cellId, dir);
1382 if(!grid().tree().isLeafCell(nghbrId)) {
1383 // check the two childs!
1384 for(MInt child = 0; child < ipow(2, nDim); child++) {
1385 if(!childCode[dir][child]) continue;
1386 if(grid().tree().child(nghbrId, child) < 0) {
1387 if(grid().noHaloLayers() % 2 == 0) {
1388 // check if the neighboring child could be outside:
1389 MBool outside = true;
1390 const MFloat cellHalfLength = F1B4 * grid().cellLengthAtCell(cellId);
1391 MFloat corner[3] = {0, 0, 0};
1392
1393 // computing the coordinates the neighbor should be having!
1394 MFloat coords[nDim];
1395 for(MInt k = 0; k < nDim; k++) {
1396 coords[k] = grid().tree().coordinate(nghbrId, k);
1397 }
1398 for(MInt k = 0; k < nDim; k++) {
1399 coords[k] += cornerIndices[child][k] * cellHalfLength;
1400 }
1401
1402 for(MInt node = 0; node < ipow(2, nDim); node++) {
1403 for(MInt dim = 0; dim < nDim; dim++) {
1404 corner[dim] = coords[dim] + cornerIndices[node][dim] * F1B2 * cellHalfLength;
1405 }
1406 IF_CONSTEXPR(nDim == 2) {
1407 if(!solver().geometry().pointIsInside(corner)) outside = false;
1408 }
1409 else {
1410 if(!solver().geometry().pointIsInside2(corner)) outside = false;
1411 }
1412 }
1413
1414 if(outside) continue;
1415
1416 uncorrectLayerCount = true;
1417 std::cerr << "Incorrect halo-Layer count: 1 " << cellId << " " << solver().domainId() << " "
1418 << grid().tree().globalId(cellId) << " " << solver().solverId() << " " << dir << std::endl;
1419 std::cerr << grid().tree().coordinate(cellId, 0) << " " << grid().tree().coordinate(cellId, 1) << " "
1420 << grid().tree().coordinate(cellId, nDim - 1) << std::endl;
1421 break;
1422 }
1423 } else {
1424 if(nghbrId1 < 0) {
1425 nghbrId1 = grid().tree().child(nghbrId, child);
1426 } /*else {
1427 ASSERT(nghbrId2 < 0, "");
1428 nghbrId2 = grid().tree().child( nghbrId , child );
1429 }
1430 */
1431 }
1432 }
1433 }
1434 } else if(grid().tree().hasParent(cellId) && grid().tree().parent(cellId) < grid().tree().size()) {
1435 nghbrId = grid().tree().neighbor(grid().tree().parent(cellId), dir);
1436 }
1437
1438 // check if the cell is at the domain-bndry
1439 // meaning that the potential neighbor would have to be outside the geometry!
1440 if(nghbrId < 0 || nghbrId > grid().tree().size()) {
1441 // prepare to call pointIsInside in the geometry class:
1442 MBool outside = true;
1443 const MFloat cellHalfLength = F1B2 * grid().cellLengthAtCell(cellId);
1444 MFloat corner[3] = {0, 0, 0};
1445
1446 // computing the coordinates the neighbor should be having!
1447 MFloat coords[nDim];
1448 for(MInt k = 0; k < nDim; k++) {
1449 coords[k] = grid().tree().coordinate(cellId, k);
1450 }
1451 coords[dir / 2] += ((dir % 2 == 0) ? -F1 : F1) * grid().cellLengthAtCell(cellId);
1452
1453 for(MInt node = 0; node < ipow(2, nDim); node++) {
1454 for(MInt dim = 0; dim < nDim; dim++) {
1455 corner[dim] = coords[dim] + cornerIndices[node][dim] * cellHalfLength;
1456 }
1457 IF_CONSTEXPR(nDim == 2) {
1458 if(!solver().geometry().pointIsInside(corner)) outside = false;
1459 }
1460 else {
1461 if(!solver().geometry().pointIsInside2(corner)) outside = false;
1462 }
1463 }
1464
1465 if(outside) {
1466 break;
1467 } else {
1468 // check if the cell is at the cutOff
1469 if(grid().hasCutOff() && grid().tree().cutOff(cellId) > -1) {
1470 break;
1471 }
1472 }
1473 }
1474
1475 // if the neighbor is not present at all
1476 if(nghbrId < 0) {
1477 uncorrectLayerCount = true;
1478 std::cerr << "Incorrect halo-Layer count: 2 " << solver().domainId() << " " << solver().solverId() << " "
1479 << layer << " " << cellId << "/" << grid().tree().globalId(cellId) << " " << dir << " "
1480 << grid().tree().coordinate(cellId, 0) << " " << grid().tree().coordinate(cellId, 1) << " "
1481 << grid().tree().coordinate(cellId, nDim - 1) << " " << grid().isPeriodic(cellId) << std::endl;
1482 break;
1483 }
1484
1485 /*
1486 //check if the corresponding neighbor is a window cell,
1487 //meaning that the direction changed due to a domain-corner!
1488 if(nghbrId1 < 0 ) {
1489 if(!solver().a_isHalo(nghbrId)) {
1490 ASSERT(grid().tree().hasProperty(nghbrId, Cell::IsWindow), "");
1491 break;
1492 }
1493 } else if ( nghbrId1 > -1 ) {
1494 if(!solver().a_isHalo(nghbrId1)) {
1495 //neighbor appears to be a window cell again, meaning, that we went back into the domain.
1496 ASSERT(grid().tree().hasProperty(nghbrId1, Cell::IsWindow), "");
1497 break;
1498 }
1499 }
1500 if(nghbrId2 > -1 ) {
1501 //check ngbrId1
1502 if(!solver().a_isHalo(nghbrId2)) {
1503 ASSERT(grid().tree().hasProperty(nghbrId2, Cell::IsWindow), "");
1504 break;
1505 }
1506 }
1507
1508 // check that the ngbrId is also a haloCell and is a halo for the correct domain!
1509 MBool cellFound = false;
1510 MBool cellFound2 = false;
1511 if(nghbrId1 < 0 ) {
1512 //check ngbrId
1513 for(MInt j = 0; j < grid().noHaloCells(nghbrDomain); j++){
1514 const MInt haloCell = grid().haloCell(nghbrDomain,j);
1515 if(haloCell == nghbrId) {
1516 cellFound = true;
1517 break;
1518 }
1519 }
1520 } else if ( nghbrId1 > -1 ) {
1521 //check ngbrId1
1522 for(MInt j = 0; j < grid().noHaloCells(nghbrDomain); j++){
1523 const MInt haloCell = grid().haloCell(nghbrDomain,j);
1524 if(haloCell == nghbrId1) {
1525 cellFound = true;
1526 break;
1527 }
1528 }
1529 }
1530 if(nghbrId2 > -1 ) {
1531 //check ngbrId1
1532 for(MInt j = 0; j < grid().noHaloCells(nghbrDomain); j++){
1533 const MInt haloCell = grid().haloCell(nghbrDomain,j);
1534 if(haloCell == nghbrId2) {
1535 cellFound2 = true;
1536 break;
1537 }
1538 }
1539 }
1540
1541 if(!cellFound) {
1542 std::cerr << "Incorrect halo-Layer count: 3 " << cellId << " " << solver().domainId() << " "
1543 << grid().tree().globalId(cellId) << " " << solver().solverId() << " " << dir << " "
1544 << layer << " " << nghbrId << " " << solver().a_isHalo(nghbrId) << std::endl;
1545 uncorrectLayerCount = true;
1546 break;
1547 }
1548
1549 if(nghbrId2 > -1 && !cellFound2) {
1550 std::cerr << "Incorrect halo-Layer count: 4 " << cellId << " " << solver().domainId() << " "
1551 << grid().tree().globalId(cellId) << " " << solver().solverId() << " " << dir << " "
1552 << layer << " " << nghbrId << std::endl;
1553 uncorrectLayerCount = true;
1554 break;
1555 }
1556 */
1557
1558 // continue towards the next layer
1559 if(nghbrId1 < 0) {
1560 cellId = nghbrId;
1561 } else {
1562 // two remaining childs to iterate forward
1563 // for now just the first child is used...
1564 cellId = nghbrId1;
1565 }
1566 }
1567 }
1568
1569
1570 if(uncorrectLayerCount) {
1571 mTerm(1, AT_, "Incorrect number of halo-layers!");
1572 }
1573
1574
1575 m_log << "done" << std::endl;
1576
1577#endif
1578}
1579
1580// Check if point is located in cell
1581template <MInt nDim, class SolverType>
1582MInt CartesianSolver<nDim, SolverType>::inCell(const MInt cellId, MFloat* point, MFloat fac) {
1583 TRACE();
1584
1585 MFloat xmin, xmax;
1586 MFloat ymin, ymax;
1587 MFloat zmin, zmax;
1588 MFloat cellHalfLength = F1B2 * grid().cellLengthAtCell(cellId);
1589
1590 xmin = grid().tree().coordinate(cellId, 0) - (cellHalfLength * fac);
1591 xmax = grid().tree().coordinate(cellId, 0) + (cellHalfLength * fac);
1592 ymin = grid().tree().coordinate(cellId, 1) - (cellHalfLength * fac);
1593 ymax = grid().tree().coordinate(cellId, 1) + (cellHalfLength * fac);
1594
1595 if(nDim == 2) {
1596 if(point[0] < xmax && point[0] >= xmin && point[1] < ymax && point[1] >= ymin) {
1597 return true;
1598 } else {
1599 return false;
1600 }
1601 } else {
1602 zmin = grid().tree().coordinate(cellId, 2) - (cellHalfLength * fac);
1603 zmax = grid().tree().coordinate(cellId, 2) + (cellHalfLength * fac);
1604
1605 if(point[0] < xmax && point[0] >= xmin && point[1] < ymax && point[1] >= ymin && point[2] < zmax
1606 && point[2] >= zmin) {
1607 return true;
1608 } else {
1609 return false;
1610 }
1611 }
1612}
1613
1614
1621
1622template <MInt nDim, class SolverType>
1623MInt CartesianSolver<nDim, SolverType>::setUpInterpolationStencil(const MInt cellId,
1624 MInt* interpolationCells,
1625 const MFloat* point,
1626 std::function<MBool(MInt, MInt)>
1627 neighborCheck,
1628 MBool allowIncompleteStencil) {
1629 TRACE();
1630
1631 MInt position = 0;
1632 const MInt maxPosition = IPOW2(nDim) - 1;
1633 MBool invalidStencil = false;
1634
1635 //------------------------------
1636
1637 if(point[0] > grid().tree().coordinate(cellId, 0)) {
1638 position += 1;
1639 }
1640 if(point[1] > grid().tree().coordinate(cellId, 1)) {
1641 position += 2;
1642 }
1643
1644 IF_CONSTEXPR(nDim == 3) {
1645 if(point[2] > grid().tree().coordinate(cellId, 2)) {
1646 position += 4;
1647 }
1648 }
1649
1650 position = maxPosition - position;
1651
1652 // check if cell is boundary cell and located at position at the boundary in the stencil. if yes, take other stencil.
1653
1654 MInt xNghbrDir = (position + 1) % 2;
1655 MInt posIncrementX = -1 + 2 * xNghbrDir;
1656 MInt yNghbrDir = (position / 2 + 1) % 2;
1657 MInt posIncrementY = -2 + 4 * yNghbrDir;
1658 yNghbrDir += 2;
1659 MInt zNghbrDir = -1;
1660 MInt posIncrementZ = -1;
1661
1662
1663 if(!neighborCheck(cellId, xNghbrDir) || !grid().tree().hasNeighbor(cellId, xNghbrDir)) {
1664 position += posIncrementX;
1665 }
1666 if(!neighborCheck(cellId, yNghbrDir) || !grid().tree().hasNeighbor(cellId, yNghbrDir)) {
1667 position += posIncrementY;
1668 }
1669
1670 IF_CONSTEXPR(nDim == 3) {
1671 zNghbrDir = (position / 4 + 1) % 2;
1672 posIncrementZ = -4 + 8 * zNghbrDir;
1673 zNghbrDir += 4;
1674
1675 if(!neighborCheck(cellId, zNghbrDir) || !grid().tree().hasNeighbor(cellId, zNghbrDir)) {
1676 position += posIncrementZ;
1677 }
1678 }
1679
1680 interpolationCells[position] = cellId;
1681
1682 IF_CONSTEXPR(nDim == 2) {
1683 MInt nghbrX = -1;
1684 MInt nghbrY = -1;
1685 if(position % 2 == 0) {
1686 xNghbrDir = 1;
1687 nghbrX = neighborCheck(cellId, xNghbrDir) ? grid().tree().neighbor(cellId, xNghbrDir) : -1;
1688 posIncrementX = 1;
1689 } else {
1690 xNghbrDir = 0;
1691 nghbrX = neighborCheck(cellId, xNghbrDir) ? grid().tree().neighbor(cellId, xNghbrDir) : -1;
1692 posIncrementX = -1;
1693 }
1694 if(((MInt)(position / 2)) % 2 == 0) {
1695 yNghbrDir = 3;
1696 nghbrY = neighborCheck(cellId, yNghbrDir) ? grid().tree().neighbor(cellId, yNghbrDir) : -1;
1697 posIncrementY = 2;
1698 } else {
1699 yNghbrDir = 2;
1700 nghbrY = neighborCheck(cellId, yNghbrDir) ? grid().tree().neighbor(cellId, yNghbrDir) : -1;
1701 posIncrementY = -2;
1702 }
1703 interpolationCells[position + posIncrementX] = nghbrX;
1704 interpolationCells[position + posIncrementY] = nghbrY;
1705 if(nghbrX != -1)
1706 interpolationCells[position + posIncrementX + posIncrementY] =
1707 neighborCheck(nghbrX, yNghbrDir) ? grid().tree().neighbor(nghbrX, yNghbrDir) : -1;
1708 else if(nghbrY != -1)
1709 interpolationCells[position + posIncrementX + posIncrementY] =
1710 neighborCheck(nghbrY, xNghbrDir) ? grid().tree().neighbor(nghbrY, xNghbrDir) : -1;
1711 else
1712 interpolationCells[position + posIncrementX + posIncrementY] = -1;
1713 }
1714 else {
1715 MInt nghbrX = -1;
1716 MInt nghbrY = -1;
1717 MInt nghbrZ = -1;
1718 if(position % 2 == 0) {
1719 xNghbrDir = 1;
1720 nghbrX = neighborCheck(cellId, xNghbrDir) ? grid().tree().neighbor(cellId, xNghbrDir) : -1;
1721 posIncrementX = 1;
1722 } else {
1723 xNghbrDir = 0;
1724 nghbrX = neighborCheck(cellId, xNghbrDir) ? grid().tree().neighbor(cellId, xNghbrDir) : -1;
1725 posIncrementX = -1;
1726 }
1727 if(((MInt)(position / 2)) % 2 == 0) {
1728 yNghbrDir = 3;
1729 nghbrY = neighborCheck(cellId, yNghbrDir) ? grid().tree().neighbor(cellId, yNghbrDir) : -1;
1730 posIncrementY = 2;
1731 } else {
1732 yNghbrDir = 2;
1733 nghbrY = neighborCheck(cellId, yNghbrDir) ? grid().tree().neighbor(cellId, yNghbrDir) : -1;
1734 posIncrementY = -2;
1735 }
1736 if((MInt)(position / 4) == 0) {
1737 zNghbrDir = 5;
1738 nghbrZ = neighborCheck(cellId, zNghbrDir) ? grid().tree().neighbor(cellId, zNghbrDir) : -1;
1739 posIncrementZ = 4;
1740 } else {
1741 zNghbrDir = 4;
1742 nghbrZ = neighborCheck(cellId, zNghbrDir) ? grid().tree().neighbor(cellId, zNghbrDir) : -1;
1743 posIncrementZ = -4;
1744 }
1745 interpolationCells[position + posIncrementX] = nghbrX;
1746 interpolationCells[position + posIncrementY] = nghbrY;
1747 interpolationCells[position + posIncrementZ] = nghbrZ;
1748 if(nghbrX > -1) {
1749 interpolationCells[position + posIncrementX + posIncrementZ] =
1750 neighborCheck(nghbrX, zNghbrDir) ? grid().tree().neighbor(nghbrX, zNghbrDir) : -1;
1751 interpolationCells[position + posIncrementX + posIncrementY] =
1752 neighborCheck(nghbrX, yNghbrDir) ? grid().tree().neighbor(nghbrX, yNghbrDir) : -1;
1753 } else {
1754 interpolationCells[position + posIncrementX + posIncrementZ] = -1;
1755 interpolationCells[position + posIncrementX + posIncrementY] = -1;
1756 }
1757 if(nghbrY != -1) {
1758 interpolationCells[position + posIncrementY + posIncrementZ] =
1759 neighborCheck(nghbrY, zNghbrDir) ? grid().tree().neighbor(nghbrY, zNghbrDir) : -1;
1760 } else {
1761 interpolationCells[position + posIncrementY + posIncrementZ] = -1;
1762 }
1763 const MInt nghbrXY = interpolationCells[position + posIncrementX + posIncrementY];
1764 if(nghbrX > -1 && nghbrXY > -1)
1765 interpolationCells[position + posIncrementX + posIncrementY + posIncrementZ] =
1766 neighborCheck(nghbrXY, zNghbrDir) ? grid().tree().neighbor(nghbrXY, zNghbrDir) : -1;
1767 else
1768 interpolationCells[position + posIncrementX + posIncrementY + posIncrementZ] = -1;
1769 }
1770
1771 for(MInt p = 0; p <= maxPosition; p++) {
1772 if(interpolationCells[p] == -1) {
1773 invalidStencil = true;
1774 break;
1775 }
1776 }
1777 if(invalidStencil && !allowIncompleteStencil) {
1778 for(MInt p = 0; p <= maxPosition; p++) {
1779 interpolationCells[p] = cellId;
1780 }
1781 position = -1;
1782 }
1783
1784 return position;
1785}
1786
1787
1790// mode: cartesian interpolation = true:
1791// -> cartesian coordinates (c_coordinate/a_coordinate of non-boundary-cells!)
1792// based on pure cartesian cooordinates (i.e. known reconstruction constants!)
1793// NOTE: requires ordering of the cells in the cartesin-directions
1794// cartesian interpolation = false:
1797template <MInt nDim, class SolverType>
1798template <MBool cartesianInterpolation>
1799MFloat CartesianSolver<nDim, SolverType>::interpolateFieldData(MInt* interpolationCells, MFloat* point, MInt varId,
1800 std::function<MFloat(MInt, MInt)> scalarField,
1801 std::function<MFloat(MInt, MInt)> coordinate) {
1802 TRACE();
1803
1804 IF_CONSTEXPR(cartesianInterpolation) {
1805 const MFloat xMin = coordinate(interpolationCells[0], 0);
1806 const MFloat xMax = coordinate(interpolationCells[1], 0);
1807 const MFloat yMin = coordinate(interpolationCells[0], 1);
1808 const MFloat yMax = coordinate(interpolationCells[2], 1);
1809 const MFloat deltaX_Minus = point[0] - xMin;
1810 const MFloat deltaX_Plus = xMax - point[0];
1811 const MFloat deltaY_Minus = point[1] - yMin;
1812 const MFloat deltaY_Plus = yMax - point[1];
1813 const MFloat Delta = deltaX_Minus + deltaX_Plus;
1814
1815 IF_CONSTEXPR(nDim == 2) {
1816 const MFloat phi = (scalarField(interpolationCells[0], varId) * deltaX_Plus * deltaY_Plus
1817 + scalarField(interpolationCells[1], varId) * deltaX_Minus * deltaY_Plus
1818 + scalarField(interpolationCells[2], varId) * deltaX_Plus * deltaY_Minus
1819 + scalarField(interpolationCells[3], varId) * deltaX_Minus * deltaY_Minus)
1820 / (Delta * Delta);
1821 return phi;
1822 }
1823 else { // nDim == 3
1824 const MFloat zMin = coordinate(interpolationCells[0], 2);
1825 const MFloat zMax = coordinate(interpolationCells[4], 2);
1826 const MFloat deltaZ_Minus = point[2] - zMin;
1827 const MFloat deltaZ_Plus = zMax - point[2];
1828 const MFloat phi = (scalarField(interpolationCells[0], varId) * deltaX_Plus * deltaY_Plus * deltaZ_Plus
1829 + scalarField(interpolationCells[1], varId) * deltaX_Minus * deltaY_Plus * deltaZ_Plus
1830 + scalarField(interpolationCells[2], varId) * deltaX_Plus * deltaY_Minus * deltaZ_Plus
1831 + scalarField(interpolationCells[3], varId) * deltaX_Minus * deltaY_Minus * deltaZ_Plus
1832 + scalarField(interpolationCells[4], varId) * deltaX_Plus * deltaY_Plus * deltaZ_Minus
1833 + scalarField(interpolationCells[5], varId) * deltaX_Minus * deltaY_Plus * deltaZ_Minus
1834 + scalarField(interpolationCells[6], varId) * deltaX_Plus * deltaY_Minus * deltaZ_Minus
1835 + scalarField(interpolationCells[7], varId) * deltaX_Minus * deltaY_Minus * deltaZ_Minus)
1836 / (Delta * Delta * Delta);
1837 return phi;
1838 }
1839 }
1840 else {
1841 MFloat phi = 0;
1842 MFloat sumDist = 0;
1843 MFloat dist = -1;
1844 for(MInt i = 0; i < IPOW2(nDim); i++) {
1845 if(interpolationCells[i] < 0) continue;
1846 IF_CONSTEXPR(nDim == 2) {
1847 dist = std::sqrt(POW2(point[0] - coordinate(interpolationCells[i], 0))
1848 + POW2(point[1] - coordinate(interpolationCells[i], 1)));
1849 }
1850 else {
1851 dist = std::sqrt(POW2(point[0] - coordinate(interpolationCells[i], 0))
1852 + POW2(point[1] - coordinate(interpolationCells[i], 1))
1853 + POW2(point[2] - coordinate(interpolationCells[i], 2)));
1854 }
1855 sumDist += dist;
1856 phi += scalarField(interpolationCells[i], varId) * dist;
1857 }
1858 return phi / sumDist;
1859 }
1860}
1861
1867template <MInt nDim, class SolverType>
1868MFloat CartesianSolver<nDim, SolverType>::leastSquaresInterpolation(MInt* interpolationCells, MFloat* point, MInt varId,
1869 std::function<MFloat(MInt, MInt)> scalarField,
1870 std::function<MFloat(MInt, MInt)> coordinate) {
1871 std::array<MFloat, nDim + 1> b;
1872 std::array<std::array<MFloat, (nDim + 1)>, (nDim + 1)> A;
1873 std::array<MFloat, nDim + 1> adjA;
1874 b.fill(0.0);
1875 for(MInt i = 0; i < nDim + 1; i++) {
1876 for(MInt j = 0; j < nDim + 1; j++) {
1877 A[i][j] = 0.0;
1878 }
1879 }
1880 for(MInt n = 0; n < IPOW2(nDim); n++) {
1881 const MInt cellId = interpolationCells[n];
1882 b[0] += scalarField(cellId, varId);
1883 MFloat deltax = point[0] - coordinate(cellId, 0);
1884 MFloat deltay = point[1] - coordinate(cellId, 1);
1885 b[1] += scalarField(cellId, varId) * deltax;
1886 b[2] += scalarField(cellId, varId) * deltay;
1887
1888 A[0][0]++;
1889 A[1][1] += POW2(deltax);
1890 A[2][2] += POW2(deltay);
1891 A[1][0] += deltax;
1892 A[2][0] += deltay;
1893 A[2][1] += deltax * deltay;
1894 IF_CONSTEXPR(nDim == 3) {
1895 MFloat deltaz = point[2] - coordinate(cellId, 2);
1896 b[3] += scalarField(cellId, varId) * deltaz;
1897 A[3][0] += deltaz;
1898 A[3][1] += deltaz * deltax;
1899 A[3][2] += deltaz * deltay;
1900 A[3][3] += POW2(deltaz);
1901 }
1902 }
1903
1904 A[0][1] = A[1][0];
1905 A[0][2] = A[2][0];
1906 A[1][2] = A[2][1];
1907 IF_CONSTEXPR(nDim == 3) {
1908 A[0][3] = A[3][0];
1909 A[1][3] = A[3][1];
1910 A[2][3] = A[3][2];
1911 }
1912
1913 MFloat det = maia::math::determinant(A);
1914 ASSERT(fabs(det) > MFloatEps, "Poor least-squares stencil!");
1915
1916 maia::math::adjoint1stRow(A, adjA);
1917
1918 MFloat sum = 0;
1919 for(MInt i = 0; i < nDim + 1; i++) {
1920 sum += adjA[i] * b[i];
1921 }
1922
1923 return sum / det;
1924}
1925
1926
1930
1931template <MInt nDim, class SolverType>
1932void CartesianSolver<nDim, SolverType>::mapInterpolationCells(std::map<MInt, MInt>& cellMap) {
1933 TRACE();
1934
1935 ParallelIo srcGrid("src/grid.Netcdf", maia::parallel_io::PIO_READ, solver().mpiComm());
1936
1937 MInt noMinCells = solver().m_localMinCellsOffsets[1] - solver().m_localMinCellsOffsets[0] + 1; // see cartesiangrid
1938 srcGrid.setOffset(noMinCells, solver().m_localMinCellsOffsets[0]);
1939
1940 // memory
1941 MIntScratchSpace srcMinCellsGlobalIds(noMinCells, AT_, "srcMinCellsGlobalIds");
1942 MIntScratchSpace srcMinCellsNoOffsprings(noMinCells, AT_, "srcMinCellsNoOffsprings");
1943 // read
1944 srcGrid.readArray(srcMinCellsGlobalIds.getPointer(), "minCellsId");
1945 srcGrid.readArray(srcMinCellsNoOffsprings.getPointer(), "minCellsNoOffsprings");
1946
1947 // sum offsprings
1948 MInt srcFirstMinCellGlobalId = srcMinCellsGlobalIds(0);
1949 MInt srcNoCells = 0;
1950 for(MInt i = 0; i < noMinCells; i++) {
1951 srcNoCells += srcMinCellsNoOffsprings(i);
1952 srcMinCellsGlobalIds(i) -= srcFirstMinCellGlobalId;
1953 }
1954 // new offsets
1955 srcGrid.setOffset(srcNoCells, srcFirstMinCellGlobalId);
1956
1957
1958 // more memory
1959 // grid
1960 MFloatScratchSpace srcCoords(srcNoCells, nDim, AT_, "srcCoords");
1961 MIntScratchSpace srcNoChildIds(srcNoCells, AT_, "srcNoChildIds");
1962 MIntScratchSpace srcLevel(srcNoCells, AT_, "srcLevel");
1963 MIntScratchSpace srcChildIds(srcNoCells, IPOW2(nDim), AT_, "srcChildIds");
1964
1965 // load coords
1966 std::vector<MString> varNames;
1967 for(MInt j = 0; j < nDim; ++j) {
1968 std::stringstream ss;
1969 ss << "coordinates_" << j;
1970 varNames.push_back(ss.str());
1971 ss.clear();
1972 ss.str("");
1973 }
1974 srcGrid.readArray(&srcCoords(0, 0), varNames, nDim);
1975 varNames.clear();
1976 // load nmbr child
1977 srcGrid.readArray(srcNoChildIds.begin(), "noChildIds");
1978 // level
1979 srcGrid.readArray(srcLevel.begin(), "level_0");
1980 // load childIds
1981 for(MInt j = 0; j < IPOW2(nDim); ++j) {
1982 std::stringstream ss;
1983 ss << "childIds_" << j;
1984 varNames.push_back(ss.str());
1985 ss.clear();
1986 ss.str("");
1987 }
1988 srcGrid.readArray(&srcChildIds(0, 0), varNames, IPOW2(nDim));
1989 varNames.clear();
1990 // correct Ids
1991 for(MInt i = 0; i < srcNoCells; i++) {
1992 for(MInt j = 0; j < IPOW2(nDim); j++) {
1993 srcChildIds(i, j) -= srcFirstMinCellGlobalId;
1994 }
1995 }
1996
1997 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
1998 if(solver().c_noChildren(cellId)) {
1999 continue;
2000 }
2001 auto* coord = (MFloat*)(&(solver().a_coordinate(cellId, 0)));
2002 MInt srcId = -1;
2003 // find min Cell
2004 for(MInt srcMinCellId = 0; srcMinCellId < noMinCells; srcMinCellId++) {
2005 MInt srcMinCellGlobalId = srcMinCellsGlobalIds(srcMinCellId);
2006 MFloat* srcCoord = &srcCoords(srcMinCellGlobalId, 0);
2007 MFloat hdc = solver().c_cellLengthAtLevel(srcLevel(srcMinCellGlobalId) + 1);
2008 MInt insideCnt = 0;
2009 for(MInt dim = 0; dim < nDim; dim++) {
2010 if(coord[dim] < srcCoord[dim] + hdc && coord[dim] >= srcCoord[dim] - hdc) {
2011 insideCnt++;
2012 }
2013 }
2014
2015 if(insideCnt == nDim) {
2016 srcId = srcMinCellGlobalId;
2017 break;
2018 }
2019 if(srcMinCellId == (noMinCells - 1)) {
2020 mTerm(1, " should not happen, loc cell not in src region");
2021 }
2022 }
2023 // loop up (or down)
2024 while(srcNoChildIds(srcId) != 0) {
2025 for(MInt cc = 0; cc < IPOW2(nDim); cc++) {
2026 MInt tcid = srcChildIds(srcId, cc);
2027 if(tcid >= 0) {
2028 MFloat* srcCoord = &srcCoords(tcid, 0);
2029 MFloat hdc = solver().c_cellLengthAtLevel(srcLevel(tcid) + 1);
2030 // if the cell is inside the child break loop
2031 MInt insideCnt = 0;
2032 for(MInt dim = 0; dim < nDim; dim++) {
2033 if(coord[dim] < srcCoord[dim] + hdc && coord[dim] >= srcCoord[dim] - hdc) {
2034 insideCnt++;
2035 }
2036 }
2037 if(insideCnt == nDim) {
2038 srcId = tcid;
2039 break;
2040 }
2041 }
2042 if(cc == (IPOW2(nDim) - 1)) {
2043 mTerm(1, " should not happen, loc cell not in any child");
2044 }
2045 }
2046 }
2047
2048 cellMap[cellId] = srcId;
2049 }
2050}
2051
2052
2053template <MInt nDim, class SolverType>
2054void CartesianSolver<nDim, SolverType>::patchRefinement(std::vector<std::vector<MFloat>>& sensors,
2055 std::vector<std::bitset<64>>& sensorCellFlag,
2056 std::vector<MFloat>& sensorWeight, MInt sensorOffset,
2057 MInt sen) {
2058 if(solver().m_testPatch && globalTimeStep < solver().m_patchStartTimeStep
2059 && globalTimeStep > solver().m_patchStopTimeStep) {
2060 return;
2061 }
2062
2063 m_log << " - Sensor preparation for the patch sensor" << std::endl;
2064
2065 MFloat coordinates[nDim]{};
2066 MFloat patchDimensions[nDim * 2]{};
2067
2068 for(MInt lvl = 0; lvl < solver().m_patchRefinement->noLocalPatchRfnLvls(); lvl++) {
2069 const MInt level = lvl + solver().maxUniformRefinementLevel();
2070 for(MInt patch = 0; patch < solver().m_patchRefinement->noPatchesPerLevel(lvl); patch++) {
2071 MString patchStr = solver().m_patchRefinement->localRfnLevelMethods[lvl].substr(patch, 1);
2072 const MInt pos = solver().m_patchRefinement->localRfnLevelPropertiesOffset[lvl] + patch;
2073 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2074 if(solver().grid().tree().level(cellId) > level) continue;
2075 for(MInt i = 0; i < nDim; i++) {
2076 coordinates[i] = solver().grid().tree().coordinate(cellId, i);
2077 patchDimensions[i] = solver().m_patchRefinement->localRfnPatchProperties[pos][i];
2078 if(patchStr == "B") {
2079 patchDimensions[nDim + i] = solver().m_patchRefinement->localRfnPatchProperties[pos][nDim + i];
2080 }
2081 }
2082 MBool keepCell = false;
2083
2084 if(patchStr == "B") {
2085 if(maia::geom::isPointInsideBox<MFloat, nDim>(&coordinates[0], &patchDimensions[0])) {
2086 keepCell = true;
2087 }
2088 } else if(patchStr == "R") {
2089 if(maia::geom::isPointInsideSphere<nDim>(&coordinates[0], &patchDimensions[0],
2090 solver().m_patchRefinement->localRfnPatchProperties[pos][nDim])) {
2091 keepCell = true;
2092 }
2093 } else if(patchStr == "H") {
2094 // checks if point lies within a ring segment, requires center, r_min, r_max, phi_min, phi_max, axis
2095
2096 std::array<MFloat, nDim> center{};
2097 for(MInt dim = 0; dim < nDim; dim++) {
2098 center[dim] = solver().m_patchRefinement->localRfnPatchProperties[pos][dim];
2099 }
2100
2101 std::array<MFloat, 2> R{};
2102 R[0] = solver().m_patchRefinement->localRfnPatchProperties[pos][nDim];
2103 R[1] = solver().m_patchRefinement->localRfnPatchProperties[pos][nDim + 1];
2104
2105 std::array<MFloat, 2> phi{};
2106 phi[0] = solver().m_patchRefinement->localRfnPatchProperties[pos][nDim + 2];
2107 phi[1] = solver().m_patchRefinement->localRfnPatchProperties[pos][nDim + 3];
2108
2109 constexpr MInt axis = 2; // currently hardcoded
2110
2111 if(maia::geom::isPointInsideRingSegment<nDim>(&coordinates[0], &center[0], R.data(), phi.data(), axis)) {
2112 keepCell = true;
2113 }
2114 }
2115
2116 MBool refineCell = false;
2117 if(keepCell && solver().grid().tree().level(cellId) <= level) {
2118 refineCell = true;
2119 keepCell = false;
2120 }
2121
2122 if(refineCell) {
2123 const MInt gridCellId = grid().tree().solver2grid(cellId);
2124 sensors[sensorOffset + sen][gridCellId] = 1.0;
2125 sensorCellFlag[gridCellId][sensorOffset + sen] = true;
2126 } else if(keepCell) {
2127 const MInt gridCellId = grid().tree().solver2grid(cellId);
2128 sensors[sensorOffset + sen][gridCellId] = 0.0;
2129 sensorCellFlag[gridCellId][sensorOffset + sen] = false;
2130 }
2131 }
2132 }
2133 }
2134
2135 sensorWeight[sensorOffset + sen] = solver().m_sensorWeight[sen];
2136}
2137
2140template <MInt nDim, class SolverType>
2141void CartesianSolver<nDim, SolverType>::setHaloCellsOnInactiveRanks() {
2142 TRACE();
2143
2144 ASSERT(solver().m_cells.size() == 0, "");
2145 ASSERT(!solver().isActive(), "");
2146
2147 for(MInt cellId = 0; cellId < solver().grid().tree().size(); cellId++) {
2148 MInt solverCellId = solver().m_cells.size();
2149 solver().m_cells.append();
2150 solver().a_resetPropertiesSolver(solverCellId);
2151 solver().a_isHalo(solverCellId) = true;
2152
2153 ASSERT(solver().grid().tree().grid2solver(solver().grid().tree().solver2grid(solverCellId)) == solverCellId, "");
2154 }
2155}
2156
2160template <MInt nDim, class SolverType>
2161void CartesianSolver<nDim, SolverType>::readPatchProperties() {
2162 m_log << " + reading patch refinement properties for solver " << solver().solverId() << std::endl;
2163
2164 mAlloc(m_patchRefinement, 1, "m_patchRefinement", AT_);
2165
2166 MString localRfnLevelMethods = Context::getSolverProperty<MString>("localRfnLvlMethods", solver().solverId(), AT_);
2167
2168 // sanity check
2169 if(localRfnLevelMethods.front() == '-' || localRfnLevelMethods.back() == '-')
2170 mTerm(1, AT_, "ERROR: localRfnLvlMethods begins or ends with hyphen!");
2171
2172 if(localRfnLevelMethods.find(" ") != MString::npos) mTerm(1, AT_, "ERROR: localRfnLvlMethods contains space!");
2173
2174 MString::size_type prev_pos = MString::npos;
2175 MString::size_type nextMarker = 0;
2176
2177 // save substrings marked with '-' as the patch-refinement methods
2178 while(nextMarker != MString::npos) {
2179 prev_pos++;
2180 nextMarker = localRfnLevelMethods.find("-", prev_pos);
2181 MString substring(localRfnLevelMethods.substr(prev_pos, nextMarker - prev_pos));
2182 m_patchRefinement->localRfnLevelMethods.push_back(substring);
2183 prev_pos = nextMarker;
2184 }
2185
2186 // how many patches do we have
2187 MInt numPatches = 0;
2188 for(MInt l = 0; l < m_patchRefinement->noLocalPatchRfnLvls(); l++) {
2189 numPatches += m_patchRefinement->noPatchesPerLevel(l);
2190 }
2191
2192 m_log << " - " << numPatches << " Total-Patches. Ranging from level " << solver().maxUniformRefinementLevel()
2193 << " - " << solver().maxUniformRefinementLevel() + m_patchRefinement->noLocalPatchRfnLvls() << std::endl;
2194
2195 // allocate vector containing the offset to each rfnLvl in the patchProperties matrix
2196 mAlloc(m_patchRefinement->localRfnLevelPropertiesOffset, m_patchRefinement->noLocalPatchRfnLvls() + 1,
2197 "localRfnLevelPropertiesOffset", 0, AT_);
2198
2199 // allocate vector containing the number of properties for each patch
2200 mAlloc(m_patchRefinement->noLocalRfnPatchProperties, numPatches, "noLocalRfnPatchProperties", 0, AT_);
2201
2202 // calculate number of properties required for each patchRfnLvl and set offsets
2203 m_patchRefinement->localRfnLevelPropertiesOffset[0] = 0;
2204 MInt count_req = 0;
2205 MInt s = 0;
2206 for(MInt i = 0; i < m_patchRefinement->noLocalPatchRfnLvls(); i++) {
2207 const MString lvlStr = m_patchRefinement->localRfnLevelMethods[i];
2208 for(MString::size_type j = 0; j < lvlStr.size(); j++) {
2209 const MString patchStr = lvlStr.substr(j, 1);
2210 if(patchStr == "B") {
2211 m_patchRefinement->noLocalRfnPatchProperties[s] = 2 * nDim;
2212 } else if(patchStr == "R") {
2213 m_patchRefinement->noLocalRfnPatchProperties[s] = nDim + 1;
2214 } else if(patchStr == "H") {
2215 m_patchRefinement->noLocalRfnPatchProperties[s] = nDim + 4;
2216 } else {
2217 mTerm(1, AT_, "Not yet implemented, please do so...");
2218 }
2219 count_req += m_patchRefinement->noLocalRfnPatchProperties[s];
2220 s++;
2221 }
2222 m_patchRefinement->localRfnLevelPropertiesOffset[i + 1] = s;
2223 }
2224
2225 // allocate matrix containing the properties for each patch
2226 mAlloc(m_patchRefinement->localRfnPatchProperties, numPatches, m_patchRefinement->noLocalRfnPatchProperties,
2227 "localRfnPatchProperties", AT_);
2228
2229 MInt count_prov = Context::propertyLength("localRfnLevelProperties", solver().solverId());
2230
2231 // sanity check
2232 if(count_prov < count_req)
2233 mTerm(1, AT_, "ERROR: number of localRfnLevelProperties does not match the requested value!");
2234
2235 // load properties of patch refinement methods
2236 MInt lvl = 0;
2237 MInt j = 0;
2238 for(MInt i = 0; i < count_req; i++) {
2239 m_patchRefinement->localRfnPatchProperties[lvl][j] =
2240 Context::getSolverProperty<MFloat>("localRfnLevelProperties", solver().solverId(), AT_, i);
2241 j++;
2242 if(j == m_patchRefinement->noLocalRfnPatchProperties[lvl]) {
2243 j = 0;
2244 lvl++;
2245 }
2246 }
2247
2248 // properties to test patch build-up and reduction
2249 m_testPatch = Context::getSolverProperty<MBool>("testPatchRefinement", solver().solverId(), AT_, &m_testPatch);
2250
2251 if(m_testPatch) {
2252 m_patchStartTimeStep =
2253 Context::getSolverProperty<MInt>("patchStart", solver().solverId(), AT_, &m_patchStartTimeStep);
2254 m_patchStopTimeStep = Context::getSolverProperty<MInt>("patchStop", solver().solverId(), AT_, &m_patchStopTimeStep);
2255 }
2256}
2257
2258template <MInt nDim, class SolverType>
2259void CartesianSolver<nDim, SolverType>::initAdaptation() {
2272 m_adaptation = false;
2273 m_adaptation = Context::getSolverProperty<MBool>("adaptation", solver().solverId(), AT_, &m_adaptation);
2274
2275 m_adapts = true;
2276 m_adapts = Context::getSolverProperty<MBool>("solverAdapts", solver().solverId(), AT_, &m_adapts);
2277
2278 m_resTriggeredAdapt = false;
2279 m_resTriggeredAdapt =
2280 Context::getSolverProperty<MBool>("resTriggeredAdaptation", solver().solverId(), AT_, &m_resTriggeredAdapt);
2281
2282 m_noSensors = 0;
2283 m_noInitialSensors = 0;
2284 m_rfnBandWidth = nullptr;
2285 m_sensorInterface = false;
2286 m_sensorParticle = false;
2287
2288 if(!m_adaptation) return;
2289
2290 const MInt maxLevel = Context::getSolverProperty<MInt>("maxRfnmntLvl", solver().solverId(), AT_);
2291
2292 m_log << "##################################################################" << std::endl;
2293 m_log << "##################### Adaptation is active #######################" << std::endl;
2294 m_log << "##################################################################" << std::endl << std::endl;
2295
2296 solver().m_singleAdaptation = false;
2297 solver().m_singleAdaptation =
2298 Context::getSolverProperty<MBool>("singleAdaptationLoop", solver().solverId(), AT_, &solver().m_singleAdaptation);
2299
2300 ASSERT(Context::propertyLength("sensorType") == Context::propertyLength("sensorWeight"),
2301 "Lenght of sensorType and sensorWeight not equal.");
2302
2303 if(Context::propertyExists("sensorType", solver().solverId())) {
2304 m_noSensors = Context::propertyLength("sensorType", solver().solverId());
2305 }
2306
2307 if(m_noSensors > 0 && !Context::propertyExists("sensorWeight", solver().solverId())) {
2308 TERMM(-1, "sensorWeight property not set!");
2309 }
2310
2311 ASSERT(m_noSensors == Context::propertyLength("sensorWeight", solver().solverId()),
2312 "Length of sensorType " + std::to_string(m_noSensors) + " and sensorWeight "
2313 + std::to_string(Context::propertyLength("sensorWeight", solver().solverId())) + " not equal.");
2314
2315 m_saveSensorData = false;
2316 if(Context::propertyExists("saveSensorData", solver().solverId())) {
2317 m_saveSensorData = Context::getSolverProperty<MBool>("saveSensorData", solver().solverId(), AT_);
2318 }
2319
2320 m_log << " * Sensors for adaptive mesh refinement for solver " << solver().solverId() << " are active:" << std::endl;
2321
2322 MInt der = 0;
2323 for(MInt s = 0; s < m_noSensors; s++) {
2324 const MString sensor = Context::getSolverProperty<MString>("sensorType", solver().solverId(), AT_, s);
2325 const MFloat weight = Context::getSolverProperty<MFloat>("sensorWeight", solver().solverId(), AT_, s);
2326 const MInt level = Context::getSolverProperty<MInt>("maxSensorLevel", solver().solverId(), AT_, &maxLevel, s);
2327
2328 m_sensorType.push_back(sensor);
2329 m_sensorWeight.push_back(weight);
2330 m_maxSensorRefinementLevel.push_back(level);
2331
2332 m_log << " - sensor: " << sensor << std::endl;
2333 m_log << " - Weight: " << weight << std::endl;
2334 m_log << " - maxLevel: " << level << std::endl;
2335
2336 switch(string2enum(sensor)) {
2337 case DERIVATIVE: {
2338 ASSERT(Context::propertyLength("sensorType", solver().solverId()) >= der,
2339 "sensorDerivativeVariables not stated for this sensor");
2340 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorDerivative);
2341 MInt derivative = Context::getSolverProperty<MInt>("sensorDerivativeVariables", solver().solverId(), AT_, der);
2342 m_sensorDerivativeVariables.push_back(derivative);
2343 der++;
2344 break;
2345 }
2346 case DIVERGENCE:
2347 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorDivergence);
2348 break;
2349 case TOTALPRESSURE:
2350 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorTotalPressure);
2351 break;
2352 case ENTROPY_GRADIENT:
2353 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorEntropyGrad);
2354 break;
2355 case ENTROPY_QUOTIENT:
2356 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorEntropyQuot);
2357 break;
2358 case MEANSTRESS:
2359 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorMeanStress);
2360 break;
2361 case VORTICITY:
2362 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorVorticity);
2363 break;
2364 case INTERFACE:
2365 ASSERT(s == 0, "Interface sensor has to be the first sensor");
2366 m_sensorInterface = true;
2367 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorInterface);
2368 ASSERT((MInt)m_sensorWeight[s] == -1, "Must be discrete sensor!");
2369 m_noInitialSensors++;
2370 break;
2371 case PARTICLE:
2372 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorParticle);
2373 m_sensorParticle = true;
2374 ASSERT((MInt)m_sensorWeight[s] == -1, "Must be discrete sensor!");
2375 m_noInitialSensors++;
2376 break;
2377 case SPECIES:
2378 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorSpecies);
2379 break;
2380 case PATCH:
2381 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorPatch);
2382 readPatchProperties();
2383 // Patch sensor must be second or first, as its also not solution dependand!
2384 ASSERT(s == 0 || (m_sensorInterface && s == 1), "");
2385 ASSERT((MInt)m_sensorWeight[s] == -1, "Must be discrete sensor!");
2386 m_noInitialSensors++;
2387 break;
2388 case CUTOFF:
2389 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorCutOff);
2390 m_noInitialSensors++;
2391 break;
2392 case SMOOTH:
2393 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorSmooth);
2394 m_noInitialSensors++;
2395 m_noSmoothingLayers = Context::getSolverProperty<MInt>("smoothingLayers", solver().solverId(), AT_);
2396 break;
2397 case BAND:
2398 m_sensorFnPtr.push_back(&CartesianSolver<nDim, SolverType>::sensorBand);
2399 m_sensorBandAdditionalLayers = 10;
2400 m_sensorBandAdditionalLayers = Context::getSolverProperty<MInt>("sensorBandAdditionalLayers", m_solverId, AT_,
2401 &m_sensorBandAdditionalLayers);
2402 ASSERT(s == (m_noSensors - 1), "Please ensure this is the last sensor to be called.");
2403 break;
2404 default:
2405 m_sensorDerivativeVariables.push_back(-1);
2406 mTerm(1, AT_, "Sensor not found! Check property sensorType!");
2407 }
2408 }
2409
2410 // todo labels:toenhance @Moritz: Switch to sensors and uncomment the Assert below!
2411 // ASSERT(m_noSensors > 0, "Adaptation without sensors is pointless...");
2412
2413 m_log << std::endl;
2414}
2415
2422template <MInt nDim, class SolverType>
2423void CartesianSolver<nDim, SolverType>::sensorBand(std::vector<std::vector<MFloat>>& sensors,
2424 std::vector<std::bitset<64>>& sensorCellFlag,
2425 std::vector<MFloat>& sensorWeight,
2426 MInt sensorOffset,
2427 MInt sen) {
2428 TRACE();
2429
2430 if(this->m_adaptationStep < (maxRefinementLevel() - minLevel())) return;
2431
2432 MIntScratchSpace markedCells(solver().a_noCells(), AT_, "markedCells");
2433 MIntScratchSpace bandWidth(this->m_maxSensorRefinementLevel[sen], AT_, "bandWidth");
2434 markedCells.fill(0);
2435 bandWidth.fill(0);
2436
2437 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2438 if(c_level(cellId) != m_maxSensorRefinementLevel[sen]) continue;
2439 markedCells[cellId] = 1;
2440 MInt parentId = c_parentId(cellId);
2441 while(parentId > -1 && parentId < solver().c_noCells()) {
2442 markedCells[parentId] = 1;
2443 parentId = c_parentId(parentId);
2444 }
2445 }
2446
2447 bandWidth[this->m_maxSensorRefinementLevel[sen] - 1] =
2448 m_bandWidth[this->m_maxSensorRefinementLevel[sen] - 1] + this->m_sensorBandAdditionalLayers;
2449
2450 for(MInt i = this->m_maxSensorRefinementLevel[sen] - 2; i >= 0; i--) {
2451 bandWidth[i] = (bandWidth[i + 1] / 2) + this->m_sensorBandAdditionalLayers;
2452 }
2453
2454 for(MInt level = minLevel(); level < this->m_maxSensorRefinementLevel[sen]; level++) {
2455 this->markSurrndCells(markedCells, bandWidth[level], level, true);
2456 }
2457
2458 MInt refinedCells = 0;
2459 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2460 if(markedCells(cellId) > 0) {
2461 const MInt gridCellId = grid().tree().solver2grid(cellId);
2462 refinedCells++;
2463 sensors[sensorOffset + sen][gridCellId] = 1;
2464 sensorCellFlag[gridCellId][sensorOffset + sen] = true;
2465 MInt parent = c_parentId(cellId);
2466 if(parent > -1 && parent < solver().a_noCells()) {
2467 MInt parentGridCellId = grid().tree().solver2grid(parent);
2468 if(parentGridCellId > -1 && parentGridCellId < grid().raw().m_noInternalCells) {
2469 sensors[sensorOffset + sen][parentGridCellId] = 1;
2470 sensorCellFlag[parentGridCellId][sensorOffset + sen] = true;
2471 }
2472 }
2473 }
2474 }
2475 sensorWeight[sensorOffset + sen] = this->m_sensorWeight[sen];
2476}
2477
2491template <MInt nDim, class SolverType>
2492void CartesianSolver<nDim, SolverType>::sensorLimit(std::vector<std::vector<MFloat>>& sensors,
2493 std::vector<std::bitset<64>>& sensorCellFlag,
2494 std::vector<MFloat>& sensorWeight, MInt sensorOffset, MInt sen,
2495 std::function<MFloat(MInt)> value, const MFloat limit,
2496 const MInt* bandWidth, const MBool refineDiagonals,
2497 const MBool allowCoarsening) {
2498 TRACE();
2499
2500 std::ignore = sensorWeight[sensorOffset];
2501
2502 MIntScratchSpace markedCells(solver().a_noCells(), AT_, "markedCells");
2503 markedCells.fill(0);
2504
2505 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2506 if(value(cellId) >= limit) {
2507 markedCells[cellId] = 1;
2508 MInt parentId = c_parentId(cellId);
2509 while(parentId > -1 && parentId < solver().a_noCells()) {
2510 markedCells[parentId] = 1;
2511 parentId = c_parentId(parentId);
2512 }
2513 }
2514 }
2515
2516 exchangeData(&markedCells[0], 1);
2517
2518 for(MInt level = minLevel(); level < m_maxSensorRefinementLevel[sen]; level++) {
2519 this->markSurrndCells(markedCells, bandWidth[level], level, refineDiagonals);
2520 }
2521
2522 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2523 if(markedCells(cellId) == 0) { // coarsen
2524 if(solver().grid().tree().level(cellId) == minLevel()) continue;
2525 if(!solver().c_isLeafCell(cellId)) continue;
2526 if(markedCells[solver().c_parentId(cellId)]) continue;
2527 if(!allowCoarsening) continue;
2528
2529 const MInt gridCellId = grid().tree().solver2grid(cellId);
2530 sensors[sensorOffset + sen][gridCellId] = -1.0;
2531 sensorCellFlag[gridCellId][sensorOffset + sen] = true;
2532 //}
2533 } else {
2534 if(solver().c_level(cellId) < m_maxSensorRefinementLevel[sen]) { // refine cell
2535 // if(solver().c_noChildren(cellId) > 0) continue;
2536
2537 const MInt gridCellId = grid().tree().solver2grid(cellId);
2538 sensors[sensorOffset + sen][gridCellId] = 1.0;
2539 sensorCellFlag[gridCellId][sensorOffset + sen] = true;
2540
2541 MInt parent = solver().c_parentId(cellId);
2542 if(parent > -1 && parent < solver().c_noCells()) {
2543 MInt parentGridCellId = grid().tree().solver2grid(parent);
2544 if(parentGridCellId > -1 && parentGridCellId < solver().grid().raw().m_noInternalCells) {
2545 sensors[sensorOffset + sen][parentGridCellId] = 1.0;
2546 sensorCellFlag[parentGridCellId][sensorOffset + sen] = true;
2547 }
2548 }
2549 }
2550 }
2551 }
2552}
2553
2565template <MInt nDim, class SolverType>
2566void CartesianSolver<nDim, SolverType>::saveSensorData(const std::vector<std::vector<MFloat>>& sensors,
2567 const MInt& level, const MString& gridFileName,
2568 const MInt* const recalcIdsTree) {
2569 using namespace maia::parallel_io;
2570 ASSERT(nDim == 2 || nDim == 3, "wrong number of dimensions!");
2571
2572 // Update recalc ids
2573 std::vector<MInt> recalcCellIdsSolver(0);
2574 MInt noCells;
2575 MInt noInternalCellIds;
2576 std::vector<MInt> reorderedCellIds(0);
2577 this->calcRecalcCellIdsSolver(recalcIdsTree, noCells, noInternalCellIds, recalcCellIdsSolver, reorderedCellIds);
2578
2579 MLong totalNoInternalCells = -1;
2580 const MLong longNoInternalCells = noInternalCellIds;
2581 MPI_Allreduce(&longNoInternalCells, &totalNoInternalCells, 1, MPI_LONG, MPI_SUM, mpiComm(), AT_,
2582 "longNoInternalCells", "totalNoInternalCells");
2583
2584 // Open file
2585 std::stringstream fileName;
2586 fileName << outputDir() << "sensorData_" << solverId() << "_" << std::to_string(level) << "_" << globalTimeStep
2587 << ParallelIo::fileExt();
2588 ParallelIo parallelIo(fileName.str(), PIO_REPLACE, mpiComm());
2589 // Write header information
2590 parallelIo.defineScalar(PIO_INT, "noCells");
2591 parallelIo.setAttribute(solverId(), "solverId");
2592 for(MUint counter = 0; counter < sensors.size(); counter++) {
2593 const MString arrayName = "variables" + std::to_string(counter);
2594 parallelIo.defineArray(PIO_FLOAT, arrayName, totalNoInternalCells);
2595 const MString varName = "sensor_" + m_sensorType[counter];
2596 parallelIo.setAttribute(varName, "name", arrayName);
2597 }
2598 // Write global attributes
2599 parallelIo.setAttribute(gridFileName, "gridFile", "");
2600 parallelIo.setAttribute(solver().time(), "time");
2601 parallelIo.setAttribute(globalTimeStep, "globalTimeStep");
2602 // Write scalars
2603 parallelIo.writeScalar(totalNoInternalCells, "noCells");
2604 // Set file offsets for field data
2605 ParallelIo::size_type offset = 0;
2606 MPI_Exscan(&longNoInternalCells, &offset, 1, MPI_LONG, MPI_SUM, mpiComm(), AT_, "longNoInternalCells", "offset");
2607 const ParallelIo::size_type noInternalCells = longNoInternalCells;
2608 parallelIo.setOffset(noInternalCells, offset);
2609 // Write sensor values
2610 const MString name = "variables";
2611 MInt suffix = 0;
2612 for(auto&& sensor : sensors) {
2613 ScratchSpace<MFloat> buffer(noInternalCells, AT_, "buffer");
2614 for(MInt cellId = 0; cellId < noInternalCells; cellId++) {
2615 const MInt cellIdRecalc = (recalcIdsTree == nullptr) ? cellId : recalcCellIdsSolver[cellId];
2616 const MInt cellIdGrid = grid().tree().solver2grid(cellIdRecalc);
2617 buffer[cellId] = sensor[cellIdGrid];
2618 }
2619 parallelIo.writeArray(buffer.data(), name + std::to_string(suffix));
2620 suffix++;
2621 }
2622}
2623
2624
2631template <MInt nDim, class SolverType>
2632void CartesianSolver<nDim, SolverType>::sensorSmooth(std::vector<std::vector<MFloat>>& sensors,
2633 std::vector<std::bitset<64>>& sensorCellFlag,
2634 std::vector<MFloat>& sensorWeight,
2635 MInt sensorOffset,
2636 MInt sen) {
2637 TRACE();
2638
2639 static constexpr MInt noDirs = 2 * nDim;
2640 ASSERT(m_noSmoothingLayers > 0, "No-smoothing layers not specified!");
2641 static constexpr std::array<MInt, 4> diagDirs{3, 2, 0, 1};
2642
2643 MFloatScratchSpace markedCells(solver().c_noCells(), AT_, "markedCells");
2644 MIntScratchSpace markedLevel(solver().c_noCells(), AT_, "markedCells");
2645
2646 // 1) loop over all refinement levels and set sensor for smooth level transition
2647 // top-down:
2648 for(MInt lvl = m_maxSensorRefinementLevel[sen]; lvl > maxUniformRefinementLevel(); lvl--) {
2649 markedCells.fill(-1.0);
2650
2651 for(MInt cellId = 0; cellId < c_noCells(); cellId++) {
2652 markedLevel[cellId] = c_level(cellId);
2653 }
2654
2655 // 2) mark cells at matching level jump
2656 // const MInt parentLvl = lvl - 1;
2657 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2658 if(solver().c_level(cellId) != lvl) continue;
2659 MBool atLvlJump = false;
2660 for(MInt dir = 0; dir < noDirs; dir++) {
2661 if(solver().a_hasNeighbor(cellId, dir)) continue;
2662 // note: cell must have a valid parent as its level is above the minLevel!
2663 if(solver().a_hasNeighbor(c_parentId(cellId), dir)) {
2664 atLvlJump = true;
2665 }
2666 }
2667 if(atLvlJump) {
2668 markedCells(c_parentId(cellId)) = 1.0;
2669 // markedLevel(c_parentId(cellId)) = parentLvl;
2670 }
2671 }
2672
2673 exchangeData(&markedCells[0], 1);
2674
2675 // 3) mark band around the level jump on the parent level to ensure smooth transition
2676 // note: similar to markSurrndCells cells, but different:
2677 // in this case the refinement is done top-down instead of bottom-up level-wise!
2678 // thus marking across level-jumps must be considered
2679 for(MInt loopMarker = 1; loopMarker <= m_noSmoothingLayers + 1; loopMarker++) {
2680 for(MInt cellId = 0; cellId < solver().c_noCells(); cellId++) {
2681 if((MInt)floor(markedCells(cellId)) != loopMarker) continue;
2682 const MInt orgLvl = markedLevel[cellId];
2683
2684 MFloat cellDist = -1;
2685 if(solver().c_level(cellId) == orgLvl) {
2686 cellDist = loopMarker + 1;
2687 } else {
2688 // note: only mark half as many cells when going down-wards
2689 ASSERT(solver().c_level(cellId) < orgLvl, "");
2690 cellDist = loopMarker + 2 * (orgLvl - solver().c_level(cellId));
2691 }
2692
2693 // direct neighbors +x-x+y-y+z-z
2694 for(MInt mainDir = 0; mainDir < noDirs; mainDir++) {
2695 MInt nghbrId = solver().c_neighborId(cellId, mainDir);
2696 if(nghbrId < 0) {
2697 if(c_parentId(cellId) > -1 && solver().a_hasNeighbor(c_parentId(cellId), mainDir)) {
2698 nghbrId = solver().c_neighborId(c_parentId(cellId), mainDir);
2699 } else {
2700 continue;
2701 }
2702 }
2703 const MFloat nghbDist = markedCells(nghbrId);
2704 if(nghbDist < cellDist) {
2705 markedCells(nghbrId) = cellDist;
2706 markedLevel(nghbrId) = orgLvl;
2707 }
2708 if(mainDir < 4 && loopMarker < m_noSmoothingLayers) {
2709 // add diagonal neighbors in x-y plane
2710 MInt nghbrId2 = solver().c_neighborId(nghbrId, diagDirs.at(mainDir));
2711 if(nghbrId2 < 0) {
2712 if(c_parentId(nghbrId) > -1 && solver().a_hasNeighbor(c_parentId(nghbrId), diagDirs.at(mainDir))) {
2713 nghbrId2 = solver().c_neighborId(c_parentId(nghbrId), diagDirs.at(mainDir));
2714 }
2715 }
2716 if(nghbrId2 > -1) {
2717 const MFloat nghbDist2 = markedCells(nghbrId2);
2718 if(nghbDist2 < cellDist) {
2719 markedCells(nghbrId2) = cellDist;
2720 markedLevel(nghbrId2) = orgLvl;
2721 }
2722 }
2723 IF_CONSTEXPR(nDim == 3) {
2724 // add both diagonal neighbors for the direct neighbors
2725 for(MInt dirZ = 4; dirZ < 6; dirZ++) {
2726 MInt nghbrIdZ = solver().c_neighborId(nghbrId, dirZ);
2727 if(nghbrIdZ < 0) {
2728 if(c_parentId(nghbrId) > -1 && solver().a_hasNeighbor(c_parentId(nghbrId), dirZ)) {
2729 nghbrIdZ = solver().c_neighborId(c_parentId(nghbrId), dirZ);
2730 }
2731 }
2732 if(nghbrIdZ > -1) {
2733 const MFloat nghbDistZ = markedCells(nghbrIdZ);
2734 if(nghbDistZ < cellDist) {
2735 markedCells(nghbrIdZ) = cellDist;
2736 markedLevel(nghbrIdZ) = orgLvl;
2737 }
2738 }
2739 }
2740 }
2741 }
2742 }
2743 }
2744 exchangeData(&markedCells[0], 1);
2745 }
2746
2747 // 4) marks cells for refinement
2748 for(MInt cellId = 0; cellId < solver().noInternalCells(); cellId++) {
2749 if(markedCells(cellId) < 0) continue;
2750 if(solver().c_level(cellId) < m_maxSensorRefinementLevel[sen]) {
2751 const MInt gridCellId = grid().tree().solver2grid(cellId);
2752
2753 // refine
2754 if(markedLevel[cellId] > solver().c_level(cellId)) {
2755 sensors[sensorOffset + sen][gridCellId] = 1.0;
2756 sensorCellFlag[gridCellId][sensorOffset + sen] = true;
2757
2758 MInt parent = solver().c_parentId(cellId);
2759 if(parent > -1 && parent < solver().c_noCells()) {
2760 MInt parentGridCellId = grid().tree().solver2grid(parent);
2761 if(parentGridCellId > -1 && parentGridCellId < solver().grid().raw().m_noInternalCells) {
2762 sensors[sensorOffset + sen][parentGridCellId] = 1.0;
2763 sensorCellFlag[parentGridCellId][sensorOffset + sen] = true;
2764 }
2765 }
2766 }
2767 }
2768 }
2769 }
2770 sensorWeight[sensorOffset + sen] = m_sensorWeight[sen];
2771}
2772
2778template <MInt nDim, class SolverType>
2779void CartesianSolver<nDim, SolverType>::reOrderCellIds(std::vector<MInt>& reOrderedCells) {
2780 TRACE();
2781
2782 if(grid().newMinLevel() < 0) {
2783 return;
2784 }
2785
2786 // see cartesiangrid storeMinLevelCells
2787 std::map<MLong, MInt> minLevelCells;
2788
2789 for(MInt cellId = 0; cellId < solver().a_noCells(); cellId++) {
2790 // allow partition-levelAnchestor minLevel halo-cell!
2791 if(grid().tree().hasProperty(cellId, Cell::IsHalo)
2792 && !grid().raw().a_hasProperty(grid().tree().solver2grid(cellId), Cell::IsPartLvlAncestor)) {
2793 continue;
2794 }
2795 if(c_level(cellId) == grid().newMinLevel()) {
2796 const MLong hilbertId = grid().generateHilbertIndex(cellId, grid().newMinLevel());
2797 if(minLevelCells.count(hilbertId) > 0) {
2798 mTerm(1, AT_, "duplicate hilbertId.");
2799 }
2800 minLevelCells.insert(std::make_pair(hilbertId, cellId));
2801 }
2802 }
2803 // the new order of mincell levels is stored as the block ids of the new grid block
2804 const MInt size = minLevelCells.size();
2805 std::vector<MInt> newMinCellOrder;
2806 newMinCellOrder.reserve(size);
2807 for(auto& minLevelCell : minLevelCells) {
2808 newMinCellOrder.push_back(minLevelCell.second);
2809 }
2810
2811 // now generate the full new cell order for writeout
2812 reOrderedCells.reserve(size);
2813 for(auto it = newMinCellOrder.begin(); it != newMinCellOrder.end(); it++) {
2814 reOrderedCells.push_back(*it);
2815 addChildren(reOrderedCells, *it);
2816 }
2817}
2818
2819
2825template <MInt nDim, class SolverType>
2826void CartesianSolver<nDim, SolverType>::addChildren(std::vector<MInt>& reOrderedIds, const MInt parentId) {
2827 for(MInt child = 0; child < grid().m_maxNoChilds; child++) {
2828 const MInt childId = grid().tree().child(parentId, child);
2829 // skip if there is no child
2830 if(childId < 0) continue;
2831 // add the child to the list
2832 reOrderedIds.push_back(childId);
2833 // if the cell has even more children repeat this for those children
2834 if(grid().tree().noChildren(childId) > 0) {
2835 addChildren(reOrderedIds, childId);
2836 }
2837 }
2838}
2844template <MInt nDim, class SolverType>
2845void CartesianSolver<nDim, SolverType>::recomputeGlobalIds(std::vector<MInt>& reOrderedCells,
2846 std::vector<MLong>& newGlobalIds) {
2847 std::vector<MLong> domainOffsets;
2848
2849 MInt countInternal = 0;
2850 for(MUint i = 0; i < reOrderedCells.size(); i++) {
2851 if(solver().a_isHalo(reOrderedCells[i])) continue;
2852 countInternal++;
2853 }
2854
2855 MIntScratchSpace newInternalCells(solver().grid().noDomains(), AT_, "newInternalCells");
2856 newInternalCells[solver().domainId()] = countInternal;
2857 MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &newInternalCells[0], 1, type_traits<MInt>::mpiType(),
2858 solver().mpiComm(), AT_, "MPI_IN_PLACE", "newInternalCells");
2859
2860 domainOffsets.assign(solver().grid().noDomains() + 1, -1);
2861 domainOffsets[0] = solver().grid().bitOffset();
2862 for(MInt d = 1; d < solver().grid().noDomains() + 1; d++) {
2863 domainOffsets[d] = domainOffsets[d - 1] + (MLong)newInternalCells[d - 1];
2864 }
2865
2866 newGlobalIds.clear();
2867 for(MUint i = 0; i < reOrderedCells.size(); i++) {
2868 if(solver().a_isHalo(reOrderedCells[i])) continue;
2869 newGlobalIds.push_back(domainOffsets[solver().domainId()] + i);
2870 }
2871}
2877template <MInt nDim, class SolverType>
2878template <typename T>
2879void CartesianSolver<nDim, SolverType>::collectVariables(T* variablesIn, ScratchSpace<T>& variablesOut,
2880 const std::vector<MString>& variablesNameIn,
2881 std::vector<MString>& variablesNameOut, const MInt noVars,
2882 const MInt noCells, const MBool reverseOrder) {
2883 TRACE();
2884
2885 const MInt variablesOffset = variablesNameOut.size();
2886 const MInt noTotalVars = noVars + variablesOffset;
2887 for(MInt j = variablesOffset; j < noTotalVars; j++) {
2888 const MInt k = j - variablesOffset;
2889 for(MInt i = 0; i < noCells; i++) {
2890 if(reverseOrder) {
2891 variablesOut.p[j * noCells + i] = variablesIn[k * noCells + i];
2892 } else {
2893 variablesOut.p[j * noCells + i] = variablesIn[i * noVars + k];
2894 }
2895 }
2896 variablesNameOut.push_back(variablesNameIn[k]);
2897 }
2898}
2899
2905template <MInt nDim, class SolverType>
2906template <typename T>
2907void CartesianSolver<nDim, SolverType>::collectVariables(T** variablesIn, ScratchSpace<T>& variablesOut,
2908 const std::vector<MString>& variablesNameIn,
2909 std::vector<MString>& variablesNameOut, const MInt noVars,
2910 const MInt noCells) {
2911 TRACE();
2912
2913 const MInt variablesOffset = variablesNameOut.size();
2914 const MInt noTotalVars = noVars + variablesOffset;
2915 for(MInt j = variablesOffset; j < noTotalVars; j++) {
2916 const MInt k = j - variablesOffset;
2917 for(MInt i = 0; i < noCells; i++) {
2918 variablesOut.p[j * noCells + i] = variablesIn[i][k];
2919 }
2920 variablesNameOut.push_back(variablesNameIn[k]);
2921 }
2922}
2923
2928template <MInt nDim, class SolverType>
2929template <typename T>
2930void CartesianSolver<nDim, SolverType>::collectParameters(T parametersIn, ScratchSpace<T>& parametersOut,
2931 const MChar* parametersNameIn,
2932 std::vector<MString>& parametersNameOut) {
2933 TRACE();
2934
2935 const MInt parsOffset = parametersNameOut.size();
2936 parametersOut.p[parsOffset] = parametersIn;
2937 parametersNameOut.push_back(parametersNameIn);
2938}
2939
2948template <MInt nDim, class SolverType>
2949void CartesianSolver<nDim, SolverType>::calcRecalcCellIdsSolver(const MInt* const recalcIdsTree, MInt& noCells,
2950 MInt& noInternalCellIds,
2951 std::vector<MInt>& recalcCellIdsSolver,
2952 std::vector<MInt>& reorderedCellIds) {
2953 MBool needRecalcCellIdsSolver = (recalcIdsTree != nullptr);
2954 noCells = noInternalCells();
2955 noInternalCellIds = noInternalCells();
2956
2957 if(grid().newMinLevel() > 0) {
2958 if(domainId() == 0) {
2959 std::cerr << "Increasing minLevel for solver " << m_solverId << std::endl;
2960 }
2961 this->reOrderCellIds(reorderedCellIds);
2962 MInt countInternal = 0;
2963 for(MUint i = 0; i < reorderedCellIds.size(); i++) {
2964 if(grid().tree().hasProperty(reorderedCellIds[i], Cell::IsHalo)) {
2965 continue;
2966 }
2967 countInternal++;
2968 }
2969 needRecalcCellIdsSolver = false;
2970 noCells = reorderedCellIds.size();
2971 noInternalCellIds = countInternal;
2972 }
2973
2974 if(needRecalcCellIdsSolver) {
2975 recalcCellIdsSolver.resize(grid().tree().size());
2976 // for multisolver recalc size needs to be raw grid size
2977 MInt recalcCounter = 0;
2978 for(MInt cellIdGrid = 0; cellIdGrid < grid().raw().noInternalCells(); cellIdGrid++) {
2979 if(grid().raw().a_hasProperty(cellIdGrid, Cell::IsHalo)) {
2980 continue;
2981 }
2982 const MInt cellIdSolver = grid().tree().grid2solver(recalcIdsTree[cellIdGrid]);
2983 if(cellIdSolver > -1) {
2984 recalcCellIdsSolver[recalcCounter] = cellIdSolver;
2985 recalcCounter++;
2986 ASSERT(grid().solverFlag(recalcIdsTree[cellIdGrid], m_solverId), "");
2987 }
2988 }
2989 ASSERT(recalcCounter == grid().noInternalCells(), "recalc ids size is wrong");
2990 }
2991
2992#ifdef LS_DEBUG
2993 MInt internalGCells = 0;
2994 for(MInt k = 0; k < a_noCells(); ++k) {
2995 if(grid().tree().hasProperty(k, Cell::IsHalo)) continue;
2996 ++internalGCells;
2997 }
2998 ASSERT(internalGCells == noInternalCells(), "");
2999#endif
3000}
3001
3007template <MInt nDim, class SolverType>
3008void CartesianSolver<nDim, SolverType>::saveGridFlowVars(
3009 const MChar* fileName, const MChar* gridFileName, const MInt noTotalCells, const MInt noInternal,
3010 MFloatScratchSpace& dbVariables, std::vector<MString>& dbVariablesName, MInt noDbVars,
3011 MIntScratchSpace& idVariables, std::vector<MString>& idVariablesName, MInt noIdVars,
3012 MFloatScratchSpace& dbParameters, std::vector<MString>& dbParametersName, MIntScratchSpace& idParameters,
3013 std::vector<MString>& idParametersName, MInt* recalcIds, MFloat time) {
3014 TRACE();
3015
3016 MString tmpNcVariablesName = "";
3017 noDbVars = dbVariablesName.size();
3018 noIdVars = idVariablesName.size();
3019 MInt noIdPars = idParametersName.size();
3020 MInt noDbPars = dbParametersName.size();
3021 std::vector<MString> ncVariablesName;
3022 MFloat* tmpDbVariables;
3023 MInt* tmpIdVariables;
3024
3025 using namespace maia::parallel_io;
3026 ParallelIo parallelIo(fileName, PIO_REPLACE, solver().mpiComm());
3027
3028 MLong num = -1;
3029 MLong longNoInternalCells = (MLong)noInternal;
3030 MPI_Allreduce(&longNoInternalCells, &num, 1, MPI_LONG, MPI_SUM, solver().mpiComm(), AT_, "noInternalCells", "num");
3031
3032 for(MInt j = 0; j < noDbVars; j++) {
3033 tmpNcVariablesName = "variables" + std::to_string(j);
3034 ncVariablesName.push_back(tmpNcVariablesName);
3035 parallelIo.defineArray(PIO_FLOAT, ncVariablesName[j], num);
3036 m_log << dbVariablesName[j].c_str() << std::endl;
3037 parallelIo.setAttribute(dbVariablesName[j], "name", ncVariablesName[j]);
3038 }
3039 MInt k;
3040 for(MInt j = 0; j < noIdVars; j++) {
3041 k = j + noDbVars;
3042 tmpNcVariablesName = "variables" + std::to_string(k);
3043 ncVariablesName.push_back(tmpNcVariablesName);
3044 parallelIo.defineArray(PIO_INT, ncVariablesName[k], num);
3045 parallelIo.setAttribute(idVariablesName[j], "name", ncVariablesName[k]);
3046 }
3047
3048 for(MInt j = 0; j < noIdPars; j++) {
3049 parallelIo.defineScalar(PIO_INT, idParametersName[j]);
3050 }
3051
3052 for(MInt j = 0; j < noDbPars; j++) {
3053 parallelIo.defineScalar(PIO_FLOAT, dbParametersName[j]);
3054 }
3055
3056 parallelIo.setAttribute(gridFileName, "gridFile");
3057 if(g_multiSolverGrid) {
3058 parallelIo.setAttribute(solver().solverId(), "solverId");
3059 }
3060 parallelIo.setAttribute(num, "noCells");
3061 // TODO labels:IO,toenhance make universal by setting as attribute and update header files in testcases!
3062 if(time > -1) {
3063 parallelIo.setAttribute(time, "levelSetTime");
3064 }
3065
3066 ParallelIo::size_type start = 0;
3067 ParallelIo::size_type count = 1;
3068
3069 MInt noInternalCellIds = noInternal;
3070 MPI_Exscan(&noInternalCellIds, &start, 1, MPI_INT, MPI_SUM, solver().mpiComm(), AT_, "noInternalCellIds", "start");
3071 count = noInternalCellIds;
3072
3073 parallelIo.setOffset(count, start);
3074
3075 for(MInt j = 0; j < noDbVars; j++) {
3076 tmpDbVariables = &(dbVariables.p[j * noTotalCells]);
3077 MFloatScratchSpace tmpDouble(count, AT_, "tmpDouble");
3078 if(recalcIds != nullptr) {
3079 for(MInt l = 0; l < count; ++l) {
3080 tmpDouble[l] = tmpDbVariables[recalcIds[l]];
3081 }
3082 } else {
3083 for(MInt l = 0; l < count; ++l) {
3084 tmpDouble[l] = tmpDbVariables[l];
3085 }
3086 }
3087 parallelIo.writeArray(tmpDouble.begin(), ncVariablesName[j]);
3088 }
3089
3090 for(MInt j = 0; j < noIdVars; j++) {
3091 tmpIdVariables = &(idVariables.p[j * noTotalCells]);
3092 MIntScratchSpace tmpInt(count, AT_, "tmpInt");
3093 if(recalcIds != nullptr) {
3094 for(MInt l = 0; l < count; ++l) {
3095 tmpInt[l] = tmpIdVariables[recalcIds[l]];
3096 }
3097 } else {
3098 for(MInt l = 0; l < count; ++l) {
3099 tmpInt[l] = tmpIdVariables[l];
3100 }
3101 }
3102 parallelIo.writeArray(tmpInt.begin(), ncVariablesName[j + noDbVars]);
3103 }
3104
3105 for(MInt j = 0; j < noIdPars; j++) {
3106 parallelIo.writeScalar(idParameters.p[j], idParametersName[j]);
3107 }
3108
3109 for(MInt j = 0; j < noDbPars; j++) {
3110 parallelIo.writeScalar(dbParameters.p[j], dbParametersName[j]);
3111 }
3112
3113 m_log << Scratch::printSelfReport();
3114 m_log << Scratch::printSelfReport();
3115}
3116
3125template <MInt nDim, class SolverType>
3126template <typename T>
3127void CartesianSolver<nDim, SolverType>::exchangeData(T* data, const MInt dataBlockSize) {
3128 TRACE();
3129 if(noNeighborDomains() == 0) {
3130 return;
3131 }
3132
3133 maia::mpi::exchangeData(solver().grid().neighborDomains(), solver().grid().haloCells(), solver().grid().windowCells(),
3134 solver().mpiComm(), data, dataBlockSize);
3135}
3136
3146template <MInt nDim, class SolverType>
3147template <typename T>
3148void CartesianSolver<nDim, SolverType>::exchangeLeafData(std::function<T&(MInt, MInt)> data, const MInt noDat) {
3149 TRACE();
3150
3151 if(grid().noDomains() < 2) {
3152 return;
3153 }
3154
3155 const MInt tag = 613;
3156 auto DTYPE = type_traits<T>::mpiType();
3157 ScratchSpace<T> receiveBuffer(noDat * grid().leafRecSize(), AT_, "windowBuffer");
3158 ScratchSpace<T> sendBuffer(noDat * grid().leafSendSize(), AT_, "windowBuffer");
3159
3160 ScratchSpace<MPI_Request> recvRequests(grid().noLeafRecvNeighborDomains(), AT_, "recvRequests");
3161 ScratchSpace<MPI_Request> sendRequests(grid().noLeafSendNeighborDomains(), AT_, "sendRequests");
3162
3163 // 1) start receiving
3164 MInt receiveCount = 0;
3165 for(MInt n = 0; n < grid().noLeafRecvNeighborDomains(); n++) {
3166 const MInt d = grid().leafRecvNeighborDomain(n);
3167 MPI_Irecv(&(receiveBuffer[receiveCount]), noDat * grid().noLeafHaloCells(d), DTYPE, grid().neighborDomain(d), tag,
3168 grid().mpiComm(), &recvRequests[n], AT_, "receiveBuffer");
3169 receiveCount += noDat * grid().noLeafHaloCells(d);
3170 }
3171
3172 // 2) fill send buffer
3173 MInt sendCount = 0;
3174 for(MInt n = 0; n < grid().noLeafSendNeighborDomains(); n++) {
3175 const MInt d = grid().leafSendNeighborDomain(n);
3176 for(MInt j = 0; j < grid().noLeafWindowCells(d); j++) {
3177 const MInt cellId = grid().leafWindowCell(d, j);
3178 for(MInt k = 0; k < noDat; k++) {
3179 sendBuffer[sendCount] = data(cellId, k);
3180 sendCount++;
3181 }
3182 }
3183 }
3184
3185 // 3) start sending
3186 sendCount = 0;
3187 for(MInt n = 0; n < grid().noLeafSendNeighborDomains(); n++) {
3188 const MInt d = grid().leafSendNeighborDomain(n);
3189 MPI_Isend(&sendBuffer[sendCount], noDat * grid().noLeafWindowCells(d), DTYPE, grid().neighborDomain(d), tag,
3190 grid().mpiComm(), &sendRequests[n], AT_, "&sendBuffer[sendCount]");
3191 sendCount += noDat * grid().noLeafWindowCells(d);
3192 }
3193
3194 // 4) wait for all send and receive requests to finish
3195 MPI_Waitall(grid().noLeafRecvNeighborDomains(), &recvRequests[0], MPI_STATUSES_IGNORE, AT_);
3196
3197 // 5) scatter date from receive buffers
3198 MInt recvCount = 0;
3199 for(MInt n = 0; n < grid().noLeafRecvNeighborDomains(); n++) {
3200 const MInt d = grid().leafRecvNeighborDomain(n);
3201 for(MInt j = 0; j < grid().noLeafHaloCells(d); j++) {
3202 const MInt cellId = grid().leafHaloCell(d, j);
3203 for(MInt k = 0; k < noDat; k++) {
3204 data(cellId, k) = receiveBuffer(recvCount);
3205 recvCount++;
3206 }
3207 }
3208 }
3209
3210 MPI_Waitall(grid().noLeafSendNeighborDomains(), &sendRequests[0], MPI_STATUSES_IGNORE, AT_);
3211}
3212
3227template <MInt nDim, class SolverType>
3228template <class G, class S, class M>
3229void CartesianSolver<nDim, SolverType>::exchangeSparseLeafValues(G getData, S setData, const MInt dataSize,
3230 M cellMapping) {
3231 TRACE();
3232
3233 auto* mpi_send_req = new MPI_Request[grid().noNeighborDomains()];
3234 auto* mpi_recv_req = new MPI_Request[grid().noNeighborDomains()];
3235
3236 MIntScratchSpace sendBufferSize(grid().noNeighborDomains(), AT_, "sendBufferSize");
3237 MIntScratchSpace receiveBufferSize(grid().noNeighborDomains(), AT_, "receiveBufferSize");
3238
3239 MIntScratchSpace sendBuffersInt(grid().noNeighborDomains(), AT_, "sendBuffersInt");
3240 MIntScratchSpace receiveBuffersInt(grid().noNeighborDomains(), AT_, "receiveBuffersInt");
3241
3242 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3243 mpi_send_req[i] = MPI_REQUEST_NULL;
3244 mpi_recv_req[i] = MPI_REQUEST_NULL;
3245 }
3246
3247 // 0. gather information about the number of values to be exchanged:
3248 MInt sendBufferCounter;
3249 MInt sendBufferCounterOverall = 0;
3250 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3251 sendBufferCounter = 0;
3252 sendBufferCounter++;
3253 for(MInt j = 0; j < grid().noLeafWindowCells(i); j++) {
3254 const MInt cellId = grid().leafWindowCell(i, j);
3255 const MInt index = cellMapping(cellId, 0);
3256 if(index < 0) continue;
3257 for(MInt d = 0; d < dataSize; d++) {
3258 ASSERT(!std::isnan(getData(index, d)), grid().tree().globalId(cellId));
3259 }
3260
3261 //+check
3262 sendBufferCounter++;
3263
3264 //+data
3265 sendBufferCounter += dataSize;
3266 }
3267 sendBufferSize[i] = sendBufferCounter;
3268 sendBuffersInt[i] = sendBufferCounter;
3269 sendBufferCounterOverall += sendBufferCounter;
3270 }
3271
3272 // 0.a. exchange number data solvers to be exchanged:
3273 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3274 MPI_Irecv(&receiveBuffersInt[i], 1, MPI_INT, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_recv_req[i], AT_,
3275 "receiveBuffers[i]");
3276 }
3277 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3278 MPI_Isend(&sendBuffersInt[i], 1, MPI_INT, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_send_req[i], AT_,
3279 "sendBuffersInt[i]");
3280 }
3281 MPI_Waitall(grid().noNeighborDomains(), mpi_recv_req, MPI_STATUSES_IGNORE, AT_);
3282 MPI_Waitall(grid().noNeighborDomains(), mpi_send_req, MPI_STATUSES_IGNORE, AT_);
3283
3284 // 0.b setup real exchange framework:
3285 MInt receiveBufferCounterOverall = 0;
3286 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3287 receiveBufferCounterOverall += receiveBuffersInt[i];
3288 receiveBufferSize[i] = receiveBuffersInt[i];
3289 }
3290
3291 MFloatScratchSpace sendBuffersOverall(sendBufferCounterOverall, AT_, "sendBuffersOverall");
3292 MFloatScratchSpace receiveBuffersOverall(receiveBufferCounterOverall, AT_, "receiveBuffersOverall");
3293 MFloatPointerScratchSpace sendBuffers(grid().noNeighborDomains(), AT_, "sendBuffers");
3294 MFloatPointerScratchSpace receiveBuffers(grid().noNeighborDomains(), AT_, "receiveBuffers");
3295 MInt counterSend = 0;
3296 MInt counterReceive = 0;
3297 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3298 sendBuffers.p[i] = &sendBuffersOverall.p[counterSend];
3299 counterSend += sendBufferSize[i];
3300 receiveBuffers.p[i] = &receiveBuffersOverall.p[counterReceive];
3301 counterReceive += receiveBufferSize[i];
3302 }
3303
3304 // sicherheitshalber zuruecksetzen
3305 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3306 mpi_send_req[i] = MPI_REQUEST_NULL;
3307 mpi_recv_req[i] = MPI_REQUEST_NULL;
3308 }
3309
3310 // 1. gather relevant information
3311 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3312 sendBufferCounter = 0;
3313 sendBuffers[i][0] = (MFloat)(-1);
3314 sendBufferCounter++;
3315 for(MInt j = 0; j < grid().noLeafWindowCells(i); j++) {
3316 const MInt cellId = grid().leafWindowCell(i, j);
3317 const MInt index = cellMapping(cellId, 0);
3318 if(index < 0) continue;
3319 // check
3320 sendBuffers[i][sendBufferCounter] = (MFloat)j;
3321 sendBufferCounter++;
3322 // data
3323 for(MInt d = 0; d < dataSize; d++) {
3324 ASSERT(!std::isnan(getData(index, d)), "");
3325 sendBuffers[i][sendBufferCounter] = getData(index, d);
3326 sendBufferCounter++;
3327 }
3328 }
3329 sendBufferSize[i] = sendBufferCounter;
3330 sendBuffers[i][0] = (MFloat)(sendBufferCounter);
3331 }
3332
3333
3334 // 2. receive data
3335 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3336 MInt bufSize = receiveBufferSize[i];
3337 MPI_Irecv(receiveBuffers[i], bufSize, MPI_DOUBLE, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_recv_req[i],
3338 AT_, "receiveBuffers[i]");
3339 }
3340
3341 // 3. send data
3342 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3343 MInt bufSize = sendBufferSize[i];
3344 MPI_Isend(sendBuffers[i], bufSize, MPI_DOUBLE, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_send_req[i], AT_,
3345 "sendBuffers[i]");
3346 }
3347 MPI_Waitall(grid().noNeighborDomains(), mpi_recv_req, MPI_STATUSES_IGNORE, AT_);
3348 MPI_Waitall(grid().noNeighborDomains(), mpi_send_req, MPI_STATUSES_IGNORE, AT_);
3349
3350 // 4. store recieved data
3351 MInt receiveBufferCounter = 0;
3352 MInt j = -1;
3353 for(MInt i = 0; i < grid().noNeighborDomains(); i++) {
3354 receiveBufferCounter = 0;
3355 if(receiveBufferSize[i] != receiveBuffersInt[i]) {
3356 mTerm(1, AT_, " receiveBufferSize doesn't match expected size from previous communication! ");
3357 }
3358 if(receiveBufferSize[i] == 0) {
3359 m_log << "Warning: empty message from rank " << grid().neighborDomain(i) << std::endl;
3360 }
3361 receiveBufferCounter++;
3362 while(receiveBufferCounter < receiveBufferSize[i]) {
3363 j = (MInt)receiveBuffers[i][receiveBufferCounter];
3364 receiveBufferCounter++;
3365 const MInt cellId = grid().leafHaloCell(i, j);
3366 MBool skip = false;
3367
3368 if(cellId > grid().tree().size()) skip = true;
3369 if(!skip) {
3370 // add halo-cutCandidate
3371 // const MInt candId = candidates.size();
3372 // candidates.emplace_back();
3373 // candidates[candId].cellId = cellId;
3374
3375 const MInt index = cellMapping(cellId, 0);
3376
3377 ASSERT(index > -1, "No corresponding halo cell found!");
3378
3379 // add infomation from computeNodalvalues:
3380 for(MInt d = 0; d < dataSize; d++) {
3381 setData(index, d) = receiveBuffers[i][receiveBufferCounter];
3382 receiveBufferCounter++;
3383 }
3384
3385 } else {
3386 receiveBufferCounter += dataSize;
3387 }
3388 }
3389 }
3390
3391 delete[] mpi_send_req;
3392 delete[] mpi_recv_req;
3393}
3394
3406template <MInt nDim, class SolverType>
3407template <typename T>
3408void CartesianSolver<nDim, SolverType>::exchangeAzimuthalPer(T* data, MInt dataBlockSize, MInt firstBlock) {
3409 TRACE();
3410
3411 if(grid().noAzimuthalNeighborDomains() == 0) {
3412 return;
3413 }
3414
3415 MUint sndSize = maia::mpi::getBufferSize(grid().azimuthalWindowCells());
3416 ScratchSpace<T> windowData(sndSize * dataBlockSize, AT_, "windowData");
3417 windowData.fill(0);
3418 MUint rcvSize = maia::mpi::getBufferSize(grid().azimuthalHaloCells());
3419 ScratchSpace<T> haloData(rcvSize * dataBlockSize, AT_, "haloData");
3420 haloData.fill(0);
3421
3422 MInt sndCnt = 0;
3423 if(std::is_same<T, MInt>::value || std::is_same<T, MBool>::value || std::is_same<T, MLong>::value) {
3424 for(MInt i = 0; i < grid().noAzimuthalNeighborDomains(); i++) {
3425 for(MInt j = 0; j < grid().noAzimuthalWindowCells(i); j++) {
3426 MInt windowId = grid().azimuthalWindowCell(i, j);
3427 for(MInt b = firstBlock; b < firstBlock + dataBlockSize; b++) {
3428 windowData[sndCnt] = data[windowId * dataBlockSize + b];
3429 sndCnt++;
3430 }
3431 }
3432 }
3433 } else if(std::is_same<T, MFloat>::value) {
3434 std::function<MBool(const MInt, const MInt)> neighborCheck = [&](const MInt cell, const MInt id) {
3435 return static_cast<MBool>(grid().tree().hasNeighbor(cell, id));
3436 };
3437 std::function<MFloat(const MInt, const MInt)> coordinate = [&](const MInt cell, const MInt id) {
3438 return static_cast<MFloat>(grid().tree().coordinate(cell, id));
3439 };
3440 std::function<MFloat(const MInt, const MInt)> scalarField = [&](const MInt cell, const MInt id) {
3441 return data[cell * dataBlockSize + id];
3442 };
3443 MInt cellCnt = 0;
3444 for(MInt i = 0; i < grid().noAzimuthalNeighborDomains(); i++) {
3445 for(MInt j = 0; j < grid().noAzimuthalWindowCells(i); j++) {
3446 MInt windowId = grid().azimuthalWindowCell(i, j);
3447 std::array<MFloat, nDim> recCoord;
3448 std::copy_n(&(solver().m_azimuthalCartRecCoord[cellCnt * nDim]), nDim, &recCoord[0]);
3449
3450 MInt position = -1;
3451 const MInt magic_number = 8; // pow(2, nDim);
3452 std::array<MInt, magic_number> interpolationCells = {0, 0, 0, 0, 0, 0, 0, 0};
3453 position =
3454 setUpInterpolationStencil(windowId, interpolationCells.data(), recCoord.data(), neighborCheck, false);
3455 for(MInt b = firstBlock; b < firstBlock + dataBlockSize; b++) {
3456 if(position > -1) {
3457 windowData[sndCnt] =
3458 interpolateFieldData<true>(&interpolationCells[0], &recCoord[0], b, scalarField, coordinate);
3459 } else {
3460 windowData[sndCnt] = scalarField(windowId, b);
3461 }
3462 sndCnt++;
3463 }
3464 cellCnt++;
3465 }
3466 }
3467 } else {
3468 mTerm(1, AT_, "Non implemented data type.");
3469 }
3470
3471 // Exchange
3472 maia::mpi::exchangeBuffer(grid().azimuthalNeighborDomains(), grid().azimuthalHaloCells(),
3473 grid().azimuthalWindowCells(), mpiComm(), windowData.getPointer(), haloData.getPointer(),
3474 dataBlockSize);
3475
3476 // Extract
3477 MInt rcvCnt = 0;
3478 for(MInt i = 0; i < grid().noAzimuthalNeighborDomains(); i++) {
3479 for(MInt j = 0; j < grid().noAzimuthalHaloCells(i); j++) {
3480 MInt haloId = grid().azimuthalHaloCell(i, j);
3481
3482 for(MInt b = firstBlock; b < firstBlock + dataBlockSize; b++) {
3483 data[haloId * dataBlockSize + b] = haloData[rcvCnt * dataBlockSize + b];
3484 }
3485 rcvCnt++;
3486 }
3487 }
3488
3489 // Unmapped Halos
3490 MBool valueSet;
3491 const MInt magic_number = 126; // Should be the highest possible count
3492 MIntScratchSpace nghbrList(magic_number, AT_, "nghbrList");
3493 for(MInt i = 0; i < grid().noAzimuthalUnmappedHaloCells(); i++) {
3494 MInt haloId = grid().azimuthalUnmappedHaloCell(i);
3495 valueSet = false;
3496 MInt counter = grid().getAdjacentGridCells(haloId, 2, nghbrList, grid().tree().level(haloId), 0);
3497 for(MInt n = 0; n < counter; n++) {
3498 MInt nghbrId = nghbrList[n];
3499 if(nghbrId < 0) {
3500 continue;
3501 }
3502 for(MInt b = firstBlock; b < firstBlock + dataBlockSize; b++) {
3503 data[haloId * dataBlockSize + b] = data[nghbrId * dataBlockSize + b];
3504 }
3505 valueSet = true;
3506 break;
3507 }
3508 if(!valueSet) {
3509 std::cerr << "Unmapped not set:" << domainId() << " " << haloId << " " << grid().tree().coordinate(haloId, 0)
3510 << " " << grid().tree().coordinate(haloId, 1) << " " << grid().tree().coordinate(haloId, 2)
3511 << std::endl;
3512 }
3513 ASSERT(valueSet, "No value set!");
3514 }
3515}
3516
3523template <MInt nDim, class Solver>
3524void CartesianSolver<nDim, Solver>::extractPointIdsFromGrid(Collector<PointBasedCell<nDim>>*& extractedCells,
3525 Collector<CartesianGridPoint<nDim>>*& gridPoints,
3526 const MBool extractHaloCells,
3527 const std::map<MInt, MInt>& splitChildToSplitCell,
3528 MInt levelThreshold, MFloat* bBox, MBool levelSetMb) const {
3529 TRACE();
3530
3531 const MBool isFV =
3532 (string2enum(solver().solverType()) == MAIA_FINITE_VOLUME) || (string2enum(solver().solverType()) == MAIA_FV_MB);
3533
3534#ifndef NDEBUG
3535 if(isFV) {
3536 for(MInt cellId = 0; cellId < solver().c_noCells(); cellId++) {
3537 ASSERT(solver().a_hasProperty(cellId, FvCell::IsSplitChild)
3538 ? splitChildToSplitCell.count(cellId) == 1
3539 && splitChildToSplitCell.find(cellId)->second < solver().c_noCells()
3540 && splitChildToSplitCell.find(cellId)->second > -1
3541 : true,
3542 "associated BndryCell for splitChild is missing.");
3543 }
3544 }
3545#endif
3546
3547
3548 const MInt noCells = solver().c_noCells();
3549 const MInt maxNoGridPoints = noCells + (noCells / 2);
3550 const MInt noPoints = IPOW2(nDim);
3551 const MInt sign[8][3] = {{-1, -1, -1}, {1, -1, -1}, {-1, 1, -1}, {1, 1, -1},
3552 {-1, -1, 1}, {1, -1, 1}, {-1, 1, 1}, {1, 1, 1}};
3553 const MFloat DX = solver().c_cellLengthAtLevel(maxRefinementLevel());
3554 const MFloat EPS = 0.1 * DX;
3555 ScratchSpace<MFloat> cellLength(maxRefinementLevel() + 1, AT_, "cellLength");
3556 ScratchSpace<MBool> cellIsActive(noCells, AT_, "isActive");
3557 for(MInt l = 0; l <= maxRefinementLevel(); l++) {
3558 cellLength.p[l] = solver().c_cellLengthAtLevel(l);
3559 }
3560 MFloat bbox_mem[6];
3561 if(bBox == nullptr) {
3562 bBox = &bbox_mem[0];
3563 solver().geometry().getBoundingBox(bBox);
3564 for(MInt i = 0; i < nDim; i++) {
3565 bBox[i] = bBox[i] - F1B2 * cellLength(maxUniformRefinementLevel());
3566 bBox[nDim + i] = bBox[nDim + i] + F1B2 * cellLength(maxUniformRefinementLevel());
3567 }
3568 }
3569 MFloat bbox2[6];
3570 for(MInt i = 0; i < nDim; i++) {
3571 bBox[i] -= EPS;
3572 bBox[nDim + i] += EPS;
3573 // bbox2[i] = bBox[i] - cellLength(maxRefinementLevel());// - EPS;
3574 // bbox2[nDim+i] = bBox[nDim+i] + cellLength(maxRefinementLevel());// + EPS;
3575 bbox2[i] = bBox[i] - cellLength(maxUniformRefinementLevel());
3576 bbox2[nDim + i] = bBox[nDim + i] + cellLength(maxUniformRefinementLevel());
3577 // bBox[i] = bbox2[i];
3578 // bBox[nDim+i] = bbox2[nDim+i];
3579 }
3580 if(levelThreshold < minLevel()) {
3581 if(domainId() == 0)
3582 std::cerr << "level threshold reset to minimum refinement level, since it was below." << std::endl;
3583 levelThreshold = minLevel();
3584 }
3585 if(levelSetMb) {
3586 cellIsActive.fill(false);
3587 for(MInt cellId = 0; cellId < noCells; cellId++) {
3588 if(solver().a_isBndryGhostCell(cellId)) continue;
3589 if(!solver().a_hasProperty(cellId, FvCell::IsOnCurrentMGLevel)) continue;
3590 if(!solver().a_hasProperty(cellId, FvCell::IsSplitChild) && solver().c_noChildren(cellId) > 0) continue;
3591 cellIsActive(cellId) = true;
3592 MInt parentId = solver().a_hasProperty(cellId, FvCell::IsSplitChild) ? splitChildToSplitCell.find(cellId)->second
3593 : c_parentId(cellId);
3594 while(parentId > -1) {
3595 cellIsActive(parentId) = true;
3596 parentId = c_parentId(parentId);
3597 }
3598 }
3599 }
3600 MUlong N[3] = {1, 1, 1};
3601 for(MInt i = 0; i < nDim; i++)
3602 N[i] = 1 + (size_t)((bbox2[nDim + i] - bbox2[i]) / DX);
3603 if(N[0] * N[1] * N[2] > std::numeric_limits<MUlong>::max()) {
3604 std::cerr << "Warning: MUlong exceeded by hash function!" << std::endl;
3605 }
3606 std::unordered_map<size_t, MInt> table;
3607
3608 //------
3609 // 0. create collectors
3610 if(extractedCells) {
3611 std::cerr << "Warning: extractedCells is not a nullptr pointer as expected." << std::endl;
3612 delete extractedCells;
3613 extractedCells = nullptr;
3614 }
3615 if(gridPoints) {
3616 std::cerr << "Warning: gridPoints is not a nullptr pointer as expected." << std::endl;
3617 delete gridPoints;
3618 gridPoints = nullptr;
3619 }
3620 extractedCells = new Collector<PointBasedCell<nDim>>(noCells, nDim, 0, 0);
3621 gridPoints = new Collector<CartesianGridPoint<nDim>>(maxNoGridPoints, nDim, 0);
3622 if(!extractedCells) {
3623 mTerm(1, AT_, "Allocation of extractedCells failed.");
3624 }
3625 if(!gridPoints) {
3626 mTerm(1, AT_, "Allocation of gridPoints failed.");
3627 }
3628
3629 //------
3630 // 1. determine extracted cells
3631 MInt noExtractedCells = 0;
3632 extractedCells->resetSize(0);
3633 for(MInt cellId = 0; cellId < noCells; cellId++) {
3634 if(a_isBndryGhostCell(cellId)) continue;
3635 if(isFV) {
3636 if((solver().a_hasProperty(cellId, FvCell::IsSplitChild) || solver().c_noChildren(cellId) == 0)
3637 && !solver().a_hasProperty(cellId, FvCell::IsOnCurrentMGLevel))
3638 continue;
3639 if(!solver().a_hasProperty(cellId, FvCell::IsSplitChild) && solver().c_noChildren(cellId) > 0
3640 && solver().a_level(cellId) < levelThreshold)
3641 continue;
3642 if(solver().a_level(
3643 solver().a_hasProperty(cellId, FvCell::IsSplitChild) ? splitChildToSplitCell.find(cellId)->second : cellId)
3644 > levelThreshold)
3645 continue;
3646 if(levelSetMb && !cellIsActive(cellId)) continue;
3647 if(!(/*(grid().azimuthalPeriodicity() && solver().a_isPeriodic(cellId)) ||*/ extractHaloCells)
3648 && solver().a_isHalo(cellId)) {
3649 continue;
3650 }
3651 if(solver().a_level(
3652 solver().a_hasProperty(cellId, FvCell::IsSplitChild) ? splitChildToSplitCell.find(cellId)->second : cellId)
3653 < mMin(maxUniformRefinementLevel(), levelThreshold))
3654 continue;
3655 // if ( a_hasProperty( cellId , SolverCell::IsSplitChild ) ) continue;
3656 if(solver().a_hasProperty(cellId, FvCell::IsSplitCell)) continue;
3657 }
3658 // if ( a_isPeriodic(cellId) ) continue;
3659 MBool outside = false;
3660 for(MInt i = 0; i < nDim; i++) {
3661 if(solver().a_coordinate(cellId, i) < bBox[i] || solver().a_coordinate(cellId, i) > bBox[nDim + i])
3662 outside = true;
3663 }
3664 if(isFV) {
3665 if(grid().azimuthalPeriodicity() && solver().a_isPeriodic(cellId)) {
3666 outside = false;
3667 }
3668 }
3669 if(outside) continue;
3670 extractedCells->append();
3671 extractedCells->a[noExtractedCells].m_cellId = cellId;
3672 for(MInt p = 0; p < noPoints; p++) {
3673 extractedCells->a[noExtractedCells].m_pointIds[p] = -1;
3674 }
3675 noExtractedCells++;
3676 }
3677
3678 //------
3679 // 2. determine grid points
3680 gridPoints->resetSize(0);
3681 for(MInt c = 0; c < noExtractedCells; c++) {
3682 const MInt cellId = extractedCells->a[c].m_cellId;
3683 for(MInt p = 0; p < noPoints; p++) {
3684 MInt gridPointId = -1;
3685 MFloat coords[3] = {std::numeric_limits<MFloat>::max(), std::numeric_limits<MFloat>::max(),
3686 std::numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized
3687 size_t index[3] = {0, 0, 0};
3688 if(isFV) {
3689 for(MInt i = 0; i < nDim; i++) {
3690 coords[i] = solver().a_coordinate(cellId, i)
3691 + sign[p][i] * F1B2
3692 * cellLength.p[solver().a_level(solver().a_hasProperty(cellId, FvCell::IsSplitChild)
3693 ? splitChildToSplitCell.find(cellId)->second
3694 : cellId)];
3695 ASSERT(coords[i] + EPS > bbox2[i] && coords[i] - EPS < bbox2[nDim + i],
3696 std::to_string(cellId) << " "
3697 << std::to_string(
3698 solver().a_level(solver().a_hasProperty(cellId, FvCell::IsSplitChild)
3699 ? splitChildToSplitCell.find(cellId)->second
3700 : cellId))
3701 << " " << std::to_string(i) << " " << std::to_string(EPS) << " "
3702 << std::to_string(coords[i]) << " " << std::to_string(bbox2[i]) << " "
3703 << std::to_string(bbox2[nDim + i]));
3704 index[i] = (size_t)((coords[i] - bbox2[i] + EPS) / DX);
3705 }
3706 } else {
3707 for(MInt i = 0; i < nDim; i++) {
3708 coords[i] = solver().a_coordinate(cellId, i) + sign[p][i] * F1B2 * cellLength.p[solver().a_level(cellId)];
3709 ASSERT(coords[i] + EPS > bbox2[i] && coords[i] - EPS < bbox2[nDim + i],
3710 std::to_string(cellId) << " " << std::to_string(solver().a_level(cellId)) << " " << std::to_string(i)
3711 << " " << std::to_string(EPS) << " " << std::to_string(coords[i]) << " "
3712 << std::to_string(bbox2[i]) << " " << std::to_string(bbox2[nDim + i]));
3713 index[i] = (size_t)((coords[i] - bbox2[i] + EPS) / DX);
3714 }
3715 }
3716 size_t key = index[0] + N[0] * index[1] + N[0] * N[1] * index[2];
3717 std::pair<typename std::unordered_map<size_t, MInt>::iterator, MBool> ret = table.insert(std::make_pair(key, -1));
3718 if(ret.second) {
3719 gridPointId = gridPoints->size();
3720 ASSERT(gridPointId < maxNoGridPoints, "");
3721 gridPoints->append();
3722 gridPoints->a[gridPointId].m_noAdjacentCells = 0;
3723 for(MInt i = 0; i < noPoints; i++)
3724 gridPoints->a[gridPointId].m_cellIds[i] = -1;
3725 for(MInt i = 0; i < nDim; i++)
3726 gridPoints->a[gridPointId].m_coordinates[i] = coords[i];
3727 ret.first->second = gridPointId;
3728 } else {
3729 ASSERT(gridPointId < maxNoGridPoints, "");
3730 gridPointId = ret.first->second;
3731 }
3732 ASSERT(gridPointId > -1 && gridPointId < maxNoGridPoints, "");
3733 extractedCells->a[c].m_pointIds[p] = gridPointId;
3734 MInt rootId =
3735 (solver().a_hasProperty(cellId, FvCell::IsSplitChild)) ? splitChildToSplitCell.find(cellId)->second : cellId;
3736 MBool found = false;
3737 for(MInt n = 0; n < gridPoints->a[gridPointId].m_noAdjacentCells; n++) {
3738 if(gridPoints->a[gridPointId].m_cellIds[n] == rootId) found = true;
3739 }
3740 if(!found) {
3741 // ASSERT( gridPoints->a[ gridPointId ].m_noAdjacentCells > -1 && gridPoints->a[ gridPointId ].m_noAdjacentCells
3742 // < IPOW2(nDim), "");
3743 if(gridPoints->a[gridPointId].m_noAdjacentCells < IPOW2(nDim)) {
3744 gridPoints->a[gridPointId].m_cellIds[gridPoints->a[gridPointId].m_noAdjacentCells] = rootId;
3745 gridPoints->a[gridPointId].m_noAdjacentCells++;
3746 } else
3747 std::cerr << "Warning: grid point with more than " << IPOW2(nDim)
3748 << " neighbor cells: " << gridPoints->a[gridPointId].m_noAdjacentCells << std::endl;
3749 }
3750 }
3751 for(MInt p = 0; p < noPoints; p++) {
3752 if(extractedCells->a[c].m_pointIds[p] < 0) {
3753 std::cerr << "Warning: no point for cell " << cellId << " " << p << std::endl;
3754 }
3755 }
3756 }
3757
3758 //------
3759 // 3. determine grid point connectivity at domain interfaces
3760 if(!extractHaloCells) {
3761 for(MInt d = 0; d < noNeighborDomains(); d++) {
3762 for(MInt c = 0; c < noHaloCells(d); c++) {
3763 MInt cellId = haloCellId(d, c);
3764 if(isFV) {
3765 if((solver().a_hasProperty(cellId, FvCell::IsSplitChild) || solver().c_noChildren(cellId) == 0)
3766 && !solver().a_hasProperty(cellId, FvCell::IsOnCurrentMGLevel))
3767 continue;
3768 if((!solver().a_hasProperty(cellId, FvCell::IsSplitChild) || solver().c_noChildren(cellId) > 0)
3769 && solver().a_level(cellId) < levelThreshold)
3770 continue;
3771 if(levelSetMb && !cellIsActive(cellId)) continue;
3772 if(solver().a_hasProperty(cellId, FvCell::IsNotGradient)) continue;
3773 // if ( solver().a_hasProperty( cellId , FvCell::IsSplitChild ) ) continue;
3774 if(solver().a_hasProperty(cellId, FvCell::IsSplitCell)) continue;
3775 if(solver().a_level(solver().a_hasProperty(cellId, FvCell::IsSplitChild)
3776 ? splitChildToSplitCell.find(cellId)->second
3777 : cellId)
3778 > levelThreshold)
3779 continue;
3780 }
3781 MBool outside = false;
3782 for(MInt i = 0; i < nDim; i++)
3783 if(solver().a_coordinate(cellId, i) < bBox[i] || solver().a_coordinate(cellId, i) > bBox[nDim + i])
3784 outside = true;
3785 if(outside) continue;
3786 for(MInt p = 0; p < noPoints; p++) {
3787 MInt gridPointId = -1;
3788 MFloat coords[3];
3789 size_t index[3] = {0, 0, 0};
3790 if(isFV) {
3791 for(MInt i = 0; i < nDim; i++) {
3792 coords[i] = solver().a_coordinate(cellId, i)
3793 + sign[p][i] * F1B2
3794 * cellLength.p[solver().a_level(solver().a_hasProperty(cellId, FvCell::IsSplitChild)
3795 ? splitChildToSplitCell.find(cellId)->second
3796 : cellId)];
3797 ASSERT(coords[i] + EPS > bbox2[i] && coords[i] - EPS < bbox2[nDim + i], "");
3798 index[i] = (size_t)((coords[i] - bbox2[i] + EPS) / DX);
3799 }
3800 } else {
3801 for(MInt i = 0; i < nDim; i++) {
3802 coords[i] = solver().a_coordinate(cellId, i) + sign[p][i] * F1B2 * cellLength.p[solver().a_level(cellId)];
3803 ASSERT(coords[i] + EPS > bbox2[i] && coords[i] - EPS < bbox2[nDim + i], "");
3804 index[i] = (size_t)((coords[i] - bbox2[i] + EPS) / DX);
3805 }
3806 }
3807 size_t key = index[0] + N[0] * index[1] + N[0] * N[1] * index[2];
3808 std::pair<typename std::unordered_map<size_t, MInt>::iterator, MBool> ret =
3809 table.insert(std::make_pair(key, -1));
3810 if(ret.second)
3811 continue;
3812 else {
3813 ASSERT(gridPointId < maxNoGridPoints, "");
3814 gridPointId = ret.first->second;
3815 }
3816 if(gridPointId < 0) continue;
3817 ASSERT(gridPointId > -1 && gridPointId < maxNoGridPoints, std::to_string(gridPointId));
3818 ASSERT(gridPoints->a[gridPointId].m_noAdjacentCells > -1
3819 && gridPoints->a[gridPointId].m_noAdjacentCells < IPOW2(nDim),
3820 std::to_string(gridPoints->a[gridPointId].m_noAdjacentCells));
3821 gridPoints->a[gridPointId].m_cellIds[gridPoints->a[gridPointId].m_noAdjacentCells] = cellId;
3822 gridPoints->a[gridPointId].m_noAdjacentCells++;
3823 }
3824 }
3825 }
3826 }
3827}
3828
3829
3830} // namespace maia
3831
3832#endif // ifndef CARTESIANSOLVER_H_
mAlloc
void mAlloc(T *&a, const MLong N, const MString &objectName, MString function)
allocates memory for one-dimensional array 'a' of size N
Definition: alloc.h:173
GridCell
GridCell
Grid cell Property Labels.
Definition: cartesiangridcellproperties.h:14
CartesianGridPoint
This class defines a grid point, grid cells consist of.
Definition: cartesiangridpoint.h:18
Collector::resetSize
MInt resetSize(MInt inputSize)
Definition: collector.h:45
Collector::append
void append()
Definition: collector.h:153
Collector::size
MInt size()
Definition: collector.h:29
Context::propertyLength
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
Definition: context.cpp:538
Context::propertyExists
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494
Geometry::getIntersectionElements
virtual MInt getIntersectionElements(MFloat *, std::vector< MInt > &)
Definition: geometry.h:63
PointBasedCell
This class defines a grid cell which is a internal data structure for cartesian cells needed for an e...
Definition: cartesianpointbasedcell.h:21
Scratch::printSelfReport
static MString printSelfReport()
Returns a shortened string summing up the scratch space state information.
Definition: scratch.cpp:129
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
ScratchSpace::size
size_type size() const
Definition: scratch.h:302
ScratchSpace::data
pointer data()
Definition: scratch.h:289
ScratchSpace::getPointer
T * getPointer() const
Deprecated: use begin() instead!
Definition: scratch.h:316
ScratchSpace::fill
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
ScratchSpace::p
pointer p
Deprecated: use [] instead!
Definition: scratch.h:315
ScratchSpace::empty
MBool empty() const
Definition: scratch.h:306
ScratchSpace::begin
iterator begin()
Definition: scratch.h:273
Solver
Parent class of all solvers This class is the base for all solvers. I.e. all solver class (e....
Definition: solver.h:29
Solver::time
virtual MFloat time() const =0
Return the time.
Solver::domainId
virtual MInt domainId() const
Return the domainId (rank)
Definition: solver.h:383
Solver::solverId
MInt solverId() const
Return the solverId.
Definition: solver.h:425
maia::CartesianSolver::exchangeLeafData
void exchangeLeafData(std::function< T &(MInt, MInt)> data, const MInt noDat=1)
Blocking exchange memory in 'data' assuming a solver size of 'dataBlockSize' per cell NOTE: exchange ...
Definition: cartesiansolver.h:3148
maia::CartesianSolver::sensorCutOff
virtual void sensorCutOff(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:150
maia::CartesianSolver::sensorBand
void sensorBand(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
This sensor generates a max refinement band around the cells with max refinement level....
Definition: cartesiansolver.h:2423
maia::CartesianSolver::sensorDerivative
virtual void sensorDerivative(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:97
maia::CartesianSolver::sensorVorticity
virtual void sensorVorticity(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:118
maia::CartesianSolver::c_parentId
MLong c_parentId(const MInt cellId) const
Returns the grid parent id of the cell cellId.
Definition: cartesiansolver.h:161
maia::CartesianSolver::setUpInterpolationStencil
MInt setUpInterpolationStencil(const MInt cellId, MInt *, const MFloat *, std::function< MBool(MInt, MInt)>, MBool allowIncompleteStencil)
Definition: cartesiansolver.h:1623
maia::CartesianSolver::gridInputFileName
MString gridInputFileName() const
Definition: cartesiansolver.h:78
maia::CartesianSolver::removeCellId
void removeCellId(const MInt cellId)
Definition: cartesiansolver.h:912
maia::CartesianSolver::noWindowCells
MInt noWindowCells(const MInt domainId) const
Definition: cartesiansolver.h:76
maia::CartesianSolver::solver
constexpr const SolverType & solver() const
Definition: cartesiansolver.h:256
maia::CartesianSolver::minLevel
MInt minLevel() const
Read-only accessors for grid data.
Definition: cartesiansolver.h:65
maia::CartesianSolver::grid
constexpr GridProxy & grid() const
Definition: cartesiansolver.h:94
maia::CartesianSolver::m_azimuthalCartRecCoord
std::vector< MFloat > m_azimuthalCartRecCoord
Definition: cartesiansolver.h:288
maia::CartesianSolver::domainOffset
MLong domainOffset(const MInt id) const
Definition: cartesiansolver.h:72
maia::CartesianSolver::GridProxy
typename maia::grid::Proxy< nDim > GridProxy
Definition: cartesiansolver.h:54
maia::CartesianSolver::saveGridFlowVars
void saveGridFlowVars(const MChar *fileName, const MChar *gridFileName, const MInt noTotalCells, const MInt noInternal, MFloatScratchSpace &dbVariables, std::vector< MString > &dbVariablesName, MInt noDbVars, MIntScratchSpace &idVariables, std::vector< MString > &idVariablesName, MInt noIdVars, MFloatScratchSpace &dbParameters, std::vector< MString > &dbParametersName, MIntScratchSpace &idParameters, std::vector< MString > &idParametersName, MInt *recalcIds, MFloat time)
This function writes the parallel Netcdf cartesian grid cell based solution/restart file currently us...
Definition: cartesiansolver.h:3008
maia::CartesianSolver::assertValidGridCellId
void assertValidGridCellId(const MInt) const
Definition: cartesiansolver.h:159
maia::CartesianSolver::collectVariables
void collectVariables(T *variablesIn, ScratchSpace< T > &variablesOut, const std::vector< MString > &variablesNameIn, std::vector< MString > &variablesNameOut, const MInt noVars, const MInt noCells, const MBool reverseOrder=false)
generalised helper function for writing restart files! This is necessary for example if the minLevel ...
Definition: cartesiansolver.h:2879
maia::CartesianSolver::maxRefinementLevel
MInt maxRefinementLevel() const
Definition: cartesiansolver.h:68
maia::CartesianSolver::extractPointIdsFromGrid
void extractPointIdsFromGrid(Collector< PointBasedCell< nDim > > *&, Collector< CartesianGridPoint< nDim > > *&, const MBool, const std::map< MInt, MInt > &, MInt levelThreshold=999999, MFloat *bBox=nullptr, MBool levelSetMb=false) const
Creates a list of unique corner points for all cells using a hash map levelThreshold optionally speci...
Definition: cartesiansolver.h:3524
maia::CartesianSolver::receiveWindowTriangles
void receiveWindowTriangles()
Receives triangles from neighbors contained in their window cells and inserts them locally.
Definition: cartesiansolver.h:447
maia::CartesianSolver::collectParameters
void collectParameters(T, ScratchSpace< T > &, const MChar *, std::vector< MString > &)
This function collects a single parameters for the massivley parallel IO functions.
Definition: cartesiansolver.h:2930
maia::CartesianSolver::noHaloCells
MInt noHaloCells(const MInt domainId) const
Definition: cartesiansolver.h:74
maia::CartesianSolver::leastSquaresInterpolation
MFloat leastSquaresInterpolation(MInt *, MFloat *, MInt varId, std::function< MFloat(MInt, MInt)> scalarField, std::function< MFloat(MInt, MInt)> coordinate)
Definition: cartesiansolver.h:1868
maia::CartesianSolver::setBoundaryDistance
void setBoundaryDistance(const MBool *const interfaceCell, const MFloat *const outerBandWidth, MFloatScratchSpace &distance)
transverses over all neighboring cells for a specified length
Definition: cartesiansolver.h:940
maia::CartesianSolver::m_sensorWeight
std::vector< MFloat > m_sensorWeight
Definition: cartesiansolver.h:264
maia::CartesianSolver::interpolateFieldData
MFloat interpolateFieldData(MInt *, MFloat *, MInt varId, std::function< MFloat(MInt, MInt)> scalarField, std::function< MFloat(MInt, MInt)> coordinate)
Definition: cartesiansolver.h:1799
maia::CartesianSolver::c_neighborId
MLong c_neighborId(const MInt cellId, const MInt dir) const
Returns the grid neighbor id of the grid cell cellId dir.
Definition: cartesiansolver.h:163
maia::CartesianSolver::exchangeAzimuthalPer
void exchangeAzimuthalPer(T *data, MInt dataBlockSize=1, MInt firstBlock=0)
Exchange of sparse data structures on max Level.
Definition: cartesiansolver.h:3408
maia::CartesianSolver::collectVariables
void collectVariables(T **variablesIn, ScratchSpace< T > &variablesOut, const std::vector< MString > &variablesNameIn, std::vector< MString > &variablesNameOut, const MInt noVars, const MInt noCells)
generalised helper function for writing restart files! This is necessary for example if the minLevel ...
Definition: cartesiansolver.h:2907
maia::CartesianSolver::cellIsOnGeometry
MBool cellIsOnGeometry(MInt cellId, Geometry< nDim > *geom)
checks whether a cell lies on a certain geometry copied the essential part from identifyBoundaryCells
Definition: cartesiansolver.h:399
maia::CartesianSolver::compactCells
void compactCells()
Removes all holes in the cell collector and moves halo cells to the back of the collector.
Definition: cartesiansolver.h:723
maia::CartesianSolver::saveSensorData
void saveSensorData(const std::vector< std::vector< MFloat > > &sensors, const MInt &level, const MString &gridFileName, const MInt *const recalcIds) override
Saves all sensor values for debug/tunig purposes.
Definition: cartesiansolver.h:2566
maia::CartesianSolver::windowCell
const MInt & windowCell(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:89
maia::CartesianSolver::m_sensorFnPtr
std::vector< fun > m_sensorFnPtr
Definition: cartesiansolver.h:282
maia::CartesianSolver::mapInterpolationCells
void mapInterpolationCells(std::map< MInt, MInt > &cellMap)
Definition: cartesiansolver.h:1932
maia::CartesianSolver::isActive
MBool isActive() const override
Definition: cartesiansolver.h:91
maia::CartesianSolver::windowCellId
MInt windowCellId(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:77
maia::CartesianSolver::sensorSpecies
virtual void sensorSpecies(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:142
maia::CartesianSolver::sensorDivergence
virtual void sensorDivergence(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:101
maia::CartesianSolver::m_sensorType
std::vector< MString > m_sensorType
Definition: cartesiansolver.h:276
maia::CartesianSolver::neighborList
MInt neighborList(const MInt cellId, const MInt dir) const
Definition: cartesiansolver.h:81
maia::CartesianSolver::reOrderCellIds
void reOrderCellIds(std::vector< MInt > &reOrderedCells)
reOrder cellIds before writing the restart file! This is necessary for example if the minLevel shall ...
Definition: cartesiansolver.h:2779
maia::CartesianSolver::centerOfGravity
MFloat centerOfGravity(const MInt dir) const
Definition: cartesiansolver.h:80
maia::CartesianSolver::haloCell
const MInt & haloCell(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:88
maia::CartesianSolver::noNeighborDomains
MInt noNeighborDomains() const
Definition: cartesiansolver.h:70
maia::CartesianSolver::sensorInterface
virtual void sensorInterface(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:122
maia::CartesianSolver::checkNoHaloLayers
void checkNoHaloLayers()
check that each solver has the required number of haloLayers for leaf cells!!! TODO labels:toenhance ...
Definition: cartesiansolver.h:1137
maia::CartesianSolver::createCellId
MInt createCellId(const MInt gridCellId)
Definition: cartesiansolver.h:880
maia::CartesianSolver::sensorLimit
void sensorLimit(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt, std::function< MFloat(MInt)>, const MFloat, const MInt *, const MBool, const MBool allowCoarsening=true)
simple sensor to apply a limit for a value
Definition: cartesiansolver.h:2492
maia::CartesianSolver::sensorMeanStress
virtual void sensorMeanStress(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:134
maia::CartesianSolver::reductionFactor
MFloat reductionFactor() const
Definition: cartesiansolver.h:79
maia::CartesianSolver::haloCellId
MInt haloCellId(const MInt domainId, const MInt cellId) const
Definition: cartesiansolver.h:75
maia::CartesianSolver::localPartitionCellOffsets
MLong localPartitionCellOffsets(const MInt index) const
Definition: cartesiansolver.h:85
maia::CartesianSolver::maxUniformRefinementLevel
MInt maxUniformRefinementLevel() const
Definition: cartesiansolver.h:69
maia::CartesianSolver::sensorEntropyQuot
virtual void sensorEntropyQuot(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:114
maia::CartesianSolver::m_sensorDerivativeVariables
std::vector< MFloat > m_sensorDerivativeVariables
Definition: cartesiansolver.h:265
maia::CartesianSolver::neighborDomain
MInt neighborDomain(const MInt id) const
Definition: cartesiansolver.h:71
maia::CartesianSolver::calcRecalcCellIdsSolver
void calcRecalcCellIdsSolver(const MInt *const recalcIdsTree, MInt &noCells, MInt &noInternalCellIds, std::vector< MInt > &recalcCellIdsSolver, std::vector< MInt > &reorderedCellIds)
Derive recalc cell ids of the solver and reordered cell ids.
Definition: cartesiansolver.h:2949
maia::CartesianSolver::m_maxSensorRefinementLevel
std::vector< MInt > m_maxSensorRefinementLevel
Definition: cartesiansolver.h:263
maia::CartesianSolver::fun
void(CartesianSolver< nDim, SolverType >::*)(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt) fun
Definition: cartesiansolver.h:281
maia::CartesianSolver::sensorTotalPressure
virtual void sensorTotalPressure(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:105
maia::CartesianSolver::markSurrndCells
void markSurrndCells(MIntScratchSpace &inList, const MInt bandWidth, const MInt level, const MBool refineDiagonals=true)
Definition: cartesiansolver.h:1057
maia::CartesianSolver::sensorEntropyGrad
virtual void sensorEntropyGrad(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:110
maia::CartesianSolver::identifyBoundaryCells
void identifyBoundaryCells(MBool *const isInterface, const std::vector< MInt > &bndCndIds=std::vector< MInt >())
Definition: cartesiansolver.h:326
maia::CartesianSolver::patchRefinement
void patchRefinement(std::vector< std::vector< MFloat > > &sensors, std::vector< std::bitset< 64 > > &sensorCellFlag, std::vector< MFloat > &sensorWeight, MInt sensorOffset, MInt sen)
Definition: cartesiansolver.h:2054
maia::CartesianSolver::recomputeGlobalIds
void recomputeGlobalIds(std::vector< MInt > &, std::vector< MLong > &)
reOrder cellIds before writing the restart file! This is necessary for example if the minLevel shall ...
Definition: cartesiansolver.h:2845
maia::CartesianSolver::maxNoGridCells
MInt maxNoGridCells() const
Definition: cartesiansolver.h:67
maia::CartesianSolver::CartesianSolver
CartesianSolver(const MInt solverId, GridProxy &gridProxy_, const MPI_Comm comm, const MBool checkActive=false)
Definition: cartesiansolver.h:311
maia::CartesianSolver::localPartitionCellGlobalIds
const MLong & localPartitionCellGlobalIds(const MInt cellId) const
Definition: cartesiansolver.h:82
maia::CartesianSolver::c_level
MInt c_level(const MInt cellId) const
Definition: cartesiansolver.h:165
maia::CartesianSolver::exchangeData
void exchangeData(T *data, const MInt dataBlockSize=1)
Exchange memory in 'data' assuming a solver size of 'dataBlockSize' per cell.
Definition: cartesiansolver.h:3127
maia::CartesianSolver::sensorPatch
virtual void sensorPatch(std::vector< std::vector< MFloat > > &sensor, std::vector< std::bitset< 64 > > &sensorCellFlag, std::vector< MFloat > &sensorWeight, MInt sensorOffset, MInt sen)
Definition: cartesiansolver.h:146
maia::CartesianSolver::noHaloLayers
MInt noHaloLayers() const
Definition: cartesiansolver.h:73
maia::CartesianSolver::inCell
MInt inCell(const MInt cellId, MFloat *point, MFloat fac=F1)
Definition: cartesiansolver.h:1582
maia::CartesianSolver::minCell
MInt minCell(const MInt id) const
Definition: cartesiansolver.h:87
maia::CartesianSolver::exchangeSparseLeafValues
void exchangeSparseLeafValues(G getData, S setData, const MInt dataSize, M cellMapping)
Exchange of sparse data structures on max Level.
Definition: cartesiansolver.h:3229
maia::CartesianSolver::sensorParticle
virtual void sensorParticle(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
Definition: cartesiansolver.h:138
maia::CartesianSolver::addChildren
void addChildren(std::vector< MInt > &reOrderedIds, const MInt parentId)
add childs to reOrdered cellIds This is necessary for example if the minLevel shall be raised at the ...
Definition: cartesiansolver.h:2826
maia::CartesianSolver::c_globalGridId
MLong c_globalGridId(const MInt cellId)
Definition: cartesiansolver.h:168
maia::CartesianSolver::m_patchRefinement
PatchRefinement * m_patchRefinement
Definition: cartesiansolver.h:304
maia::CartesianSolver::sensorSmooth
void sensorSmooth(std::vector< std::vector< MFloat > > &, std::vector< std::bitset< 64 > > &, std::vector< MFloat > &, MInt, MInt)
sensor to smooth level jumps NOTE: only refines additional cells to ensure a smooth level transition ...
Definition: cartesiansolver.h:2632
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
SPECIES
@ SPECIES
Definition: enums.h:112
CUTOFF
@ CUTOFF
Definition: enums.h:116
SMOOTH
@ SMOOTH
Definition: enums.h:118
DIVERGENCE
@ DIVERGENCE
Definition: enums.h:115
INTERFACE
@ INTERFACE
Definition: enums.h:110
ENTROPY_QUOTIENT
@ ENTROPY_QUOTIENT
Definition: enums.h:108
PATCH
@ PATCH
Definition: enums.h:113
DERIVATIVE
@ DERIVATIVE
Definition: enums.h:106
ENTROPY_GRADIENT
@ ENTROPY_GRADIENT
Definition: enums.h:107
PARTICLE
@ PARTICLE
Definition: enums.h:111
MEANSTRESS
@ MEANSTRESS
Definition: enums.h:117
BAND
@ BAND
Definition: enums.h:119
VORTICITY
@ VORTICITY
Definition: enums.h:109
TOTALPRESSURE
@ TOTALPRESSURE
Definition: enums.h:114
SolverType
SolverType
Definition: enums.h:22
MAIA_FINITE_VOLUME
@ MAIA_FINITE_VOLUME
Definition: enums.h:23
MAIA_FV_MB
@ MAIA_FV_MB
Definition: enums.h:27
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
ipow
MInt ipow(MInt base, MInt exp)
Integer exponent function for non-negative exponents.
Definition: functions.h:317
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
m_revDir
const MInt m_revDir[6]
Definition: fvmbcartesiansolverxd.h:45
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
g_multiSolverGrid
MBool g_multiSolverGrid
Definition: globalvariables.cpp:37
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
IPOW2
constexpr MLong IPOW2(MInt x)
Definition: maiaconstants.h:154
MInt
int32_t MInt
Definition: maiatypes.h:62
MUint
uint32_t MUint
Definition: maiatypes.h:63
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
noSolvers
MInt noSolvers
Definition: maiatypes.h:73
MFloat
double MFloat
Definition: maiatypes.h:52
MLong
int64_t MLong
Definition: maiatypes.h:64
MBool
bool MBool
Definition: maiatypes.h:58
MChar
char MChar
Definition: maiatypes.h:56
id
MInt id
Definition: maiatypes.h:71
MUlong
uint64_t MUlong
Definition: maiatypes.h:65
MPI_Recv
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status, const MString &name, const MString &varname)
same as MPI_Recv
Definition: mpioverride.cpp:565
MPI_Isend
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Isend
Definition: mpioverride.cpp:481
MPI_Irecv
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Irecv
Definition: mpioverride.cpp:606
MPI_Issend
int MPI_Issend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Issend
Definition: mpioverride.cpp:523
MPI_Wait
int MPI_Wait(MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Wait
Definition: mpioverride.cpp:732
MPI_Exscan
int MPI_Exscan(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Exscan
Definition: mpioverride.cpp:1439
MPI_Waitall
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
Definition: mpioverride.cpp:766
MPI_Allreduce
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
MPI_Allgather
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
Definition: mpioverride.cpp:1274
maia::grid::tree::Cell
GridCell Cell
Underlying enum type for property access.
Definition: cartesiangridtree.h:101
maia::math::adjoint1stRow
void adjoint1stRow(std::array< std::array< T, N >, N > &m, std::array< T, N > &A)
Definition: maiamath.cpp:114
maia::math::determinant
MFloat determinant(std::array< T, N > &m)
Definition: maiamath.cpp:130
maia::mpi::getBufferSize
MUint getBufferSize(const std::vector< std::vector< MInt > > &exchangeCells)
Generic exchange of data.
Definition: mpiexchange.h:681
maia::mpi::exchangeData
void exchangeData(const MInt noNghbrDomains, const MInt *const nghbrDomains, const MInt *const noHaloCells, const MInt **const, const MInt *const noWindowCells, const MInt **const windowCells, const MPI_Comm comm, const U *const data, U *const haloBuffer, const MInt noDat=1)
Generic exchange of data.
Definition: mpiexchange.h:295
maia::mpi::exchangeBuffer
void exchangeBuffer(const MInt noExDomains, const MInt *const exDomainId, const MInt *const recvSize, const MInt *const sendSize, const MPI_Comm comm, U *const receiveBuffer, const U *const sendBuffer, const MInt noDat=1)
Generic exchange of data.
Definition: mpiexchange.h:44
maia::parallel_io
Definition: parallelio.h:26
maia::parallel_io::PIO_READ
const MInt PIO_READ
Definition: parallelio.h:40
maia
Namespace for auxiliary functions/classes.
Definition: acaobserverdatacollector.h:76
ParallelIo
PARALLELIO_DEFAULT_BACKEND ParallelIo
Definition: parallelio.h:292
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54
maia::PatchRefinement::noPatchesPerLevel
MInt noPatchesPerLevel(MInt addLevel)
Definition: cartesiansolver.h:40
maia::PatchRefinement::localRfnPatchProperties
MFloat ** localRfnPatchProperties
Definition: cartesiansolver.h:46
maia::PatchRefinement::localRfnLevelPropertiesOffset
MInt * localRfnLevelPropertiesOffset
Definition: cartesiansolver.h:44
maia::PatchRefinement::localRfnLevelMethods
std::vector< MString > localRfnLevelMethods
Definition: cartesiansolver.h:42
maia::type_traits
Definition: typetraits.h:22