diffuse_tensor_Vsync solves with viscosity 1.0, plus RHS-scaling and C-F BC errors

[WeiqunZhang]

Severity: high/medium/low · Category: correctness, physics-numerics · Fix order: 1 of 21 — fix this first.

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

Locations: Source/Diffusion.cpp:1121, Source/Diffusion.cpp:926, Source/Diffusion.cpp:1048, Source/Diffusion.cpp:1147, Source/Diffusion.cpp:518, Source/Diffusion.cpp:898

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

Nine reports, five distinct defects: F015 and F052 duplicate F002, F056 duplicates F021, and F061 duplicates F022 (same defect anchored at line 898 rather than 926). The remaining two are F020 (line 1048) and F062 (line 518). Three of the five — F002, F020 and F021 — live in diffuse_tensor_Vsync and touch the same solve, so they are best fixed in one pass.

The defect

Source/Diffusion.cpp:1121 — diffuse_tensor_Vsync ignores its beta (viscosity) argument entirely and performs the implicit viscous sync solve with shear viscosity hard-coded to 1.0 on all faces (and 1.0 for EB), even though the caller passes the real viscosity and diffuse_Vsync asserts on it. _{Reported independently at this same line by F015, F052; each reviewer's own wording and evidence is under Verification evidence below.}

Source/Diffusion.cpp:926 — diffuse_tensor_velocity computes rhsscale via computeAlpha (line 898) and applies it to the operator scalars (setScalars at line 901), but never multiplies Rhs by rhsscale before mlmg.solve, so the solved velocity is wrong by a factor 1/rhsscale. _{Reported independently at this same line by F061; each reviewer's own wording and evidence is under Verification evidence below.}

Source/Diffusion.cpp:1048 — In diffuse_tensor_Vsync with rho_flag==3 (do_mom_diff=1) the RHS is multiplied by old-time density while the operator acoef uses new-time density (line 1107) and the caller normalized Vsync by new-time density, so the momentum sync is rescaled by rho^n/rho^{n+1}.

Source/Diffusion.cpp:1147 — diffuse_tensor_Vsync scales only component 0 of the AMREX_SPACEDIM-component RHS by rhsscale ('Rhs.mult(rhsscale,0,1)'), while the operator scalars are scaled by rhsscale for all components, so the y (and z) Vsync solutions are off by a factor 1/rhsscale. _{Reported independently at this same line by F056; each reviewer's own wording and evidence is under Verification evidence below.}

Source/Diffusion.cpp:518 — diffuse_scalar never calls setCoarseFineBC when has_coarse_data is false but the solve is on level>0 (the mac_sync scalar path passes nlev=1), so MLMG falls back to its default assumed coarse/fine ratio of 2 for the homogeneous Dirichlet C-F boundary.

Why it matters

F002: Any multilevel run (max_level>=1, do_reflux on) with viscous velocity (ns.vel_visc_coef>0): mac_sync calls diffuse_Vsync(...,loc_viscn,...); the (rho - thetadtdiv tau) solve for Vsync uses mu=1 instead of mu, so the coarse-level velocity sync correction and the viscflux_reg FineAdd fluxes (line 1177) are wrong for any mu != 1.

F022: Any run with implicit velocity diffusion and diffuse.scale_abec=1 and non-unit density: rhsscale=1/max(rhoalpha)!=1, the solve is (rhsscale(rho - b div beta grad))U = Rhs, so U_new = max(rho)*U_correct every time step -- blatantly wrong momentum solution. diffuse_scalar (line 563) and diffuse_Ssync (line 1301) both scale their RHS; this near-twin routine omits it.

F020: do_mom_diff=1, variable density, multilevel with do_reflux: NavierStokes::mac_sync divides Vsync by S_new Density (NavierStokes.cpp:1543), then diffuse_tensor_Vsync multiplies by get_old_data Density; in the inviscid limit the solve returns (rho^n/rho^{n+1})*Vsync instead of Vsync — an O(delta rho) error in the coarse velocity sync.

F021: Multi-level run (Vsync only exists with more than one AMR level) with diffuse.scale_abec=1 and non-unit density: computeAlpha returns rhsscale=1/max(rho)!=1; the x-velocity sync is consistent but the y/z velocity sync corrections come out multiplied by max(rho), corrupting the coarse velocity after mac_sync. Note the correct pattern in diffuse_scalar line 563: Rhs.mult(rhsscale,0,nComp).

F062: amr.max_level>=2 with amr.ref_ratio=4, diffusive scalars, do_reflux: mac_sync on level 1 (NavierStokes.cpp:1626) solves the Ssync diffusion with the C-F Dirichlet value interpolated as if the coarse cell were at ratio-2 distance (1.5 fine cells) instead of ratio-4 (2.5), giving a wrong stencil at every coarse-fine boundary; diffuse_Ssync and diffuse_tensor_Vsync handle this case explicitly.

How to reach it

• Tutorials/Bubble/inputs.2d.bubble: amr.max_level=2, ns.vel_visc_coef=0.001, defaults do_mac_proj=1/do_reflux=1; post_timestep→mac_sync→diffuse_Vsync every coarse step. Also Exec/eb_run2d/regtest.2d.flow_past_cylinder-x (visc=0.1, max_level=1).
• Exec/run2d/regtest.2d.hotspot: amr.max_level=2, ns.vel_visc_coef=0.001, defaults do_reflux=1, do_mac_proj=1. post_timestep -> mac_sync -> diffuse_Vsync -> diffuse_tensor_Vsync, unconditionally. RZ+viscosity aborts, Cartesian unaffected.

Suggested fix

F002/F015/F052 — restore the beta/betaComp parameters and call setViscosity(tensorop, beta, betaComp), as diffuse_tensor_velocity:909 and the pre-MLMG getTensorOp(a,b,rho,beta,betaComp) (52d2e6d, dropped in 1f37175) did. Judgement call: MLEBTensorOp also needs a cell-centered eta for setEBShearViscosity, so either thread viscn_cc through diffuse_Vsync the way diffuse_velocity threads betanCC (mac_sync uses prevTime viscosity), or average faces to centers — the former changes the public signature.

F020 — use get_new_data Density for rho_flag==3, matching acoef (1107) and mac_sync's division by S_new density; the same correction PR #160 (17cb508) made in diffuse_tensor_velocity.

F021/F056, F022/F061 — scale all AMREX_SPACEDIM components, and add the missing Rhs.mult in diffuse_tensor_velocity (hoisting rhsscale out of its block). Flux registers are safe either way: MLCellABecLap::getFluxes divides the b scalar back out.

F062 — mirror diffuse_Ssync:1261 with else if (level>0) setCoarseFineBC(nullptr,cratio[0]) on opn and opnp1; otherwise needsCoarseDataForBC() is false, br_ref_ratio becomes 1, and MLMGBndry puts the homogeneous C-F Dirichlet at 0.5*dx instead of 0.5*ratio*dx.

For Source/Diffusion.cpp:1121 (F002):

--- a/Source/Diffusion.cpp
+++ b/Source/Diffusion.cpp
@@ -1017,8 +1017,8 @@
                                  Real                   be_cn_theta,
                                  const MultiFab&        rho_half,
                                  int                    rho_flag,
-                                 const MultiFab* const* /*beta*/,
-                                 int                    /*betaComp*/,
+                                 const MultiFab* const* beta,
+                                 int                    betaComp,
                                  bool                   update_fluxreg)
 {
     AMREX_ASSERT(rho_flag == 1 || rho_flag == 3);
@@ -1113,22 +1113,11 @@
       tensorop.setACoeffs(0, acoef);
     }
 
-    {
-      FluxBoxes  fb_bcoef;
-      MultiFab** face_bcoef = nullptr;
-      face_bcoef = fb_bcoef.define(navier_stokes);
-      for (int dir=0; dir<AMREX_SPACEDIM; dir++) {
-         face_bcoef[dir]->setVal(1.0);
-      }
-
 #ifdef AMREX_USE_EB
-      MultiFab bcoefCC(grids,dmap,1,0,MFInfo(),navier_stokes->Factory());
-      bcoefCC.setVal(1.0);
-      setViscosity(tensorop, face_bcoef, 0, bcoefCC);
+    setViscosity(tensorop, beta, betaComp, *navier_stokes->viscn_cc);
 #else
-      setViscosity(tensorop, face_bcoef, 0);
+    setViscosity(tensorop, beta, betaComp);
 #endif
-    }
 
     MLMG mlmg(tensorop);
     if (max_iter > 0) {

_{Restores use of the passed viscosity via the existing setViscosity helper. EB branch passes *navier_stokes->viscn_cc (old-time cc viscosity backing the caller's loc_viscn; Diffusion is a friend of NavierStokesBase) to MLEBTensorOp::setEBShearViscosity, mirroring diffuse_tensor_velocity's betanCC. Maintainer may instead prefer threading a betaCC parameter through diffuse_Vsync.}

For Source/Diffusion.cpp:926 (F022):

--- a/Source/Diffusion.cpp
+++ b/Source/Diffusion.cpp
@@ -889,18 +889,19 @@
          tensorop.setLevelBC(0, &Soln);
       }
 
+      Real rhsscale = 1.0;
       {
          MultiFab acoef;
          std::pair<Real,Real> scalars;
-         Real rhsscale = 1.0;
          const MultiFab& rho = (rho_flag == 1) ? rho_half : navier_stokes->get_new_data(State_Type);
          const int rho_comp = (rho_flag == 1) ? 0 : Density;
          computeAlpha(acoef, scalars, a, b,
                       &rhsscale, nullptr, 0,
                       rho_flag, &rho, rho_comp);
          tensorop.setScalars(scalars.first, scalars.second);
          tensorop.setACoeffs(0, acoef);
       }
+      Rhs.mult(rhsscale,0,AMREX_SPACEDIM);
 
 #ifdef AMREX_USE_EB
       setViscosity(tensorop, betanp1, betaComp, *betanp1CC);

_{Hoists rhsscale out of the acoef scope and multiplies Rhs by it before mlmg.solve, matching diffuse_scalar (line 563) and diffuse_Ssync (line 1301). Only active with diffuse.scale_abec=1. Minor residue: tol_abs at line 850 is still computed from the unscaled Rhs; move get_scaled_abs_tol after the scaling for exact parity if desired.}

For Source/Diffusion.cpp:1048 (F020):

--- a/Source/Diffusion.cpp
+++ b/Source/Diffusion.cpp
@@ -1045,7 +1045,7 @@
     {
        const Box& bx = mfi.tilebox();
        auto const& rhs = Rhs.array(mfi);
-       auto const& rho = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_old_data(State_Type).array(mfi,Density);
+       auto const& rho = (rho_flag == 1) ? rho_half.array(mfi) : navier_stokes->get_new_data(State_Type).array(mfi,Density);
 
        amrex::ParallelFor(bx, [rhs, rho]
        AMREX_GPU_DEVICE (int i, int j, int k) noexcept

_{Uses new-time density for the rho_flag==3 RHS multiply, matching the acoef built from get_new_data (line 1107), the caller's Vsync /= rho^{n+1} normalization (NavierStokes.cpp:1543), and the identical line in twin diffuse_tensor_velocity (line 819, fixed the same way by PR #160).}

For Source/Diffusion.cpp:1147 (F021):

--- a/Source/Diffusion.cpp
+++ b/Source/Diffusion.cpp
@@ -1144,7 +1144,7 @@
     mlmg.setMaxFmgIter(max_fmg_iter);
     mlmg.setVerbose(verbose);
 
-    Rhs.mult(rhsscale,0,1);
+    Rhs.mult(rhsscale,0,AMREX_SPACEDIM);
 
     mlmg.setFinalFillBC(true);
     mlmg.solve({&Soln}, {&Rhs}, tol_rel, tol_abs);

_{Scales all AMREX_SPACEDIM components of Rhs by rhsscale, matching setScalars which scales the operator for all components; follows the correct twin pattern Rhs.mult(rhsscale,0,nComp) in diffuse_scalar (line 563). Only active with diffuse.scale_abec=1.}

For Source/Diffusion.cpp:518 (F062):

--- a/Source/Diffusion.cpp
+++ b/Source/Diffusion.cpp
@@ -394,6 +394,9 @@
                 }
                 opn.setCoarseFineBC(Solnc.get(), cratio[0]);
             }
+            else if (level > 0) {
+                opn.setCoarseFineBC(nullptr, cratio[0]);
+            }
             MultiFab::Copy(Soln,*S_old[0],sigma,0,nComp,ng);
             if (rho_flag == 2) {
 #ifdef AMREX_USE_OMP
@@ -516,7 +519,10 @@
             });
         }
         opnp1.setCoarseFineBC(Solnc.get(), cratio[0]);
     }
+    else if (level > 0) {
+        opnp1.setCoarseFineBC(nullptr, cratio[0]);
+    }
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif

_{Mirrors diffuse_Ssync (line 1261) and diffuse_tensor_Vsync (line 1100): when no coarse data is supplied on a level>0 solve, register the homogeneous C-F Dirichlet BC with the true refinement ratio via setCoarseFineBC(nullptr, cratio[0]) (amrex MLLinOp::setCoarseFineBC accepts null coarse data). Applied to both opn and opnp1 for consistency; the failing mac_sync path exercises opnp1.}

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

F002 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): Params are '/beta/','/betaComp/' (1020-21); 'face_bcoef[dir]->setVal(1.0)' (1121) and 'bcoefCC.setVal(1.0)' feed setViscosity. Caller NavierStokes.cpp:1562 passes loc_viscn from getViscosity; diffuse_Vsync asserts beta[d]->min>=0. Pre-MLMG code used getTensorOp(a,b,rho,beta,betaComp) (52d2e6d); commit 1f37175 (2020 'WIP Diffusion clean up') introduced the 1.0 placeholder. viscflux_reg->FineAdd at 1177 uses these unit-mu fluxes. Not related to issue #162.

