Undefined tmf in post_timestep_particle; ParticleDerive reads device data on host

[WeiqunZhang]

Severity: high/medium · Category: gpu-parallel, memory-ub · Fix order: 4 of 21 — fix this 4th.

_{Filenames are numbered in reverse fix order: 001 = fix last, 021 = fix first. This file is 018.}

Locations: Source/NavierStokesBase.cpp:3950, Source/NavierStokesBase.cpp:4035, Source/NavierStokesBase.cpp:3921

Based on commit 9bf664bf (line numbers refer to that tree).

The five findings collapse into two edits: the tmf define in post_timestep_particle, and the fine-to-coarse accumulation loop in ParticleDerive("total_particle_count").

The defect

Source/NavierStokesBase.cpp:3950 — post_timestep_particle passes a default-constructed (never defined) MultiFab tmf to AmrTracerParticleContainer::Timestamp whenever particles.timestamp_indices is not supplied, because tmf.define() is guarded by if (tindices.size() > 0) while the Timestamp call is not.

Source/NavierStokesBase.cpp:4035 — ParticleDerive("total_particle_count") iterates cells on the HOST via BoxIterator, dereferencing ffab(p)/cfab(...) whose data live in device memory on GPU builds (temp_dat/ctemp_dat use the default device arena), the exact class fixed in 92f0e72/#147. _{Reported independently at this same line by F029; each reviewer's own wording and evidence is under Verification evidence below.}

Source/NavierStokesBase.cpp:3921 — tmf is built on level lev's BoxArray/DistributionMap but is handed Factory(), the FabFactory of the current level level; in EB builds the EB cell-flag FabArray belongs to a different BoxArray/DM, so EBFArrayBoxFactory::create indexes it with foreign box indices. _{Reported independently at this same line by F050; each reviewer's own wording and evidence is under Verification evidence below.}

Why it matters

F005: AMREX_PARTICLES build, defaults do_nspc=true and timestamp_dir="Timestamps", user sets particles.particle_init_file but no particles.timestamp_indices. Timestamp then executes ba[grid] on an empty BoxArray and mf[grid] with localindex()==-1 (AMReX_TracerParticles.cpp:278-282) -> out-of-bounds read / segfault on the first step with particles.

F011: GPU build (USE_CUDA/HIP) with AMREX_PARTICLES, tracer particles active, multi-level run (finest_level>level), and amr.derive_plot_vars containing total_particle_count: writing a plotfile executes the host loop at lines 4033-4038 and segfaults because the_arena_is_managed=false by default. Additionally an OpenMP CPU race exists when the refinement-ratio product exceeds the tile size.

F028: AMREX_USE_EB + AMREX_PARTICLES, >=2 AMR levels, particles.timestamp_indices set. For lev>level, EBFArrayBoxFactory::create does getMultiEBCellFlagFab()[box_index] with an index from lev's BoxArray -> assertion ebcellflag.box().contains(enclosedCells(bx)) fires in debug, or m_fabs_v[-1] UB / wrong flags in release.

How to reach it

• Build Exec/run_2d_particles (ships USE_PARTICLES=TRUE); run its regtest.inputs with the particles.timestamp_indices line deleted, keeping particles.particle_init_file=particle_file. First post_timestep calls Timestamp with the undefined tmf: segfault (release) or AMREX_ASSERT abort (debug).
• Exec/run_2d_particles/regtest.inputs (amr.max_level=1, plot_int=50, particles.particle_init_file) built with USE_PARTICLES=TRUE USE_CUDA=TRUE, plus amr.derive_plot_vars=total_particle_count; first plotfile write segfaults at level 0.

Suggested fix

In post_timestep_particle, hoist tmf.define out of if (tindices.size() > 0) and pass amr_level.Factory() — the factory of the level being timestamped, not this level's. Prefer that to skipping Timestamp: particles.timestamp_dir alone is a request for trajectory output, and AMReX's Timestamp still writes id/position/time/velocity, yet it unconditionally evaluates ba[grid]/mf[grid] (dtoh-copying that fab on GPU), so it needs a real BoxArray and one allocated component. In ParticleDerive, make the same factory fix at line 4009 and replace the BoxIterator loop with ParallelFor + Gpu::Atomic::AddNoRet; ParticleContainer::Increment fills temp_dat from a device kernel, and this is the bug class already fixed for Projection::computeRhoG in 92f0e72 (#147). Drop the assert at 4031 in the same change — it compares a coarsened grid box against a coarsened tilebox and fires under tiling in debug. Maintainer's call: whether temp_dat needs an EB factory at all, since ctemp_dat already uses the default.

For Source/NavierStokesBase.cpp:3950 (F005):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -3946,6 +3946,15 @@
                         timestamp_add_extras(lev, curr_time, tmf);
                     }
                 }
