Skip to content

Commit d1eecf6

Browse files
debogclaude
debog
and
claude
committed
Move particle-to-mesh deposition to ERFPC base class, remove duplicate m_zlevels_d
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
1 parent 3d3fa5c commit d1eecf6

6 files changed

Lines changed: 145 additions & 154 deletions

File tree

‎Source/Particles/ERFPC.H‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -333,6 +333,17 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps,
333333
// the following functions should ideally be private or protected, but need to be
334334
// public due to CUDA extended lambda capture rules
335335

336+
/*! \brief Helper for ParticleToMesh operations with non-uniform z support
337+
*
338+
* Generic CIC deposition from particles to a MultiFab. The caller provides
339+
* a device lambda `value_func(ptd, i)` that returns the per-particle scalar
340+
* to deposit (e.g. mass, number, flux). Handles both uniform and non-uniform
341+
* vertical grids. Implementation in ERFPCParticleToMesh.H.
342+
*/
343+
template<typename ValueFunc>
344+
void ERFPCParticleToMesh(amrex::MultiFab& a_mf, int a_lev, int a_comp,
345+
ValueFunc&& value_func) const;
346+
336347
/*! Default particle initialization */
337348
void initializeParticlesUniformDistributionInBox (const std::unique_ptr<amrex::MultiFab>& a_ptr,
338349
const amrex::RealBox& particle_box);

‎Source/Particles/ERFPCParticleToMesh.H‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
/*! \file ERFPCParticleToMesh.H
2+
* \brief Template implementation of ERFPC::ERFPCParticleToMesh
3+
*
4+
* Separated from ERFPC.H to keep the header lightweight.
5+
* Include this in .cpp files that call ERFPCParticleToMesh().
6+
*/
7+
8+
#ifndef ERFPC_PARTICLE_TO_MESH_H_
9+
#define ERFPC_PARTICLE_TO_MESH_H_
10+
11+
#include <ERFPC.H>
12+
#include <AMReX_ParticleInterpolators.H>
13+
14+
template<typename ValueFunc>
15+
void ERFPC::ERFPCParticleToMesh(amrex::MultiFab& a_mf, int a_lev, int a_comp,
16+
ValueFunc&& value_func) const
17+
{
18+
AMREX_ASSERT(OK());
19+
AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0);
20+
21+
const auto& geom = Geom(a_lev);
22+
const auto plo = geom.ProbLoArray();
23+
const auto dxi = geom.InvCellSizeArray();
24+
25+
const amrex::Real* zlevels = m_zlevels_d.empty() ? nullptr : m_zlevels_d.data();
26+
const int nz_levels = static_cast<int>(m_zlevels_d.size());
27+
28+
a_mf.setVal(0.0);
29+
30+
amrex::ParticleToMesh(*this, a_mf, a_lev,
31+
[=] AMREX_GPU_DEVICE (const typename ERFPC::ParticleTileType::ConstParticleTileDataType& ptd,
32+
int i, amrex::Array4<amrex::Real> const& rho)
33+
{
34+
auto p = ptd.m_aos[i];
35+
auto particle_value = value_func(ptd, i);
36+
37+
if (zlevels) {
38+
amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + amrex::Real(0.5);
39+
int ix = static_cast<int>(amrex::Math::floor(lx)) - 1;
40+
amrex::Real wx1 = lx - static_cast<amrex::Real>(ix + 1);
41+
amrex::Real wx0 = amrex::Real(1.0) - wx1;
42+
43+
amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + amrex::Real(0.5);
44+
int iy = static_cast<int>(amrex::Math::floor(ly)) - 1;
45+
amrex::Real wy1 = ly - static_cast<amrex::Real>(iy + 1);
46+
amrex::Real wy0 = amrex::Real(1.0) - wy1;
47+
48+
int kz = p.idata(ERFParticlesIntIdxAoS::k);
49+
int nz = nz_levels - 1;
50+
51+
amrex::Real dz_k = zlevels[kz+1] - zlevels[kz];
52+
amrex::Real z_center_k = amrex::Real(0.5) * (zlevels[kz] + zlevels[kz+1]);
53+
54+
int kz_lo, kz_hi;
55+
amrex::Real wz_lo, wz_hi;
56+
amrex::Real dz_lo, dz_hi;
57+
58+
if (p.pos(2) >= z_center_k) {
59+
kz_lo = kz;
60+
kz_hi = kz + 1;
61+
dz_lo = dz_k;
62+
if (kz_hi < nz) {
63+
dz_hi = zlevels[kz_hi+1] - zlevels[kz_hi];
64+
amrex::Real z_center_hi = amrex::Real(0.5) * (zlevels[kz_hi] + zlevels[kz_hi+1]);
65+
amrex::Real span = z_center_hi - z_center_k;
66+
wz_hi = (p.pos(2) - z_center_k) / span;
67+
wz_lo = amrex::Real(1.0) - wz_hi;
68+
} else {
69+
kz_hi = kz;
70+
dz_hi = dz_k;
71+
wz_lo = amrex::Real(1.0);
72+
wz_hi = amrex::Real(0.0);
73+
}
74+
} else {
75+
kz_lo = kz - 1;
76+
kz_hi = kz;
77+
dz_hi = dz_k;
78+
if (kz_lo >= 0) {
79+
dz_lo = zlevels[kz_lo+1] - zlevels[kz_lo];
80+
amrex::Real z_center_lo = amrex::Real(0.5) * (zlevels[kz_lo] + zlevels[kz_lo+1]);
81+
amrex::Real span = z_center_k - z_center_lo;
82+
wz_lo = (z_center_k - p.pos(2)) / span;
83+
wz_hi = amrex::Real(1.0) - wz_lo;
84+
} else {
85+
kz_lo = kz;
86+
dz_lo = dz_k;
87+
wz_lo = amrex::Real(0.0);
88+
wz_hi = amrex::Real(1.0);
89+
}
90+
}
91+
92+
amrex::Real inv_vol_lo = dxi[0] * dxi[1] / dz_lo;
93+
amrex::Real inv_vol_hi = dxi[0] * dxi[1] / dz_hi;
94+
95+
amrex::Real pval_lo = particle_value * wz_lo * inv_vol_lo;
96+
amrex::Real pval_hi = particle_value * wz_hi * inv_vol_hi;
97+
98+
amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy , kz_lo, a_comp), wx0*wy0*pval_lo);
99+
amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy , kz_lo, a_comp), wx1*wy0*pval_lo);
100+
amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy+1, kz_lo, a_comp), wx0*wy1*pval_lo);
101+
amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy+1, kz_lo, a_comp), wx1*wy1*pval_lo);
102+
amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy , kz_hi, a_comp), wx0*wy0*pval_hi);
103+
amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy , kz_hi, a_comp), wx1*wy0*pval_hi);
104+
amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy+1, kz_hi, a_comp), wx0*wy1*pval_hi);
105+
amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy+1, kz_hi, a_comp), wx1*wy1*pval_hi);
106+
} else {
107+
amrex::ParticleInterpolator::Linear interp(p, plo, dxi);
108+
amrex::ParticleReal inv_cell_volume = dxi[0] * dxi[1] * dxi[2];
109+
interp.ParticleToMesh(p, rho, 0, a_comp, 1,
110+
[=] AMREX_GPU_DEVICE (const ERFPC::ParticleType&, int) {
111+
return particle_value * inv_cell_volume;
112+
});
113+
}
114+
});
115+
}
116+
117+
#endif

