MLALaplacian/MLEBABecLap singular-flag tests drifted from MLABecLaplacian (missing std::abs, #3926 CF check)

[WeiqunZhang]

Severity: medium/low · Category: correctness · Subsystem: LinearSolvers (blast radius 4/5)

Locations: Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:300, Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:304, Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:830

Based on commit 248fcfb7a5 (branch ai_audit; line numbers refer to that tree).

MLABecLaplacian::update_singular_flags is the reference implementation of the cell-centered singularity test; two fixes landed there (std::abs on the a-sum in #2477/300985a3e5, the Neumann coarse/fine-interface check in #3926) but were never mirrored into MLALaplacian::updateSingularFlag or MLEBABecLap::prepareForSolve/update.

The defect

Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:300 — updateSingularFlag tests asum <= amax*1e-12 without std::abs, unlike the sibling MLABecLaplacian (AMReX_MLABecLaplacian.H:820), so any negative sum of a-coefficients wrongly flags the level singular.

Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:304 — updateSingularFlag lacks the Neumann coarse/fine-interface singularity check that PR #3926 (a89b465) added to MLABecLaplacian::update_singular_flags and MLPoisson::prepareForSolve; MLALaplacian was never updated.

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:830 — Singularity test uses asum <= amax*1.e-12 without std::abs, unlike sibling MLABecLaplacian::update_singular_flags, so a negative-sum a coefficient wrongly marks the level singular; same drift duplicated in update() at line 1469.

Why it matters

F209: Any DIM/backend/precision. App: MLALaplacian on grids covering the domain, all Neumann/periodic domain BCs, setScalars(-1.0, 1.0), setACoeffs with a<0 (so alphaa>0, operator positive definite), MLMG::solve. asum<0 satisfies asum<=amax1e-12 -> isSingular true -> MLMG's makeSolvable subtracts the RHS mean -> silently wrong solution to a well-posed problem.

F210: Any DIM/backend. App: MLALaplacian whose coarsest solve level does not cover the domain, setCoarseFineBC(&crse, ratio, LinOpBCType::Neumann) (generic MLLinOp/MLMGBndry machinery from #3926 accepts it), all Neumann/periodic domain BCs, setScalars(0.0, 1.0). Problem is singular, but m_domain_covered[0] is false so m_is_singular[0] stays false -> MLMG skips singular handling -> residual stalls / wrong answer.

F405: EB build, single fully-covering AMR level, all-Neumann/periodic domain BCs, no EB Dirichlet, setScalars(a!=0,b), setACoeffs with negative a (e.g. implicit negative reaction coefficient): asum<0 <= amax*1e-12 -> m_is_singular=true although the operator is nonsingular; MLMG's makeSolvable then subtracts the RHS mean, silently shifting the converged solution. MLABecLaplacian (AMReX_MLABecLaplacian.H:820 uses std::abs(asum)) gives the correct answer for the identical problem.

How to reach it

• F209: config: Any AMReX_SPACEDIM/PRECISION, AMReX_LINEAR_SOLVERS=ON, CPU or GPU, MPI optional. Compiled by all CI jobs (MLMG.cpp instantiation); no CI test executes MLALaplacian::updateSingularFlag.
• F209: API path: MLALaplacian op({geom},{ba},{dm}); // grids cover domain op.setDomainBC(all Neumann or Periodic); op.setLevelBC(0,nullptr); op.setScalars(1.0,-1.0); acoef.setVal(-c); op.setACoeffs(0,acoef); // Lap phi - c phi convention, positive-definite up to sign MLMG mlmg(op); mlmg.solve({&sol},{&rhs},tol,0); // updateSingularFlag: asum=-cN <= amax1e-12 -> singular -> makeSolvable subtracts mean(rhs) -> wrong converged answer Public API, no in-tree caller with negative a.
• F209: test coverage: none - no Tests/ program instantiates MLALaplacian; MLPoisson::makeNLinOp uses alpha in {0,1e30*dx^-2} (asum>0) and is unaffected.
• F210: config: Any AMReX_SPACEDIM/PRECISION, AMReX_LINEAR_SOLVERS=ON, CPU or GPU, MPI optional. Compiled by every CI job; no CI test runs MLALaplacian with a coarse/fine BC.
• F210: API path: // level-0 grids strictly inside the domain (m_domain_covered[0]==false) MLALaplacian op({geom},{ba_sub},{dm}); op.setDomainBC(all Neumann/Periodic); op.setCoarseFineBC(&crse, 2, LinOpBCType::Neumann); op.setLevelBC(0,&bcdata); op.setScalars(0.0,1.0); MLMG mlmg(op); mlmg.solve({&sol},{&rhs},tol,0); // isSingular(0)==false for a pure-Neumann Poisson -> no makeSolvable -> stall/'MLMG failed.' for rhs with nonzero mean Public API; no in-tree caller.
• F210: test coverage: none - no Tests/ program uses MLALaplacian; MLPoisson::makeNLinOp (only Src/ caller) sets no coarse/fine BC type and has a!=0.
• F405: config: AMReX_EB=ON, SPACEDIM 2 or 3 (EB builds in macos.yml/cuda.yml). Also reachable with an all-regular EB factory.
• F405: API path: MLEBABecLap op(geom,grids,dmap,info,factory) on a single domain-covering level; op.setDomainBC(all Neumann/Periodic); op.setLevelBC(0,nullptr); op.setScalars(1.0,b); op.setACoeffs(0,-1.0); MLMG(op).solve -> prepareForSolve marks level singular -> makeSolvable shifts RHS. Public API; no in-tree caller passes negative a.
• F405: test coverage: none — Tests/LinearSolvers/CellEB2 uses positive a and Dirichlet BCs; no test asserts on m_is_singular for negative a.

Suggested fix

Bring the three copies back in line with MLABecLaplacianT<MF>::update_singular_flags (AMReX_MLABecLaplacian.H:806-855), which is the path that has been kept correct: 1. std::abs(asum) <= amax * 1e-12 in MLALaplacian.H:300 and in both MLEBABecLap.cpp copies (prepareForSolve, line 830; update, line 1469). The test is meant to ask "is a effectively zero", and only the twin says so; a negative-sum a is a legitimate, nonsingular operator that today gets its RHS mean subtracted by MLMG::makeSolvable. These are one-line, self-contained changes and can land first, in one PR. 2. Port the #3926 block (ABecLap lines 826-855) into MLALaplacian::updateSingularFlag so a level-0 grid that does not cover the domain but has a Neumann coarse/fine BC is recognised as singular. MLLinOp/MLCellLinOp already accept LinOpBCType::Neumann from setCoarseFineBC for any cell-centered operator and getSolvabilityOffset already handles it, so only the flag is missing. Maintainer decisions: - MLALaplacian never consults m_overset_mask in its domain-covered loop, so the sketch omits the twin's AMREX_ASSERT(m_overset_mask[0][0]==nullptr) and the !m_overset_mask guard. If MLALaplacian is meant to support overset masks, add both as in ABecLap. - MLEBABecLap also has no #3926 Neumann coarse/fine block; that was not verified here (EB Dirichlet interplay), but worth checking in the same pass. - Open PRs #4930/#4922 touch MLEBABecLap.cpp; the verifier did not confirm they address this, but check for overlap before landing.

For Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:300 (F209):

--- a/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H
@@ -297,7 +297,7 @@ MLALaplacianT<MF>::updateSingularFlag ()
                     // are similar.
                     RT asum = m_a_coeffs[alev].back().sum(0,IntVect(0));
                     RT amax = m_a_coeffs[alev].back().norminf(0,1,IntVect(0));
-                    m_is_singular[alev] = (asum <= amax * RT(1.e-12));
+                    m_is_singular[alev] = (std::abs(asum) <= amax * RT(1.e-12));
                 }
             }
         }

