DiffusionScalarOp: EB Dirichlet BCs inconsistent across solve/apply ops (ntrac>=2 abort)

[WeiqunZhang]

Severity: medium · Category: correctness, physics-numerics · Fix order: 11 of 29 — fix this 11th.

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

Locations: src/diffusion/DiffusionScalarOp.cpp:230, src/diffusion/DiffusionScalarOp.cpp:554, src/diffusion/DiffusionScalarOp.cpp:555, src/diffusion/DiffusionScalarOp.cpp:639

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

All five findings reduce to one rule applied inconsistently in this file: each MLEBABecLap here is built with ncomp=1, so EB Dirichlet phi and beta must be per-component aliases, and each apply op must carry the same EB BC as its solve twin.

The defect

src/diffusion/DiffusionScalarOp.cpp:230 — diffuse_scalar passes the full m_ntrac-component eta MultiFab as the beta argument of setEBDirichlet on a 1-component MLEBABecLap, instead of the comp-th component being solved.

src/diffusion/DiffusionScalarOp.cpp:554 — compute_laps never sets an EB Dirichlet BC on m_eb_scal_apply_op (defaults to homogeneous Neumann), even though diffuse_scalar applies setEBDirichlet from tracer_eb/temperature_eb for the same fields, so the explicit diffusion term uses a different EB BC than the implicit solve.

src/diffusion/DiffusionScalarOp.cpp:555 — compute_laps applies m_eb_scal_apply_op without ever calling setEBDirichlet, so the explicit tracer/temperature diffusion term sees homogeneous-Neumann EB, while diffuse_scalar (line 228-231) applies the configured tracer_eb/temperature_eb Dirichlet values in the implicit solve.

src/diffusion/DiffusionScalarOp.cpp:639 — DiffusionScalarOp::compute_divtau always calls setEBHomogDirichlet, ignoring hasEBFlow(), while both diffuse_vel_components (lines 382-387) and DiffusionTensorOp::compute_divtau (lines 262-266) use the EB inflow velocity via setEBDirichlet/setEBShearViscosityWithInflow.

src/diffusion/DiffusionScalarOp.cpp:639 — DiffusionScalarOp::compute_divtau always calls setEBHomogDirichlet (u=0 on EB), ignoring EB inflow velocity, while the matching implicit solve diffuse_vel_components (line 382-387) and both DiffusionTensorOp ops use setEBDirichlet/setEBShearViscosityWithInflow(velocity_eb) when hasEBFlow().

Why it matters

F029: EB build, incflo.ntrac>=2, eb_flow.tracer specified (mandatory whenever eb_flow is enabled with tracer advection), diff_type implicit or CN: MLEBABecLap::setEBDirichlet hits AMREX_ALWAYS_ASSERT(beta_ncomp == 1 || beta_ncomp == ncomp) (1==getNComp(), beta_ncomp==m_ntrac) and aborts. With unequal mu_s it would otherwise use tracer-0 diffusivity for every tracer's EB flux.

F031: EB build with eb_flow.tracer (or eb_flow.temperature) specified and diff_type=Explicit: tracer/temperature diffusion is laps-only, so the user-specified EB Dirichlet value is silently ignored and no diffusive flux crosses the EB (e.g., a heated EB never heats the fluid). With CN, predictor laps_o uses Neumann while the half-implicit solve uses Dirichlet, giving an inconsistent O(1) flux error at the EB.

F045: EB build with eb_flow.tracer (or temperature) set (intended 'Dirichlet on EB', incflo_arrays.cpp:43) and mu_s>0: with diffusion_type=0 the EB Dirichlet value never enters diffusion (zero EB flux); with Crank-Nicolson the two half-operators use different EB BCs, so the tracer/temperature boundary layer at the EB is wrong.

