MacProj: known-edgestate sync overload segfaults, sync divu omits dsdt, umac test on host

[WeiqunZhang]

Severity: medium/low · Category: gpu-parallel, memory-ub, physics-numerics · Fix order: 11 of 21 — fix this 11th.

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

Locations: Source/MacProj.cpp:778, Source/MacProj.cpp:562, Source/MacProj.cpp:1068

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

The three defects are independent and can be applied separately; only the first also touches NavierStokesBase::ComputeAofs.

The defect

Source/MacProj.cpp:778 — The known-edgestate mac_sync_compute overload passes a default-constructed temporary MultiFab() as the state, but NavierStokesBase::ComputeAofs unconditionally evaluates S.const_array(mfi, S_comp) (NavierStokesBase.cpp:4665) for every grid, indexing the empty FabArray's internal m_fabs_v vector out of bounds.

Source/MacProj.cpp:562 — mac_sync_compute builds the divu constraint at prev_time only, omitting the 0.5dtdsdt half-time correction that both velocity_advection (NavierStokesBase.cpp:3423-3424) and scalar_advection (NavierStokes.cpp:732-733) apply before Godunov/BDS edge-state prediction, so the sync reconstructs edge states with a different divu than the advance used.

Source/MacProj.cpp:1068 — test_umac_periodic performs host-side FArrayBox operations (mfcd.FillFab into 'diff', diff.minusRunOn::Host, diff.normRunOn::Host) on data allocated in device memory ('diff' uses the default arena, u_mac fabs are device-resident), the #147/#109 host-touches-device class.

Why it matters

F024: Any invocation of this public API (kept and documented in MacProj.cpp:733-739 for derived physics codes that supply precomputed edge states, e.g. div(rho U h) syncs) fails on the first grid: AMREX_ASSERT(mfi.LocalIndex() < indexArray.size()) fires in debug builds; in release builds m_fabs_v[li] on an empty vector yields a garbage/null FAB pointer that is immediately dereferenced, segfaulting. The overload cannot work in any configuration.

F065: Multilevel (max_level>=1, do_reflux) run with advection_scheme Godunov or BDS and nonzero, time-varying divu (e.g. do_temp=1 with diffusive temperature so have_divu and have_dsdt are set): during mac_sync the recomputed edge states differ from the advance's by the O(0.5dtdsdt) source-term difference, so Vsync/Ssync and the sync flux-register increments are systematically inconsistent with the advance at coarse-fine boundaries, violating the design requirement stated in fix 82d3106 (#118) that the sync rebuild edge states exactly as the advance did.

F077: GPU build with a periodic direction and mac_proj.check_umac_periodicity=1 in the inputs (the option is user-settable; only its default is 0 on GPU): after each mac projection, FillFab writes and minusRunOn::Host reads device pointers on the host, segfaulting since the_arena_is_managed=false by default.

Suggested fix

F024: keep the documented overload, but stop materializing an undefined state. AMReX-Hydro's ComputeFluxesOnBoxFromState never reads q when knownFaceState is true (an empty Array4 is the accepted sentinel there, exactly as for divu/fq), so ComputeAofs should build S_arr only where it is used, and shape the EB StateRedist rstate_tmp from advc rather than S.boxArray()/S.nGrow() — both must change together or the EB sync still dereferences the empty state. Having MacProj pass a defined alias instead is the maintainer's call.

F065: mirror the 0.5dtdsdt Saxpy the advance applies, so the sync rebuilds edge states with the advance's divu as #118 requires. The pre-port Fortran had the same asymmetry, so this is a long-standing gap, not a migration regression; getDsdt returns a zeroed MultiFab when have_divu/have_dsdt are unset, so runs without divu are unchanged.

