CD_EBParticleMesh.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_EBPARTICLEMESH_H
#define CD_EBPARTICLEMESH_H
 
// Std includes
#include <cstddef>
#include <type_traits>
#include <utility>
 
// Chombo includes
#include <CH_Timer.H>
#include <ProblemDomain.H>
#include <EBISBox.H>
#include <EBCellFAB.H>
#include <RealVect.H>
#include <Box.H>
#include <IntVect.H>
 
// Our includes
#include <CD_DepositionType.H>
#include <CD_ParticleOps.H>
#include <CD_BoxLoops.H>
#include <CD_ParticleSoA.H>
#include <CD_NamespaceHeader.H>
 
class EBParticleMesh
{
public:
  EBParticleMesh() : m_verbose(false)
  {}
 
  EBParticleMesh(const ProblemDomain& a_domain,
                 const Box&           a_region,
                 const EBISBox&       a_ebisbox,
                 const RealVect&      a_dx,
                 const RealVect&      a_probLo)
  {
    this->define(a_domain, a_region, a_ebisbox, a_dx, a_probLo);
  }
 
  void
  define(const ProblemDomain& a_domain,
         const Box&           a_region,
         const EBISBox&       a_ebisbox,
         const RealVect&      a_dx,
         const RealVect&      a_probLo)
  {
    CH_TIME("EBParticleMesh::define");
 
    m_verbose = false;
    m_domain  = a_domain;
    m_region  = a_region;
    m_ebisbox = &a_ebisbox; // alias, NOT a copy -- see m_ebisbox declaration
    m_probLo  = a_probLo;
    m_dx      = a_dx;
  }
 
  template <typename P, typename Traits>
  void
  depositWeight(EBCellFAB&                    a_meshData,
                const int                     a_comp,
                const ParticleSoA<P, Traits>& a_particles,
                const DepositionType          a_depositionType,
                const Real                    a_widthScale,
                const bool                    a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::depositWeight");
 
    this->depositCore(a_meshData,
                      a_comp,
                      a_particles,
                      1,
                      a_depositionType,
                      a_widthScale,
                      a_forceIrregNGP,
                      [&](const std::size_t a_i, Real* a_out) {
                        a_out[0] = a_particles.weight(a_i);
                      });
  }
 
  template <typename P, typename Traits>
  void
  depositWeight(EBCellFAB&                    a_meshData,
                const ParticleSoA<P, Traits>& a_particles,
                const DepositionType          a_depositionType,
                const Real                    a_widthScale,
                const bool                    a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::depositWeight");
 
    this->depositWeight(a_meshData, 0, a_particles, a_depositionType, a_widthScale, a_forceIrregNGP);
  }
 
  template <auto... Members, typename P, typename Traits>
  void
  deposit(EBCellFAB&                    a_meshData,
          const int                     a_comp,
          const ParticleSoA<P, Traits>& a_particles,
          const DepositionType          a_depositionType,
          const Real                    a_widthScale,
          const bool                    a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::deposit");
 
    static_assert(sizeof...(Members) == 1 || sizeof...(Members) == SpaceDim,
                  "deposit selects either one scalar column or SpaceDim vector-component columns");
 
    this->depositCore(a_meshData,
                      a_comp,
                      a_particles,
                      static_cast<int>(sizeof...(Members)),
                      a_depositionType,
                      a_widthScale,
                      a_forceIrregNGP,
                      [&](const std::size_t a_i, Real* a_out) {
                        std::size_t k = 0;
                        ((a_out[k++] = static_cast<Real>(a_particles.template get<Members>(a_i))), ...);
                      });
  }
 
  template <auto... Members, typename P, typename Traits>
  void
  deposit(EBCellFAB&                    a_meshData,
          const ParticleSoA<P, Traits>& a_particles,
          const DepositionType          a_depositionType,
          const Real                    a_widthScale,
          const bool                    a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::deposit");
 
    this->deposit<Members...>(a_meshData, 0, a_particles, a_depositionType, a_widthScale, a_forceIrregNGP);
  }
 
  template <typename P, typename Traits, typename GatherFunc>
  void
  depositGathered(EBCellFAB&                    a_meshData,
                  const int                     a_comp,
                  const ParticleSoA<P, Traits>& a_particles,
                  const DepositionType          a_depositionType,
                  const Real                    a_widthScale,
                  const bool                    a_forceIrregNGP,
                  GatherFunc&&                  a_gather) const
  {
    CH_TIME("EBParticleMesh::depositGathered");
 
    this->depositCore(a_meshData,
                      a_comp,
                      a_particles,
                      1,
                      a_depositionType,
                      a_widthScale,
                      a_forceIrregNGP,
                      std::forward<GatherFunc>(a_gather));
  }
 
