ERFPCParticleToMesh.H

Go to the documentation of this file.

/*! \file ERFPCParticleToMesh.H
 *  \brief Template implementation of ERFPC::ERFPCParticleToMesh
 *
 *  Separated from ERFPC.H to keep the header lightweight.
 *  Include this in .cpp files that call ERFPCParticleToMesh().
 */
 
#ifndef ERFPC_PARTICLE_TO_MESH_H_
#define ERFPC_PARTICLE_TO_MESH_H_
 
#include <ERFPC.H>
#include <ERF_TerrainConversion.H>
#include <AMReX_MultiFabUtil.H>
 
template<typename ValueFunc>
void ERFPC::ERFPCParticleToMesh (amrex::MultiFab& a_mf,
                                 const amrex::MultiFab& a_z_phys_nd,
                                 int a_lev, int a_comp,
                                 ValueFunc&& value_func,
                                 bool composite) const
{
    using namespace amrex;
 
    AMREX_ASSERT(OK());
    AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0);
 
    // CIC deposition writes one cell outside the particle's own cell
    AMREX_ALWAYS_ASSERT(a_mf.nGrowVect().allGE(amrex::IntVect(1)));
 
    const auto& geom = Geom(a_lev);
    const auto plo = geom.ProbLoArray();
    const auto dxi = geom.InvCellSizeArray();
    const Real dx  = geom.CellSize(0);
    const Real dy  = geom.CellSize(1);
    const Box  domain = geom.Domain();
    const int  k_max  = domain.bigEnd(AMREX_SPACEDIM-1);
    const int  i_lo   = domain.smallEnd(0), i_hi = domain.bigEnd(0);
    const int  j_lo   = domain.smallEnd(1), j_hi = domain.bigEnd(1);
    const int  k_lo   = domain.smallEnd(AMREX_SPACEDIM-1);
    const bool per_x  = geom.isPeriodic(0);
    const bool per_y  = geom.isPeriodic(1);
 
    iMultiFab owned;
    if (composite) {
        owned.define(a_mf.boxArray(), a_mf.DistributionMap(), 1, IntVect(1));
        owned.setVal(0);
        if (a_lev < finestLevel()) {
            iMultiFab exposed = makeFineMask(a_mf, ParticleBoxArray(a_lev+1),
                                             GetParGDB()->refRatio(a_lev), 1, 0);
            iMultiFab::Copy(owned, exposed, 0, 0, 1, 0);
        } else {
            owned.setVal(1, 0, 1, 0);
        }
        owned.FillBoundary(geom.periodicity());
    }
 
    a_mf.setVal(0.0);
 
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (ParConstIterType pti(*this, a_lev); pti.isValid(); ++pti)
    {
        const int grid = pti.index();
        const auto& ptile = ParticlesAt(a_lev, pti);
        const auto& aos = ptile.GetArrayOfStructs();
        const int np = aos.numParticles();
        if (np == 0) { continue; }
 
        auto rho     = a_mf[grid].array();
        auto zheight = a_z_phys_nd[grid].array();
        const Array4<const int> valid = composite ? owned.const_array(grid) : Array4<const int>{};
 
        const auto ptd = ptile.getConstParticleTileData();
 
        ParallelFor(np, [=] AMREX_GPU_DEVICE (int i)
        {
            const auto& p = ptd.m_aos[i];
            if (p.id() <= 0) { return; }
 
            const Real pval = static_cast<Real>(value_func(ptd, i));
 
            const Real lx = static_cast<Real>((p.pos(0) - plo[0]) * dxi[0]) + Real(0.5);
            const int  ix = static_cast<int>(Math::floor(lx)) - 1;
            const Real wx1 = lx - static_cast<Real>(ix + 1);
            const Real wx0 = Real(1.0) - wx1;
 
            const Real ly = static_cast<Real>((p.pos(1) - plo[1]) * dxi[1]) + Real(0.5);
            const int  iy = static_cast<int>(Math::floor(ly)) - 1;
            const Real wy1 = ly - static_cast<Real>(iy + 1);
            const Real wy0 = Real(1.0) - wy1;
 
            const int kz = static_cast<int>(amrex::Math::floor((p.pos(AMREX_SPACEDIM-1) - plo[AMREX_SPACEDIM-1])
                                                               * dxi[AMREX_SPACEDIM-1]));
            const int ix_p = static_cast<int>(Math::floor((p.pos(0) - plo[0]) * dxi[0]));
            const int iy_p = static_cast<int>(Math::floor((p.pos(1) - plo[1]) * dxi[1]));
            const Real fx  = static_cast<Real>((p.pos(0) - plo[0]) * dxi[0]) - static_cast<Real>(ix_p);
            const Real fy  = static_cast<Real>((p.pos(1) - plo[1]) * dxi[1]) - static_cast<Real>(iy_p);
 
            auto zface = [=] AMREX_GPU_DEVICE (int k_face) -> Real {
                return zheight(ix_p  , iy_p  , k_face) * (Real(1.0)-fx) * (Real(1.0)-fy)
                     + zheight(ix_p+1, iy_p  , k_face) *            fx  * (Real(1.0)-fy)
                     + zheight(ix_p  , iy_p+1, k_face) * (Real(1.0)-fx) *            fy
                     + zheight(ix_p+1, iy_p+1, k_face) *            fx  *            fy;
            };
 
            // Effective k bounds for safe zface access: must respect both the
            // geometry domain (k_max) and the local fab's z extent. With AMR
            // partial-z fabs (especially with many small boxes from large rank
            // counts), a particle at the top/bottom of its fab can have kz
            // such that kz+2 / kz-1 read outside the fab. Use the fallback
            // (one-sided) span formula in those cases.
            const int kz_hi_safe = amrex::min(k_max, zheight.end[2] - 2);
            const int kz_lo_safe = amrex::max(0,     zheight.begin[2] + 1);
 
            const Real z_lo_k = zface(kz);
            const Real z_hi_k = zface(kz + 1);
            const Real z_ctr_k = Real(0.5) * (z_lo_k + z_hi_k);
 
            // Particle pos(2) is computational zeta; convert to physical z so
            // the cell-center comparison and weight span are in the same units
            // as zface() (which is built from the height array).
            const Real p_z = static_cast<Real>(ERF::ParticlePos::z_from_zeta(
                p.pos(0), p.pos(1), p.pos(AMREX_SPACEDIM-1),
                plo, dxi, zheight));
 
            int  kz_lo, kz_hi;
            Real wz_lo, wz_hi;
            if (p_z >= z_ctr_k) {
                kz_lo = kz;
                kz_hi = kz + 1;
                Real span;
                if (kz < kz_hi_safe) {
                    const Real z_hi2    = zface(kz + 2);
                    const Real z_ctr_hi = Real(0.5) * (z_hi_k + z_hi2);
                    span = z_ctr_hi - z_ctr_k;
                } else {
                    span = z_hi_k - z_lo_k;
                }
                wz_hi = (p_z - z_ctr_k) / span;
                wz_lo = Real(1.0) - wz_hi;
            } else {
                kz_hi = kz;
                kz_lo = kz - 1;
                Real span;
                if (kz > kz_lo_safe) {
                    const Real z_lo2    = zface(kz - 1);
                    const Real z_ctr_lo = Real(0.5) * (z_lo2 + z_lo_k);
                    span = z_ctr_k - z_ctr_lo;
                } else {
                    span = z_hi_k - z_lo_k;
                }
                wz_lo = (z_ctr_k - p_z) / span;
                wz_hi = Real(1.0) - wz_lo;
            }
 
            // Per-corner node-averaged cell thickness for volume normalization.
            auto dz_cell = [=] AMREX_GPU_DEVICE (int ic, int jc, int kc) -> Real {
                const Real dz = Real(0.25) * (
                      zheight(ic  , jc  , kc+1) - zheight(ic  , jc  , kc)
                    + zheight(ic+1, jc  , kc+1) - zheight(ic+1, jc  , kc)
                    + zheight(ic  , jc+1, kc+1) - zheight(ic  , jc+1, kc)
                    + zheight(ic+1, jc+1, kc+1) - zheight(ic+1, jc+1, kc) );
                return dz;
            };
 
            const Real inv_dxdy = Real(1.0) / (dx * dy);
 
            // Restrict appended-species CIC to the owning level's composite
            // cells. Renormalization preserves each particle's mass at walls
            // and coarse/fine edges without creating another grid reservoir.
            Real weight = Real(1);
            if (composite) {
                weight = Real(0);
                for (int iz = 0; iz < 2; ++iz) {
                    const int kc = iz ? kz_hi : kz_lo;
                    const Real wz = iz ? wz_hi : wz_lo;
                    for (int jy = 0; jy < 2; ++jy) {
                        for (int jx = 0; jx < 2; ++jx) {
                            if (valid(ix+jx, iy+jy, kc)) {
                                weight += (jx ? wx1 : wx0) * (jy ? wy1 : wy0) * wz;
                            }
                        }
                    }
                }
                AMREX_ALWAYS_ASSERT(weight > Real(0));
            }
 
            auto deposit = [=] AMREX_GPU_DEVICE (int ic, int jc, int kc, Real w_horiz, Real w_z)
            {
                // A corner with zero CIC weight deposits nothing.  Skip it before
                // touching dz_cell: when a particle sits exactly on a cell-center
                // plane the outermost corner lands in the ghost region, where the
                // node-averaged thickness can come out exactly zero (z_phys_nd
                // ghost nodes outside the domain are not filled), and 0/0 raises
                // an invalid-operation FPE.
                if (w_horiz == Real(0.0) || w_z == Real(0.0)) { return; }
                if (composite && !valid(ic,jc,kc)) { return; }
 
                // Clamp only for the dz lookup; still write to the ghost rho cell. In a
                // periodic direction the exterior nodes hold the wrapped image, so the
                // unclamped stencil already gives the thickness of the cell SumBoundary
                // folds into; only the non-periodic exterior nodes need clamping.
                const int ic_dz = per_x ? ic : amrex::max(i_lo, amrex::min(i_hi, ic));
                const int jc_dz = per_y ? jc : amrex::max(j_lo, amrex::min(j_hi, jc));
                const int kc_dz = amrex::max(k_lo, amrex::min(k_max, kc));
                const Real dz = dz_cell(ic_dz, jc_dz, kc_dz);
                AMREX_ALWAYS_ASSERT(dz > Real(0.0));
                const Real contrib = pval * w_horiz * w_z * inv_dxdy / (dz * weight);
                Gpu::Atomic::AddNoRet(&rho(ic, jc, kc, a_comp), contrib);
            };
 
            deposit(ix  , iy  , kz_lo, wx0*wy0, wz_lo);
            deposit(ix+1, iy  , kz_lo, wx1*wy0, wz_lo);
            deposit(ix  , iy+1, kz_lo, wx0*wy1, wz_lo);
            deposit(ix+1, iy+1, kz_lo, wx1*wy1, wz_lo);
            deposit(ix  , iy  , kz_hi, wx0*wy0, wz_hi);
            deposit(ix+1, iy  , kz_hi, wx1*wy0, wz_hi);
            deposit(ix  , iy+1, kz_hi, wx0*wy1, wz_hi);
            deposit(ix+1, iy+1, kz_hi, wx1*wy1, wz_hi);
        });
    }
 
    a_mf.SumBoundary(geom.periodicity());
}
 
#endif