incflo_PCInit 2D cylinder cull: ALWAYS_ASSERT rejects cylinder.direction=2 (the only value its disk formula fits) and culls a disk for directions 0/1 where AdvectWithFlow and EB2::CylinderIF use slabs (regression from PR #213) #243

Severity: medium · Category: correctness
Locations: src/particles/incflo_PCInit.cpp:243-245, src/particles/incflo_PCInit.cpp:266-269, src/particles/incflo_PCEvolve.cpp:169-171, src/particles/incflo_PCEvolve.cpp:195-210, src/embedded_boundaries/eb_cylinder.cpp:20-34, src/embedded_boundaries/eb_cylinder.cpp:53
Based on commit 46de3367 (line numbers refer to that tree).

The defect

PR #213 (fix for #171, "cylinder cull ignores cylinder.direction") rewrote the init-time cull in
initializeParticlesUniformDistributionInBox:

        int cyl_direction;
        pp.get("direction",cyl_direction);
        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction < AMREX_SPACEDIM);   // :245
        ...
#if (AMREX_SPACEDIM == 3)
                Real r = (cyl_direction == 0) ? std::sqrt(y*y + z*z) : ... ;             // :263-265
#else
                // In 2D the only cylinder axis is the out-of-plane one
                Real r =  std::sqrt(x*x + y*y);                                          // :268
#endif

In a 2D build this is self-contradictory. The assert admits only direction = 0 or 1, but the 2D
formula is the distance from the out-of-plane axis, i.e. the direction = 2 disk, and it ignores the
direction it just validated. The step-time cull that PR #215 wrote in AdvectWithFlow uses the opposite
conventions (incflo_PCEvolve.cpp:171, :195-210):

        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= 2);
        ...
                if (cyl_direction == 2) {  r = std::sqrt(x*x + y*y);          // disk
#else
                } else if (cyl_direction == 1) { r = std::abs(x);            // slab about the y axis
                } else {                         r = std::abs(y);            // slab about the x axis

which matches amrex::EB2::CylinderIF in 2D (AMReX_EB2_IF_Cylinder.H:51-81: case 0 -> d2 = y*y,
case 1 -> d2 = x*x, default -> x*x + y*y) and the EB builder make_eb_cylinder
(eb_cylinder.cpp:53, default direction = 0). So, in 2D:

• cylinder.direction = 2 (a disk, the configuration AdvectWithFlow and the docs' "2 = z" support):
  the run aborts in particle initialisation with the AMREX_ALWAYS_ASSERT at :245.
• cylinder.direction = 0 or 1 (a slab/channel, what make_eb_cylinder builds): initialisation keeps
  only the particles inside the disk x^2+y^2 <= R^2 while every later step keeps those inside the slab
  |y| <= R (or |x| <= R). Particles seeded in the channel outside the disk are silently deleted at
  t = 0, so the tracer field never covers the flow domain it was meant to trace.

The 3D branch is consistent between the two functions.

Why it matters

With a 2D USE_PARTICLES build and any cylinder.* block (which is also what incflo.geometry = cylinder
reads), the only accepted directions produce an initial particle set that does not match the region the
solver subsequently confines particles to, and the natural 2D disk setting is rejected outright. The
inconsistency was introduced by the fix for #171 and is exactly the class of bug that fix set out to remove.

How to reach it

DIM = 2, USE_PARTICLES = TRUE (-DINCFLO_PARTICLES=ON), EB or not, CPU or GPU.
Take test_no_eb_2d/benchmark.taylor_green_vortices and add

incflo.use_tracer_particles = 1
cylinder.radius    = 0.3
cylinder.center    = 0.5 0.5 0.0
cylinder.direction = 2

-> abort at incflo_PCInit.cpp:245. Change cylinder.direction = 0: the run proceeds, "Initialized N tracer
particles" reports only the particles inside the disk of radius 0.3 (approx. 0.28 of the domain), whereas
AdvectWithFlow will from then on keep everything with |y-0.5| <= 0.3 (0.6 of the domain).
For the EB twin use test_2d/ with incflo.geometry = cylinder, cylinder.internal_flow = true,
cylinder.direction = 0: the fluid region is the channel, the initial tracer set is the disk.

Suggested fix

Use the same direction convention as AdvectWithFlow (and EB2::CylinderIF) at init time. Diff checked
with git apply --check (audit/notes/T8-scratch/028.diff):

--- a/src/particles/incflo_PCInit.cpp
+++ b/src/particles/incflo_PCInit.cpp
@@ -238,11 +238,12 @@
         Real z_ctr = cyl_center[2];
 #endif
 
-        // The cylinder axis matters here exactly as it does in AdvectWithFlow;
-        // assuming a z-parallel axis culls against the wrong axis for direction 0/1.
+        // The cylinder axis matters here exactly as it does in AdvectWithFlow,
+        // which accepts 0, 1 and 2 in both 2D and 3D; the culling below must
+        // use the same distance for each direction.
         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 cylinder
         for (ParIterType pti(*this, lev); pti.isValid(); ++pti)
@@ -264,8 +265,11 @@
                        : (cyl_direction == 1) ? std::sqrt(x*x + z*z)
                                               : std::sqrt(x*x + y*y);
 #else
-                // In 2D the only cylinder axis is the out-of-plane one
-                Real r =  std::sqrt(x*x + y*y);
+                // Same convention as AdvectWithFlow and amrex::EB2::CylinderIF:
+                // direction 2 is a disk, directions 0 and 1 are slabs.
+                Real r = (cyl_direction == 2) ? std::sqrt(x*x + y*y)
+                       : (cyl_direction == 1) ? std::abs(x)
+                                              : std::abs(y);
 #endif
 
                 if (r > cyl_radius) {

A follow-up cleanup would be to share one cylinder_distance(x,y,z,direction) helper between the two
functions so the two culls cannot drift apart again.