cartesiangrid.h

Go to the documentation of this file.

// Copyright (C) 2024 The m-AIA AUTHORS
//
// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
//
// SPDX-License-Identifier: LGPL-3.0-only
 
#ifndef CARTESIANGRID_H
#define CARTESIANGRID_H
//#define IMPAIR_SUBDIVISION
 
#include <algorithm>
#include <array>
#include <cstring>
#include <functional>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <unordered_map>
#include <vector>
 
#include "COMM/mpioverride.h"
#include "INCLUDE/maiatypes.h"
#include "MEMORY/alloc.h"
#include "MEMORY/scratch.h"
#include "UTIL/debug.h"
#include "UTIL/functions.h"
#include "UTIL/timer.h"
#include "cartesiangridtree.h"
 
template <class T>
class Collector;
 
// Forward declaration of cartesiangridio.h:
namespace maia {
namespace grid {
template <typename Grid>
class IO;
template <MInt nDim>
class Controller;
} // namespace grid
} // namespace maia
 
template <MInt nDim>
class CartesianGrid {
  template <MInt nDim_, class SysEqn>
  friend class FvCartesianSolverXD;
  template <MInt nDim_>
  friend class LsCartesianSolver;
  template <MInt nDim_, class SysEqn>
  friend class FvCartesianSolverXD;
  template <MInt nDim_, class SysEqn>
  friend class FvBndryCnd2D;
  template <MInt nDim_, class SysEqn>
  friend class FvBndryCnd3D;
  template <MInt nDim_, class SysEqn>
  friend class FvMbCartesianSolverXD;
  template <MInt nDim_>
  friend class LbSolver;
#if not defined(MAIA_MS_COMPILER)
  template <MInt nDim_, class SysEqn>
  friend class DgCartesianSolver;
#endif
  friend class maia::grid::Controller<nDim>;
 
 public:
  // cartesiangridio.h API:
  friend class maia::grid::IO<CartesianGrid<nDim>>;
  static constexpr MInt nDim_() { return nDim; }
  using MInt_dyn_array = MInt*;
  using Tree = maia::grid::tree::Tree<nDim>;
 
  // Type for cell properties
  using Cell = typename Tree::Cell;
 
  MBool hasCutOff() const { return m_cutOff; }
 
  // NOTE: the cartesiangrid requires parent/child/neighbor relations as Longs
  //      during initialisation and balance of the grid,
  //      localToGlobalIds and changeGlobalToLocalIds switches between local and global storage!
  MLong& a_parentId(const MInt cellId) { return m_tree.parent(cellId); }
 
  MLong a_parentId(const MInt cellId) const { return m_tree.parent(cellId); }
 
  MLong a_parentId(const MInt cellId, const MInt solverId) const {
    return (m_tree.parent(cellId) > -1 && m_tree.solver(m_tree.parent(cellId), solverId)) ? m_tree.parent(cellId) : -1;
  }
 
  MBool a_hasParent(const MInt cellId) const { return m_tree.hasParent(cellId); }
 
  MBool a_hasParent(const MInt cellId, const MInt solverId) const {
    return m_tree.hasParent(cellId) && m_tree.solver(m_tree.parent(cellId), solverId);
  }
 
  MLong& a_childId(const MInt cellId, const MInt pos) { return m_tree.child(cellId, pos); }
 
  MLong a_childId(const MInt cellId, const MInt pos) const { return m_tree.child(cellId, pos); }
 
  MLong a_childId(const MInt cellId, const MInt pos, const MInt solverId) const {
    return (m_tree.child(cellId, pos) > -1 && m_tree.solver(m_tree.child(cellId, pos), solverId))
               ? m_tree.child(cellId, pos)
               : -1;
  }
 
  MLong noPartitionCellsGlobal() const { return m_noPartitionCellsGlobal; }
 
  MInt& a_level(const MInt cellId) { return m_tree.level(cellId); }
 
  MInt a_level(const MInt cellId) const { return m_tree.level(cellId); }
 
  MFloat a_workload(const MInt cellId) const { return m_tree.workload(cellId); }
 
  MFloat& a_workload(const MInt cellId) { return m_tree.workload(cellId); }
 
  MInt a_noOffsprings(const MInt cellId) const { return m_tree.noOffsprings(cellId); }
 
  MInt& a_noOffsprings(const MInt cellId) { return m_tree.noOffsprings(cellId); }
 
  MLong& a_globalId(const MInt cellId) { return m_tree.globalId(cellId); }
 
  MLong a_globalId(const MInt cellId) const { return m_tree.globalId(cellId); }
 
  MFloat& a_weight(const MInt cellId) { return m_tree.weight(cellId); }
 
  MFloat a_weight(const MInt cellId) const { return m_tree.weight(cellId); }
 
  MInt a_noChildren(const MInt cellId) const { return m_tree.noChildren(cellId); }
 
  MInt a_hasChildren(const MInt cellId) const { return m_tree.hasChildren(cellId); }
 
  MInt a_hasChildren(const MInt cellId, const MInt solverId) const { return m_tree.hasChildren(cellId, solverId); }
 
  MInt a_hasChild(const MInt cellId, const MInt pos) const { return m_tree.hasChild(cellId, pos); }
 
  MInt a_hasChild(const MInt cellId, const MInt pos, const MInt solverId) const {
    return m_tree.hasChild(cellId, pos) && m_tree.solver(m_tree.child(cellId, pos), solverId);
  }
 
  maia::grid::cell::BitsetType::reference a_isToDelete(const MInt cellId) {
    return m_tree.hasProperty(cellId, Cell::IsToDelete);
  }
 
  MBool a_isToDelete(const MInt cellId) const { return m_tree.hasProperty(cellId, Cell::IsToDelete); }
 
  MLong& a_neighborId(const MInt cellId, const MInt dir) { return m_tree.neighbor(cellId, dir); }
 
  MLong a_neighborId(const MInt cellId, const MInt dir) const { return m_tree.neighbor(cellId, dir); }
 
  MLong a_neighborId(const MInt cellId, const MInt dir, const MInt solverId) const {
    return (m_tree.neighbor(cellId, dir) > -1 && m_tree.solver(m_tree.neighbor(cellId, dir), solverId))
               ? m_tree.neighbor(cellId, dir)
               : -1;
  }
 
