CD_ParticleContainer.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_PARTICLECONTAINER_H
#define CD_PARTICLECONTAINER_H
 
// Std includes
#include <cmath>
#include <cstdint>
#include <map>
#include <string>
#include <utility>
#include <vector>
 
// Chombo includes
#include <REAL.H>
#include <SPACE.H>
#include <RealVect.H>
#include <IntVect.H>
#include <Box.H>
#include <ProblemDomain.H>
#include <DisjointBoxLayout.H>
#include <BoxLayout.H>
#include <DataIterator.H>
#include <LayoutData.H>
#include <LevelData.H>
#include <BaseFab.H>
#include <RefCountedPtr.H>
#include <Vector.H>
#include <SPMD.H>
#include <CH_Timer.H>
#include <CH_assert.H>
 
// Our includes
#include <CD_LevelTiles.H>
#include <CD_ParticleGhostMask.H>
#include <CD_ParallelOps.H>
#include <CD_ParticleSoA.H>
#include <CD_NamespaceHeader.H>
 
template <typename P, typename Traits>
using AMRParticlesSoA = Vector<RefCountedPtr<LayoutData<ParticleSoA<P, Traits>>>>;
 
template <typename P = NoPayload, typename Traits = ParticleTraits<P>>
class ParticleContainer
{
public:
  using Leaf = ParticleSoA<P, Traits>;
 
  using LevelParticles = LayoutData<Leaf>;
 
 
  ParticleContainer()
    : m_probLo(RealVect::Zero),
      m_minBlockSize(0),
      m_finestLevel(-1),
      m_numOutcastLocal(0),
      m_isDefined(false),
      m_isOrganizedByCell(false)
  {}
 
  ~ParticleContainer() = default;
 
  ParticleContainer(const ParticleContainer&) = delete;
 
  ParticleContainer&
  operator=(const ParticleContainer&) = delete;
 
  void
  define(const Vector<DisjointBoxLayout>&                       a_grids,
         const Vector<ProblemDomain>&                           a_domains,
         const Vector<Real>&                                    a_dx,
         const Vector<int>&                                     a_refRat,
         const RealVect&                                        a_probLo,
         const int                                              a_minBlockSize,
         const Vector<RefCountedPtr<LevelTiles>>&               a_levelTiles,
         const int                                              a_finestLevel,
         const std::string&                                     a_realm,
         const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>* a_validCells)
  {
    CH_TIME("ParticleContainer::define");
 
    CH_assert(a_grids.size() > a_finestLevel);
    CH_assert(a_domains.size() > a_finestLevel);
    CH_assert(a_dx.size() > a_finestLevel);
    CH_assert(a_refRat.size() > a_finestLevel);
    CH_assert(a_levelTiles.size() > a_finestLevel);
 
    m_grids           = a_grids;
    m_domains         = a_domains;
    m_refRat          = a_refRat;
    m_probLo          = a_probLo;
    m_minBlockSize    = a_minBlockSize;
    m_finestLevel     = a_finestLevel;
    m_realm           = a_realm;
    m_validCells      = a_validCells;
    m_numOutcastLocal = 0;
 
    m_dx.resize(a_finestLevel + 1);
    for (int lvl = 0; lvl <= a_finestLevel; lvl++) {
      m_dx[lvl] = a_dx[lvl] * RealVect::Unit;
    }
 
    // Alias the Realm's per-level tile->box maps (the single source of truth); never build our own.
    m_levelTiles = a_levelTiles;
 
    m_particles.resize(a_finestLevel + 1);
    m_maskParticles.resize(a_finestLevel + 1);
    for (int lvl = 0; lvl <= a_finestLevel; lvl++) {
      m_particles[lvl]     = RefCountedPtr<LevelParticles>(new LevelParticles(a_grids[lvl]));
      m_maskParticles[lvl] = RefCountedPtr<LevelParticles>(new LevelParticles(a_grids[lvl]));
    }
 
    // Buffer particles live on the grown grids (boxes grown by the refinement factor on finer levels).
    this->setupGrownGrids();
    m_bufferParticles.resize(a_finestLevel + 1);
    for (int lvl = 0; lvl <= a_finestLevel; lvl++) {
      m_bufferParticles[lvl] = RefCountedPtr<LevelParticles>(new LevelParticles(m_grownGrids[lvl]));
    }
 
    m_isOrganizedByCell = false;
    m_isDefined         = true;
  }
 
 
 
  const Vector<DisjointBoxLayout>&
  getGrids() const
  {
    return m_grids;
  }
 
  Vector<RealVect>
  getDx() const
  {
    return m_dx;
  }
 
  RealVect
  getProbLo() const
  {
    return m_probLo;
  }
 
  int
  getFinestLevel() const
  {
    return m_finestLevel;
  }
 
  std::string
  getRealm() const
  {
    return m_realm;
  }
 
 
 
  AMRParticlesSoA<P, Traits>&
  getParticles()
  {
    return m_particles;
  }
 
  const AMRParticlesSoA<P, Traits>&
  getParticles() const
  {
    return m_particles;
  }
 
