compute_divtau use_tensor_correction: size-1 divtau_scal indexed for all levels, leaks MultiFab, corrects level 0 only

[WeiqunZhang]

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

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

Locations: src/diffusion/incflo_diffusion.cpp:24, src/diffusion/incflo_diffusion.cpp:20

Based on commit 7307d872 (line numbers refer to that tree).

All six findings trace to the same few lines in the use_tensor_correction branch of incflo::compute_divtau: a scratch vector built for level 0 only, held by raw owning pointers, and post-processed only on level 0.

The defect

src/diffusion/incflo_diffusion.cpp:24 — In the use_tensor_correction branch, divtau_scal holds only level 0, but DiffusionScalarOp::compute_divtau loops lev=0..finestLevel(), indexing a_divtau[lev] out of bounds on multilevel runs; the Saxpy at line 40 also corrects only level 0.

src/diffusion/incflo_diffusion.cpp:24 — With use_tensor_correction, compute_divtau builds divtau_scal containing only a level-0 MultiFab but passes it to DiffusionScalarOp::compute_divtau, which indexes a_divtau[lev] for all levels 0..finest_level; the subsequent EB_set_covered and Saxpy also touch only level 0.

src/diffusion/incflo_diffusion.cpp:20 — In compute_divtau's use_tensor_correction path, divtau_scal is a Vector<MultiFab*> holding one leaked 'new MultiFab' for level 0 only, but DiffusionScalarOp::compute_divtau indexes a_divtau[lev] for all levels, and the Saxpy correction at line 40 is applied only on level 0.

src/diffusion/incflo_diffusion.cpp:20 — divtau_scal.push_back(new MultiFab(...)) is never deleted, leaking a full AMREX_SPACEDIM-component MultiFab on every compute_divtau call when use_tensor_correction is on.

src/diffusion/incflo_diffusion.cpp:24 — The use_tensor_correction branch of incflo::compute_divtau builds divtau_scal for level 0 only (grids[0], dmap[0]) and passes the size-1 vector to DiffusionScalarOp::compute_divtau, which loops lev=0..finest_level and indexes a_divtau[lev]; the Saxpy on line 40 also corrects only level 0.

src/diffusion/incflo_diffusion.cpp:20 — divtau_scal[0] is allocated with new in the use_tensor_correction branch of incflo::compute_divtau and never deleted, leaking one AMREX_SPACEDIM-component MultiFab per call.

Why it matters

F007: incflo.use_tensor_correction=1 (with implicit diffusion, as required) and amr.max_level>=1: first ApplyPredictor calls compute_divtau; DiffusionScalarOp::compute_divtau dereferences a_divtau[1] of a size-1 std::vector -> UB/segfault. Even if it survived, levels>0 would get the uncorrected tensor divtau.

F012: incflo.use_tensor_correction=1 (with implicit diffusion) and amr.max_level>=1: DiffusionScalarOp::compute_divtau dereferences a_divtau[1] past the end of a size-1 amrex::Vector, causing out-of-bounds access/crash in release builds; even absent a crash, divtau on levels >=1 never receives the tensor-minus-scalar correction.

F022: incflo.use_tensor_correction=true (with diffusion_type Implicit) and amr.max_level>=1: every ApplyPredictor (incflo_apply_predictor.cpp:148) calls compute_divtau; DiffusionScalarOp::compute_divtau dereferences a_divtau[1] on a size-1 vector -> out-of-bounds/segfault. Even single-level, one 3-component EB MultiFab leaks per call, every step (twice with MOL corrector).

F028: incflo.use_tensor_correction=1: compute_divtau runs every predictor (and corrector for MOL) step; each call leaks one velocity-sized MultiFab (Arena memory never freed). Long runs, especially on GPU where arena memory is scarce, grow until allocation failure. Same class as leak fixed in #147.

F055: incflo.use_tensor_correction=1 (allowed with Implicit diffusion, nothing gates AMR) with amr.max_level>=1: a_divtau[1] indexes past the end of the size-1 Vector inside DiffusionScalarOp::compute_divtau -> undefined behavior / crash on the first predictor step.

F072: incflo.use_tensor_correction=1: compute_divtau runs every predictor step (need_divtau() is true), leaking a full velocity-sized MultiFab each step; memory grows without bound over long runs.

How to reach it