F030: EB build with eb_flow.velocity set plus incflo.use_tensor_solve=false and explicit/CN diffusion (or use_tensor_correction=1): the explicit divtau is computed with u=0 on the EB instead of the inflow velocity, so the viscous term near the EB inflow is wrong and inconsistent with the implicit solve of the same step; wrong shear layers near the EB.

F044: EB inflow (eb_flow.velocity set) with use_tensor_solve=false (or use_tensor_correction) and explicit/Crank-Nicolson diffusion or godunov_include_diff_in_forcing: explicit divtau near the EB is computed as if the EB were a stationary no-slip wall, inconsistent with the implicit half/tensor path, giving wrong viscous stresses in cut cells at the inflowing EB.

Suggested fix

In diffuse_scalar, alias component comp of eta exactly as line 229 already does for phi: MLEBABecLap::setEBDirichlet asserts beta_ncomp == 1 || beta_ncomp == ncomp, so the full tra_eta aborts for ntrac>=2 (and would silently use tracer-0 diffusivity otherwise). In compute_laps, take an eb_dirichlet vector like diffuse_scalar's, plumb get_tracer_eb()/get_temperature_eb() through incflo::compute_laps/compute_laps_T, and call setEBDirichlet per component when non-empty — keeping the default homogeneous Neumann otherwise, so apply matches the solve at lines 228-231. In compute_divtau, mirror diffuse_vel_components lines 382-387 (hasEBFlow() ? setEBDirichlet(velocity_eb comp alias) : setEBHomogDirichlet), moved inside the comp loop since velocity_eb differs per component. Maintainer decision: the scalar solve/apply ops are shared between tracer and temperature passes, and MLEBABecLap has no way to revert EB Dirichlet back to Neumann; if only one of tracer_eb/temperature_eb is defined, the later pass would see stale EB phi — decide whether to require both together or reset/rebuild the op.

For src/diffusion/DiffusionScalarOp.cpp:230 (F029):

--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -227,5 +227,6 @@
 
                 if (!eb_dirichlet[lev]->empty()) {
                     MultiFab phi(*eb_dirichlet[lev], amrex::make_alias, comp, 1);
-                  m_eb_scal_solve_op->setEBDirichlet(lev, phi, *eta[lev]);
+                    MultiFab beta(*eta[lev], amrex::make_alias, comp, 1);
+                    m_eb_scal_solve_op->setEBDirichlet(lev, phi, beta);
                 } // else use default homogeneous Neumann on EB

_{Aliases the comp-th component of eta as the 1-component EB beta, matching the existing per-component phi alias, so MLEBABecLap::setEBDirichlet's beta_ncomp assert passes for ntrac>=2 and each tracer's EB flux uses its own diffusivity. The MultiFab alias constructor accepts const sources (already used at line 535).}

For src/diffusion/DiffusionScalarOp.cpp:554 (F031):

--- a/src/diffusion/DiffusionScalarOp.H
+++ b/src/diffusion/DiffusionScalarOp.H
@@ -29,10 +29,11 @@
                          amrex::Vector<amrex::MultiFab const*> const& eta,
                          amrex::Real dt);
 
     void compute_laps (amrex::Vector<amrex::MultiFab*> const& laps,
                        amrex::Vector<amrex::MultiFab const*> const& a_scalar,
                        amrex::Vector<amrex::MultiFab const*> const& eta,
+                       amrex::Vector<amrex::MultiFab*> const& eb_dirichlet,
                        amrex::Vector<amrex::BCRec> bcrec);
 
     void compute_divtau (amrex::Vector<amrex::MultiFab*> const& a_divtau,
                          amrex::Vector<amrex::MultiFab const*> const& a_vel,
--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -485,9 +485,10 @@
     }
 }
 
 void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
                                       Vector<MultiFab const*> const& a_scalar,
                                       Vector<MultiFab const*> const& a_eta,
