Tracer advection: corrector base mis-reflected, double advection at C-F, uninitialized r

[WeiqunZhang]

Severity: high/medium · Category: correctness, memory-ub · Fix order: 9 of 29 — fix this 9th.

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

Locations: src/particles/incflo_PCEvolve.cpp:99, src/particles/incflo_PCEvolve.cpp:21, src/particles/incflo_PCEvolve.cpp:180

Based on commit 7307d872 (line numbers refer to that tree).

All three defects are in the tracer-particle advection path (incflo_PC::EvolveParticles / AdvectWithFlow in src/particles/incflo_PCEvolve.cpp); they are independent, but one sweep of the file fixes them.

The defect

src/particles/incflo_PCEvolve.cpp:99 — In ipass 0 of AdvectWithFlow, v_ptr holds the particle's ORIGINAL POSITION (line 90), yet boundary reflection does v_ptr[dim][i] *= -1.0 (lines 99, 107); negation only equals reflection about coordinate 0, so the corrector base position is corrupted.

src/particles/incflo_PCEvolve.cpp:21 — EvolveParticles calls a full Redistribute() after advecting each level, but evolveTracerParticles (incflo_Tracers.cpp:74-78) calls it inside an ascending level loop, so a particle advected at level L and redistributed into level L+1 grids is advected again at L+1 the same step.

src/particles/incflo_PCEvolve.cpp:180 — Real r is uninitialized and the assert at line 160 (cyl_direction >= 0 && cyl_direction <= AMREX_SPACEDIM) admits unhandled values: in 2D only direction==2 sets r (0 and 1 compile out), and in 3D direction==3 passes the assert; r is then compared at line 191.

Why it matters

F006: Tracer particles, non-periodic domain not anchored at 0 (e.g. prob_lo=-1, prob_hi=1): particle at x=0.99 whose half-step crosses the high boundary gets corrector base -0.99; ipass 1 places it near x=-0.94, teleporting it across the domain. Even with plo=0, high-boundary reflections land at x0-dtv instead of 2phi-x0-dt*v.

F026: amr.max_level>=1 + incflo.use_tracer_particles=1: any particle crossing from the coarse region into a fine-grid region during the step is moved twice (about 2*dt) that step, systematically distorting tracer trajectories near coarse-fine boundaries.

F025: 2D build + tracer particles with test_2d/benchmark.poiseuille_cylinder_bingham (cylinder.direction=1): every step each particle reads uninitialized stack/register garbage for r, nondeterministically deleting particles (id=-1) whenever garbage > radius. Same for 3D with cylinder.direction=3.

How to reach it

• Build -DINCFLO_PARTICLES=ON; inputs: incflo.use_tracer_particles=1, geometry.prob_lo=-1 -1 -1, prob_hi=1 1 1, non-periodic with xhi.type="po" (outflow): particles crossing the high boundary during the half-step hit lines 104-108.

Suggested fix

F006: in the predictor pass (ipass 0) v_ptr caches the particle's pre-step position, so a boundary bounce must mirror that cached base about the wall actually crossed — 2plo[d]-x or 2phi[d]-x — exactly as lines 98/106 already do for p.pos. The bare *= -1.0 is correct only in the corrector pass, where v_ptr holds velocity.

F026: redistribute once per step, after every level has advected. evolveTracerParticles loops levels ascending, and AMReX's Redistribute assigns each particle to the finest level containing it, so the full per-level Redistribute() inside EvolveParticles hands coarse-to-fine crossers to L+1 before L+1 is advected. Move the call after the level loop in incflo::evolveTracerParticles; EvolveParticles has no other caller today, but decide whether the virtual API should still promise a redistribute.

F025: initialize r and tighten the assert to only handled cylinder.direction values (3D: 0 <= dir < AMREX_SPACEDIM; 2D: dir == 2) so bad inputs fail at setup instead of reading indeterminate stack in the device lambda. Maintainer call: the shipped test_2d/benchmark.poiseuille_cylinder_bingham sets cylinder.direction=1 and would now be rejected — update the input or make 2D ignore direction.

For src/particles/incflo_PCEvolve.cpp:99 (F006):

--- a/src/particles/incflo_PCEvolve.cpp
+++ b/src/particles/incflo_PCEvolve.cpp
@@ -96,15 +96,15 @@
                         if (!is_periodic[dim] && p.pos(dim) < plo[dim])
                         {
                             p.pos(dim) = 2.0*plo[dim] - p.pos(dim);
-                            v_ptr[dim][i] *= -1.0;
+                            v_ptr[dim][i] = 2.0*plo[dim] - v_ptr[dim][i];
                         }
                         //
                         // Reflect off high domain boundaries if not periodic
                         //
                         if (!is_periodic[dim] && p.pos(dim) > phi[dim])
                         {
                             p.pos(dim) = 2.0*phi[dim] - p.pos(dim);
-                            v_ptr[dim][i] *= -1.0;
+                            v_ptr[dim][i] = 2.0*phi[dim] - v_ptr[dim][i];
                         }
                     }
                 } else {

_{In ipass 0, v_ptr holds the saved base position (line 90), so a boundary reflection must mirror it about the actual boundary plane (2plo-x, 2phi-x), matching how p.pos is reflected one line above. The ipass 1 branch is untouched: there v_ptr holds velocity, for which *= -1.0 is correct.}

For src/particles/incflo_PCEvolve.cpp:21 (F026):

--- a/src/particles/incflo_PCEvolve.cpp
+++ b/src/particles/incflo_PCEvolve.cpp
@@ -18,7 +18,6 @@
         AdvectWithFlow(a_lev, a_dt_lev, AMREX_D_DECL(a_umac, a_vmac, a_wmac));
     }
 