• incflo.use_tensor_solve=0, incflo.use_tensor_correction=1, default incflo.diffusion_type=2, amr.max_level=1 with tagging that creates level 1 (any multilevel deck); first ApplyPredictor calls compute_divtau and goes out of bounds.
• amr.max_level=1, incflo.diffusion_type=2 (default), incflo.use_tensor_correction=true, any refining problem — e.g. test_3d/benchmark.upipe (already max_level=1) with use_tensor_solve=false replaced by use_tensor_correction=true.

Suggested fix

Make the branch multilevel and value-owned, mirroring the plain branches that already pass full finest_level+1 vectors into the ops. Replace the Vector<MultiFab*> holding a new MultiFab with a Vector<MultiFab> sized finestLevel()+1, defining each element on grids[lev], dmap[lev], *m_factory[lev] with divtau[lev]->nComp()/nGrow(), and pass GetVecOfPtrs(divtau_scal). This is required because both DiffusionTensorOp::compute_divtau and DiffusionScalarOp::compute_divtau unconditionally loop lev = 0..finestLevel() over their argument vectors (amrex::Vector::operator[] is unchecked std::vector indexing). Then loop the EB_set_covered calls — including the one on divtau[lev] after the tensor apply, which is also level-0-only today — and the MultiFab::Saxpy over all levels so every level gets the tensor-minus-scalar correction. Value semantics also fixes the per-call leak (same class as the leak fixed in #147). One judgement call for the maintainer: this makes multilevel + use_tensor_correction actually run for the first time, so its results are untested; if that combination is not meant to be supported, an amrex::Abort when finestLevel() > 0 would be the conservative alternative — but the leak fix is needed either way.

For src/diffusion/incflo_diffusion.cpp:24 (F007):

--- a/src/diffusion/incflo_diffusion.cpp
+++ b/src/diffusion/incflo_diffusion.cpp
@@ -12,19 +12,25 @@
     if (use_tensor_correction) {
 
         get_diffusion_tensor_op()->compute_divtau(divtau, vel, density, eta);
 #ifdef AMREX_USE_EB
-        EB_set_covered(*divtau[0]     , 0.0);
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            EB_set_covered(*divtau[lev], 0.0);
+        }
 #endif
 
-        Vector<MultiFab*> divtau_scal;
-        divtau_scal.push_back(new MultiFab(grids[0], dmap[0], divtau[0]->nComp(),
-                                           divtau[0]->nGrow(),MFInfo(),*m_factory[0]));
-        divtau_scal[0]->setVal(0.);
+        Vector<MultiFab> divtau_scal(finest_level+1);
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            divtau_scal[lev].define(grids[lev], dmap[lev], divtau[lev]->nComp(),
+                                    divtau[lev]->nGrow(), MFInfo(), *m_factory[lev]);
+            divtau_scal[lev].setVal(0.);
+        }
 
-        get_diffusion_scalar_op()->compute_divtau({divtau_scal}, vel, density, eta);
+        get_diffusion_scalar_op()->compute_divtau(GetVecOfPtrs(divtau_scal), vel, density, eta);
 #ifdef AMREX_USE_EB
-        EB_set_covered(*divtau_scal[0], 0.0);
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            EB_set_covered(divtau_scal[lev], 0.0);
+        }
 #endif
 
         // Define divtau to be (divtau_full - divtau_separate)
         if (m_verbose > 0)
@@ -37,7 +43,9 @@
         // amrex::Print() << "Z-comp: Norm of tensor apply vs scalar apply " <<
         //                    divtau[0]->norm0(2) << " " << divtau_scal[0]->norm0(2) << "\n";
 
-        MultiFab::Saxpy(*divtau[0], -1.0, *divtau_scal[0], 0, 0, AMREX_SPACEDIM, 0);
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            MultiFab::Saxpy(*divtau[lev], -1.0, divtau_scal[lev], 0, 0, AMREX_SPACEDIM, 0);
+        }
 
         // amrex::Print() << "X-comp: Norm of difference of tensor apply vs scalar apply " <<
         //                    divtau[0]->norm0(0) << "\n";

_{Allocates divtau_scal on every level 0..finest_level and loops EB_set_covered and Saxpy over all levels, matching DiffusionScalarOp::compute_divtau (DiffusionScalarOp.cpp:611-636), which indexes a_divtau[lev] for all levels. Same patch also resolves F022 (apply once). GetVecOfPtrs matches the Vector<MultiFab*> const& signature in DiffusionScalarOp.H:37.}

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