+                                      Vector<MultiFab*> const& eb_dirichlet,
                                       amrex::Vector<amrex::BCRec> bcrec)
 {
     BL_PROFILE("DiffusionScalarOp::compute_laps");
@@ -533,6 +534,12 @@
             for (int lev = 0; lev <= finest_level; ++lev) {
                 laps_comp.emplace_back(laps_tmp[lev],amrex::make_alias,comp,1);
                 scalar_comp.emplace_back(*a_scalar[lev],amrex::make_alias,comp,1);
+
+                if (!eb_dirichlet[lev]->empty()) {
+                    MultiFab phi_eb(*eb_dirichlet[lev], amrex::make_alias, comp, 1);
+                    MultiFab beta_eb(*a_eta[lev], amrex::make_alias, eta_comp, 1);
+                    m_eb_scal_apply_op->setEBDirichlet(lev, phi_eb, beta_eb);
+                } // else use default homogeneous Neumann on EB
 
                 Array<MultiFab,AMREX_SPACEDIM>
                     b = m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]);
@@ -572,6 +579,7 @@
     else
 #endif
     {
+        amrex::ignore_unused(eb_dirichlet);
         // We want to return div (mu grad)) phi
         m_reg_scal_apply_op->setScalars(0.0, -1.0);
 
--- a/src/diffusion/incflo_diffusion.cpp
+++ b/src/diffusion/incflo_diffusion.cpp
@@ -59,16 +59,17 @@
                      Vector<MultiFab const*> const& scalar,
                      Vector<MultiFab const*> const& eta)
 {
-    get_diffusion_scalar_op()->compute_laps(laps, scalar, eta,
+    get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, get_tracer_eb(),
                                             get_tracer_bcrec());
 
 }
 
 void
 incflo::compute_laps_T(Vector<MultiFab      *> const& laps,
                        Vector<MultiFab const*> const& scalar,
                        Vector<MultiFab const*> const& eta)
 {
     get_diffusion_scalar_op()->compute_laps(laps, scalar, eta,
+                                            get_temperature_eb(),
                                             get_temperature_bcrec());
 }