F077: follow 92f0e72 (#147) — pinned FAB plus dtoh_memcpy_async/streamSynchronize before minus/norm; hard-disabling the option on GPU is an acceptable alternative.

For Source/MacProj.cpp:778 (F024):

--- a/Source/MacProj.cpp
+++ b/Source/MacProj.cpp
@@ -770,12 +770,20 @@
     bool do_crse_add = false;
     bool do_fine_add = update_fluxreg;
 
+    //
+    // ComputeAofs always builds an Array4 from the state MultiFab, even though
+    // it does not use the values when the edge states are known. So it must be
+    // given a defined MultiFab; an alias of Sync has the correct BoxArray and
+    // DistributionMap.
+    //
+    MultiFab Smf(Sync, amrex::make_alias, Sync_indx, ncomp);
+
     //
     // Compute the mac sync correction.
     //
     ns_level.ComputeAofs(Sync, /*Ssync_comp*/ Sync_indx, state_comp,
                          ncomp,
-                         /*State*/ MultiFab(), /*S_comp*/ int(),//not used when known_edgestates
+                         /*State*/ Smf, /*S_comp*/ 0,           //not used when known_edgestates
                          /*forcing*/ nullptr, /*f_comp*/ int(), //not used when known_edgestates
                          /*constraint divU*/ nullptr,           //not used when known_edgestates
                          fluxes, /*flux_comp*/ 0,

_{Aliases Sync (AMReX MultiFab alias ctor) so ComputeAofs's unconditional S.const_array(mfi,S_comp) is valid; values are unused once knownFaceState is true (hydro_compute_edgestate_and_flux.cpp:428). In EB builds this also gives the non-StateRedist rstate=&S path the sync data, matching the StateRedist branch's documented choice. Alternative: guard S_arr inside ComputeAofs.}

For Source/MacProj.cpp:562 (F065):

--- a/Source/MacProj.cpp
+++ b/Source/MacProj.cpp
@@ -561,6 +561,13 @@
         forcing_term = std::make_unique<MultiFab>(grids, dmap, num_state_comps, NavierStokesBase::nghost_force());
         divu_fp.reset(ns_level.getDivCond(NavierStokesBase::nghost_force(),prev_time));
 
+        // Get divu to time n+1/2, exactly as velocity_advection and
+        // scalar_advection do before predicting the edge states.
+        {
+            std::unique_ptr<MultiFab> dsdt(ns_level.getDsdt(NavierStokesBase::nghost_force(),prev_time));
+            MultiFab::Saxpy(*divu_fp, 0.5*dt, *dsdt, 0, 0, 1, NavierStokesBase::nghost_force());
+        }
+
         MultiFab& Gp = ns_level.get_old_data(Gradp_Type);
 
         visc_terms.setVal(0.0); // Initialize to make calls below safe

_{Copies the half-time divu update from the twins at NavierStokesBase.cpp:3423 and NavierStokes.cpp:732. getDsdt (public, non-const ns_level) returns a zeroed MF unless have_dsdt&&have_divu, so runs without dsdt are unchanged. dt here is parent->dtLevel(level) (NavierStokes.cpp:1454), the same dt the advance used.}

For Source/MacProj.cpp:1068 (F077):

--- a/Source/MacProj.cpp
+++ b/Source/MacProj.cpp
@@ -989,13 +989,25 @@
     std::vector< std::pair<int,Box> > isects;
 
 
+    //
+    // MultiFabCopyDescriptor and the FArrayBox comparison below both work on
+    // the host, so they must be given host-accessible copies of u_mac.
+    //
+    Array<MultiFab,AMREX_SPACEDIM> h_umac;
+
     for (int dim = 0; dim < AMREX_SPACEDIM; dim++)
     {
         if (geom.isPeriodic(dim))
         {
             Box eDomain = amrex::surroundingNodes(geom.Domain(),dim);
 
-            mfid[dim] = mfcd.RegisterMultiFab(&u_mac[dim]);
+            h_umac[dim].define(u_mac[dim].boxArray(), u_mac[dim].DistributionMap(),
+                               u_mac[dim].nComp(), u_mac[dim].nGrowVect(),
+                               MFInfo().SetArena(The_Pinned_Arena()));
+            amrex::dtoh_memcpy(h_umac[dim], u_mac[dim]);
+            Gpu::streamSynchronize();
+
+            mfid[dim] = mfcd.RegisterMultiFab(&h_umac[dim]);
 
             // How to combine pirm into one global pirm?
             // don't think std::vector::push_back() is thread safe
@@ -1061,10 +1073,10 @@
         AMREX_ASSERT(pirm_i.m_srcBox.sameSize(pirm_i.m_dstBox));
         AMREX_ASSERT(u_mac[dim].DistributionMap()[pirm_i.m_idx] == ParallelDescriptor::MyProc());
 
-        diff.resize(pirm_i.m_srcBox, 1);
+        diff.resize(pirm_i.m_srcBox, 1, The_Pinned_Arena());
 
         mfcd.FillFab(mfid[dim], pirm_i.m_fbid, diff);
 
-        diff.minus<RunOn::Host>(u_mac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);
+        diff.minus<RunOn::Host>(h_umac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);
 
         const Real max_norm = diff.norm<RunOn::Host>(0);

_{Registers pinned host copies with the descriptor (so FillFab/copyToMem read host memory) and puts diff in The_Pinned_Arena. amrex::dtoh_memcpy (AMReX_FabArrayUtility.H:1722) falls back to Copy on CPU builds. Caveat: for remote fills AMReX's FabArrayCopyDescriptor still allocates its cache FAB and MPI buffers from The_Arena, so GPU+MPI needs an AMReX-side fix too.}

_{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

F024 — confirmed (one verifier lens)

Lens 1 (refutation attempt): MacProj.cpp:778 passes /*State*/ MultiFab() to ComputeAofs; NavierStokesBase.cpp:4665 runs const auto& S_arr = S.const_array(mfi, S_comp); unconditionally per grid before known_edge_state is consulted. AMReX fabPtr: AMREX_ASSERT(mfi.LocalIndex() < indexArray.size()); return m_fabs_v[li]; — both empty for MultiFab(), then ->const_array() dereferences the garbage pointer. No in-repo caller exists, so it is latent within stock IAMR, but any invocation of this documented public API crashes as claimed.

F065 — confirmed (one verifier lens)

Lens 1 (refutation attempt): MacProj.cpp:562 divu_fp.reset(ns_level.getDivCond(nghost_force(),prev_time)); with no dsdt Saxpy, vs NSB.cpp:3423-24 and NS.cpp:732-33 MultiFab::Saxpy(*divu_fp, 0.5*dt, *dsdt,...) before edge prediction; divu feeds Godunov::/BDS::ComputeEdgeState via ComputeFluxesOnBoxFromState. Commit 82d3106: "we need to reconstruct the edge states exactly as we made them during the advection step". Reachable: do_temp=1 sets have_divu/have_dsdt (NS_setup.cpp:372-391); mac_sync gated only by do_reflux (default 1) && level<finest_level.

F077 — confirmed (one verifier lens)

Lens 1 (refutation attempt): MacProj.cpp:1064-1070: diff.resize(...) uses default DataAllocator -> The_Arena() (device; AMReX_Arena.cpp:59 the_arena_is_managed = false); FACopyDescriptor.H:852 FillFab does destFab.template copy<RunOn::Host>(*fcdp->localFabSource,...) where localFabSource aliases the device u_mac FAB; then diff.minus<RunOn::Host>(u_mac[dim][idx],...). Reachable on GPU: default 0 (lines 56-60) but pp.query("check_umac_periodicity",...) (line 69) allows opt-in, and mac_project calls it when any direction is periodic.

---

_{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: F024, F065, F077. Reviewer unit(s): MacProj, 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.}