  MFloat& a_coordinate(const MInt cellId, const MInt dir) { return m_tree.coordinate(cellId, dir); }
 
  MFloat a_coordinate(const MInt cellId, const MInt dir) const { return m_tree.coordinate(cellId, dir); }
 
  template <class U>
  MFloat a_coordinate(Collector<U>* NotUsed(cells_), const MInt cellId, const MInt dir) const {
    return a_coordinate(cellId, dir);
  }
 
 
  MBool a_isLeafCell(const MInt cellId) const { return m_tree.isLeafCell(cellId); }
 
  maia::grid::tree::SolverBitsetType::reference a_isLeafCell(const MInt cellId, const MInt solverId) {
    return m_tree.isLeafCell(cellId, solverId);
  }
 
  MBool a_isLeafCell(const MInt cellId, const MInt solverId) const { return m_tree.isLeafCell(cellId, solverId); }
 
  MBool a_hasNeighbor(const MInt cellId, const MInt dir) const { return m_tree.hasNeighbor(cellId, dir); }
 
  MBool a_hasNeighbor(const MInt cellId, const MInt dir, const MInt solverId) const {
    if(treeb().noSolvers() == 1 && m_tree.hasNeighbor(cellId, dir)
       && !m_tree.solver(m_tree.neighbor(cellId, dir), solverId)) {
      std::cerr << "mis " << cellId << " " << dir << " " << m_tree.neighbor(cellId, dir) << " " << solverId << " "
                << m_tree.solver(m_tree.neighbor(cellId, dir), solverId) << std::endl;
      ASSERT(false, "");
    }
    return m_tree.hasNeighbor(cellId, dir) && m_tree.solver(m_tree.neighbor(cellId, dir), solverId);
  }
 
  void a_resetProperties(const MInt cellId) { m_tree.resetProperties(cellId); }
 
  void a_copyProperties(const MInt fromCellId, const MInt toCellId) {
    m_tree.properties(toCellId) = m_tree.properties(fromCellId);
  }
 
 
  maia::grid::cell::BitsetType::reference a_hasProperty(const MInt cellId, const Cell p) {
    return m_tree.hasProperty(cellId, p);
  }
 
  MBool a_hasProperty(const MInt cellId, const Cell p) const { return m_tree.hasProperty(cellId, p); }
 
  MString a_propertiesToString(const MInt cellId) { return m_tree.propertiesToString(cellId); }
 
  MBool a_isHalo(const MInt cellId) const { return m_tree.hasProperty(cellId, Cell::IsHalo); }
 
  maia::grid::tree::SolverBitsetType::reference a_solver(const MInt cellId, const MInt solverId) {
    return m_tree.solver(cellId, solverId);
  }
 
  void setSolver(const MInt cellId, const MInt solverId, MBool flag) { m_tree.solver(cellId, solverId) = flag; }
 
  MBool addSolverToGrid() { return m_addSolverToGrid; }
 
  MBool wasAdapted() const { return m_wasAdapted; }
  MBool wasBalanced() const { return m_wasBalanced; }
  MBool wasBalancedAtLeastOnce() const { return m_wasBalancedAtLeastOnce; }
 
  MBool m_wasAdapted = false;
  MBool m_wasBalanced = false;
  MBool m_wasBalancedAtLeastOnce = false;
 
  MInt m_noPeriodicCartesianDirs = 0;
  std::array<MInt, 3> m_periodicCartesianDir{};
  MFloat m_periodicCartesianLength[3]{};
  MBool m_azimuthalPer = false;
  std::array<MFloat, 3> m_azimuthalPerCenter{};
  std::array<MInt, 2> m_azimuthalPeriodicDir{};
  MInt m_azimuthalAxialDir = -1;
  MFloat m_azimuthalAngle = F0;
  struct azimuthalBbox {
    MFloat azimuthalAngle = F0;
    MFloat minCoord[nDim];
    MFloat boxLength[nDim];
    MInt nBins[nDim];
    MFloat* minPer3 = nullptr;
    MFloat* maxPer3 = nullptr;
    MFloat center = F0;
    MInt axDir = -1;
    MInt perDir1 = -1;
    MInt perDir2 = -1;
    MInt noHaloLayers = -1;
 
    void init(MFloat angle, MInt periodicDir1, MInt periodicDir2, MInt noLayers) {
      azimuthalAngle = angle;
      perDir1 = periodicDir1;
      perDir2 = periodicDir2;
      noHaloLayers = noLayers;
      MInt perSize1 = nBins[perDir1];
      MInt perSize2 = nBins[perDir2];
      mDeallocate(minPer3);
      mDeallocate(maxPer3);
      mAlloc(minPer3, perSize1 * perSize2, "minPer3", std::numeric_limits<MFloat>::max(), AT_);
      mAlloc(maxPer3, perSize1 * perSize2, "maxPer3", -std::numeric_limits<MFloat>::max(), AT_);
    }
    MInt azimuthalSide(const MFloat phi) {
      MFloat mid = center;
      if(mid > phi) {
        return -1;
      } else {
        return 1;
      }
    }
    MFloat azimuthalCenter() { return center; }
    MFloat azimuthalBoundary(const MFloat* coords, MInt side) {
      MInt perSize2 = nBins[perDir2];
      MInt ind1 = index(coords, perDir1);
      MInt ind2 = index(coords, perDir2);
      if(side < 0) {
        MFloat minPhi = center - F1B2 * azimuthalAngle; // std::numeric_limits<MFloat>::max();
        for(MInt i = -(noHaloLayers + 1); i < (noHaloLayers + 2); i++) {
          MInt ind1_ = ind1 + i;
          ind1_ = mMax(0, mMin(ind1_, nBins[perDir1] - 1));
          MInt ind2_ = ind2 + i;
          ind2_ = mMax(0, mMin(ind2_, nBins[perDir2] - 1));
          minPhi = mMin(minPhi, minPer3[ind1 * perSize2 + ind2_]);
          minPhi = mMin(minPhi, minPer3[ind1_ * perSize2 + ind2]);
        }
        return minPhi;
      } else if(side > 0) {
        MFloat maxPhi = center + F1B2 * azimuthalAngle; //-std::numeric_limits<MFloat>::max();
        for(MInt i = -(noHaloLayers + 1); i < (noHaloLayers + 2); i++) {
          MInt ind1_ = ind1 + i;
          ind1_ = mMax(0, mMin(ind1_, nBins[perDir1] - 1));
          MInt ind2_ = ind2 + i;
          ind2_ = mMax(0, mMin(ind2_, nBins[perDir2] - 1));
          maxPhi = mMax(maxPhi, maxPer3[ind1 * perSize2 + ind2_]);
          maxPhi = mMax(maxPhi, maxPer3[ind1_ * perSize2 + ind2]);
        }
        return maxPhi;
      } else {
        mTerm(1, AT_, "side = 0?");
      }
    }
    MInt index(const MFloat* coord, MInt dir) {
      MInt ind = -1;
      ind = (MInt)((coord[dir] - minCoord[dir]) * nBins[dir] / boxLength[dir]);
      if(ind < 0)
        ind = 0;
      else if(ind >= nBins[dir])
        ind = nBins[dir] - 1;
 
      return ind;
    }
  };
  azimuthalBbox m_azimuthalBbox;
 