Lens 2 (reachability/intent): Diffusion.cpp:1020-21 params are /beta/,/betaComp/; :1121 face_bcoef setVal(1.0), :1126 bcoefCC=1.0 -> setViscosity. Caller NavierStokes.cpp:1562 passes real loc_viscn; wrapper asserts on beta it never uses (:975-980). Pre-MLMG code (f590444) passed beta to getTensorOp; WIP rewrite abc6488 dropped it, 3e6224a just silenced the unused warning — no comment documents intent. Result feeds Vsync (:1155), level_sync SyncInterp, and viscflux_reg->FineAdd (:1177); diffuse_Ssync uses its beta (:1280).

F015 — confirmed (two independent verifier lenses)

Reported as: diffuse_tensor_Vsync ignores the viscosity passed in (signature comments out beta/betaComp at line 1020) and solves the viscous velocity-sync system with unit shear viscosity (face_bcoef=1.0), unlike diffuse_tensor_velocity (line 909) which sets the real viscosity, and unlike the pre-MLMG code which used getTensorOp(a,b,rho,beta,betaComp).

Failure scenario: Any multilevel run with viscous velocity (ns.vel_visc_coef>0, do_reflux default): mac_sync's Vsync solve becomes (rho + thetadtdiv grad)dV instead of (rho + thetadtdiv mu grad)dV. For mu=1e-3, dt=1e-3, dx=1/64 the sync correction's fine-scale content is damped ~20x too much; the FineAdd'ed viscflux_reg fluxes (line 1177) are likewise computed with eta=1, corrupting the next-coarser viscous reflux.