-    Redistribute();
     return;
 }
 
--- a/src/particles/incflo_Tracers.cpp
+++ b/src/particles/incflo_Tracers.cpp
@@ -75,7 +75,8 @@
         {
             particleData[incfloParticleNames::tracers]->EvolveParticles(lev, m_dt,
                                                                         AMREX_D_DECL(u_mac[lev],v_mac[lev],w_mac[lev]));
         }
+        particleData.Redistribute();
     }
 }
 #endif

_{Two-file fix: drops the all-level Redistribute() from per-level EvolveParticles and calls ParticleData::Redistribute() (ParticleData.H:118, loops all species' full Redistribute) once after the ascending level loop, so a coarse-to-fine crosser is not advected twice per step. Only in-repo caller is evolveTracerParticles; maintainer should confirm EvolveParticles has no external standalone users relying on its redistribute.}

For src/particles/incflo_PCEvolve.cpp:180 (F025):

--- a/src/particles/incflo_PCEvolve.cpp
+++ b/src/particles/incflo_PCEvolve.cpp
@@ -157,7 +157,7 @@
 
         int cyl_direction;
         pp.get("direction",cyl_direction);
-        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= AMREX_SPACEDIM);
+        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= 2);
 
         // Remove particles that are outside of the cylindner
         for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti)
@@ -177,16 +177,21 @@
                 Real z =  p.pos(2) - z_ctr;
 #endif
 
-                Real r;
+                Real r = 0.;
                 if (cyl_direction == 2) {
                     r =  std::sqrt(x*x + y*y);
 #if (AMREX_SPACEDIM == 3)
                 } else if (cyl_direction == 1) {
                     r =  std::sqrt(x*x + z*z);
                 } else if (cyl_direction == 0) {
                     r =  std::sqrt(y*y + z*z);
+#else
+                } else if (cyl_direction == 1) {
+                    r =  std::abs(x);
+                } else if (cyl_direction == 0) {
+                    r =  std::abs(y);
 #endif
                 }
 
                 if (r > cyl_radius) {
                     p.id() = -1;

_{Initializes r, tightens the assert to reject direction 3 in 3D, and handles 2D directions 0/1 following amrex::EB2::CylinderIF 2D semantics (AMReX_EB2_IF_Cylinder.H: dir 1 -> distance |x|, dir 0 -> |y|), matching the shipped test_2d/benchmark.poiseuille_cylinder_bingham (direction=1). Maintainer must confirm 2D slab removal is intended; alternative is asserting cyl_direction==2 in 2D.}

_{Diff(s) are against 7307d872, 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

F006 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): ipass 0: line 90 v_ptr[dim][i] = p.pos(dim); stores original position x0; lines 99/107 v_ptr[dim][i] *= -1.0; on reflection; ipass 1 line 113 p.pos(dim) = v_ptr[dim][i] + a_dt*v[dim];. Negating x0 equals reflection 2*b-x0 only for b=0; high-boundary or plo!=0 corrupts the corrector base exactly as claimed. git log -S shows no fix.

Lens 2 (reachability/intent): incflo_PCEvolve.cpp:90 stores base position in v_ptr; :99/:107 apply v_ptr[dim][i] *= -1.0 while p.pos gets true reflection 2*plo/phi - p.pos; :113 uses v_ptr as corrector base. Negation equals reflection only about 0. In ipass 1 v_ptr holds velocity (:114), where the sign flip is correct — the operation was pasted into ipass 0 with wrong semantics. Introduced undocumented in 011b26a (#99); still present; corrupted base consumed immediately, Redistribute keeps the in-domain teleported particle.

F026 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 21 Redistribute(); (full, all-level) after advecting only a_lev; incflo_Tracers.cpp:74-78 loops lev 0..finest_level ascending with same m_dt. AMReX AmrAssignGrid (AMReX_ParticleLocator.H) iterates for (int lev = lev_max; lev >= lev_min; --lev) assigning particles to the finest containing level, so a coarse->fine crosser is advected again at L+1 the same step.

F025 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 180 Real r; uninitialized; setters for direction 1/0 are inside #if (AMREX_SPACEDIM == 3) (lines 183-188); assert line 160 cyl_direction >= 0 && cyl_direction <= AMREX_SPACEDIM admits 2D dir 0/1 and 3D dir 3; line 191 reads if (r > cyl_radius). Shipped test_2d/benchmark.poiseuille_cylinder_bingham sets cylinder.radius=1, cylinder.direction=1.

---

_{Based on commit 7307d872, which is also the tree the audit verified against. From an automated audit of src/. Audit finding ids: F006, F026, F025. Reviewer unit(s): Derive+Particles. 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.}