  MLong m_noCellsGlobal{};
 
  // Marks the beginning of halo cells in the collector
  MInt m_noInternalCells{};
 
  MInt m_noDomains{};
 
  MInt m_maxUniformRefinementLevel = -1;
  MInt m_maxRfnmntLvl = -1;
  MInt m_maxLevel = -1;
  MInt m_minLevel = -1;
  MInt m_newMinLevel = -1;
 
  MBool m_allowInterfaceRefinement = false;
  MInt m_cutOff = 0;
  MBool m_lowMemAdaptation = true;
 
  // Grid cell volume
  MFloat* m_gridCellVolume = nullptr;
 
  // File processing
  MString m_gridInputFileName = "";
 
  template <class F>
  MInt reduceToLevel(const MInt reductionLevel, F interpolateToParentCells);
 
  MFloat cellLengthAtLevel(const MInt level) const { return m_lengthLevel0 * FFPOW2(level); }
  MFloat cellLengthAtCell(const MInt cellId) const { return cellLengthAtLevel(a_level(cellId)); }
  MFloat halfCellLength(const MInt cellId) const { return F1B2 * cellLengthAtCell(cellId); }
  MFloat cellVolumeAtLevel(const MInt level) const {
    return (nDim == 3) ? POW3(cellLengthAtLevel(level)) : POW2(cellLengthAtLevel(level));
  }
 
  void boundingBox(MFloat* const box) const { std::copy_n(&m_boundingBox[0], 2 * nDim, box); }
 
  const MFloat* globalBoundingBox() const { return &m_boundingBox[0]; }
 
  void centerOfGravity(MFloat* const center) const { std::copy_n(m_centerOfGravity, nDim, center); }
 
  MFloat lengthLevel0() const { return m_lengthLevel0; }
 
  MInt noInternalCells() const { return m_noInternalCells; }
 
  MInt noNeighborDomains() const { return (signed)m_nghbrDomains.size(); }
 
  const MInt& neighborDomain(const MInt id) const { return m_nghbrDomains[id]; }
 
  const std::vector<MInt>& neighborDomains() const { return m_nghbrDomains; }
 
  const MLong& domainOffset(const MInt id) const { return m_domainOffsets[id]; }
 
  MInt noHaloCells(const MInt domainId) const { return m_haloCells[domainId].size(); }
  constexpr MBool paraViewPlugin() const { return m_paraViewPlugin; }
 
  MInt haloCell(const MInt domainId, const MInt cellId) const { return m_haloCells[domainId][cellId]; }
 
  // Returns vector of halo cells
  const std::vector<std::vector<MInt>>& haloCells() const { return m_haloCells; };
 
  MInt noWindowCells(const MInt domainId) const { return m_windowCells[domainId].size(); }
 
  MInt windowCell(const MInt domainId, const MInt cellId) const { return m_windowCells[domainId][cellId]; }
 
  // Returns vector of window cells
  const std::vector<std::vector<MInt>>& windowCells() const { return m_windowCells; };
 
  //
  MInt windowLayer(const MInt domainId, const MInt cellId, const MInt solverId) const {
    return m_windowLayer_[domainId].at(cellId).get(solverId);
  }
 
  MBool isSolverWindowCell(const MInt domainId, const MInt cellId, const MInt solverId) const {
#ifndef NDEBUG
    if(m_windowLayer_[domainId].find(cellId) != m_windowLayer_[domainId].end()
       && m_windowLayer_[domainId].at(cellId).get(solverId) <= m_noSolverHaloLayers[solverId])
      TERMM_IF_COND(!m_tree.solver(cellId, solverId), "");
#endif
    return m_windowLayer_[domainId].find(cellId) != m_windowLayer_[domainId].end()
           && m_windowLayer_[domainId].at(cellId).get(solverId) <= m_noSolverHaloLayers[solverId];
  }
 
  MInt noSolverWindowCells(const MInt domainId, const MInt solverId) const {
    MInt cnt = 0;
    for(const auto& cellId : m_windowCells[domainId]) {
      ASSERT(m_windowLayer_[domainId].find(cellId) != m_windowLayer_[domainId].end(), "");
      if(m_windowLayer_[domainId].at(cellId).get(solverId) <= m_noSolverHaloLayers[solverId]) {
        ASSERT(m_tree.solver(cellId, solverId), "");
        ++cnt;
      }
    }
    return cnt;
  }
 
  MInt noSolverHaloCells(const MInt domainId, const MInt solverId, const MBool periodicCells) const {
    MInt cnt = 0;
    for(const auto& cellId : m_haloCells[domainId]) {
      if(m_tree.solver(cellId, solverId)) {
        if(periodicCells || !a_hasProperty(cellId, Cell::IsPeriodic)) ++cnt;
      }
    }
    return cnt;
  }
 
  MInt noAzimuthalNeighborDomains() const { return (signed)m_azimuthalNghbrDomains.size(); }
 
  const MInt& azimuthalNeighborDomain(const MInt id) const { return m_azimuthalNghbrDomains[id]; }
 
  const std::vector<MInt>& azimuthalNeighborDomains() const { return m_azimuthalNghbrDomains; }
 
