RZ scaleVar zeroes axial inflow ghosts; computeRhoG y-hi rho stencil wrong and OOB

[WeiqunZhang]

Severity: high/medium · Category: correctness, memory-ub, physics-numerics · Fix order: 2 of 21 — fix this 2nd.

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

Locations: Source/Projection.cpp:1312, Source/Projection.cpp:2308, Source/Projection.cpp:2323

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

F008 and F014 are the same defect: one edit fixes both, and the two attached diffs are alternate spellings of it — apply either, not both. F036 and F037 are two independent transcription slips inside the same computeRhoG y-hi block and should land together.

The defect

Source/Projection.cpp:1312 — scaleVar's RZ branch zeroes velocity ghost cells outside the domain in the axial (j) direction, destroying inflow values the nodal projector's divergence stencil must read; the old Fortran radmpyvel multiplied those ghosts by radius precisely to preserve inflow values. _{Reported independently at this same line by F014; each reviewer's own wording and evidence is under Verification evidence below.}

Source/Projection.cpp:2308 — In computeRhoG's 3D y-hi outflow-face main loop, rho_ii mixes density rows: it averages rho(i,j-1,k) with rho(i-1,j-2,k), while rho_ii must be the second interior row j-2 averaged across cells i-1 and i, as in the twin xlo/xhi/ylo blocks.

Source/Projection.cpp:2323 — In computeRhoG's 3D y-hi face, x-lo ext_dir edge branch, rho_ii = rho(i-2,j-1,k) keeps the xhi-face stencil (normal = x) instead of the y-normal stencil rho(i-1,j-2,k); with i = domlo.x, i-2 = domlo.x-2 reads outside the rho FArrayBox, which is grown only 1 cell tangentially.

Why it matters

F008: 2D RZ run with axial inflow (e.g. Exec/run2d/test_grids/inputs_2_rz_zinflow, ns.lo_bc="5 1"): every level_project/initialSyncProject call zeroes the z-inflow ghost velocity before doMLMGNodalProjection; AMReX mlndlap_divu then sees zero normal velocity at inflow boundary nodes, so the projection treats the inflow as a wall, corrupting velocity and pressure near the inflow each step (uniform inflow is not preserved).

F036: 3D build, ns.gravity != 0, outflow BC on the y-hi face, do_outflow_bcs=1 (default), variable density near that face: the extrapolated boundary density rhoExt=0.5*(3*rho1-rho2) uses a wrong rho2 at every node of the face, so the hydrostatic Dirichlet phi for the nodal projection is wrong, producing spurious pressure/velocity near the outflow. Constant-density runs mask it.

F037: 3D, ns.gravity != 0, outflow at y-hi, inflow (density ext_dir) at x-lo: the rho FAB spans x in [domlo.x-1, domhi.x+1] (state_strip grown 1 tangentially, Projection.cpp:1796-1797,1888), so rho(domlo.x-2,...) is an out-of-bounds host read -> abort with bound checking, garbage otherwise, and in any case the wrong stencil corrupts phi along that edge column of the outflow face.

How to reach it

• 2D non-EB run2d build: geometry.coord_sys=1, ns.vel_visc_coef=0 (required by RZ abort), ns.lo_bc="3 1", ns.hi_bc="4 2" (axial inflow); per Exec/run2d/test_grids/inputs_2_rz_zinflow with visc set to 0. Hit every level_project.
• 2D non-EB build; geometry.coord_sys=1; ns.lo_bc="5 1" (axial inflow); ns.vel_visc_coef=0 (RZ+viscosity aborts). Essentially Exec/run2d/test_grids/inputs_2_rz_zinflow with vel_visc_coef set to 0. Hit every step in level_project.

Note from the audit

F008: Reaching-configuration correction. The shipped deck Exec/run2d/test_grids/inputs_2_rz_zinflow has geometry.coord_sys = 1, ns.lo_bc = 5 1, ns.hi_bc = 5 2 — not the 3 1 / 4 2 quoted in the reaching configuration above, which is not from this deck. The radial BCs are incidental to the defect (it concerns the axial out-of-domain ghost strip), so either set reproduces it; quote the deck's own values when reproducing. One documented edit is required: the deck ships ns.vel_visc_coef = 0.01, and RZ with viscosity hits the amrex::Abort at Source/NavierStokesBase.cpp:253, so set ns.vel_visc_coef = 0 to run it at all.

Suggested fix

RZ (line 1312). In scaleVar's RZ velocity loop, drop the j test so every ghost with i in [domlox,domhix] is scaled by (i+0.5)*dxr, keeping the axis and hi-r inflow special cases. Pre-port radmpyvel (dbf01617^:Source/Src_2d/PROJECTION_2D.F90) did exactly that, noting the divu stencil must "include these values ... because they might contain inflow values". Zeroing out-of-domain ghosts is set_boundary_velocity's job (line 2570), and it protects only inflow faces. Mirror the change in rescaleVar — divide those ghosts instead of writing BogusValue, as fort_raddiv did — so the pair stays invertible; decide also whether rescaleVar should restore the hi-r inflow ghost it now bogus-fills.

computeRhoG (2308, 2323). Pre-port rhogbc YHI (358fdce7^:Source/Src_3d/PROJECTION_3D.F90) uses 0.5*(rho(i,j-2,k)+rho(i-1,j-2,k)), and rho(i-1,j-2,k) at the x-lo ext_dir edge: precisely these two fixes. Every other y-hi branch still matches that Fortran, so no wider sweep is needed; the stale "compute y-edges" comment betrays the x-hi copy-paste.

For Source/Projection.cpp:1312 (F008):

--- a/Source/Projection.cpp
+++ b/Source/Projection.cpp
@@ -1292,8 +1292,12 @@
           amrex::ParallelFor(bx, AMREX_SPACEDIM, [=]
           AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
           {
-            if ( i >= domlox && i <= domhix &&
-                 j >= domloy && j <= domhiy)
+            // NOTE: cells outside the domain in the z-direction get multiplied
+            // by the radius rather than zeroed because the divu routine in the
+            // nodal solver includes them in its stencil and they might contain
+            // inflow values. set_boundary_velocity() zeroes the ones that must
+            // be zero.
+            if ( i >= domlox && i <= domhix )
             {
               velarr(i,j,k,n) = (static_cast<Real>(i)+ Real(0.5))*dxr*velarr(i,j,k,n);
             }
@@ -1411,8 +1415,7 @@
           amrex::ParallelFor(bx, AMREX_SPACEDIM, [=]
           AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
           {
-            if ( i >= domlox && i <= domhix &&
-                 j >= domloy && j <= domhiy)
+            if ( i >= domlox && i <= domhix )
             {
               velarr(i,j,k,n) /= (static_cast<Real>(i)+ Real(0.5))*dxr;
             }

_{Restores deleted FORT_RADMPYVEL/fort_raddiv semantics (Src_2d/PROJECTION_2D.F90, dropped in dbf0161): scale all axial ghosts, zero/bogus only i outside [domlo.x,domhi.x]. Safe because set_boundary_velocity() runs later in doMLMGNodalProjection and zeroes every non-inflow ghost, while amrex mlndlap_divu (AMReX_MLNodeLap_2D_K.H) reads normal velocity from the inflow ghost row. Maintainer may drop the rescaleVar hunk if BogusValue ghosts are wanted.}

For Source/Projection.cpp:2308 (F036):

--- a/Source/Projection.cpp
+++ b/Source/Projection.cpp
@@ -2305,7 +2305,7 @@
 
                   for (int k = hi.z-1; k >= lo.z; k--) {
                     rho_i   = 0.5 * (rho(i,j-1,k) + rho(i-1,j-1,k));
-                    rho_ii = 0.5 * (rho(i,j-1,k) + rho(i-1,j-2,k));
+                    rho_ii = 0.5 * (rho(i,j-2,k) + rho(i-1,j-2,k));
                     add_rhog(rho_i, rho_ii, rhog, phi(i,j,k));
                   }
                 }

_{Makes rho_ii the pure second-interior row j-2 averaged over i-1 and i, matching the xhi twin (rho(i-2,j)/rho(i-2,j-1)) and the ylo twin (rho(i,j+1)/rho(i-1,j+1)). Index j-2 is in bounds (the rho strip covers both interior rows). Result changes only where density varies near the y-hi outflow face.}

For Source/Projection.cpp:2323 (F037):

--- a/Source/Projection.cpp
+++ b/Source/Projection.cpp
@@ -2320,7 +2320,7 @@
                   if ( has_extdir_lo ) {
                     for (int k = hi.z-1; k >= lo.z; k--) {
                       rho_i   = rho(i-1,j-1,k);
-                      rho_ii = rho(i-2,j-1,k);
+                      rho_ii = rho(i-1,j-2,k);
                       add_rhog(rho_i, rho_ii, rhog, phi(i,j,k));
                     }
                   } else if ( has_hoextrap_lo ) {

_{Uses the x-lo ext_dir ghost column i-1 at the second interior row j-2, mirroring the ylo/x-lo twin (rho(i-1,j)/rho(i-1,j+1)) and the yhi/x-hi twin (rho(i,j-1)/rho(i,j-2)). Also eliminates the i-2 = domlo.x-2 read, which is outside the rho FAB (state_strip grown only one cell tangentially).}

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

F008 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): scaleVar RZ branch: else { velarr(i,j,k,n) = 0.; } zeroes axial ghost vel; set_boundary_velocity (line 2432) only protects, never refills; amrex mlndlap_divu reads normal vel(i,j-1,k,1) from inflow ghosts (zero_jlo zeroes tangential only). Old Fortran FORT_RADMPYVEL multiplied all axial ghosts by r(i): 'in the divu routine ... we need to include these values because they might contain inflow values'. inputs_2_rz_zinflow (lo_bc=5 1, coord_sys=1) exercises it.

Lens 2 (reachability/intent): Projection.cpp:1312 zeroes axial ghost vel; removed Fortran radmpyvel (commit dbf0161) scaled growntilebox j-ghosts by r, commenting they "might contain inflow values" needed by the divu stencil. set_boundary_velocity:2569 only protects, never refills. AMReX mlndlap_divu (AMReX_MLNodeLap_2D_K.H:785-796) reads normal ghost vel at inflow nodes (zero_jlo kills only tangential terms), so the projection sees zero inflow — wall behavior, unmasked.

F014 — confirmed (two independent verifier lenses)

Reported as: In RZ, Projection::scaleVar zeroes velocity ghost cells outside the domain in the axial (z) direction, destroying the ext_dir inflow values that the nodal projection's divergence stencil reads, whereas the Cartesian path preserves them (scaleVar never touches vel; set_boundary_velocity explicitly protects inflow ghosts).

Failure scenario: 2D RZ run (geometry.coord_sys=1) with inflow at an axial boundary, e.g. ns.lo_bc="5 1" as in Exec/run2d/test_grids/inputs_2_rz_zinflow: level_project (every step) and initialVelocityProject scale U_new via scaleVar before doMLMGNodalProjection; mlndlap_divu reads the normal velocity from the j-1 ghost cells at the inflow node but finds 0, so the projection sees a spurious divergence sheet at the inlet and produces wrong pressure/velocity there.

Lens 1 (refutation attempt): Same defect as F008. Cartesian scaleVar branch (lines 1320-1350) touches only sig, never vel, so inflow ghosts survive there; RZ branch zeroes them after setPhysBoundaryValues filled them (level_project lines 212-218). mlndlap_divu at j==domlo.y inflow node reads vel(i,j-1,k,1)=0 -> spurious divergence at inlet. No fix since dbf0161 (git log -S 'set other vals outside the domain to 0').

Lens 2 (reachability/intent): Projection.cpp RZ scaleVar zeroes vel ghosts with j outside [domloy,domhiy] (only n==0,i>domhix spared); set_boundary_velocity (called after, in doMLMGNodalProjection) only protects, never refills. AMReX MLNodeLaplacian::compRHS grows bx_vel across inflow faces, so mlndlap_divu reads the zeroed ghost as inflow flux. The Fortran radmpyvel this replaced (deleted in dbf0161) scaled all j-rows by r, commenting ghosts "might contain inflow values" and must stay in the stencil — the C++ rewrite regressed that.

F036 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 2308: rho_ii = 0.5 * (rho(i,j-1,k) + rho(i-1,j-2,k)); mixes rows j-1/j-2; twins are pure second-row: xhi line 2139 0.5*(rho(i-2,j,k)+rho(i-2,j-1,k)), ylo line 2241 0.5*(rho(i,j+1,k)+rho(i-1,j+1,k)). Should be rho(i,j-2,k). In-bounds (strip y=[domhi.y-1,domhi.y]) so wrong-stencil only; masked when density constant.

F037 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Lines 2322-2323 (yhi face, has_extdir_lo): rho_i = rho(i-1,j-1,k); rho_ii = rho(i-2,j-1,k); with i=lo.x=domlo.x (phi strip starts at domlo.x). rho FAB = state_strip grown 1 tangentially (lines 1796-1797, 1888), x-range [domlo.x-1, domhi.x+1], so rho(domlo.x-2,..) is OOB by 1. Twins move in y: yhi/x-hi line 2349 rho_ii = rho(i,j-2,k), ylo/x-lo line 2256 rho_ii = rho(i-1,j+1,k); correct is rho(i-1,j-2,k). Reachable: 3D, y-hi outflow, x-lo inflow (density ext_dir per NS_BC.H), gravity != 0.

---

_{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: F008, F014, F036, F037. Reviewer unit(s): Projection-1, Projection-2, theme:rz-metrics. 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.}