ScalMinMax limiters seed max with numeric_limits::min(); do_denminmax limiter is a no-op

[WeiqunZhang]

Severity: medium/low · Category: correctness · Fix order: 10 of 21 — fix this 10th.

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

Locations: Source/NavierStokesBase.cpp:2785, Source/NavierStokesBase.cpp:4281, Source/NavierStokesBase.cpp:4339

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

The six findings collapse to two edits best landed as one patch: the shared seed line in both ScalMinMax kernels, and the do_denminmax call site — which only begins to do any work once the seeds are right.

The defect

Source/NavierStokesBase.cpp:2785 — The ns.do_denminmax density limiter is an exact no-op: ConservativeScalMinMax is called with Density as both the scalar and the normalizing density, so it clamps rho/rho==1 into the neighborhood range of rhoo/rhoo==1 and writes rho back unchanged.

Source/NavierStokesBase.cpp:4281 — ConservativeScalMinMax has the same defect: Real smx = std::numeric_limits<Real>::min(); is the smallest positive normal, not the lowest representable value, so the upper bound on the limited quantity s/rho is wrong when the whole neighborhood of s/rho is negative. _{Reported independently at this same line by F072, F076; each reviewer's own wording and evidence is under Verification evidence below.}

Source/NavierStokesBase.cpp:4339 — ConvectiveScalMinMax initializes the running maximum with std::numeric_limits<Real>::min(), which for floating point is the smallest positive normal (~2.2e-308), not the most negative value; the upper clip bound is therefore wrong for any stencil whose values are all negative. _{Reported independently at this same line by F073; each reviewer's own wording and evidence is under Verification evidence below.}

Why it matters

F047: Any variable-density run (EB or not) with ns.do_denminmax=1: in the kernel sn and rhon alias S_new(Density) and so/rhoo alias old density, so smn=smx=1 and sn=min(max(1,1),1)*rhon=rhon. Density overshoots at sharp interfaces are never corrected; results are bit-identical to do_denminmax=0, silently defeating the advertised limiter.

F031: ns.do_scalminmax=1 with ns.do_cons_trac=1 (conservative tracer) and a tracer whose specific value rho*q/rho is negative locally: smx collapses to ~2.2e-308, so sn is clipped to ~0 instead of the true local max of so/rhoo, corrupting the conservatively advected scalar.

F030: ns.do_scalminmax=1 with a non-conservative tracer/temperature that is negative in a region (e.g. tracer initialized to -1). For a cell whose 3x3(x3) old-state neighborhood max is -0.8, smx becomes 2.2e-308 instead of -0.8, so sn = min(max(sn,smn),smx) clips an overshoot to ~0 instead of -0.8 — an actively wrong value, and legitimate overshoots below 0 are never limited.

Note from the audit

F031: Twin of F030 (Source/NavierStokesBase.cpp:4339) — both lines must be fixed together. The two are the same one-line seed bug (std::numeric_limits<Real>::min(), the smallest positive normal, where lowest() is meant) in the sibling functions ConservativeScalMinMax (line 4281) and ConvectiveScalMinMax (line 4339); the folded F076 covers both lines in a single claim. Fixing one leaves the other wrong.

F030: Twin of F031 (Source/NavierStokesBase.cpp:4281, ConservativeScalMinMax) — the identical min()-vs-lowest() seed bug in the sibling function. Apply both hunks together.

Suggested fix

Seed both accumulators with std::numeric_limits<Real>::lowest() — lines 4281 and 4339 must land together, since they are the same bug in sibling routines. The pre-port kernels consscalminmax/convscalminmax (Source/Src_2d/GODUNOV_2D.F90, deleted in cfe6d44) took an explicit 9-point min/max with no sentinel, so signed fields limited correctly; the sentinel arrived with the C++ port d14328d, making this a migration regression.

For do_denminmax the quantity to bound is rho, not rho/rho, so the raw-value limiter is right: ConvectiveScalMinMax(S_new, Density, Smf, 0). That aliasing predates the port; it is not a regression.

Your calls: once lowest() binds, an EB cell with a fully covered stencil takes -1.8e308 (times rhon = COVERED_VAL, so -inf), so also guard with if (vfrac(i,j,k) == 0.) return;, preserving the invariant set at NavierStokes.cpp:439; and Exec/eb_run2d/regtest.2d.hotspot plus eb_run3d/regtest.3d.hotspot enable both flags, so their benchmarks change.

For Source/NavierStokesBase.cpp:2785 (F047):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -2771,20 +2771,18 @@
         if (do_denminmax)
         {
             //
             // Must do FillPatch here instead of MF iterator because we need the
             // boundary values in the old data (especially at inflow)
             //
-            const int index_new_s   = Density;
-            const int index_new_rho = Density;
-            const int index_old_s   = index_new_s   - Density;
-            const int index_old_rho = index_new_rho - Density;
-
             FillPatchIterator S_fpi(*this,S_old,1,prev_time,State_Type,Density,1);
             MultiFab& Smf=S_fpi.get_mf();
 
-            ConservativeScalMinMax(S_new, index_new_s, index_new_rho,
-                                   Smf,   index_old_s, index_old_rho);
-
+            //
+            // Clamp the new density to the min/max of the old-time density in
+            // the neighborhood. (The conservative s/rho form with s = rho is
+            // an identity, i.e. a no-op.)
+            //
+            ConvectiveScalMinMax(S_new, Density, Smf, 0);
         }
 
         ++sComp;

_{Clamps density directly to the 3x3(x3) min/max of old density via ConvectiveScalMinMax (Smf has the 1 ghost cell the stencil needs); the conservative s/rho form with s==rho is provably an identity, and was in the original Fortran call too. Note: makes ns.do_denminmax=1 change results for the first time (previously bit-exact no-op).}

For Source/NavierStokesBase.cpp:4281 (F031):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -4278,7 +4278,7 @@
         AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real smn = std::numeric_limits<Real>::max();
-            Real smx = std::numeric_limits<Real>::min();
+            Real smx = std::numeric_limits<Real>::lowest();
 
 #if (AMREX_SPACEDIM==3)
             int ks = -1;

_{In ConservativeScalMinMax, std::numeric_limits::min() is the smallest positive normal; lowest() is the correct identity for a running maximum, so an all-negative s/rho neighborhood bounds correctly.}

For Source/NavierStokesBase.cpp:4339 (F030):

--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -4336,7 +4336,7 @@
         AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real smn = std::numeric_limits<Real>::max();
-            Real smx = std::numeric_limits<Real>::min();
+            Real smx = std::numeric_limits<Real>::lowest();
 
 #if (AMREX_SPACEDIM==3)
             int ks = -1;

_{Same fix in ConvectiveScalMinMax: initialize the running maximum with std::numeric_limits::lowest() instead of ::min(), so all-negative old-state neighborhoods produce the true (negative) upper clip bound.}

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

F047 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Lines 2777-2786: index_new_s=index_new_rho=Density, Smf FillPatched with 1 comp (Density), index_old_s=index_old_rho=0. In the kernel sn/rhon and so/rhoo alias the same arrays, so so/rhoo==1 exactly, smn=smx=1, and sn = min(max(1,1),1)*rhon = rhon bit-exactly (4301-4307). ns.do_denminmax=1 is a provable no-op.

F031 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 4281: 'Real smx = std::numeric_limits::min();' — smallest positive normal, not lowest(). For an all-negative s/rho neighborhood smx stays ~2.2e-308, so min(max(sn/rhon,smn),smx) at 4307 fails to bound at the true (negative) local max; a positive overshoot is clipped to ~0 instead. Reachable: ns.do_scalminmax=1 + ns.do_cons_trac=1 with negative tracer (call at 2927).

F030 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 4339 in ConvectiveScalMinMax: 'Real smx = std::numeric_limits::min();'. All-negative neighborhoods leave smx~+2.2e-308, so overshoots between the true max and 0 pass unclipped and positive overshoots clip to ~0 instead of the negative max (4365). Reachable via ns.do_scalminmax=1 with default NonConservative tracer/temperature (call at 2932). Claim's 'clips -0.5 to ~0' is only exact for positive sn, but the defect and both failure modes are real.

F072 — confirmed (one verifier lens)

Reported as: ConservativeScalMinMax seeds the neighborhood-maximum accumulator with std::numeric_limits::min() (smallest positive double ~2.2e-308) instead of lowest(), so the computed local max is wrong whenever all neighborhood values of s/rho are negative.

Failure scenario: ns.do_scalminmax=1 with a conservative scalar that is locally negative (signed tracer/temperature anomaly), EB or non-EB: in an all-negative 3^D neighborhood smx stays +2.2e-308, so min(max(sn/rhon,smn),smx) never clamps overshoots above the true local maximum -- the limiter silently fails exactly where redistribution/advection overshoots occur.

Lens 1 (refutation attempt): Same confirmed seed bug at 4281: numeric_limits::min() instead of lowest() in ConservativeScalMinMax. In an all-negative (vfrac>0) neighborhood smx stays +2.2e-308 so the upper clamp at 4307 never binds at the true local max; EB and non-EB paths identical. Reachable with ns.do_scalminmax=1 and a locally negative conservative scalar.

F076 — confirmed (one verifier lens)

Reported as: ConservativeScalMinMax (and ConvectiveScalMinMax at line 4339) initialize the stencil maximum accumulator with std::numeric_limits::min() — the smallest POSITIVE double (~2.2e-308) — instead of lowest(), so the computed local max is clamped below at ~0 and is wrong whenever all stencil values are negative.

Failure scenario: Run with ns.do_scalminmax=1 and a tracer (or temperature) field containing negative values, e.g. tracer in [-1,0]: after advection, an overshoot above the local 3x3(x3) max (say sn=-0.5 where the true stencil max is -1) is not clipped because min(sn, smx~2.2e-308) leaves it unchanged; the limiter silently fails to prevent overshoots.

Lens 1 (refutation attempt): Confirmed at both 4281 (ConservativeScalMinMax) and 4339 (ConvectiveScalMinMax): max accumulator seeded with std::numeric_limits::min() (~2.2e-308), not lowest(). For tracer in [-1,0], sn=-0.5 with true stencil max -1: min(max(-0.5,smn), 2.2e-308) = -0.5, overshoot not clipped — exactly the alleged silent limiter failure.

F073 — confirmed (one verifier lens)

Reported as: ConvectiveScalMinMax has the same defect: the max accumulator smx is initialized with std::numeric_limits::min() (tiny positive) instead of lowest(), giving a wrong neighborhood maximum for negative-valued scalars.

Failure scenario: ns.do_scalminmax=1 with a NonConservative scalar (default advectionType for tracer/temp) that takes negative values: for all-negative neighborhoods the upper clamp is never applied, so new local maxima created by advection overshoot are left uncorrected while undershoots are clamped, an asymmetric, silently broken limiter.

Lens 1 (refutation attempt): Line 4339: same min()-vs-lowest() seed in ConvectiveScalMinMax, used for NonConservative scalars (default tracer/temp advectionType) via 2930-2932 when ns.do_scalminmax=1. All-negative neighborhoods get smx~+2.2e-308: undershoots still clamped by smn but new maxima are never limited — the asymmetric failure as claimed.

---

_{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: F047, F031, F030, F072, F076, F073. Reviewer unit(s): NSB-3 syncinterp/vel-advance/particles, theme:eb-smallcell, 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.}