  template <auto... Members, typename P, typename Traits>
  void
  interpolate(ParticleSoA<P, Traits>& a_particles,
              const EBCellFAB&        a_meshData,
              const int               a_comp,
              const DepositionType    a_interpType,
              const bool              a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::interpolate");
 
    static_assert(sizeof...(Members) == 1 || sizeof...(Members) == SpaceDim,
                  "interpolate targets either one scalar column or SpaceDim vector-component columns");
 
    this->interpolateCore(
      a_particles,
      a_meshData,
      a_comp,
      static_cast<int>(sizeof...(Members)),
      a_interpType,
      a_forceIrregNGP,
      [&](const std::size_t a_i, const Real* a_in) {
        std::size_t k = 0;
        ((a_particles.template get<Members>(
            a_i) = static_cast<std::decay_t<decltype(a_particles.template get<Members>(a_i))>>(a_in[k++])),
         ...);
      });
  }
 
  template <auto... Members, typename P, typename Traits>
  void
  interpolate(ParticleSoA<P, Traits>& a_particles,
              const EBCellFAB&        a_meshData,
              const DepositionType    a_interpType,
              const bool              a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::interpolate");
 
    this->interpolate<Members...>(a_particles, a_meshData, 0, a_interpType, a_forceIrregNGP);
  }
 
  template <typename P, typename Traits>
  void
  interpolateWeight(ParticleSoA<P, Traits>& a_particles,
                    const EBCellFAB&        a_meshData,
                    const int               a_comp,
                    const DepositionType    a_interpType,
                    const bool              a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::interpolateWeight");
 
    this->interpolateCore(a_particles,
                          a_meshData,
                          a_comp,
                          1,
                          a_interpType,
                          a_forceIrregNGP,
                          [&](const std::size_t a_i, const Real* a_in) {
                            a_particles.weight(a_i) = a_in[0];
                          });
  }
 
  template <typename P, typename Traits>
  void
  interpolateWeight(ParticleSoA<P, Traits>& a_particles,
                    const EBCellFAB&        a_meshData,
                    const DepositionType    a_interpType,
                    const bool              a_forceIrregNGP) const
  {
    CH_TIME("EBParticleMesh::interpolateWeight");
 
    this->interpolateWeight(a_particles, a_meshData, 0, a_interpType, a_forceIrregNGP);
  }
 
protected:
  ProblemDomain m_domain;
 
  bool m_verbose;
 
  Box m_region;
 
  const EBISBox* m_ebisbox = nullptr;
 
  RealVect m_dx;
 
  RealVect m_probLo;
 
  template <typename P, typename Traits, typename StrengthFn>
  void
  depositCore(EBCellFAB&                    a_meshData,
              const int                     a_comp,
              const ParticleSoA<P, Traits>& a_particles,
              const int                     a_numComp,
              const DepositionType          a_depositionType,
              const Real                    a_widthScale,
              const bool                    a_forceIrregNGP,
              StrengthFn                    a_gather) const
  {
    CH_TIME("EBParticleMesh::depositCore");
 
    CH_assert(a_meshData.nComp() >= a_comp + a_numComp);
 
    Real invVol = 1.0;
    for (int dir = 0; dir < SpaceDim; dir++) {
      invVol /= m_dx[dir];
    }
 
    const RealVect cicWidth = 1 * a_widthScale * RealVect::Unit;
    const RealVect tscWidth = 2 * a_widthScale * RealVect::Unit;
 
    Real ngpFactor = invVol;
    Real cicFactor = invVol;
    Real tscFactor = invVol;
    for (int dir = 0; dir < SpaceDim; dir++) {
      cicFactor *= 1.0 / cicWidth[dir];
      tscFactor *= 2.0 / tscWidth[dir];
    }
 
    const std::size_t numParticles = a_particles.size();
 
    Real strength[SpaceDim];
 
    switch (a_depositionType) {
    case DepositionType::NGP: {
      for (std::size_t i = 0; i < numParticles; i++) {
        a_gather(i, strength);
        this->depositParticleNGP(a_meshData, a_comp, a_particles.position(i), ngpFactor, strength, a_numComp);
      }
      break;
    }
    case DepositionType::CIC: {
      for (std::size_t i = 0; i < numParticles; i++) {
        a_gather(i, strength);
        this->depositParticleCIC(a_meshData,
                                 a_comp,
                                 a_particles.position(i),
                                 cicWidth,
                                 cicFactor,
                                 strength,
                                 a_numComp,
                                 a_forceIrregNGP);
      }
      break;
    }
    case DepositionType::TSC: {
      for (std::size_t i = 0; i < numParticles; i++) {
        a_gather(i, strength);
        this->depositParticleTSC(a_meshData,
                                 a_comp,
                                 a_particles.position(i),
                                 tscWidth,
                                 tscFactor,
                                 strength,
                                 a_numComp,
                                 a_forceIrregNGP);
      }
      break;
    }
    default: {
      MayDay::Abort("EBParticleMesh::depositCore -- logic bust");
      break;
    }
    }
  }
 