Lens 1 (refutation attempt): Duplicate of F002. 'face_bcoef[dir]->setVal(1.0)' at 1121 vs diffuse_tensor_velocity's 'setViscosity(tensorop, betanp1, betaComp)' at 909. git log confirms pre-MLMG diffuse_tensor_Vsync built 'std::unique_ptr tensor_op(getTensorOp(a,b,rho,beta,betaComp))'; the MLMG rewrite (1f37175) dropped beta and it was never restored.

Lens 2 (reachability/intent): Diffusion.cpp:1020-1021 ignores beta; 1121 sets face_bcoef=1.0 fed to setViscosity, while diffuse_tensor_velocity uses real betanp1 (line 909) and diffuse_Ssync uses real beta. Pre-MLMG code (504dca8^ line 1512) used getTensorOp(a,b,rho,beta,betaComp); unit coefficient entered via computeBeta(face_bcoef,nullptr,0) in WIP commit 504dca8, no documenting comment. Wrapper still asserts beta[d]->min>=0. Result feeds level_sync projection (NavierStokesBase.cpp:1993-2001) and viscflux_reg FineAdd (1177).

F022 — confirmed (one verifier lens)

Lens 1 (refutation attempt): computeAlpha(acoef,scalars,a,b,&rhsscale,...) at 898-900 scales scalars by rhsscale=1/max(rho) when diffuse.scale_abec=1; rhsscale goes out of scope at 903 and Rhs is never multiplied before 'mlmg.solve({&Soln},{&Rhs},...)' at 926, unlike 'Rhs.mult(rhsscale,0,nComp)' at 563. Solve returns max(rho)*U. Real, but gated by non-default diffuse.scale_abec (default 0, line 87), hence low.