+                else
+                {
+                    //
+                    // Timestamp indexes tmf's BoxArray/DistributionMap even
+                    // when tindices is empty, so tmf must always be defined.
+                    //
+                    tmf.define(S_new.boxArray(), S_new.DistributionMap(), 1,
+                               ng, MFInfo(), amr_level.Factory());
+                }
 
                 NSPC->Timestamp(basename, tmf, lev, curr_time, tindices);
             }

_{AmrTracerParticleContainer::Timestamp (amrex AMReX_TracerParticles.cpp) does ba[grid]/mf[grid] before checking indices, so tmf must be defined even with empty tindices; a 1-component MultiFab (never read) preserves position-only timestamp output. Uses lev's factory (amr_level.Factory()). Adds an else branch, so it does not overlap F028's hunk at the existing define.}

For Source/NavierStokesBase.cpp:4035 (F011):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -4019,23 +4019,25 @@
 
                 NSPC->Increment(temp_dat,lev);
 
 #ifdef _OPENMP
-#pragma omp parallel
+#pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
-                for (MFIter mfi(temp_dat,true); mfi.isValid(); ++mfi)
+                for (MFIter mfi(temp_dat,TilingIfNotGPU()); mfi.isValid(); ++mfi)
                 {
-                    const FArrayBox& ffab =  temp_dat[mfi];
-                    FArrayBox&       cfab = ctemp_dat[mfi];
-                    const Box&       fbx  = mfi.tilebox();
-
-                    AMREX_ASSERT(cfab.box() == amrex::coarsen(fbx,trr));
+                    const Box& fbx = mfi.tilebox();
+                    auto const& ffab = temp_dat.const_array(mfi);
+                    auto const& cfab = ctemp_dat.array(mfi);
+                    const IntVect ratio = trr;
 
-                    for (IntVect p = fbx.smallEnd(); p <= fbx.bigEnd(); fbx.next(p))
+                    amrex::ParallelFor(fbx, [ffab, cfab, ratio]
+                    AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
-                        const Real val = ffab(p);
-                        if (val > 0)
-                            cfab(amrex::coarsen(p,trr)) += val;
-                    }
+                        const Real val = ffab(i,j,k);
+                        if (val > 0) {
+                            const Dim3 cp = amrex::coarsen(Dim3{i,j,k}, ratio);
+                            Gpu::Atomic::AddNoRet(&cfab(cp.x,cp.y,cp.z), val);
+                        }
+                    });
                 }
 
                 temp_dat.clear();

_{Replaces the host BoxIterator loop with amrex::ParallelFor over the fine tilebox using Gpu::Atomic::AddNoRet and amrex::coarsen(Dim3,IntVect) (host/device, AMReX_IntVect.H). Atomics also fix the OpenMP tile race; the tiling-invalid AMREX_ASSERT is dropped. OMP pragma gains if(Gpu::notInLaunchRegion()) and MFIter uses TilingIfNotGPU(), matching file idiom.}