  template <typename P, typename Traits, typename ScatterFn>
  void
  interpolateCore(ParticleSoA<P, Traits>& a_particles,
                  const EBCellFAB&        a_meshData,
                  const int               a_comp,
                  const int               a_numComp,
                  const DepositionType    a_interpType,
                  const bool              a_forceIrregNGP,
                  ScatterFn               a_scatter) const
  {
    CH_TIME("EBParticleMesh::interpolateCore");
 
    CH_assert(a_meshData.nComp() >= a_comp + a_numComp);
 
    // validBox is where the chosen scheme has enough cells; gatherBox is the read stencil.
    Box validBox;
    Box gatherBox;
    switch (a_interpType) {
    case DepositionType::NGP: {
      validBox  = m_domain.domainBox();
      gatherBox = Box(IntVect::Zero, IntVect::Zero);
      break;
    }
    case DepositionType::CIC: {
      validBox  = grow(m_domain.domainBox(), -1);
      gatherBox = Box(IntVect::Zero, IntVect::Unit);
      break;
    }
    case DepositionType::TSC: {
      validBox  = grow(m_domain.domainBox(), -2);
      gatherBox = Box(IntVect::Zero, 2 * IntVect::Unit);
      break;
    }
    default: {
      MayDay::Error("EBParticleMesh::interpolateCore - logic bust");
    }
    }
 
    const std::size_t numParticles = a_particles.size();
 
    Real field[SpaceDim];
 
    switch (a_interpType) {
    case DepositionType::NGP: {
      for (std::size_t i = 0; i < numParticles; i++) {
        this->interpolateParticleNGP(field, a_meshData, a_comp, a_particles.position(i), a_numComp);
        a_scatter(i, field);
      }
      break;
    }
    case DepositionType::CIC: {
      for (std::size_t i = 0; i < numParticles; i++) {
        this->interpolateParticleCIC(field,
                                     a_meshData,
                                     a_comp,
                                     validBox,
                                     gatherBox,
                                     a_particles.position(i),
                                     a_numComp,
                                     a_forceIrregNGP);
        a_scatter(i, field);
      }
      break;
    }
    case DepositionType::TSC: {
      for (std::size_t i = 0; i < numParticles; i++) {
        this->interpolateParticleTSC(field,
                                     a_meshData,
                                     a_comp,
                                     validBox,
                                     gatherBox,
                                     a_particles.position(i),
                                     a_numComp,
                                     a_forceIrregNGP);
        a_scatter(i, field);
      }
      break;
    }
    default: {
      MayDay::Abort("EBParticleMesh::interpolateCore -- logic bust");
      break;
    }
    }
  }
 
  // ===========================================================================
  // The per-particle EB/cut-cell kernels below are adapted from
  // Source/Particle/CD_EBParticleMeshImplem.H -- storage-agnostic, with an added
  // destination-component offset (a_comp) and the TSC partition-of-unity fix.
  // ===========================================================================
 
  inline void
  depositParticleNGP(EBCellFAB&      a_rho,
                     const int       a_comp,
                     const RealVect& a_position,
                     const Real&     a_volumeFactor,
                     const Real*     a_strength,
                     const int&      a_numComp) const noexcept
  {
    const IntVect particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
 
    CH_assert(m_region.contains(particleIV));
 
    FArrayBox& rho = a_rho.getFArrayBox();
 
    for (int comp = 0; comp < a_numComp; comp++) {
      rho(particleIV, a_comp + comp) += a_strength[comp] * a_volumeFactor;
    }
  }
 