F020 — confirmed (one verifier lens)

Lens 1 (refutation attempt): Line 1048: rho_flag==3 Rhs multiply uses 'get_old_data(State_Type).array(mfi,Density)'; acoef at 1107 uses get_new_data; caller NavierStokes.cpp:1543 does 'vsync/=rho_c' with S_new density. Net RHS=(rho^n/rho^{n+1})*momentum sync. Twin diffuse_tensor_velocity (819) uses new density for rho_flag==3. Mismatch is pre-MLMG-era (rho_ptime vs rho_ctime in 52d2e6d) but uncompensated. Low: needs non-default ns.do_mom_diff=1; O(drho) error on a correction term.

F052 — confirmed (one verifier lens)

Reported as: diffuse_tensor_Vsync ignores its viscosity argument (parameter is literally '/beta/') and sets the tensor-solve shear viscosity to 1.0, so the implicit velocity-sync solve uses unit viscosity instead of mu; the pre-MLMG implementation passed beta to getTensorOp, and the sibling routines diffuse_tensor_velocity and diffuse_Ssync both use the real coefficients.

Failure scenario: Any multi-level run with do_reflux and viscous velocity (ns.vel_visc_coef>0, e.g. mu=1e-3): mac_sync's diffuse_Vsync solves (rho - thetadtdiv 1grad)dV = rhodV* instead of (rho - thetadtdiv mu*grad)dV, drastically over-damping the velocity sync correction at coarse-fine boundaries; with do_LES the spatially varying mu_t is likewise dropped.