‎Source/Particles/ERFPCUtils.cpp‎

Lines changed: 4 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#include <ERFPC.H>
2+
#include <ERFPCParticleToMesh.H>
23
#include <ERF_Constants.H>
3-
#include <AMReX_ParticleInterpolators.H>
44
#include <AMReX_ParticleLocator.H>
55

66
#ifdef ERF_USE_PARTICLES
@@ -21,32 +21,10 @@ void ERFPC::massDensity ( MultiFab& a_mf,
2121
const int& a_comp ) const
2222
{
2323
BL_PROFILE("ERFPC::massDensity()");
24-
25-
AMREX_ASSERT(OK());
26-
AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0);
27-
28-
const auto& geom = Geom(a_lev);
29-
const auto plo = geom.ProbLoArray();
30-
const auto dxi = geom.InvCellSizeArray();
31-
32-
const Real inv_cell_volume = dxi[0]*dxi[1]*dxi[2];
33-
a_mf.setVal(0.0);
34-
35-
ParticleToMesh( *this, a_mf, a_lev,
36-
[=] AMREX_GPU_DEVICE ( const ERFPC::ParticleTileType::ConstParticleTileDataType& ptr,
37-
int i, Array4<Real> const& rho)
38-
{
39-
auto p = ptr.m_aos[i];
40-
ParticleInterpolator::Linear interp(p, plo, dxi);
41-
interp.ParticleToMesh ( p, rho, 0, a_comp, 1,
42-
[=] AMREX_GPU_DEVICE ( const ERFPC::ParticleType&, int)
43-
{
44-
auto mass = ptr.m_rdata[ERFParticlesRealIdxSoA::mass][i];
45-
return mass*inv_cell_volume;
46-
});
24+
ERFPCParticleToMesh(a_mf, a_lev, a_comp,
25+
[=] AMREX_GPU_DEVICE (const ERFPC::ParticleTileType::ConstParticleTileDataType& ptd, int i) {
26+
return ptd.m_rdata[ERFParticlesRealIdxSoA::mass][i];
4727
});
48-
49-
return;
5028
}
5129

5230
/*! \brief Fix k-indices for all particles after AMR regrid */

‎Source/Particles/ERF_SuperDropletPC.H‎

Lines changed: 0 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -276,10 +276,6 @@ class SuperDropletPC : public ERFPC
276276
amrex::Real,
277277
const int a_comp = 0) const;
278278

279-
/*! \brief Helper for ParticleToMesh operations */
280-
template<typename ValueFunc>
281-
void particleToMeshHelper(amrex::MultiFab& a_mf, int a_comp, ValueFunc&& value_func) const;
282-
283279
/*! \brief Compute the species mass flux component on a mesh */
284280
virtual void speciesMassFlux ( amrex::MultiFab&, const int, const int, const int a_comp=0 ) const;
285281

@@ -375,11 +371,6 @@ class SuperDropletPC : public ERFPC
375371
/*! save inactive particles to file */
376372
bool m_save_inactive;
377373

378-
/*! Staggered z-levels for non-uniform vertical grid spacing (from terrain_z_levels).
379-
* Size nz+1: z_levels[k] is the bottom face of cell k. Used for CIC interpolation
380-
* weights and cell volume computation in particle-to-mesh deposition. */
381-
amrex::Gpu::DeviceVector<amrex::Real> m_zlevels_d;
382-
383374
/*! Persistent device vectors for material properties */
384375
amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sp_density;
385376
amrex::Gpu::DeviceVector<int> m_sp_solubility;

0 commit comments

Comments
 (0)