_{Matches the twin MLABecLaplacianT::update_singular_flags (AMReX_MLABecLaplacian.H:820, fixed by #2477) so a negative a-coefficient sum no longer flags the level singular. std::abs is already available via the included AMReX headers. MLEBABecLap.cpp:830,1469 share the same missing abs and could be fixed identically. Cross-reference: F405 is the same missing-std::abs pattern in MLEBABecLap.cpp:830/1469; land both in one PR.}

For Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H:304 (F210):

--- a/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLALaplacian.H
@@ -302,6 +302,35 @@ MLALaplacianT<MF>::updateSingularFlag ()
             }
         }
     }
+
+    if (!m_is_singular[0] && this->m_needs_coarse_data_for_bc &&
+        this->m_coarse_fine_bc_type == BCType::Neumann)
+    {
+        bool lev0_a_is_zero = false;
+        if (m_a_scalar == RT(0.0)) {
+            lev0_a_is_zero = true;
+        } else {
+            RT asum = m_a_coeffs[0].back().sum(0,IntVect(0));
+            RT amax = m_a_coeffs[0].back().norminf(0,1,IntVect(0));
+            bool a_is_almost_zero = std::abs(asum) <= amax * RT(1.e-12);
+            if (a_is_almost_zero) { lev0_a_is_zero = true; }
+        }
+
+        if (lev0_a_is_zero) {
+            auto bbox = this->m_grids[0][0].minimalBox();
+            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+                if (this->m_lobc[0][idim] == BCType::Dirichlet) {
+                    bbox.growLo(idim,1);
+                }
+                if (this->m_hibc[0][idim] == BCType::Dirichlet) {
+                    bbox.growHi(idim,1);
+                }
+            }
+            if (this->m_geom[0][0].Domain().contains(bbox)) {
+                m_is_singular[0] = true;
+            }
+        }
+    }
 }
 
 template <typename MF>