For Source/NavierStokesBase.cpp:3921 (F028):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -3918,7 +3918,7 @@
 
                 if (tindices.size() > 0)
                 {
-                    tmf.define(S_new.boxArray(), S_new.DistributionMap(), tindices.size(), ng, MFInfo(), Factory());
+                    tmf.define(S_new.boxArray(), S_new.DistributionMap(), tindices.size(), ng, MFInfo(), amr_level.Factory());
 
                     if (n > 0)
                     {
@@ -4006,6 +4006,6 @@
             {
                 BoxArray ba = parent->boxArray(lev);
 
-                MultiFab temp_dat(ba,parent->DistributionMap(lev),1,0,MFInfo(),Factory());
+                MultiFab temp_dat(ba,parent->DistributionMap(lev),1,0,MFInfo(),parent->getLevel(lev).Factory());
 
                 trr *= parent->refRatio(lev-1);

_{Uses the target level's FabFactory: amr_level is already parent->getLevel(lev) in post_timestep_particle; ParticleDerive's total_particle_count temp_dat (second hunk) has the identical defect per EBFArrayBoxFactory::create indexing getMultiEBCellFlagFab() by box index (amrex EBFabFactory.cpp). If F005's diff is applied first, this first hunk needs a trivial context refresh.}

_{Diff(s) are against 9bf664bf, written from the current source and verified only with git apply --check — never compiled, never run, never applied to the tree. Treat them as precise intent, not tested patches.}

Verification evidence

F005 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): timestamp_dir defaults to "Timestamps" (line 213), so block always entered; tmf.define guarded by 'if (tindices.size() > 0)' (3919) but 'NSPC->Timestamp(basename, tmf, lev, curr_time, tindices)' (3950) is not. With no particles.timestamp_indices and timestamp_num_extras()==0 (NavierStokesBase.H:509), tmf is undefined. amrex TracerParticles.cpp does 'const Box& bx = ba[grid]; const FArrayBox& fab = mf[grid]' on the empty BoxArray/undefined MultiFab before any indices check -> OOB/UB when a rank has particles.

Lens 2 (reachability/intent): NavierStokesBase.cpp:3917-3950: tmf.define() guarded by tindices.size()>0 but NSPC->Timestamp(basename,tmf,...) unconditional; timestamp_dir defaults "Timestamps" (line 213), do_nspc defaults true (199), timestamp_num_extras()==0. AMReX_TracerParticles.cpp:281-282 executes ba[grid] on the empty BoxArray and mf[grid] -> fabPtr -> m_fabs_v[localindex(K)==-1] (AMReX_FabArray.H:2074-2078) before the if(M>0) guard, per grid with particles. Structure unchanged since commit 2ba529d (2016); no deliberate-asymmetry evidence.

F011 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): Lines 4033-4038: host loop 'for (IntVect p...; fbx.next(p)) { const Real val = ffab(p); ... cfab(coarsen(p,trr)) += val; }' on temp_dat/ctemp_dat allocated in The_Arena (amrex Arena.cpp:59 'the_arena_is_managed = false' default) -> host deref of device memory on GPU builds. Reachable: derive 'total_particle_count' registered (NS_setup.cpp:469) and routed via NavierStokes::derive->ParticleDerive. Same class as commit 92f0e72 (#147) which fixed only Projection.cpp.

Lens 2 (reachability/intent): NavierStokesBase.cpp:4009/4015 allocate temp_dat/ctemp_dat with default MFInfo() (The_Arena; amrex Arena.cpp:59 the_arena_is_managed=false), then lines 4033-4038 dereference ffab(p)/cfab on host. AMReX Increment is GPU-safe (ParticleToMesh), so the host loop is reached. Commit 92f0e72 (#147) fixed this exact class in Projection.cpp via pinned staging; GPU-port db271ac (#146) never touched ParticleDerive — legacy, not deliberate. OMP race real only when trr>8 in y/z (3D).

F028 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 3921: tmf.define(S_new.boxArray(), S_new.DistributionMap(), ..., Factory()) inside 'for (lev = level; lev <= finest_level; lev++)' where S_new is level lev's but Factory() is this level's. EBFArrayBoxFactory::create (amrex EBFabFactory.cpp:107) does m_ebdc->getMultiEBCellFlagFab()[box_index] with box_index from the fine BoxArray; FabArray::fabPtr asserts localindex>=0 in debug, m_fabs_v[-1]/wrong flags in release for lev>level.

F050 — confirmed (one verifier lens)

Reported as: post_timestep_particle defines tmf on the finer level's BoxArray/DistributionMap but passes THIS (coarser) level's Factory(); with EB, EBFArrayBoxFactory::create indexes its own flag FabArray with box indices from the mismatched BoxArray/DM.

Failure scenario: EB build + AMREX_PARTICLES with particles.timestamp_dir set and particles on a finer level: in the loop for (lev=level; lev<=finest_level; lev++), for lev>level tmf.define calls the coarse level's EBFArrayBoxFactory::create(box,...,box_index) where box_index refers to the fine BoxArray; FabArray::operator[] on the coarse flag MF hits a wrong/absent local index -> assertion failure or out-of-bounds access/wrong EB flags.

Lens 1 (refutation attempt): Same confirmed factory/BoxArray mismatch at 3921 (needs particles.timestamp_indices nonempty so the define executes — minor gate imprecision in the scenario, but the defect and EB reachability stand). EBFArrayBoxFactory::create indexes the coarse flag FabArray with a fine-level global box index -> assertion or OOB/wrong flags for lev>level.

F029 — confirmed (one verifier lens)

Reported as: ParticleDerive("total_particle_count") reads and writes MultiFab data directly from the host (ffab(p), cfab(coarsen(p,trr)) += val) although both MultiFabs are allocated in The_Arena(), which is device (non-managed) memory by default in GPU builds.

Failure scenario: CUDA/HIP build (amrex.the_arena_is_managed=0, the default) with AMREX_PARTICLES and amr.derive_plot_vars containing total_particle_count, finestLevel>level: the host dereference of device pointers segfaults when the plotfile is written. The same loop also races under OpenMP when trr exceeds the tile size in y/z.

Lens 1 (refutation attempt): Confirmed: temp_dat (4009, MFInfo()+Factory()) and ctemp_dat (4015, default MFInfo) allocate in The_Arena; amrex Arena.cpp:59 defaults the_arena_is_managed=false, so GPU builds hold device memory. Host loop 4033-4038 dereferences ffab(p)/cfab(...) directly -> segfault at plotfile time. OMP race also plausible: MFIter(temp_dat,true) tiles, and coarsen(p,trr) collides across tiles when trr exceeds tile extent (+= at 4037).

---

_{Based on commit 9bf664bf, which is also the tree the audit verified against. From an automated audit of Source/, Tutorials/ and Util/. Audit finding ids: F005, F011, F028, F050, F029. Reviewer unit(s): NSB-3 syncinterp/vel-advance/particles, theme:ghost-fillpatch, theme:gpu-capture. Nothing here was compiled or run — the failure scenarios are code reasoning, so the reaching configuration above is the cheapest way to confirm or refute it.}