  AMRParticlesSoA<P, Traits>&
  getMaskParticles()
  {
    return m_maskParticles;
  }
 
  const AMRParticlesSoA<P, Traits>&
  getMaskParticles() const
  {
    return m_maskParticles;
  }
 
  AMRParticlesSoA<P, Traits>&
  getBufferParticles()
  {
    return m_bufferParticles;
  }
 
  const AMRParticlesSoA<P, Traits>&
  getBufferParticles() const
  {
    return m_bufferParticles;
  }
 
  LevelParticles&
  operator[](const int a_lvl)
  {
    CH_assert(m_isDefined);
    CH_assert(a_lvl >= 0 && a_lvl <= m_finestLevel);
 
    return *m_particles[a_lvl];
  }
 
  const LevelParticles&
  operator[](const int a_lvl) const
  {
    CH_assert(m_isDefined);
    CH_assert(a_lvl >= 0 && a_lvl <= m_finestLevel);
 
    return *m_particles[a_lvl];
  }
 
 
 
  bool
  isOrganizedByCell() const
  {
    return m_isOrganizedByCell;
  }
 
  bool
  isDefined() const
  {
    return m_isDefined;
  }
 
 
 
  void
  clearParticles()
  {
    CH_TIME("ParticleContainer::clearParticles");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      LevelParticles&          level = *m_particles[lvl];
      const DisjointBoxLayout& dbl   = m_grids[lvl];
      const DataIterator&      dit   = dbl.dataIterator();
      const int                nbox  = dit.size();
 
#pragma omp parallel for schedule(runtime)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din = dit[mybox];
 
        level[din].clear();
      }
    }
  }
 
  void
  addParticlesDestructive(ParticleSoA<P, Traits>& a_particles)
  {
    CH_TIME("ParticleContainer::addParticlesDestructive");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      DataIterator dit = m_grids[lvl].dataIterator();
      if (dit.ok()) {
        (*m_particles[lvl])[dit()].catenate(a_particles); // move the buffer onto a local leaf
        break;
      }
    }
 
    this->remap(); // route every ingested particle to its owning patch/level/rank
  }
 
 
 
  unsigned long long
  getNumberOfValidParticlesLocal() const
  {
    CH_TIME("ParticleContainer::getNumberOfValidParticlesLocal");
 
    CH_assert(m_isDefined);
 
    unsigned long long n = 0;
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      const LevelParticles&    level = *m_particles[lvl];
      const DisjointBoxLayout& dbl   = m_grids[lvl];
      const DataIterator&      dit   = dbl.dataIterator();
      const int                nbox  = dit.size();
 
#pragma omp parallel for schedule(runtime) reduction(+ : n)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din  = dit[mybox];
        const Leaf&      leaf = level[din];
 
        for (std::size_t i = 0; i < leaf.size(); i++) {
          if (!leaf.isGhost(i)) {
            n++;
          }
        }
      }
    }
 
    return n;
  }
 
  unsigned long long
  getNumberOfValidParticlesGlobal() const
  {
    CH_TIME("ParticleContainer::getNumberOfValidParticlesGlobal");
 
    CH_assert(m_isDefined);
 
    return ParallelOps::sum(this->getNumberOfValidParticlesLocal());
  }
 
  unsigned long long
  getNumberOfOutcastParticlesLocal() const
  {
    return m_numOutcastLocal;
  }
 
  unsigned long long
  getNumberOfOutcastParticlesGlobal() const
  {
    return ParallelOps::sum(m_numOutcastLocal);
  }
 
 
 
  void
  remap();
 
  void
  preRegrid();
 
  void
  regrid(const Vector<DisjointBoxLayout>&         a_grids,
         const Vector<ProblemDomain>&             a_domains,
         const Vector<Real>&                      a_dx,
         const Vector<int>&                       a_refRat,
         const int                                a_minBlockSize,
         const Vector<RefCountedPtr<LevelTiles>>& a_levelTiles,
         const int                                a_newFinestLevel);
 
 
 
  void
  fillGhostParticles(const AMRParticleGhostMask& a_maskSame,
                     const AMRParticleGhostMask& a_maskCoarToFine,
                     const AMRParticleGhostMask& a_maskFineToCoar);
 
  void
  clearGhostParticles();
 
  void
  claimOwnership(const bool a_onlyValidRegion);
 
  long long
  resetParticleIDs(const bool a_onlyValidRegion, const ParticleID a_startID);
 
 
 
  void
  copyMaskParticles(const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>& a_mask);
 
  void
  copyMaskParticles(const int a_level, const LevelData<BaseFab<bool>>& a_mask);
 
  void
  transferMaskParticles(const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>& a_mask);
 
  void
  transferMaskParticles(const int a_level, const LevelData<BaseFab<bool>>& a_mask);
 
  void
  clearMaskParticles()
  {
    CH_TIME("ParticleContainer::clearMaskParticles");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      LevelParticles&          level = *m_maskParticles[lvl];
      const DisjointBoxLayout& dbl   = m_grids[lvl];
      const DataIterator&      dit   = dbl.dataIterator();
 
      const int nbox = dit.size();
 
#pragma omp parallel for schedule(runtime)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din = dit[mybox];
 
        level[din].clear();
      }
    }
  }
 
  void
  clearBufferParticles()
  {
    CH_TIME("ParticleContainer::clearBufferParticles");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      LevelParticles&     level = *m_bufferParticles[lvl];
      const BoxLayout&    dbl   = m_grownGrids[lvl];
      const DataIterator& dit   = dbl.dataIterator();
 
      const int nbox = dit.size();
 
#pragma omp parallel for schedule(runtime)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din = dit[mybox];
 
        level[din].clear();
      }
    }
  }
 
  const Vector<BoxLayout>&
  getGrownGrids() const
  {
    return m_grownGrids;
  }
 
  void
  transferParticles(AMRParticlesSoA<P, Traits>& a_source)
  {
    CH_TIME("ParticleContainer::transferParticles");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      LevelParticles& src = *a_source[lvl];
      LevelParticles& dst = *m_particles[lvl];
 
      const DisjointBoxLayout& dbl = m_grids[lvl];
      const DataIterator&      dit = dbl.dataIterator();
 
      const int nbox = dit.size();
 
#pragma omp parallel for schedule(runtime)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din = dit[mybox];
 
        dst[din].catenate(src[din]); // bulk move (id/rank preserved); empty dst -> O(1) swap
      }
    }
 
    m_isOrganizedByCell = false;
  }
 
 
 
  void
  organizeParticlesByCell()
  {
    CH_TIME("ParticleContainer::organizeParticlesByCell");
 
    CH_assert(m_isDefined);
 
    for (int lvl = 0; lvl <= m_finestLevel; lvl++) {
      LevelParticles& level = *m_particles[lvl];
 
      const DisjointBoxLayout& dbl = m_grids[lvl];
      const DataIterator&      dit = dbl.dataIterator();
 
      const int nbox = dit.size();
 
#pragma omp parallel for schedule(runtime)
      for (int mybox = 0; mybox < nbox; mybox++) {
        const DataIndex& din = dit[mybox];
 
        level[din].sortByCell(dbl[din], m_dx[lvl], m_probLo);
      }
    }
 
    m_isOrganizedByCell = true;
  }
 
  void
  organizeParticlesByPatch()
  {
    CH_TIME("ParticleContainer::organizeParticlesByPatch");
 
    m_isOrganizedByCell = false;
  }
 
 
  LevelTiles::LevelAndBox
  findDestination(const RealVect& a_pos) const;
 