  inline Box
  cloudBox(const RealVect& a_position, const RealVect& a_particleWidth) const noexcept
  {
    RealVect lo;
    RealVect hi;
    for (int dir = 0; dir < SpaceDim; dir++) {
      const Real halfWidth = 0.5 * a_particleWidth[dir] * m_dx[dir];
      lo[dir]              = a_position[dir] - halfWidth;
      hi[dir]              = a_position[dir] + halfWidth;
    }
    return Box(ParticleOps::getParticleCellIndex(lo, m_probLo, m_dx),
               ParticleOps::getParticleCellIndex(hi, m_probLo, m_dx));
  }
 
  inline void
  depositParticleCIC(EBCellFAB&      a_rho,
                     const int       a_comp,
                     const RealVect& a_position,
                     const RealVect& a_particleWidth,
                     const Real&     a_volumeFactor,
                     const Real*     a_strength,
                     const int&      a_numComp,
                     const bool      a_forceIrregNGP) const noexcept
  {
    const IntVect particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
 
    CH_assert(m_region.contains(particleIV));
 
    FArrayBox& rho = a_rho.getFArrayBox();
 
    // forceIrregNGP first so the common (false) case short-circuits before the EBISBox query.
    if (a_forceIrregNGP && m_ebisbox->isIrregular(particleIV)) {
      for (int comp = 0; comp < a_numComp; comp++) {
        rho(particleIV, a_comp + comp) += a_strength[comp] * a_volumeFactor;
      }
    }
    else {
      auto cicKernel = [&](const IntVect& iv) -> void {
        Real weight = a_volumeFactor;
 
        for (int dir = 0; dir < SpaceDim; dir++) {
          const Real a = (m_probLo[dir] - a_position[dir]) / m_dx[dir] + iv[dir];
          const Real b = a + 1.0;
          const Real L = 0.5 * a_particleWidth[dir];
 
          weight *= std::max(0.0, std::min(b, L) - std::max(a, -L));
        }
 
        for (int comp = 0; comp < a_numComp; comp++) {
          rho(iv, a_comp + comp) += weight * a_strength[comp];
        }
      };
 
      BoxLoops::loop<D_DECL(1, 1, 1)>(this->cloudBox(a_position, a_particleWidth), cicKernel);
    }
  }
 
  inline void
  depositParticleTSC(EBCellFAB&      a_rho,
                     const int       a_comp,
                     const RealVect& a_position,
                     const RealVect& a_particleWidth,
                     const Real&     a_volumeFactor,
                     const Real*     a_strength,
                     const int&      a_numComp,
                     const bool      a_forceIrregNGP) const noexcept
  {
    const IntVect particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
 
    CH_assert(m_region.contains(particleIV));
 
    FArrayBox& rho = a_rho.getFArrayBox();
 
    // forceIrregNGP first so the common (false) case short-circuits before the EBISBox query.
    if (a_forceIrregNGP && m_ebisbox->isIrregular(particleIV)) {
      for (int comp = 0; comp < a_numComp; comp++) {
        rho(particleIV, a_comp + comp) += a_strength[comp] * a_volumeFactor;
      }
    }
    else {
      auto tscKernel = [&](const IntVect& iv) -> void {
        Real weight = a_volumeFactor;
 
        for (int dir = 0; dir < SpaceDim; dir++) {
          const Real a = (m_probLo[dir] - a_position[dir]) / m_dx[dir] + iv[dir];
          const Real b = a + 1.0;
          const Real L = a_particleWidth[dir];
 
          const Real alpha  = std::max(a, -0.5 * L);
          const Real beta   = std::min(b, +0.5 * L);
          const Real factor = (alpha < beta) ? 1.0 : 0.0;
 
          // factor gates the ENTIRE overlap integral so out-of-support cells contribute zero.
          // This fixes a partition-of-unity bug in the original AoS depositParticleTSC, where
          // factor multiplied only the first term and the cloud over-deposited ~2.33x per
          // dimension.
          weight *= factor * ((beta - alpha) - (beta * std::abs(beta) - alpha * std::abs(alpha)) / L);
        }
 
        for (int comp = 0; comp < a_numComp; comp++) {
          rho(iv, a_comp + comp) += weight * a_strength[comp];
        }
      };
 
      BoxLoops::loop<D_DECL(1, 1, 1)>(this->cloudBox(a_position, a_particleWidth), tscKernel);
    }
  }
 