  MInt noAzimuthalHaloCells(const MInt domainId) const { return m_azimuthalHaloCells[domainId].size(); }
 
  MInt azimuthalHaloCell(const MInt domainId, const MInt cellId) const {
    return m_azimuthalHaloCells[domainId][cellId];
  }
 
  // Returns vector of azimuthal Halo cells
  const std::vector<std::vector<MInt>>& azimuthalHaloCells() const { return m_azimuthalHaloCells; };
 
  // Returns number of unmapped azimuthal halo cells
  MInt noAzimuthalUnmappedHaloCells() const { return m_azimuthalUnmappedHaloCells.size(); };
 
  // Returns unmapped azimuthal halo cell id
  MInt azimuthalUnmappedHaloCell(const MInt cellId) const { return m_azimuthalUnmappedHaloCells[cellId]; };
 
  // Returns unmapped azimuthal halo cell id
  MInt azimuthalUnmappedHaloDomain(const MInt cellId) const { return m_azimuthalUnmappedHaloDomains[cellId]; };
 
  MInt noAzimuthalWindowCells(const MInt domainId) const { return m_azimuthalWindowCells[domainId].size(); }
 
  MInt azimuthalWindowCell(const MInt domainId, const MInt cellId) const {
    return m_azimuthalWindowCells[domainId][cellId];
  }
 
  // Returns vector of azimuthal window cells
  const std::vector<std::vector<MInt>>& azimuthalWindowCells() const { return m_azimuthalWindowCells; };
  MString gridInputFileName() const { return m_gridInputFileName; }
 
  MInt maxUniformRefinementLevel() const { return m_maxUniformRefinementLevel; }
 
  MFloat reductionFactor() const { return m_reductionFactor; }
 
  MInt maxRefinementLevel() const { return m_maxRfnmntLvl; }
 
  MInt centerOfGravity(MInt dir) const { return m_centerOfGravity[dir]; }
 
  MInt noMinCells() const { return m_minLevelCells.size(); }
 
  MInt minCell(const MInt id) const { return m_minLevelCells[id]; }
 
  MLong noCellsGlobal() const { return m_noCellsGlobal; }
 
  MInt noHaloLayers() const { return m_noHaloLayers; }
  constexpr MInt noHaloLayers(const MInt solverId) const { return m_noSolverHaloLayers[solverId]; }
  // WH_old
  MInt haloMode() const { return m_haloMode; }
 
  MLong localPartitionCellOffsets(const MInt index) const { return m_localPartitionCellOffsets[index]; }
  MInt periodicCartesianDir(const MInt dir) const { return m_periodicCartesianDir[dir]; }
  MFloat periodicCartesianLength(const MInt dir) const { return m_periodicCartesianLength[dir]; }
  MFloat gridCellVolume(const MInt level) const { return m_gridCellVolume[level]; }
  MBool checkOutsideAzimuthalDomain(MFloat*);
  void tagAzimuthalHigherLevelExchangeCells(std::vector<MLong>&, /*std::vector<MLong>&,*/ std::vector<MLong>&, MInt);
  void tagAzimuthalUnmappedHaloCells(std::vector<std::vector<MLong>>&, std::vector<MLong>&, std::vector<MInt>&, MInt);
  void correctAzimuthalSolverBits();
 
  MInt localPartitionCellLocalIds(const MInt id) const { return m_localPartitionCellLocalIds[id]; }
  constexpr MInt noPartitionCells() const { return m_noPartitionCells; }
  constexpr MBool allowInterfaceRefinement() const { return m_allowInterfaceRefinement; }
 
  void gridSanityChecks();
  void checkWindowHaloConsistency(const MBool fullCheck = false, const MString a = "");
  void checkAzimuthalWindowHaloConsistency();
  void dumpCellData(const MString name);
 
  MBool updatedPartitionCells() const { return m_updatedPartitionCells; }
  void setUpdatedPartitionCells(const MBool flag) { m_updatedPartitionCells = flag; }
 
 private:
  MFloat m_centerOfGravity[3]{};
  MFloat m_boundingBox[6]{};
 
  // Remove the following line once transition to new grid format is complete
  // (see also setupWindowHaloCellConnectivity())
  MFloat m_boundingBoxBackup[6]{};
 
  MFloat m_localBoundingBox[6]{};
  MBool m_localBoundingBoxInit = false;
 
  MFloat m_targetGridBoundingBox[6]{};
  MFloat m_targetGridCenterOfGravity[3]{};
  MFloat m_targetGridLengthLevel0 = -1;
  MInt m_targetGridMinLevel = -1;
  MBool m_hasMultiSolverBoundingBox = false;
 
  MFloat m_lengthLevel0{};
 
  // Maximum number of grid cells
  MInt m_maxNoCells = -1;
 
  // Window/Halo related props
  MInt m_haloMode = -1;
  MInt m_noHaloLayers = -1;
  // solver number of haloLayers
  MInt* m_noSolverHaloLayers = nullptr;
 
  // Holds the ids of the adjacent domains
  std::vector<MInt> m_nghbrDomains;
 
  // Holds the ids of the adjacent azimuthal domains
  std::vector<MInt> m_azimuthalNghbrDomains;
 
  // Holds the Backup for to-be-deleted periodic-neighbor-Links
  std::vector<std::tuple<MInt, MInt, MInt>> m_neighborBackup;
 
  // Holds the reverse-direction-Information
  const MInt m_revDir[6] = {1, 0, 3, 2, 5, 4};
 
  std::vector<std::vector<MInt>> m_windowCells;
 
  // Holds the ids of the halo cells according to original adjacent domain
  std::vector<std::vector<MInt>> m_haloCells;
 
  /* m_windowLayer_[nghbrDomId][cellId][solverId] stores the window layer of this cell w.r.t.
   * solverId. The entries are initialized by the higheset possible value of 15 (corresponding
   * to all bits set to 1), which indicates that cell is no window, so effectively allowed range is
   * noHaloLayers<15. Note that this is not exactly the windowLayer, but an auxiliary value
   * required by the algorithm.
   */
  std::vector<std::unordered_map<MInt, M32X4bit<true>>> m_windowLayer_;
  static constexpr const auto WINDOWLAYER_MAX = decltype(m_windowLayer_)::value_type::mapped_type::MAX;
 