protected:
  void
  setupGrownGrids();
 
  using PoolKey = std::pair<int, unsigned int>;
 
  using MoverPool = std::vector<std::map<PoolKey, Leaf>>;
 
  MoverPool
  gatherToPool(AMRParticlesSoA<P, Traits>& a_source, const Vector<DisjointBoxLayout>& a_sourceGrids);
 
  MoverPool
  gatherMoversToPool();
 
  MoverPool
  gatherGhostsFromMasks(const AMRParticleGhostMask& a_maskSame,
                        const AMRParticleGhostMask& a_maskCoarToFine,
                        const AMRParticleGhostMask& a_maskFineToCoar);
 
  void
  distributeFromPool(MoverPool& a_pool, const bool a_setOwnerToReceiver);
 
  Vector<DisjointBoxLayout> m_grids;
 
  Vector<ProblemDomain> m_domains;
 
  Vector<RealVect> m_dx;
 
  Vector<int> m_refRat;
 
  Vector<RefCountedPtr<LevelTiles>> m_levelTiles;
 
  RealVect m_probLo;
 
  int m_minBlockSize;
 
  int m_finestLevel;
 
  unsigned long long m_numOutcastLocal;
 
  std::string m_realm;
 
  const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>* m_validCells = nullptr;
 
  bool m_isDefined;
 
  bool m_isOrganizedByCell;
 
  AMRParticlesSoA<P, Traits> m_particles;
 
  AMRParticlesSoA<P, Traits> m_cacheParticles;
 
  Vector<DisjointBoxLayout> m_cacheGrids;
 
  AMRParticlesSoA<P, Traits> m_maskParticles;
 
  AMRParticlesSoA<P, Traits> m_bufferParticles;
 
  Vector<BoxLayout> m_grownGrids;
 
  inline IntVect
  cellIndex(const int a_lvl, const RealVect& a_pos) const
  {
    IntVect iv;
    for (int dir = 0; dir < SpaceDim; dir++) {
      iv[dir] = static_cast<int>(std::floor((a_pos[dir] - m_probLo[dir]) / m_dx[a_lvl][dir]));
    }
    return iv;
  }
 
  inline bool
  inValidBox(const int a_lvl, const Box& a_box, const RealVect& a_pos) const
  {
    return a_box.contains(this->cellIndex(a_lvl, a_pos));
  }
};
 
#include <CD_NamespaceFooter.H>
 
#include <CD_ParticleContainerImplem.H>
 
#endif