_{Ports the #3926 Neumann coarse/fine-interface singularity block from MLABecLaplacianT::update_singular_flags (AMReX_MLABecLaplacian.H:826-855) into MLALaplacian, using the file's BCType alias. The overset-mask assert is omitted because MLALaplacian never checks m_overset_mask in its domain-covered loop either; maintainer should confirm MLALaplacian is not meant to support overset masks (otherwise add the assert and the !m_overset_mask check as in the twin).}

For Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:830 (F405):

--- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp
@@ -827,7 +827,7 @@
                 {
                     Real asum = m_a_coeffs[alev].back().sum();
                     Real amax = m_a_coeffs[alev].back().norm0();
-                    m_is_singular[alev] = (asum <= amax * 1.e-12);
+                    m_is_singular[alev] = (std::abs(asum) <= amax * 1.e-12);
                 }
             }
         }
@@ -1466,7 +1466,7 @@
                 {
                     Real asum = m_a_coeffs[alev].back().sum();
                     Real amax = m_a_coeffs[alev].back().norm0();
-                    m_is_singular[alev] = (asum <= amax * 1.e-12);
+                    m_is_singular[alev] = (std::abs(asum) <= amax * 1.e-12);
                 }
             }
         }

_{Matches the non-EB twin MLABecLaplacianT::update_singular_flags (AMReX_MLABecLaplacian.H:820, std::abs(asum)), fixed there by commit 300985a but missed here. Applied to both copies (prepareForSolve and update). std::abs is already available via the included AMReX_MLABecLaplacian.H. Open PRs #4930/#4922 touch this area. Cross-reference: F209 is the same pattern in MLALaplacian.H:300.}

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

F209 — confirmed (one verifier lens)

Lens 1 (refutation attempt): AMReX_MLALaplacian.H:300 m_is_singular[alev] = (asum <= amax * RT(1.e-12)); vs MLABecLaplacian.H:820 (std::abs(asum) <= amax * RT(1.e-12)). The abs was added by afe587b (#2477, message: 'Fix a bug in the detection of singularity') only in ABecLap; git log -S asum -- MLALaplacian.* shows no change since 2da11d4/57793fda87. MLEBABecLap.cpp:830,1469 share the missing abs. Consequence: MLMG.H:1519 if (linop.isSingular(0) && getEnforceSingularSolvable()) makeSolvable(); -> MLCellLinOp.H:1933 rhs.plus(-offset[c], c, 1) silently changes a well-posed problem. No assertion forbids negative a-coefficients; no fix on development, no open issue/PR.

F210 — confirmed (one verifier lens)

Lens 1 (refutation attempt): AMReX_MLALaplacian.H:278-305 updateSingularFlag only loops if (this->m_domain_covered[alev]) and ends; no block for m_needs_coarse_data_for_bc && m_coarse_fine_bc_type == LinOpBCType::Neumann, unlike MLABecLaplacian.H:826-855 and MLPoisson.H:231-247 added by a89b465 (#3926) whose --stat lists ABecLap, CellLinOp, LinOp, MLMG, MLMGBndry, Poisson but not MLALaplacian. The generic path accepts Neumann CF for any cell-centered op: MLLinOp.H:1734-1740 setCoarseFineBC stores bc_type unconditionally, MLCellLinOp.H:878 passes m_coarse_fine_bc_type to the level-0 bndry, MLCellLinOp.H:2336 getSolvabilityOffset already sizes by grid numPts for non-Dirichlet CF. MLLinOp.H:1200 m_domain_covered[0]=false when grids do not fill the domain, so m_is_singular[0] stays false. No fix on development, no open issue/PR.

F405 — confirmed (one verifier lens)

Lens 1 (refutation attempt): MLEBABecLap.cpp:830 and :1469 m_is_singular[alev] = (asum <= amax * 1.e-12); (no std::abs). Sibling MLABecLaplacian.H:820 uses std::abs(asum) <= amax*RT(1.e-12). Commit 300985a (2021-11-17, "Fix a bug in the detection of singularity") added std::abs only to AMReX_MLABecLaplacian.cpp (stat: 2 files, both MLABecLaplacian) — the EB sibling was missed; later 688ab79/#3926 kept std::abs only in MLABecLaplacian. MLMG.H:1519 if (linop.isSingular(0) && linop.getEnforceSingularSolvable()) makeSolvable(); subtracts the RHS mean, so a negative-a nonsingular operator is silently altered. No dev commit or PR repairs it.

Possibly related upstream: dedup flagged open PRs #4930 #4922 as touching this file (verifier did not match it).

---

_{Based on commit 248fcfb7a5, the tree the audit verified against. From an automated audit of AMReX Src/. Audit finding ids: F209, F210, F405. Reviewer unit(s): LS/ALap-Poisson-1, LS/EBABecLap-1. Nothing here was compiled or run except where the evidence says so — the failure scenarios are code reasoning, so the reaching configuration above is the cheapest way to confirm or refute it.}