  // Azimuthal periodicity exchange data
  std::vector<std::vector<MInt>> m_azimuthalHigherLevelConnectivity;
  std::vector<std::vector<MInt>> m_azimuthalWindowCells;
  std::vector<std::vector<MInt>> m_azimuthalHaloCells;
  std::vector<MInt> m_azimuthalUnmappedHaloCells;
  std::vector<MInt> m_azimuthalUnmappedHaloDomains;
  std::vector<MInt> m_gridBndryCells;
 
  // maximum size of a partitionCell (important for efficient distribution of cells)
  MInt m_partitionCellOffspringThreshold;
  MFloat m_partitionCellWorkloadThreshold;
 
  // offset array
  MLong* m_domainOffsets = nullptr;
 
  // save the offsets of the local partitionCells and the globalIds of the local partitionCells for later use (e.g. the
  // particle restarting) chrs 26.11
  MInt m_noPartitionCells{};
  MLong m_noPartitionCellsGlobal{};
  MInt m_noHaloPartitionLevelAncestors{};
  MLong m_localPartitionCellOffsets[3]{static_cast<MLong>(-1), static_cast<MLong>(-1), static_cast<MLong>(-1)};
  MLong* m_localPartitionCellGlobalIds = nullptr;
  MInt* m_localPartitionCellLocalIds = nullptr;
  MBool m_updatePartitionCellsOnRestart = true;
 
  // Indicates that the partition cells are currently being updated (e.g. to handle the case when
  // there are 0 local partition cells temporarily before balancing)
  MBool m_updatingPartitionCells = false;
  // Indicate if the partition cells have been updated/changed during balancing
  MBool m_updatedPartitionCells = false;
 
  MInt m_noPartitionLevelAncestors{};
 
  MLong m_noPartitionLevelAncestorsGlobal{};
  MLong m_32BitOffset = 0;
  MBool m_restart;
 
  MInt m_noMinLevelCellsGlobal{};
 
  std::vector<MInt> m_minLevelCells{};
 
  MInt m_maxPartitionLevelShift{};
 
  MBool* m_checkRefinementHoles = nullptr;
  MBool* m_diagSmoothing = nullptr;
  MInt m_noIdenticalSolvers = 0;
  MInt* m_identicalSolvers = nullptr;
  MInt* m_identicalSolverMaxLvl = nullptr;
  MInt* m_identicalSolverLvlJumps = nullptr;
 
 private:
 
  std::vector<MInt> m_partitionLevelAncestorIds{};
 
  std::vector<MLong> m_partitionLevelAncestorChildIds{};
  std::vector<MInt> m_partitionLevelAncestorNghbrDomains{};
  std::vector<std::vector<MInt>> m_partitionLevelAncestorHaloCells{};
  std::vector<std::vector<MInt>> m_partitionLevelAncestorWindowCells{};
 
 
  MInt** m_neighborList = nullptr;
 
  MFloat m_reductionFactor{};
  MInt m_decisiveDirection{};
 
  MString m_outputDir = "";
  MString m_restartDir = "";
 
  MBool m_loadGridPartition = false;
 
  MBool m_loadPartition = false;
 
  MBool m_partitionParallelSplit = false;
 
  MBool m_addSolverToGrid = false;
  MInt m_referenceSolver = -1;
 
 public:
  // Data structures for neighbor find and walk algorithms
  // Number of possible paths for 1D neighbors
  MInt m_counter1D{};
  // Number of possible paths for 2D neighbors
  MInt m_counter2D{};
  // Number of possible paths for 3D neighbors
  MInt m_counter3D{};
 
  // Holds all possible paths to 1D neighbors
  MInt* m_paths1D = nullptr;
  // Holds all possible paths to 2D neighbors
  MInt** m_paths2D = nullptr;
  // Holds all possible paths to 3D neighbors
  MInt** m_paths3D = nullptr;
 
  // 2D: 11 elements, 3D: 43 elements
  std::array<MInt, 11 + (nDim - 2) * 32> m_neighborCode{};
 
  std::map<MLong, MInt> m_globalToLocalId{};
 
  // Status flag for additional LB grid checks
  MBool m_lbGridChecks;
  MFloat m_coarseRatio = 0.2;
  MBool m_allowCoarsening = true;
 
  std::set<MInt> m_freeIndices;
 
  // Mapping from domainId to neighborDomainId
  std::vector<MInt> m_nghbrDomainIndex;
 
  // Mapping from domainId to azimuthalNeighborDomainId
  std::vector<MInt> m_azimuthalNghbrDomainIndex;
  MInt* m_localPartitionCellLocalIdsRestart = nullptr;
  MLong m_localPartitionCellOffsetsRestart[3]{static_cast<MLong>(-1), static_cast<MLong>(-1), static_cast<MLong>(-1)};
 
 private:
  // Store the MPI communicator to use for all MPI communication
  const MPI_Comm m_mpiComm;
  Tree m_tree;
  const MBool m_paraViewPlugin;
  const MInt m_maxNoNghbrs;
  static constexpr MInt m_noDirs = 2 * nDim;
  static constexpr MInt m_maxNoChilds = IPOW2(nDim);
  // Jannik TODO labels:GRID for now needed for testcase hack -> change tagActiveWindows
  MBool m_zonal;
 
 public:
  CartesianGrid(MInt maxCells, const MFloat* const bBox, const MPI_Comm comm, const MString& fileName = "");
  ~CartesianGrid();
 
  // Accessor to tree is named "treeb" to indicate that this is the tree of the base grid and not a
  // solver-specific grid tree
  Tree& treeb() { return m_tree; }
  const Tree& treeb() const { return m_tree; }
 
  static constexpr MInt m_maxNoSensors = 64;
 
  MInt maxLevel() const { return m_maxLevel; };
  MInt minLevel() const { return m_minLevel; };
 
  MInt targetGridMinLevel() const { return m_targetGridMinLevel; };
  MInt maxPartitionLevelShift() const { return m_maxPartitionLevelShift; };
 
 
  // Return maximum number of cells
  MInt maxNoCells() const { return m_maxNoCells; };
 
  MPI_Comm mpiComm() const { return m_mpiComm; }
 