_{Adds eb_dirichlet arg to compute_laps mirroring diffuse_scalar (PR #148); wrappers pass get_tracer_eb()/get_temperature_eb(). Beta is eta comp-aliased to satisfy MLEBABecLap::setEBDirichlet's beta-ncomp assert (AMReX_MLEBABecLap.H:149). Caveat: m_eb_scal_apply_op is shared by tracer/temperature and MLEBABecLap cannot unset EB Dirichlet, so an empty *_eb call after a non-empty one keeps stale values — same pre-existing limitation as the shared solve op.}

For src/diffusion/DiffusionScalarOp.cpp:639 (F030):

--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -635,9 +635,6 @@
             divtau_tmp[lev].setVal(0.0);
         }
 
-        for (int lev = 0; lev <= finest_level; ++lev) {
-            m_eb_vel_apply_op->setEBHomogDirichlet(lev, *a_eta[lev]);
-        }
         // We want to return div (mu grad)) phi
         m_eb_vel_apply_op->setScalars(0.0, -1.0);
 
@@ -664,6 +661,13 @@
             for (int lev = 0; lev <= finest_level; ++lev) {
                 divtau_single.emplace_back(divtau_tmp[lev],amrex::make_alias,comp,1);
                 vel_single.emplace_back(       vel[lev],amrex::make_alias,comp,1);
 
+                if (m_incflo->hasEBFlow()) {
+                    MultiFab phi(*m_incflo->get_velocity_eb()[lev], amrex::make_alias, comp, 1);
+                    m_eb_vel_apply_op->setEBDirichlet(lev, phi, *a_eta[lev]);
+                } else {
+                    m_eb_vel_apply_op->setEBHomogDirichlet(lev, *a_eta[lev]);
+                }
+
                 if ( m_incflo->m_has_mixedBC ) {
                     auto const robin = m_incflo->make_robinBC_MFs(lev, &vel_single[lev]);

_{Mirrors the diffuse_vel_components twin (lines 382-387): per component, alias comp of velocity_eb and setEBDirichlet when hasEBFlow(), else setEBHomogDirichlet; the hoisted homogeneous loop is removed since setEBDirichlet needs the per-comp phi. a_eta is 1-component so *a_eta[lev] as beta matches the twin; setEBDirichlet copies per call, so per-comp overwrite is safe.}

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

F029 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 230: m_eb_scal_solve_op->setEBDirichlet(lev, phi, *eta[lev]); where the op is built at lines 30-33 with default a_ncomp=1 (getNComp()==1) and eta=tra_eta has m_ntrac comps (incflo_apply_corrector.cpp:110). AMReX_MLEBABecLap.cpp:235: AMREX_ALWAYS_ASSERT(beta_ncomp == 1 || beta_ncomp == ncomp) aborts for ntrac>=2. Path gated by tracer_eb non-empty (incflo_arrays.cpp:43) and CN/Implicit diff (incflo_update_tracer.cpp:54); eb_flow.tracer is mandatory with eb_flow+tracers (init.cpp:228).

F031 — confirmed (one verifier lens)

Lens 1 (refutation attempt): compute_laps EB branch (lines 499-571) never calls setEBDirichlet/setEBHomogDirichlet on m_eb_scal_apply_op before MLMG mlmg(*m_eb_scal_apply_op); mlmg.apply(...) (554-555); AMReX defaults to homogeneous Neumann when m_eb_phi is null (MLEBABecLap.H:325 isEBDirichlet(){return m_eb_phi[0]!=nullptr;}). diffuse_scalar sets Dirichlet on the separate solve op (line 230, added by PR #148 to honor the EB Dirichlet BC). Explicit diff uses laps only (incflo_update_tracer.cpp:28-29); CN predictor laps_o (predictor:158-159) vs half-implicit solve mismatch.

F045 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Same defect as F031: no setEBDirichlet/setEBHomogDirichlet anywhere on m_eb_scal_apply_op (repo-wide grep), so mlmg.apply at line 555 sees homogeneous-Neumann EB, while diffuse_scalar lines 228-230 apply tracer_eb/temperature_eb Dirichlet in the implicit solve. incflo_arrays.cpp:42-48 defines tracer_eb/temperature_eb exactly to 'Allow for Dirichlet BC on EB', confirming intent; with diffusion_type=0 the value never enters diffusion.

F030 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 639: m_eb_vel_apply_op->setEBHomogDirichlet(lev, *a_eta[lev]); unconditional, vs diffuse_vel_components lines 382-387 if (m_incflo->hasEBFlow()) { ... setEBDirichlet(lev, phi, *eta[lev]); } and DiffusionTensorOp.cpp:262-263 setEBShearViscosityWithInflow(lev, *a_eta[lev], *(m_incflo->get_velocity_eb()[lev])). Reachable via incflo_diffusion.cpp:52 (use_tensor_solve=false) and :24 (use_tensor_correction); called whenever need_divtau() (incflo.H:817-818, predictor:146-150).

F044 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Duplicate of F030, same line 639 setEBHomogDirichlet ignoring hasEBFlow(). velocity_eb exists when eb_flow enabled (incflo_arrays.cpp:38-41) and the implicit component solve of the same step uses it (lines 382-387), so explicit divtau treats the inflowing EB as a stationary no-slip wall. need_divtau() is true even for Implicit when m_godunov_include_diff_in_forcing (default true), so the inconsistent explicit term feeds the forcing. History: hasEBFlow branch added to diffuse_vel_components (commit 8a92d3d) but compute_divtau never updated.

---

_{Based on commit 7307d872, which is also the tree the audit verified against. From an automated audit of src/. Audit finding ids: F029, F031, F045, F030, F044. Reviewer unit(s): Diffusion+Rheology, theme:eb-smallcell. 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.}