MLEBABecLap::averageDownEBPhi and MLEBTensorOp kappa average-down use MG ratio 2 across AMR levels with ref ratio 4

[WeiqunZhang]

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

Locations: Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:682, Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp:157

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

Both sites average a per-AMR-level field that exists only at MG level 0 down to the next coarser AMR level using mg_coarsen_ratio (2), but the true ratio between those two containers is the AMR refinement ratio, which MLLinOp::defineGrids explicitly allows to be 4.

The defect

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:682 — averageDownEBPhi averages inhomogeneous EB Dirichlet values between AMR levels with the fixed MG ratio 2 instead of the actual AMR refinement ratio, which can be 4.

Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp:157 — Cross-AMR-level averaging of m_kappa (and m_eb_kappa at lines 180-182) passes IntVect(mg_coarsen_ratio)=2, but m_kappa has only one MG level (kappa_num_mglevs=1), so m_kappa[amrlev].back() is at the fine level's finest resolution and the true ratio to m_kappa[amrlev-1].front() is AMRRefRatio(amrlev-1), which can be 4.

Why it matters

F217: Any build with EB: 2 AMR levels with ref ratio 4 (supported: mlmg_eb_cc_interp_r<4> exists), MLEBABecLap + setEBDirichlet(phi,beta) on both levels, MLMG::solve. m_eb_phi[1] (fine, mglev-0 resolution) is coarsened by 2 into a temp whose index space is still 2x finer than m_eb_phi[0]; crse.ParallelCopy(temp) then copies index-overlapping data from wrong physical regions (or nothing), silently corrupting the coarse EB Dirichlet values and the composite solution.

F223: 2D/3D, EB on, multi-level MLMG solve with AMR refinement ratio 4 (supported: MLMG itself uses AMRRefRatioVect for residual averaging; defineGrids builds the extra MG level for it) and spatially varying bulk viscosity. EB_average_down_faces coarsens fine kappa by 2 into an intermediate index space, then ParallelCopy into the coarse-level MultiFab overwrites coarse kappa with data from wrong physical locations -> silently wrong coarse operator.

How to reach it

• F217: config: AMReX_EB=ON, SPACEDIM 2 or 3, MPI optional. EB builds in macos.yml/cuda.yml; no CI test runs ratio-4 EB multilevel.
• F217: API path: Geometry/BoxArray with ref ratio 4 between levels 0,1; MLEBABecLap op(geom,grids,dmap,info,factories); op.setEBDirichlet(0,phi0,beta); op.setEBDirichlet(1,phi1,beta); op.setScalars/ACoeffs/BCoeffs; MLMG(op).solve(...) -> prepareForSolve -> averageDownEBPhi corrupts m_eb_phi[0]. Public API; no in-tree caller uses ratio 4 with EB Dirichlet.
• F217: test coverage: Tests/LinearSolvers/CellEB2 has ref_ratio (default 2, MyTest.H:29) and setEBDirichlet (MyTest.cpp:96) but never runs ratio 4 and asserts nothing on the coarse EB phi.
• F223: config: AMReX_SPACEDIM=2 or 3, AMReX_EB=ON, AMReX_LINEAR_SOLVERS=ON (default); AMReX_EB=ON is built in gcc.yml/clang.yml/cuda.yml/hip.yml, but no workflow runs a multi-level tensor solve.
• F223: API path: Vector geom(2) with geom[1].Domain()==refine(geom[0].Domain(),4); MLEBTensorOp op(geom,grids,dmap,LPInfo(),{fact0,fact1}); op.setDomainBC/setLevelBC(0..1); op.setShearViscosity(lev,eta); op.setBulkViscosity(lev,spatially-varying face kappa); op.setEBBulkViscosity(lev,kappa_eb); MLMG(op).solve(...) -> prepareForSolve line 157. Public API; no in-tree multi-level tensor caller (applications only).
• F223: test coverage: none. Tests/LinearSolvers/EBTensor/MyTest.cpp:51 builds a single-level MLEBTensorOp; no Tests/ program performs a multi-level tensor solve or asserts on coarse kappa.

Suggested fix

The rule that decides this: mg_coarsen_ratio is only correct across AMR levels when the fine-side operand is .back() of a coefficient vector that defineGrids has already coarsened down to twice the coarse resolution. That is why MLEBABecLap::averageDownCoeffsToCoarseAmrLevel (m_a_coeffs[flev].back(), m_b_coeffs[flev].back()) is correct as written for ratio 4. Neither m_eb_phi (one entry per AMR level, mglev 0 only) nor m_kappa/m_eb_kappa (kappa_num_mglevs = 1, so .back() is the full-resolution fine level) has that intermediate level, so the cross-AMR call must use the AMR ratio, matching averageDownSolutionRHS (AMRRefRatio(camrlev)) in the same file and MLCellLinOp's use of AMRRefRatioVect for sol/res. Change the four cross-AMR calls to AMRRefRatioVect(amrlev-1): averageDownEBPhi in MLEBABecLap.cpp, the two if (amrlev > 0) branches in MLEBTensorOp::prepareForSolve, and the same branch in the non-EB twin MLTensorOp::prepareForSolve (lines 138-140, average_down_faces has the identical index-space mismatch). The intra-AMR mglev loops keep mg_coarsen_ratio. Both pieces are independent one-liners and can land in one PR. Maintainer decision: if kappa_num_mglevs is ever raised above 1 again, the cross-AMR ratio must become AMRRefRatio / 2^(size-1); a short comment or an assert that m_kappa[amrlev].size() == 1 would guard that. Open PRs #4930 and #4922 touch MLEBABecLap.cpp and may be related (not verified). A CellEB2-style test with ref_ratio = 4 and inhomogeneous EB Dirichlet would cover the first site.

For Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp:682 (F217):

--- a/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp
@@ -680,7 +680,7 @@
     if (m_eb_phi[0]) {
         for (int amrlev = m_num_amr_levels-1; amrlev > 0; --amrlev) {
             amrex::EB_average_down_boundaries(*m_eb_phi[amrlev], *m_eb_phi[amrlev-1],
-                                              mg_coarsen_ratio, 0);
+                                              AMRRefRatio(amrlev-1), 0);
         }
     }
 }

_{Uses the AMR refinement ratio, matching averageDownSolutionRHS (line 1117, AMRRefRatio(camrlev)) in the same file. m_eb_phi exists only at mglev 0 so, unlike the coefficient average-downs that use .back(), the intermediate MG level is unavailable. The int overload forwards to the IntVect one and eb_avgdown_boundaries loops over the ratio, so 4 is handled. Open PRs #4930/#4922 touch this area.}

For Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp:157 (F223):

--- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
@@ -156,7 +156,7 @@ MLEBTensorOp::prepareForSolve ()
             if (amrlev > 0) {
                 amrex::EB_average_down_faces(GetArrOfConstPtrs(m_kappa[amrlev  ].back()),
                                              GetArrOfPtrs     (m_kappa[amrlev-1].front()),
-                                             IntVect(mg_coarsen_ratio), m_geom[amrlev-1][0]);
+                                             AMRRefRatioVect(amrlev-1), m_geom[amrlev-1][0]);
             }
         }
     } else {
@@ -179,7 +179,7 @@ MLEBTensorOp::prepareForSolve ()
             if (amrlev > 0) {
                 amrex::EB_average_down_boundaries(m_eb_kappa[amrlev  ].back(),
                                                   m_eb_kappa[amrlev-1].front(),
-                                                  IntVect(mg_coarsen_ratio), 0);
+                                                  AMRRefRatioVect(amrlev-1), 0);
             }
         }
     } else {
--- a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp
@@ -137,7 +137,7 @@ MLTensorOp::prepareForSolve ()
             if (amrlev > 0) {
                 amrex::average_down_faces(GetArrOfConstPtrs(m_kappa[amrlev  ].back()),
                                           GetArrOfPtrs     (m_kappa[amrlev-1].front()),
-                                          IntVect(mg_coarsen_ratio), m_geom[amrlev-1][0]);
+                                          AMRRefRatioVect(amrlev-1), m_geom[amrlev-1][0]);
             }
         }
     } else {

_{Cross-AMR-level averaging of m_kappa/m_eb_kappa now uses the actual AMR ratio (AMRRefRatioVect(amrlev-1), the convention MLCellLinOp uses for EB_average_down of sol/res). m_kappa has one MG level so .back() is at full fine resolution; MLEBABecLap's m_b_coeffs.back() is fine/2 so its ratio-2 call is correct and untouched. Same defect fixed in twin MLTensorOp.cpp. If kappa_num_mglevs ever changes, divide ratio by 2^(size-1).}

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

F217 — confirmed (one verifier lens)

Lens 1 (refutation attempt): averageDownEBPhi (lines 678-685): EB_average_down_boundaries(*m_eb_phi[amrlev], *m_eb_phi[amrlev-1], mg_coarsen_ratio, 0) with mg_coarsen_ratio=2 (MLLinOp.H:857). m_eb_phi is per-amrlev at mglev 0 only. MLLinOp::defineGrids (MLLinOp.H:1150-1170) accepts AMR ratio 4 by giving the fine level 2 mg levels and m_amr_ref_ratio=4; coefficient averaging correctly uses m_a_coeffs[flev].back() (already /2) but m_eb_phi is not coarsened per mg level. EB_average_down_boundaries (EBMultiFabUtil.cpp:669-676) takes the non-MFIter-safe path: cba=fine BA coarsened by 2 then crse.ParallelCopy(ctmp) — index-space mismatch copies wrong-region/zero data into coarse cut cells. Ratio 4 is supported for EB cell-centered ops: MLCellLinOp.H:1447 mlmg_eb_cc_interp_r<4>; averageDownSolutionRHS uses AMRRefRatio.

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

F223 — confirmed (one verifier lens)

Lens 1 (refutation attempt): MLEBTensorOp.cpp:47 m_kappa[amrlev].resize(std::min(kappa_num_mglevs,NMGLevels(amrlev))) with kappa_num_mglevs=1 (since 618ef09), so line 157 m_kappa[amrlev].back() is at fine resolution, yet lines 159/182 pass IntVect(mg_coarsen_ratio)=2. MLLinOp.H:1155-1175 explicitly supports AMR ratio 4 by adding one MG level (for (int i = 0; i < 2; ++i)), and MLEBABecLap.cpp:779-792 relies on that: it coarsens m_b_coeffs[flev].back() (fine/2 res) by 2. EBMultiFabUtil.cpp:604-614 builds ctmp=coarsen(fineBA,2) then crse->ParallelCopy(ctmp) into a fine/4 BoxArray: index-space mismatch overwrites wrong (incl. uncovered) coarse faces; same in EB_average_down_boundaries:670-675. No fix on development, not in KNOWN_STATE.

---

_{Based on commit 248fcfb7a5, the tree the audit verified against. From an automated audit of AMReX Src/. Audit finding ids: F217, F223. Reviewer unit(s): LS/EBABecLap-1, LS/EBTensor-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.}