  MLong bitOffset() const { return m_32BitOffset; }
 
  template <MBool t_correct = false>
  MInt findContainingPartitionCell(const MFloat* const coord, const MInt solverId = -1,
                                   std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr);
 
  // Return the cell Id of the partition cell intersecting with a sphere
  MInt intersectingWithPartitioncells(MFloat* const center, MFloat const radius);
 
  MInt intersectingWithHaloCells(MFloat* const center, MFloat const radius, MInt domainId, MBool onlyPartition);
 
  MBool intersectingWithLocalBoundingBox(MFloat* const center, MFloat const radius);
 
  MBool boxSphereIntersection(const MFloat* bMin, const MFloat* bMax, const MFloat* const sphereCenter,
                              MFloat const radius);
 
  MBool hollowBoxSphereIntersection(const MFloat* bMin, const MFloat* bMax, const MFloat* const sphereCenter,
                                    MFloat const radius);
 
  inline MInt ancestorPartitionCellId(const MInt cellId) const;
 
  inline MInt findContainingHaloPartitionCell(const MFloat* const coors, const MInt solverId = -1);
 
  template <MBool t_correct = false>
  MInt findContainingHaloCell(const MFloat* const coord, const MInt solverId, MInt domainId, MBool onlyPartitionCells,
                              std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr);
 
  template <MBool t_correct = false>
  MInt findContainingLeafCell(const MFloat* coord,
                              std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr,
                              const MInt solverId = -1);
 
  template <MBool t_correct = false>
  MInt findContainingLeafCell(const std::array<MFloat, nDim>& coord, const MInt startId,
                              std::function<MFloat*(MInt, MFloat* const)>* correctCellCoord = nullptr,
                              const MInt solverId = -1, const MBool allowNonLeafHalo = false) {
    return findContainingLeafCell<t_correct>(&coord[0], startId, correctCellCoord, solverId, allowNonLeafHalo);
  }
  template <MBool t_correct = false>
  MInt findContainingLeafCell(const MFloat* const coord, const MInt startId,
                              std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr,
                              const MInt solverId = -1, const MBool allowNonLeafHalo = false);
 
  template <MBool t_correct = false, MBool insideLimit>
  inline MBool pointWthCell(const std::array<MFloat, nDim>& coord, const MInt cellId,
                            std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr) const {
    return pointWthCell<t_correct, insideLimit>(&coord[0], cellId, correctCellCoord);
  }
 
  template <MBool t_correct = false, MBool insideLimit>
  inline MBool pointWthCell(const MFloat* const coord, const MInt cellId,
                            std::function<MFloat*(MInt, MFloat* const)> correctCellCoord = nullptr) const;
 
  inline MFloat* dummyCorrect(const MInt cellId, MFloat* coords) {
    std::copy(&a_coordinate(cellId, 0), &a_coordinate(cellId, nDim - 1), coords);
    return coords;
  }
 
  void computeLocalBoundingBox(MFloat* const bbox);
 
  MBool pointInLocalBoundingBox(const MFloat* coord);
 
  void propagateDistance(std::vector<MInt>& list, MIntScratchSpace& distMem, MInt dist);
  void propagationStep(MInt cellId, MInt dist, MInt* distMem, MInt endDist);
 
  template <typename DATATYPE>
  void exchangeNotInPlace(DATATYPE* exchangeVar, DATATYPE* recvMem);
  template <typename DATATYPE>
  void generalExchange(MInt noVars, const MInt* vars, DATATYPE** exchangeVar, MInt noDomSend, MInt* domSend,
                       const MInt* noCellsSendPerDom, const MInt* cellIdsSend, MInt noDomRecv, MInt* domRecv,
                       const MInt* noCellsRecvPerDom, const MInt* cellIdsRecv);
 
  // Create a 2D slice gird from 3D grid
  void createGridSlice(const MString& axis, const MFloat intercept, const MString& fileName);
  void createGridSlice(const MString& axis, const MFloat intercept, const MString& fileName, const MInt solverId,
                       MInt* const noSliceCellIds, MInt* const sliceCellIds, MInt* const noSliceHilbertIds = nullptr,
                       MInt* const sliceHilbertInfo = nullptr, MInt* const noSliceContHilbertIds = nullptr,
                       MInt* const sliceContiguousHilbertInfo = nullptr);
 
  void computeGlobalIds();
  void localToGlobalIds();
  void descendNoOffsprings(const MLong cellId, const MLong offset = 0);
  void storeMinLevelCells(const MBool updateMinlevel = false);
  void computeLeafLevel();
  void determineNoPartitionCellsAndOffsets(MLong* const noPartitionCells, MLong* const partitionCellOffsets);
 
  void determinePartLvlAncestorHaloWindowCells(std::vector<std::vector<MInt>>& partLvlAncestorHaloCells,
                                               std::vector<std::vector<MInt>>& partLvlAncestorWindowCells);
 
  void balance(const MInt* const noCellsToReceiveByDomain, const MInt* const noCellsToSendByDomain,
               const MInt* const sortedCellId, const MLong* const offset, const MLong* const globalIdOffsets);
 
  MBool updatePartitionCells(MInt offspringThreshold, MFloat workloadThreshold);
  MLong generateHilbertIndex(const MInt cellId, const MInt targetMinLevel) {
    return hilbertIndexGeneric(&a_coordinate(cellId, 0), m_targetGridCenterOfGravity, m_targetGridLengthLevel0,
                               targetMinLevel);
  }
 
 private:
  void createGlobalToLocalIdMapping();
  void changeGlobalToLocalIds();
  void descendStoreGlobalId(MInt cellId, MInt& localCnt);
  void refineCell(const MInt cellId, const MLong* const refineChildIds = nullptr, const MBool mayHaveChildren = false,
                  const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>& =
                      std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>(),
                  const std::bitset<maia::grid::tree::Tree<nDim>::maxNoSolvers()> refineFlag =
                      std::bitset<maia::grid::tree::Tree<nDim>::maxNoSolvers()>(1ul));
  void removeChilds(const MInt cellId);
  template <MBool removeChilds>
  void removeCell(const MInt cellId);
  MInt deleteCell(const MInt cellId);
  void createPaths();
  void setLevel();
  void updatePartitionCellInformation();
 