  inline void
  interpolateParticleNGP(Real*            a_particleField,
                         const EBCellFAB& a_meshData,
                         const int        a_comp,
                         const RealVect&  a_position,
                         const int&       a_numComp) const noexcept
  {
    const IntVect    particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
    const FArrayBox& meshData   = a_meshData.getFArrayBox();
    const Real       factor     = m_ebisbox->isCovered(particleIV) ? 0.0 : 1.0;
 
    for (int comp = 0; comp < a_numComp; comp++) {
      a_particleField[comp] = factor * meshData(particleIV, a_comp + comp);
    }
  }
 
  inline void
  interpolateParticleCIC(Real*            a_particleField,
                         const EBCellFAB& a_meshData,
                         const int        a_comp,
                         const Box&       a_validBox,
                         const Box&       a_gatherBox,
                         const RealVect&  a_position,
                         const int&       a_numComp,
                         const bool       a_forceIrregNGP) const noexcept
  {
    const IntVect particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
    const IntVect loIndex    = ParticleOps::getParticleCellIndex(a_position - 0.5 * m_dx, m_probLo, m_dx);
 
    CH_assert(m_region.contains(particleIV));
 
    const FArrayBox& meshData = a_meshData.getFArrayBox();
 
    for (int comp = 0; comp < a_numComp; comp++) {
      a_particleField[comp] = 0.0;
    }
 
    if ((m_ebisbox->isIrregular(particleIV) && a_forceIrregNGP) || !(a_validBox.contains(particleIV))) {
      for (int comp = 0; comp < a_numComp; comp++) {
        a_particleField[comp] = meshData(particleIV, a_comp + comp);
      }
    }
    else if (!(m_ebisbox->isCovered(particleIV))) {
      auto cicKernel = [&](const IntVect& iv) -> void {
        Real weight = 1.0;
 
        const RealVect L = (m_probLo - a_position) / m_dx + (RealVect(iv) + 0.5 * RealVect::Unit);
 
        for (int dir = 0; dir < SpaceDim; dir++) {
          weight *= (1. - std::abs(L[dir]));
        }
 
        for (int comp = 0; comp < a_numComp; comp++) {
          a_particleField[comp] += weight * meshData(iv, a_comp + comp);
        }
      };
 
      const Box gatherBox = a_gatherBox + loIndex;
 
      BoxLoops::loop<D_DECL(1, 1, 1)>(gatherBox, cicKernel);
    }
  }
 
  inline void
  interpolateParticleTSC(Real*            a_particleField,
                         const EBCellFAB& a_meshData,
                         const int        a_comp,
                         const Box&       a_validBox,
                         const Box&       a_gatherBox,
                         const RealVect&  a_position,
                         const int&       a_numComp,
                         const bool       a_forceIrregNGP) const noexcept
  {
    const IntVect particleIV = ParticleOps::getParticleCellIndex(a_position, m_probLo, m_dx);
    const IntVect loIndex    = ParticleOps::getParticleCellIndex(a_position - m_dx, m_probLo, m_dx);
 
    CH_assert(m_region.contains(particleIV));
 
    const FArrayBox& meshData = a_meshData.getFArrayBox();
 
    for (int comp = 0; comp < a_numComp; comp++) {
      a_particleField[comp] = 0.0;
    }
 
    if ((m_ebisbox->isIrregular(particleIV) && a_forceIrregNGP) || !(a_validBox.contains(particleIV))) {
      for (int comp = 0; comp < a_numComp; comp++) {
        a_particleField[comp] = meshData(particleIV, a_comp + comp);
      }
    }
    else if (!(m_ebisbox->isCovered(particleIV))) {
      auto tscKernel = [&](const IntVect& iv) -> void {
        Real weight = 1.0;
 
        const RealVect L = (m_probLo - a_position) / m_dx + (RealVect(iv) + 0.5 * RealVect::Unit);
 
        for (int dir = 0; dir < SpaceDim; dir++) {
          const Real& l = std::abs(L[dir]);
 
          if (l < 0.5) {
            weight *= 0.75 - l * l;
          }
          else {
            weight *= 0.5 * (1.5 - l) * (1.5 - l);
          }
        }
 
        for (int comp = 0; comp < a_numComp; comp++) {
          a_particleField[comp] += weight * meshData(iv, a_comp + comp);
        }
      };
 
      const Box gatherBox = a_gatherBox + loIndex;
 
      BoxLoops::loop<D_DECL(1, 1, 1)>(gatherBox, tscKernel);
    }
  }
};
 
#include <CD_NamespaceFooter.H>
 
#endif