F007 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): incflo_diffusion.cpp:19-24: divtau_scal gets one element (grids[0], dmap[0]) then get_diffusion_scalar_op()->compute_divtau({divtau_scal},...). DiffusionScalarOp.cpp:611-631: 'int finest_level = m_incflo->finestLevel(); ... for (int lev = 0; lev <= finest_level; ++lev) { divtau_tmp[lev].define(a_divtau[lev]->boxArray(),...)' -> a_divtau[1] OOB on size-1 vector. init.cpp:110-118 has no max_level gate; predictor line 146-148 calls it unconditionally when use_tensor_correction.

Lens 2 (reachability/intent): incflo_diffusion.cpp:19-24 passes size-1 divtau_scal to DiffusionScalarOp::compute_divtau, which loops lev=0..finestLevel() and indexes a_divtau[lev] (DiffusionScalarOp.cpp:631/665/691 EB, 724 non-EB) — OOB when finestLevel()>=1. Only guards are init.cpp:112-118 (tensor_solve conflict, implicit-only); no max_level guard. Introduced level-0-only in 6854c35 (2020), undocumented, unused by any test deck.

F012 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): Same defect: DiffusionScalarOp.cpp indexes a_divtau[lev] for lev=0..finest_level (lines 631, 665, 691 EB branch; 724 non-EB branch) while divtau_scal has size 1. EB_set_covered(*divtau_scal[0]) at line 26 and Saxpy(*divtau[0],...) at line 40 touch only level 0, so levels>=1 never get the tensor-minus-scalar correction even if the OOB read survived.

Lens 2 (reachability/intent): incflo_diffusion.cpp:19-24 builds size-1 divtau_scal, but DiffusionScalarOp.cpp:611 sets finest_level=m_incflo->finestLevel() and dereferences a_divtau[lev] for lev<=finest_level at lines 631, 724, 740-743 (both EB and non-EB branches). Only guards on use_tensor_correction (init.cpp:112-117, moving-EB abort from d55f16a) never check max_level. Level-0-only code unchanged since introduction in 6854c35 (2020); no comment/assert documents a single-level restriction.

F022 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 20: 'divtau_scal.push_back(new MultiFab(grids[0], dmap[0], ...))' into a function-local Vector<MultiFab*>; no delete anywhere in compute_divtau. Predictor (incflo_apply_predictor.cpp:148) and corrector (incflo_apply_corrector.cpp:139-142, gated on '(Explicit) || use_tensor_correction') both call it, so MOL runs leak twice per step. Multilevel OOB part also confirmed (a_divtau[lev] loop in DiffusionScalarOp.cpp:617-724).

F028 — confirmed (one verifier lens)

Lens 1 (refutation attempt): The 'new MultiFab' at line 20 is stored only in a local Vector<MultiFab*> and never deleted; the function returns at line 48 with the pointer dropped. Allocated with MFInfo() and *m_factory[0] (EB factory), AMREX_SPACEDIM components, divtau[0]->nGrow() ghosts — one such MultiFab leaks per compute_divtau call, every predictor step (and corrector step for MOL).

F055 — confirmed (one verifier lens)

Lens 1 (refutation attempt): divtau_scal is built for level 0 only (line 20) and passed at line 24; DiffusionScalarOp::compute_divtau loops 'for (int lev = 0; lev <= finest_level; ++lev)' over a_divtau (DiffusionScalarOp.cpp:630-631, 723-724). init.cpp:116-117 only requires Implicit diffusion; nothing forbids amr.max_level>=1, so a_divtau[1] is indexed past the end on the first predictor step. Saxpy at line 40 corrects only *divtau[0].

F072 — confirmed (one verifier lens)

Lens 1 (refutation attempt): 'divtau_scal.push_back(new MultiFab(grids[0], dmap[0], divtau[0]->nComp(), divtau[0]->nGrow(), MFInfo(), *m_factory[0]))' at line 20; no matching delete or smart pointer in the function. need_divtau()||use_tensor_correction at incflo_apply_predictor.cpp:146 makes it run every step when use_tensor_correction=1, leaking one velocity-sized MultiFab per call.

---

_{Based on commit 7307d872, which is also the tree the audit verified against. From an automated audit of src/. Audit finding ids: F007, F012, F022, F028, F055, F072. Reviewer unit(s): Core-2 (predictor/corrector + updates), Diffusion+Rheology, theme:gpu-capture, theme:symmetry-twins. 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.}