  void setGridInputFilename(MBool = false);
  void loadGridFile(const MInt, const MInt);
 
  MLong hilbertIndexGeneric(const MFloat* const coords, const MFloat* const center, const MFloat length0,
                            const MInt baseLevel) const;
 
  inline MLong hilbertIndexGeneric(const MFloat* const coords) const {
    return hilbertIndexGeneric(coords, m_targetGridCenterOfGravity, m_targetGridLengthLevel0, m_targetGridMinLevel);
  }
 
  MLong generateHilbertIndex(const MInt cellId) { return hilbertIndexGeneric(&a_coordinate(cellId, 0)); }
 
  template <class ITERATOR, typename U>
  MInt findIndex(ITERATOR, ITERATOR, const U&);
  template <typename TA, typename TB, typename TC>
  MInt setNeighborDomainIndex(const MInt, std::vector<TA>&, std::vector<TB>&, std::vector<TC>&);
  template <typename TA, typename TB>
  MInt setNeighborDomainIndex(const MInt, std::vector<TA>&, std::vector<TB>&);
 
  MInt setAzimuthalNeighborDomainIndex(const MInt, std::vector<std::vector<MInt>>&, std::vector<std::vector<MInt>>&);
  MInt setAzimuthalNeighborDomainIndex(const MInt, std::vector<std::vector<MInt>>&, std::vector<std::vector<MLong>>&);
 
  template <std::size_t N>
  void exchangeSolverBitset(std::bitset<N>* const data, const MBool defaultVal = false);
  void checkWindowLayer(const MString a);
 
 public:
  void saveGrid(const MChar*, const std::vector<std::vector<MInt>>&, const std::vector<std::vector<MInt>>&,
                const std::vector<std::vector<MInt>>&, const std::vector<MInt>&, const std::vector<std::vector<MInt>>&,
                MInt* recalcIdTree = nullptr);
  void saveGrid(const MChar*, MInt* recalcIdTree);
  void savePartitionFile();
  void loadPartitionFile(const MString& partitionFileName, MLong* partitionCellOffsets);
  void savePartitionCellWorkloadsGridFile();
 
  void deletePeriodicConnection(const MBool = true);
  void restorePeriodicConnection();
  MInt findNeighborDomainId(const MLong globalId) {
    return findNeighborDomainId(globalId, noDomains(), m_domainOffsets);
  };
  MInt findNeighborDomainId(const MLong globalId, const MInt noDomains, const MLong* domainOffsets);
 
  void createLeafCellMapping(const MInt donorSolverId, const MInt gridCellId, std::vector<MInt>& mappedLeafCells,
                             MBool allChilds = false);
  MInt globalIdToLocalId(const MLong& globalId, const MBool termIfNotExisting = false);
 
 private:
  void loadGrid(const MString& fileName);
  void resetCell(const MInt& cellId);
  void initGridMap();
  void createGridMap(const MString& donorGridFileName, const MString& gridMapFileName);
  void saveDonorGridPartition(const MString& gridMapFileName, const MString& gridPartitionFileName);
  void loadDonorGridPartition(const MLong* const partitionCellsId, const MInt noPartitionCells);
 
  void savePartitionFile(const MString& partitionFileNameBase, const MLong partitionCellOffset);
 
  template <typename CELLTYPE>
  void calculateNoOffspringsAndWorkload(Collector<CELLTYPE>* input_cells, MInt input_noCells);
  void partitionParallel(const MInt tmpCount, const MLong tmpOffset, const MFloat* const partitionCellsWorkload,
                         const MLong* const partitionCellsGlobalId, const MFloat totalWorkload,
                         MLong* partitionCellOffsets, MLong* globalIdOffsets, MBool computeOnlyPartition = false);
 
  void updateHaloCellCollectors();
  void tagActiveWindows(std::vector<MLong>&, MInt);
  void tagActiveWindows2_(std::vector<MLong>&, const MInt);
  void tagActiveWindowsOnLeafLvl3(
      const MInt maxLevel, const MBool duringMeshAdaptation,
      const std::vector<std::function<void(const MInt)>>& = std::vector<std::function<void(const MInt)>>());
  void tagActiveWindowsAtMinLevel(const MInt);
  MInt getAdjacentGridCells(MInt, MIntScratchSpace&, MInt, MBool diagonalNeighbors = true);
  template <MBool finer, MBool coarser>
  inline MInt getAdjacentGridCells5(const MInt, MInt* const, const MInt level = -1,
                                    const MBool diagonalNeighbors = true);
  template <MBool finer, MBool coarser>
  inline void getAdjacentGridCells1d5(std::set<MInt>&, const MInt, const MInt, const MInt dimNext = -1,
                                      const MInt dimNextNext = -1, const uint_fast8_t childCodes = (uint_fast8_t)~0);
  void setupWindowHaloCellConnectivity();
  void createMinLevelExchangeCells();
  void createHigherLevelExchangeCells(
      const MInt onlyLevel = -1, const MBool duringMeshAdaptation = false,
      const std::vector<std::function<void(const MInt)>>& = std::vector<std::function<void(const MInt)>>(),
      const std::vector<std::function<void(const MInt)>>& = std::vector<std::function<void(const MInt)>>(),
      const MBool forceLeafLvlCorrection = false);
  void createHigherLevelExchangeCells_old(
      const MInt onlyLevel = -1,
      const std::vector<std::function<void(const MInt)>>& = std::vector<std::function<void(const MInt)>>());
  MInt createAdjacentHaloCell(const MInt, const MInt, const MLong*, std::unordered_multimap<MLong, MInt>&,
                              const MFloat*, MInt* const, std::vector<std::vector<MInt>>&,
                              std::vector<std::vector<MLong>>&);
  MInt createAzimuthalHaloCell(const MInt, const MInt, const MLong*, std::unordered_multimap<MLong, MInt>&,
                               const MFloat*, MInt* const, std::vector<std::vector<MInt>>&,
                               std::vector<std::vector<MLong>>&);
  void meshAdaptation(std::vector<std::vector<MFloat>>&, std::vector<MFloat>&, std::vector<std::bitset<64>>&,
                      std::vector<MInt>&, const std::vector<std::function<void(const MInt)>>&,
                      const std::vector<std::function<void(const MInt)>>&,
                      const std::vector<std::function<void(const MInt)>>&,
                      const std::vector<std::function<void(const MInt, const MInt)>>&,
                      const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>&,
                      const std::vector<std::function<void()>>&);
  void meshAdaptationLowMem(std::vector<std::vector<MFloat>>&, std::vector<MFloat>&, std::vector<std::bitset<64>>&,
                            std::vector<MInt>&, const std::vector<std::function<void(const MInt)>>&,
                            const std::vector<std::function<void(const MInt)>>&,
                            const std::vector<std::function<void(const MInt)>>&,
                            const std::vector<std::function<void(const MInt, const MInt)>>&,
                            const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>&,
                            const std::vector<std::function<void()>>&);
 
