particleData.Redistribute() inside RemakeLevel/MakeNewLevelFromCoarse runs against the OLD grids/finest_level, so the first step after a regrid advects tracer particles through a stale grid mapping (wrong FAB or out-of-range dummy MultiFab)

[WeiqunZhang]

Severity: medium · Category: amr-coupling
Locations: src/incflo_regrid.cpp:63-65, src/incflo_regrid.cpp:124-126, src/incflo.cpp:147-154, src/incflo_apply_predictor.cpp:233, src/incflo_apply_corrector.cpp:183, src/particles/incflo_PCEvolve.cpp:57-71
Based on commit 46de3367 (line numbers refer to that tree).

The defect

PR #164 added particleData.Redistribute() at the end of MakeNewLevelFromCoarse (regrid.cpp:63-65) and the
same call sits at the end of RemakeLevel (regrid.cpp:124-126), "to ensure particles are on the finest level".
Both calls run too early to do that. AmrCore::regrid (AMReX Src/AmrCore/AMReX_AmrCore.cpp:96-135) calls the
callbacks before it installs the new layout:

RemakeLevel(lev, time, level_grids, level_dmap);      // :109
SetBoxArray(lev, level_grids);                        // :110
SetDistributionMap(lev, level_dmap);                  // :112
...
MakeNewLevelFromCoarse(lev, time, new_grids[lev], new_dmap);   // :121
SetBoxArray(lev, new_grids[lev]);                     // :122
...
finest_level = new_finest;                            // :135

The particle container was built on GetParGDB() (an AmrParGDB), whose ParticleBoxArray(lev) is
m_amrcore->boxArray(lev) and whose finestLevel() is m_amrcore->finestLevel()
(AMReX_AmrParGDB.H:177-183, 258-261). Redistribute_impl sizes m_particles/m_dummy_mf from
m_gdb->finestLevel() and rebuilds m_dummy_mf[lev] from ParticleBoxArray(lev) (AMReX_ParticleContainerI.H:1452-1464,
AMReX_ParticleContainerBase.cpp:84-100). So inside RemakeLevel(lev) the redistribution is done against the
old BoxArray/DistributionMapping of lev, and inside MakeNewLevelFromCoarse(lev) against the old
finest_level (the new level does not exist for the container). After regrid() returns nothing else
redistributes: the next call is Advance -> ApplyPredictor -> evolveTracerParticles (predictor.cpp:233 for
Godunov/BDS, corrector.cpp:183 for MOL), and only after that particleData.Redistribute() (Tracers.cpp:85,
advance.cpp:86).

incflo_PC::AdvectWithFlow therefore iterates with ParIterType pti(*this, a_lev), whose MFIter is built on
the stale m_dummy_mf[a_lev] (AMReX_ParIter.H:173), takes grid = pti.index() (an index into the OLD
BoxArray) and does (*a_umac)[grid] (PCEvolve.cpp:71) on u_mac[lev], which was defined on the NEW
grids[lev]/dmap[lev] (predictor.cpp:91). FabArray::fabPtr(int K) is m_fabs_v[localindex(K)]
(AMReX_FabArray.H:2074-2079) with localindex returning -1 for a global index this rank does not own
(AMReX_FabArrayBase.H:124-132): either an out-of-bounds vector access, or -- when the old index happens to be
owned -- the FAB of an unrelated box, so mac_interpolate reads outside the FAB for every particle. For a
regrid that adds a level, m_dummy_mf and m_particles still have old_finest+1 entries and
pti(*this, new_finest) indexes pc.m_dummy_mf[level] past the end. For a regrid that removes a level,
the particles on the removed level are simply not advected that step (the loop in evolveTracerParticles
stops at the new finest_level).

Why it matters

With tracer particles and amr.regrid_int > 0, the first step after any regrid that changes a level's grids or
the number of levels is undefined behaviour: a segfault or, silently, particle velocities interpolated from the
wrong box. Debug builds abort in AdvectWithFlow's AMREX_ASSERT(OK(...)) / a_lev < GetParticles().size().
This is exactly the configuration PR #164 set out to support.

How to reach it

Build with -DINCFLO_PARTICLES=ON (GNUmake USE_PARTICLES=TRUE), any DIM, EB or not. Take
test_no_eb_2d/benchmark.rayleigh_taylor (Godunov, amr.max_level = 2, amr.regrid_int = 2) and add
incflo.use_tracer_particles = 1. At step 2 regrid(0, m_cur_time) remakes level 1/2 (the interface moves);
the following ApplyPredictor calls AdvectWithFlow(lev=1) with pti.index() from the pre-regrid BoxArray.
On one rank with a different grid count the localindex lookup fails (crash); with the same count it reads the
wrong FAB. Any input with amr.max_level >= 1 and amr.regrid_int > 0 plus particles shows the same.

Suggested fix

Redistribute once, after regrid() has installed the new BoxArrays/DistributionMaps and updated
finest_level; drop the two in-callback calls, which can never see the new layout. (The
regrid_on_restart path is already covered by the Redistribute() in InitData at incflo.cpp:99, which runs
after ReadCheckpointFile returns.) Redistribute() grows m_particles/m_dummy_mf to the new
finest_level+1 and, for a shrink, moves particles from the removed level down.

--- a/src/incflo_regrid.cpp
+++ b/src/incflo_regrid.cpp
@@ -59,10 +59,6 @@
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,lev));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Remake an existing level using provided BoxArray and DistributionMapping and
@@ -120,10 +116,6 @@
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,finest_level));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Rebuild macproj on demand.  Note that this must not be done inside ClearLevel:
--- a/src/incflo.cpp
+++ b/src/incflo.cpp
@@ -148,6 +148,12 @@
         {
             if (m_verbose > 0) amrex::Print() << "Regridding...\n";
             regrid(0, m_cur_time);
+#ifdef INCFLO_USE_PARTICLES
+            // Must run after regrid() returns: inside RemakeLevel /
+            // MakeNewLevelFromCoarse the ParGDB still sees the old BoxArray,
+            // DistributionMapping and finest_level.
+            particleData.Redistribute();
+#endif
             if (m_verbose > 0 && ParallelDescriptor::IOProcessor()) {
                 printGridSummary(amrex::OutStream(), 0, finest_level);
             }

(audit/notes/T3-scratch/008.diff, checked with git apply --check.)

[asalmgren]

closed by #265