Lens 1 (refutation attempt): Same defect as F002/F015: beta parameter commented out (1020), shear viscosity hard-coded 1.0 (1121, EB bcoefCC 1126). LES remark verified: NavierStokes::getViscosity adds calc_mut_LES output into loc_viscn, which mac_sync passes and diffuse_tensor_Vsync ignores. Revised high as it is the identical defect filed high in F002.

F021 — confirmed (one verifier lens)

Lens 1 (refutation attempt): 'Rhs.mult(rhsscale,0,1);' at 1147 scales one component of the AMREX_SPACEDIM-component Rhs while setScalars (1112) scales the operator for all components; compare 'Rhs.mult(rhsscale,0,nComp)' at 563. With scale_abec=1, y/z sync come out *max(rho). Real but gated by non-default diffuse.scale_abec=1, hence low.

F056 — confirmed (one verifier lens)

Reported as: diffuse_tensor_Vsync scales only component 0 of the AMREX_SPACEDIM-component Rhs by rhsscale (Rhs.mult(rhsscale,0,1)), while computeAlpha scaled the operator's a/b scalars by rhsscale for all components; the twin routine diffuse_scalar correctly uses Rhs.mult(rhsscale,0,nComp).

Failure scenario: Any multilevel viscous run with diffusion.scale_abec=1 (documented runtime option, echoed by echo_settings): rhsscale = 1/max(rho) != 1, so the solved Vsync y (and z) components come out too large by a factor max(rho) while x is correct, silently corrupting the velocity sync. The sibling diffuse_tensor_velocity (line 895) computes rhsscale but never applies it to Rhs at all, mis-scaling all velocity components each step.

Lens 1 (refutation attempt): Duplicate of F021: 'Rhs.mult(rhsscale,0,1)' at 1147 vs SPACEDIM-comp Rhs (1028); diffuse_tensor_velocity indeed never applies rhsscale to Rhs (895-926). Key name in claim is slightly wrong (it is diffuse.scale_abec, not diffusion.scale_abec; echoed by echo_settings line 178). Gated by non-default option: low.

F062 — confirmed (one verifier lens)

Lens 1 (refutation attempt): setCoarseFineBC only under has_coarse_data (395, 518); mac_sync (NavierStokes.cpp:1626) passes {} S_old and size-1 S_new at levels>0 when max_level>=2. Fallback detail differs from claim: amrex MLCellLinOp.H:864 resets br_ref_ratio to m_coarse_data_crse_ratio=-1, so MLMGBndry.H:201-202 places the homogeneous C-F Dirichlet at bcloc=m_coarse_bc_loc=0 (at the face), wrong for ANY ref_ratio, not an assumed ratio 2. diffuse_Ssync (1261) and diffuse_tensor_Vsync (1100) set it explicitly. Not issue #162.

F061 — confirmed (one verifier lens)

Reported as: diffuse_tensor_velocity computes rhsscale via computeAlpha and scales the operator scalars (a,b) by it, but never multiplies Rhs by rhsscale before mlmg.solve, unlike diffuse_scalar (line 563) and diffuse_Ssync (line 1301).

Failure scenario: With diffuse.scale_abec=1 and density whose max != 1 (any variable-density or rho!=1 run), the solve becomes rhsscale*(aalpha - bdiv beta grad)u = Rhs with unscaled Rhs, so the diffused velocity (and the reflux fluxes derived from it) is uniformly wrong by the factor max(rho).

Lens 1 (refutation attempt): Duplicate of F022: computeAlpha(...,&rhsscale,...) at 898, tensorop.setScalars(scalars.first,scalars.second) at 901, no Rhs.mult(rhsscale,...) anywhere before mlmg.solve at 926; diffuse_scalar (563) and diffuse_Ssync (1301) both scale. Only manifests with diffuse.scale_abec=1 (default 0): low as filed.

---

_{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: F002, F015, F022, F020, F052, F021, F056, F062, F061. Reviewer unit(s): Diffusion-1, Diffusion-2 + LES, theme:gpu-capture, theme:rz-metrics, theme:sync-reflux-subcycling. 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.}