Severity: high/medium · Category: correctness · Fix order: 5 of 21 — fix this 5th.
Filenames are numbered in reverse fix order: 001 = fix last, 021 = fix first. This file is 017.
Locations: Source/NavierStokesBase.cpp:5047, Source/NavierStokesBase.cpp:4868, Source/NavierStokesBase.cpp:4993
Based on commit 9bf664bf (line numbers refer to that tree).
All three defects sit in the EB reflux/redistribution tail of NavierStokesBase::ComputeAofs (roughly lines 4850-5075), and each independently breaks conservation or grid-independence of the reflux/mac_sync at cut cells. They date from the EB redistribution rewrite (bf61d61, "Update IAMR to use new redist functionality") and its follow-ups, and they edit adjacent lines, so they are best fixed in one pass.
The defect
Source/NavierStokesBase.cpp:5047 — The cut-cell overloads of EBFluxRegister::CrseAdd/FineAdd are handed fluxes that AMReX-Hydro already multiplied by the EB area fraction, and the register multiplies by the area fraction again, so cut faces on a coarse/fine interface are refluxed with ap^2 instead of ap.
Source/NavierStokesBase.cpp:4868 — In the mac_sync path with EB StateRedist, rstate_tmp is copied from Vsync/Ssync whose ghost cells are still the zeros set in advance_setup, and no FillBoundary is done, so state redistribution reads bogus zero "state" in every grid's ghost region.
Source/NavierStokesBase.cpp:4993 — ApplyMLRedistribution is called with use_wts_in_divnc sitting in the fac_for_deltaR (Real) argument slot, so fac_for_deltaR is always 1.0 and dm_as_fine is scaled by +dt even during the mac_sync, where the accompanying fluxes are handed to FineAdd with sync_factor*dt = -dt.
Why it matters
F006: AMREX_USE_EB build, >=2 AMR levels with do_reflux=1 and an embedded boundary crossing a coarse/fine interface. HydroUtils::EB_ComputeFluxes stores fx = xedumacapxdydz; eb_flux_reg_crseadd_va then forms dt/dVapx/vfracfx (and eb_flux_reg_fineadd_va_* forms sum f*a), so the reflux correction on a cut c/f face with apx=0.5 is half its correct magnitude, breaking conservation across the interface.
F032: AMREX_USE_EB, redistribution_type=StateRedist (the default), >=2 AMR levels, more than one grid per level. mac_sync -> ComputeAofs(is_sync=true) sets rstate=&rstate_tmp; StateRedistribute builds Qhat = U_in + dt*dUdt_in out to grow(bx,3), reading zeros instead of the neighbouring grid's Vsync/Ssync, so the sync correction near every box boundary is grid-decomposition dependent.
F033: AMREX_USE_EB with ns.redistribution_type=FluxRedist, >=3 AMR levels, do_reflux=1. In mac_sync at level 1, advflux_reg->FineAdd (line 5069) adds the sync fluxes with -dt but adds dm_as_fine unscaled (+dt*drho, see EBFluxRegister::FineAdd), so the re-redistributed mass entering the level-0/1 register has the wrong sign and level 0's reflux is corrupted at cut cells near the interface.
How to reach it
- Exec/eb_run2d/regtest.2d.shock_past_cylinder: EB build, amr.max_level=1, ns.refine_cutcells=0, default ns.do_reflux=1; box tagging leaves the cylinder EB crossing the level-0/1 boundary, so fractional-ap c/f faces reflux every step.
Suggested fix
Hand the cut-cell CrseAdd/FineAdd fluxes that are not ap-weighted: AMReX's EBFluxRegister header states "the flux is not scaled", and eb_flux_reg_crseadd_va/fineadd_va apply ax themselves (Tests/EB_CNS gives one unscaled fab to both overloads). Divide by ap into per-direction temporaries (zero where ap==0) rather than rescaling cfluxes in place, since EB_ComputeDivergence needs the area-weighted values; keep dxDp=dV.
rstate_tmp needs the FillBoundary update_MF already gets (4853) — InitialRedistribution says internal ghosts must be filled before redistribution; physical-boundary ghosts staying zero is right for a correction.
Hoist sync_factor (5060) above the redistribution call and pass it as fac_for_deltaR so the bool reaches use_wts_in_divnc; those slots were swapped in d270067 (#139) to fix the build after amrex inserted fac_for_deltaR. Decide: use_wts_in_divnc=true now takes effect for FluxRedist (matches 5000), a result change to regression-test; and whether as_crse/as_fine should be gated on do_crse_add/do_fine_add.
For Source/NavierStokesBase.cpp:5047 (F006):
--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -5034,6 +5034,45 @@
AMREX_D_TERM(FArrayBox fx_fr_fab(fx_fab,amrex::make_alias,flux_comp,ncomp);,
FArrayBox fy_fr_fab(fy_fab,amrex::make_alias,flux_comp,ncomp);,
FArrayBox fz_fr_fab(fz_fab,amrex::make_alias,flux_comp,ncomp););
+
+ //
+ // The cut-cell CrseAdd/FineAdd overloads below multiply the
+ // fluxes by the area fraction, but the fluxes from AMReX-Hydro
+ // are already area-fraction weighted. Un-weight them here so
+ // cut faces are not refluxed with ap^2.
+ //
+ AMREX_D_TERM(FArrayBox fx_unwtd_fab;,
+ FArrayBox fy_unwtd_fab;,
+ FArrayBox fz_unwtd_fab;);
+ if ( do_reflux &&
+ ( (do_crse_add && level < parent->finestLevel()) ||
+ (do_fine_add && level > 0) ) &&
+ flagfab.getType(amrex::grow(bx,1)) != FabType::regular )
+ {
+ AMREX_D_TERM(auto const& apx = areafrac[0]->const_array(mfi);,
+ auto const& apy = areafrac[1]->const_array(mfi);,
+ auto const& apz = areafrac[2]->const_array(mfi););
+ AMREX_D_TERM(fx_unwtd_fab.resize(fx_fr_fab.box(),ncomp,The_Async_Arena());,
+ fy_unwtd_fab.resize(fy_fr_fab.box(),ncomp,The_Async_Arena());,
+ fz_unwtd_fab.resize(fz_fr_fab.box(),ncomp,The_Async_Arena()););
+ AMREX_D_TERM(auto const& fxu = fx_unwtd_fab.array();,
+ auto const& fyu = fy_unwtd_fab.array();,
+ auto const& fzu = fz_unwtd_fab.array(););
+ AMREX_D_TERM(auto const& fxw = fx_fr_fab.const_array();,
+ auto const& fyw = fy_fr_fab.const_array();,
+ auto const& fzw = fz_fr_fab.const_array(););
+ amrex::ParallelFor(fx_unwtd_fab.box(), ncomp, [fxu,fxw,apx]
+ AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+ { fxu(i,j,k,n) = (apx(i,j,k) > 0.) ? fxw(i,j,k,n)/apx(i,j,k) : 0.; });
+ amrex::ParallelFor(fy_unwtd_fab.box(), ncomp, [fyu,fyw,apy]
+ AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+ { fyu(i,j,k,n) = (apy(i,j,k) > 0.) ? fyw(i,j,k,n)/apy(i,j,k) : 0.; });
+#if (AMREX_SPACEDIM == 3)
+ amrex::ParallelFor(fz_unwtd_fab.box(), ncomp, [fzu,fzw,apz]
+ AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+ { fzu(i,j,k,n) = (apz(i,j,k) > 0.) ? fzw(i,j,k,n)/apz(i,j,k) : 0.; });
+#endif
+ }
// Now update the flux registers (inside test on AMREX_USE_EB)
if ( do_reflux && do_crse_add && (level < parent->finestLevel()) ) {
@@ -5045,7 +5084,7 @@
} else if (flagfab.getType(bx) != FabType::covered ) {
getAdvFluxReg(level + 1).CrseAdd(mfi,
- {AMREX_D_DECL(&fx_fr_fab,&fy_fr_fab,&fz_fr_fab)},
+ {AMREX_D_DECL(&fx_unwtd_fab,&fy_unwtd_fab,&fz_unwtd_fab)},
dxDp, dt, (*volfrac)[mfi],
{AMREX_D_DECL(&(*areafrac[0])[mfi], &(*areafrac[1])[mfi], &(*areafrac[2])[mfi])},
0, state_indx, ncomp, amrex::RunOn::Device);
@@ -5067,7 +5106,7 @@
dxDp, sync_factor*dt, 0, state_indx, ncomp, amrex::RunOn::Device);
} else if (flagfab.getType(bx) != FabType::covered ) {
advflux_reg->FineAdd(mfi,
- {AMREX_D_DECL(&fx_fr_fab,&fy_fr_fab,&fz_fr_fab)},
+ {AMREX_D_DECL(&fx_unwtd_fab,&fy_unwtd_fab,&fz_unwtd_fab)},
dxDp, sync_factor*dt, (*volfrac)[mfi],
{AMREX_D_DECL(&(*areafrac[0])[mfi], &(*areafrac[1])[mfi], &(*areafrac[2])[mfi])},
dm_as_fine, 0, state_indx, ncomp, amrex::RunOn::Device);amrex eb_flux_reg_crseadd_va/fineadd_va kernels (AMReX_EBFluxRegister_*_C.H) multiply flux by area fraction, but HydroUtils::EB_ComputeFluxes already ap-weighted them. Builds Async-arena copies divided by ap (0 where ap==0) and passes those to the cut-cell CrseAdd/FineAdd only; regular branches and the dxDp=dV hack unchanged. Maintainer may prefer un-weighting upstream instead.
For Source/NavierStokesBase.cpp:4868 (F032):
--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -4866,6 +4866,11 @@
rstate_tmp.define(S.boxArray(),S.DistributionMap(),ncomp,S.nGrow(),
MFInfo(),ebfact);
MultiFab::Copy(rstate_tmp,advc,a_comp,0,ncomp,S.nGrow());
+ //
+ // Vsync/Ssync only hold valid-region data; fill ghost cells from
+ // neighboring grids before using this as the redistribution "state".
+ //
+ rstate_tmp.FillBoundary(geom.periodicity());
}
MultiFab const* rstate = (is_sync && redistribution_type == "StateRedist")
? &rstate_tmp : &S;Vsync/Ssync ghosts are the zeros set in advance_setup (reflux writes valid cells only), and amrex ApplyRedistribution builds Qhat = U_in + dt*dUdt_in out to grow(bx,3); FillBoundary(geom.periodicity()) fills grid-overlap ghosts so the SRD state is decomposition-independent. Physical/coarse-fine ghosts remain zero (best available).
For Source/NavierStokesBase.cpp:4993 (F033):
--- a/Source/NavierStokesBase.cpp
+++ b/Source/NavierStokesBase.cpp
@@ -4990,7 +4990,9 @@
geom, dt, redistribution_type,
as_crse, p_drho_as_crse->array(), p_rrflag_as_crse->array(),
as_fine, dm_as_fine.array(), coarse_fine_mask->const_array(mfi),
- level_mask_notcovered, use_wts_in_divnc);
+ level_mask_notcovered,
+ /*fac_for_deltaR*/ do_crse_add ? Real(1.0) : Real(-1.0),
+ use_wts_in_divnc);
} else {
bool use_wts_in_divnc = true;
ApplyRedistribution( bx, ncomp, redist_arr, update_arr,Per amrex AMReX_EB_Redistribution.H, ApplyMLRedistribution takes Real fac_for_deltaR before bool use_wts_in_divnc; passing do_crse_add ? 1.0 : -1.0 makes dm_as_fine/drho_as_crse carry the same sign as the sync_factor*dt used in FineAdd (mac_sync calls always have do_crse_add=false), and use_wts_in_divnc now actually binds true.
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
F006 — confirmed (two independent verifier lenses)
Lens 1 (refutation attempt): IAMR passes cfluxes from HydroUtils with fluxes_are_area_weighted=true; EB_ComputeFluxes (AMReX-Hydro hydro_utils.cpp:474) sets fx = xedumacapxarea. amrex eb_flux_reg_crseadd_va then forms tmp = dtdxax(i,j,k)volinv * fx (already apx-weighted) and eb_flux_reg_fineadd_va_ sums f*a — area fraction applied twice, so cut c/f faces reflux with ap^2. No compensating rescale of fx_fr_fab (5034-5051, 5069-5073).
Lens 2 (reachability/intent): Hydro's EB_ComputeFluxes (hydro_utils.cpp:474) stores fx=xedumacapxarea; IAMR passes these to cut-cell CrseAdd/FineAdd (NavierStokesBase.cpp:5047,5069) whose kernels multiply by areafrac again (AMReX_EBFluxRegister_2D_C.H:20 dtdxaxvolinv; fineadd fa), though AMReX_EBFluxRegister.H requires unscaled centroid flux. The dxD hack (line 4879) compensates only dy*dz, not ap. Default ns.refine_cutcells=1 keeps EB off c/f faces (ap=1, no-op); refine_cutcells=0 only warns (line 1333), never aborts.
F032 — confirmed (one verifier lens)
Lens 1 (refutation attempt): 4866-4868: rstate_tmp copied from advc (Vsync/Ssync) with S.nGrow() ghosts; only update_MF gets FillBoundary (4853), never advc. Vsync/Ssync ghosts are the zeros from advance_setup:649-650 (reflux() writes valid cells only), and amrex ApplyRedistribution builds scratch = U_in + dt*dUdt_in over Box(scratch)=grow(bx,3) (EB_RedistributionApply.cpp:172-179), reading zeros where the neighbor grid holds data -> decomposition-dependent sync near box boundaries.
F033 — confirmed (one verifier lens)
Lens 1 (refutation attempt): Confirmed against current amrex AMReX_EB_Redistribution.H:135-166: after levmsk come 'int level_mask_not_covered, Real fac_for_deltaR=1.0, bool use_wts_in_divnc=false'. IAMR (4986-4993) passes '..., level_mask_notcovered, use_wts_in_divnc)' so the bool binds to fac_for_deltaR (=1.0) and use_wts_in_divnc silently defaults false despite 'bool use_wts_in_divnc = true' at 4985. dm/drho scale by fac_for_deltaRdt=+dt (EB_RedistributionApply.cpp:106) while sync FineAdd fluxes get sync_factordt=-dt (5060,5071) — sign mismatch as alleged.
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: F006, F032, F033. Reviewer unit(s): NSB-4 predict/ComputeAofs/redistribution. 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.
Severity: high/medium · Category: correctness · Fix order: 5 of 21 — fix this 5th.
Filenames are numbered in reverse fix order:
001= fix last,021= fix first. This file is017.Locations:
Source/NavierStokesBase.cpp:5047,Source/NavierStokesBase.cpp:4868,Source/NavierStokesBase.cpp:4993Based on commit
9bf664bf(line numbers refer to that tree).All three defects sit in the EB reflux/redistribution tail of
NavierStokesBase::ComputeAofs(roughly lines 4850-5075), and each independently breaks conservation or grid-independence of the reflux/mac_sync at cut cells. They date from the EB redistribution rewrite (bf61d61, "Update IAMR to use new redist functionality") and its follow-ups, and they edit adjacent lines, so they are best fixed in one pass.The defect
Source/NavierStokesBase.cpp:5047— The cut-cell overloads of EBFluxRegister::CrseAdd/FineAdd are handed fluxes that AMReX-Hydro already multiplied by the EB area fraction, and the register multiplies by the area fraction again, so cut faces on a coarse/fine interface are refluxed with ap^2 instead of ap.Source/NavierStokesBase.cpp:4868— In the mac_sync path with EB StateRedist, rstate_tmp is copied from Vsync/Ssync whose ghost cells are still the zeros set in advance_setup, and no FillBoundary is done, so state redistribution reads bogus zero "state" in every grid's ghost region.Source/NavierStokesBase.cpp:4993— ApplyMLRedistribution is called with use_wts_in_divnc sitting in the fac_for_deltaR (Real) argument slot, so fac_for_deltaR is always 1.0 and dm_as_fine is scaled by +dt even during the mac_sync, where the accompanying fluxes are handed to FineAdd with sync_factor*dt = -dt.Why it matters
F006: AMREX_USE_EB build, >=2 AMR levels with do_reflux=1 and an embedded boundary crossing a coarse/fine interface. HydroUtils::EB_ComputeFluxes stores fx = xedumacapxdydz; eb_flux_reg_crseadd_va then forms dt/dVapx/vfracfx (and eb_flux_reg_fineadd_va_* forms sum f*a), so the reflux correction on a cut c/f face with apx=0.5 is half its correct magnitude, breaking conservation across the interface.
F032: AMREX_USE_EB, redistribution_type=StateRedist (the default), >=2 AMR levels, more than one grid per level. mac_sync -> ComputeAofs(is_sync=true) sets rstate=&rstate_tmp; StateRedistribute builds Qhat = U_in + dt*dUdt_in out to grow(bx,3), reading zeros instead of the neighbouring grid's Vsync/Ssync, so the sync correction near every box boundary is grid-decomposition dependent.
F033: AMREX_USE_EB with ns.redistribution_type=FluxRedist, >=3 AMR levels, do_reflux=1. In mac_sync at level 1, advflux_reg->FineAdd (line 5069) adds the sync fluxes with -dt but adds dm_as_fine unscaled (+dt*drho, see EBFluxRegister::FineAdd), so the re-redistributed mass entering the level-0/1 register has the wrong sign and level 0's reflux is corrupted at cut cells near the interface.
How to reach it
Suggested fix
Hand the cut-cell CrseAdd/FineAdd fluxes that are not ap-weighted: AMReX's EBFluxRegister header states "the flux is not scaled", and eb_flux_reg_crseadd_va/fineadd_va apply ax themselves (Tests/EB_CNS gives one unscaled fab to both overloads). Divide by ap into per-direction temporaries (zero where ap==0) rather than rescaling cfluxes in place, since EB_ComputeDivergence needs the area-weighted values; keep dxDp=dV.
rstate_tmp needs the FillBoundary update_MF already gets (4853) — InitialRedistribution says internal ghosts must be filled before redistribution; physical-boundary ghosts staying zero is right for a correction.
Hoist sync_factor (5060) above the redistribution call and pass it as fac_for_deltaR so the bool reaches use_wts_in_divnc; those slots were swapped in d270067 (#139) to fix the build after amrex inserted fac_for_deltaR. Decide: use_wts_in_divnc=true now takes effect for FluxRedist (matches 5000), a result change to regression-test; and whether as_crse/as_fine should be gated on do_crse_add/do_fine_add.
For
Source/NavierStokesBase.cpp:5047(F006):amrex eb_flux_reg_crseadd_va/fineadd_va kernels (AMReX_EBFluxRegister_*_C.H) multiply flux by area fraction, but HydroUtils::EB_ComputeFluxes already ap-weighted them. Builds Async-arena copies divided by ap (0 where ap==0) and passes those to the cut-cell CrseAdd/FineAdd only; regular branches and the dxDp=dV hack unchanged. Maintainer may prefer un-weighting upstream instead.
For
Source/NavierStokesBase.cpp:4868(F032):Vsync/Ssync ghosts are the zeros set in advance_setup (reflux writes valid cells only), and amrex ApplyRedistribution builds Qhat = U_in + dt*dUdt_in out to grow(bx,3); FillBoundary(geom.periodicity()) fills grid-overlap ghosts so the SRD state is decomposition-independent. Physical/coarse-fine ghosts remain zero (best available).
For
Source/NavierStokesBase.cpp:4993(F033):Per amrex AMReX_EB_Redistribution.H, ApplyMLRedistribution takes Real fac_for_deltaR before bool use_wts_in_divnc; passing do_crse_add ? 1.0 : -1.0 makes dm_as_fine/drho_as_crse carry the same sign as the sync_factor*dt used in FineAdd (mac_sync calls always have do_crse_add=false), and use_wts_in_divnc now actually binds true.
Diff(s) are against
9bf664bf, written from the current source and verified only withgit apply --check— never compiled, never run, never applied to the tree. Treat them as precise intent, not tested patches.Verification evidence
F006— confirmed (two independent verifier lenses)Lens 1 (refutation attempt): IAMR passes cfluxes from HydroUtils with fluxes_are_area_weighted=true; EB_ComputeFluxes (AMReX-Hydro hydro_utils.cpp:474) sets fx = xedumacapxarea. amrex eb_flux_reg_crseadd_va then forms tmp = dtdxax(i,j,k)volinv * fx (already apx-weighted) and eb_flux_reg_fineadd_va_ sums f*a — area fraction applied twice, so cut c/f faces reflux with ap^2. No compensating rescale of fx_fr_fab (5034-5051, 5069-5073).
Lens 2 (reachability/intent): Hydro's EB_ComputeFluxes (hydro_utils.cpp:474) stores fx=xedumacapxarea; IAMR passes these to cut-cell CrseAdd/FineAdd (NavierStokesBase.cpp:5047,5069) whose kernels multiply by areafrac again (AMReX_EBFluxRegister_2D_C.H:20 dtdxaxvolinv; fineadd fa), though AMReX_EBFluxRegister.H requires unscaled centroid flux. The dxD hack (line 4879) compensates only dy*dz, not ap. Default ns.refine_cutcells=1 keeps EB off c/f faces (ap=1, no-op); refine_cutcells=0 only warns (line 1333), never aborts.
F032— confirmed (one verifier lens)Lens 1 (refutation attempt): 4866-4868: rstate_tmp copied from advc (Vsync/Ssync) with S.nGrow() ghosts; only update_MF gets FillBoundary (4853), never advc. Vsync/Ssync ghosts are the zeros from advance_setup:649-650 (reflux() writes valid cells only), and amrex ApplyRedistribution builds scratch = U_in + dt*dUdt_in over Box(scratch)=grow(bx,3) (EB_RedistributionApply.cpp:172-179), reading zeros where the neighbor grid holds data -> decomposition-dependent sync near box boundaries.
F033— confirmed (one verifier lens)Lens 1 (refutation attempt): Confirmed against current amrex AMReX_EB_Redistribution.H:135-166: after levmsk come 'int level_mask_not_covered, Real fac_for_deltaR=1.0, bool use_wts_in_divnc=false'. IAMR (4986-4993) passes '..., level_mask_notcovered, use_wts_in_divnc)' so the bool binds to fac_for_deltaR (=1.0) and use_wts_in_divnc silently defaults false despite 'bool use_wts_in_divnc = true' at 4985. dm/drho scale by fac_for_deltaRdt=+dt (EB_RedistributionApply.cpp:106) while sync FineAdd fluxes get sync_factordt=-dt (5060,5071) — sign mismatch as alleged.
Based on commit
9bf664bf, which is also the tree the audit verified against. From an automated audit ofSource/,Tutorials/andUtil/. Audit finding ids: F006, F032, F033. Reviewer unit(s): NSB-4 predict/ComputeAofs/redistribution. 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.