  void meshAdaptationDefault(std::vector<std::vector<MFloat>>&, std::vector<MFloat>&, std::vector<std::bitset<64>>&,
                             std::vector<MInt>&, const std::vector<std::function<void(const MInt)>>&,
                             const std::vector<std::function<void(const MInt)>>&,
                             const std::vector<std::function<void(const MInt)>>&,
                             const std::vector<std::function<void(const MInt, const MInt)>>&,
                             const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>&,
                             const std::vector<std::function<void()>>&);
 
 
  void compactCells(const std::vector<std::function<void(const MInt, const MInt)>>& =
                        std::vector<std::function<void(const MInt, const MInt)>>());
  void swapCells(const MInt cellId, const MInt otherId);
  MBool refineCheck(const MInt cellId, const MInt solverId,
                    const std::vector<std::function<MInt(const MFloat*, const MInt, const MInt)>>& cellOutsideSolver);
  MBool coarsenCheck(const MInt cellId, const MInt solverId);
  void exchangeProperties();
  void setChildSolverFlag(const MInt cellId, const MInt solver,
                          const std::function<MInt(const MFloat*, const MInt, const MInt)>& cellOutsideSolver);
 
 public:
  MInt noDomains() const { return m_noDomains; }
 
  MInt domainId() const { return m_domainId; }
 
  void cartesianToCylindric(const MFloat*, MFloat*);
  void rotateCartesianCoordinates(MFloat*, MFloat);
 
 private:
  // holds the domain id
  MInt m_domainId{};
 
  MBool isMpiRoot() const { return domainId() == 0; }
};
 
 
template <MInt nDim>
template <class F>
MInt CartesianGrid<nDim>::reduceToLevel(const MInt reductionLevel, F interpolateToParentCells) {
  TRACE();
  // 1. before we delete, generate a copy of the global ids for the domains so
  // that we can identify them later (as the global id will not be updated
  // during delete cell)
  std::vector<std::vector<MInt>> old_globalhaloCells, old_globalwindowCells;
 
  for(MInt i = 0; i < noNeighborDomains(); i++) {
    std::vector<MInt> tmp_halo;
    std::vector<MInt> tmp_window;
    for(MInt j = 0; j < (signed)m_haloCells[i].size(); j++) {
      tmp_halo.push_back(a_globalId(m_haloCells[i][j]));
    }
    for(MInt j = 0; j < (signed)m_windowCells[i].size(); j++) {
      tmp_window.push_back(a_globalId(m_windowCells[i][j]));
    }
    old_globalhaloCells.push_back(tmp_halo);
    old_globalwindowCells.push_back(tmp_window);
  }
 
  // 2. we need to realloc m_haloCells and m_windowCells so that we
  // suit the requirements of saveGridDomainPar
  m_haloCells.clear();
  m_windowCells.clear();
  m_haloCells.resize(noDomains());
  m_windowCells.resize(noDomains());
 
  // 3. Delete all cells that are above a certain level, interpolate all
  // finer cells to coarser cells
  for(MInt level = m_maxLevel; level > reductionLevel; level--) {
    interpolateToParentCells(level - 1);
    for(MInt c = m_tree.size() - 1; c >= 0; c--) {
      if(a_level(c) == level) {
        deleteCell(c);
      }
    }
  }
  m_noInternalCells = m_tree.size();
 
  // 4. Update window and halo cell arrays
 
  for(MInt c = 0; c < m_tree.size(); c++) {
    if(a_hasProperty(c, Cell::IsHalo)) {
      MInt nghr = -1;
      for(MInt i = 0; i < noNeighborDomains(); i++) {
        for(std::size_t j = 0; j < old_globalhaloCells[i].size(); j++)
          if(a_globalId(c) == old_globalhaloCells[i][j]) {
            nghr = i;
            break;
          }
        if(nghr != -1) {
          break;
        }
      }
      m_haloCells[m_nghbrDomains[nghr]].push_back(c);
 
    } else if(a_hasProperty(c, Cell::IsWindow)) {
      MInt nghr = -1;
      std::vector<MInt> t;
      for(MInt i = 0; i < noNeighborDomains(); i++) {
        for(std::size_t j = 0; j < old_globalwindowCells[i].size(); j++)
          if(a_globalId(c) == old_globalwindowCells[i][j]) {
            nghr = i;
            t.push_back(i);
            break;
          }
        if(nghr != -1) {
          continue;
        }
      }
      for(MInt& i : t) {
        m_windowCells[m_nghbrDomains[i]].push_back(c);
      }
    }
  }
 
 
  // Update the number of internal cells
  for(MInt i = 0; i < noDomains(); i++) {
    m_noInternalCells -= (signed)m_haloCells[i].size();
  }
 
  // 5. Update m_domainOffsets
  MIntScratchSpace rcvBufRed(noDomains(), AT_, "rcvBufRed");
 
  // fill buffer and exchange
  MPI_Allgather(&(m_noInternalCells), 1, MPI_INT, rcvBufRed.begin(), 1, MPI_INT, mpiComm(), AT_, "(m_noInternalCells)",
                "rcvBufRed.begin()");
 
  // put the result into the array
  MInt offset = 0;
  for(MInt i = 0; i < noDomains(); i++) {
    m_domainOffsets[i] = offset;
    offset += rcvBufRed[i];
  }
  // the last element holds the total number
  m_domainOffsets[noDomains()] = offset;
 
  return rcvBufRed[domainId()];
}
 
#endif