Repository navigation
DiffusionScalarOp: the single EB scalar op is shared by the tracer and temperature passes, so an EB Dirichlet BC set by one pass (tracer_eb or temperature_eb) is silently applied to the other, with the other field's values and coefficient #241
Description
Severity: medium · Category: eb
Locations: src/diffusion/DiffusionScalarOp.cpp:30-34, src/diffusion/DiffusionScalarOp.cpp:49-54, src/diffusion/DiffusionScalarOp.cpp:228-234, src/diffusion/DiffusionScalarOp.cpp:541-552, src/diffusion/incflo_diffusion.cpp:72-73, src/diffusion/incflo_diffusion.cpp:82-83, src/diffusion/incflo_diffusion.cpp:92-94, src/diffusion/incflo_diffusion.cpp:103-106, src/setup/incflo_arrays.cpp:42-48
Based on commit 46de3367 (line numbers refer to that tree).
The defect
DiffusionScalarOp owns one MLEBABecLap solve op (m_eb_scal_solve_op, ctor lines 30-34) and one
apply op (m_eb_scal_apply_op, lines 49-54), and incflo::diffuse_scalar / diffuse_temperature
(incflo_diffusion.cpp:92-106) and compute_laps / compute_laps_T (72-83) drive both the tracer and the
temperature through the same object, passing get_tracer_eb() or get_temperature_eb() as the
eb_dirichlet argument.
Inside the op the EB BC is applied like this (solve op, lines 228-234; apply op has the identical pattern at
541-552):
if (!eb_dirichlet[lev]->empty()) {
MultiFab phi (*eb_dirichlet[lev], amrex::make_alias, comp, 1);
MultiFab beta(*eta[lev] , amrex::make_alias, comp, 1);
m_eb_scal_solve_op->setEBDirichlet(lev, phi, beta);
} // else use default homogeneous Neumann on EBThe else comment is only true for an op on which setEBDirichlet has never been called. In AMReX
(AMReX_MLEBABecLap.cpp:251-256, .H:315-325) setEBDirichlet allocates m_eb_phi[amrlev] and
m_eb_b_coeffs, and from then on isEBDirichlet() is true and Fapply/ebflux (line 1555) use the stored
m_eb_phi unconditionally; there is no public call that reverts the op to the homogeneous-Neumann EB
(the only "reset" is setEBHomogDirichlet, which is a different BC). The code even says so in the apply
op's comment (DiffusionScalarOp.cpp:543-547, added by PR #214) and then leaves it.
tracer_eb and temperature_eb are defined independently (incflo_arrays.cpp:42-48): tracer_eb whenever
advect_tracer and eb_flow.tracer is given, temperature_eb whenever use_temperature and
eb_flow.temperature is given; ReadParameters only forces both when eb_flow.enabled
(init.cpp:234-246) -- the "Dirichlet on the EB without flow through it" use case explicitly allowed by
incflo_arrays.cpp:42 does not. So with exactly one of the two lists given:
eb_flow.traceronly: the tracer pass callssetEBDirichlet(phi = tracer_eb, beta = mu_s). The
temperature pass that follows in the sameApplyPredictor(update_tracer, then update_temperature)
findstemperature_ebempty, sets nothing, and solves with the EB still Dirichlet: the EB value is the
tracer's EB value and the EB conductivity is the tracer diffusivitymu_s, notmu_T.eb_flow.temperatureonly: step 1 tracer pass is Neumann (correct once), the temperature pass sets
Dirichlet, and from step 2 on every tracer solve has Dirichlet EB with the temperature's EB value and
mu_Tas coefficient.
The same happens in compute_laps for the explicit/Crank-Nicolson term (apply op), and there the leak also
crosses the predictor's compute_laps (tracer) -> compute_laps_T (temperature) order within one step.
Why it matters
Silently wrong physics at the EB for whichever of tracer/temperature has no EB value of its own: a
homogeneous-Neumann (insulated) EB becomes a Dirichlet EB pinned to an unrelated field's value with an
unrelated coefficient. No warning, no abort. It affects both the implicit solve and the explicit term, so
switching diffusion_type does not avoid it.
How to reach it
USE_EB=TRUE, DIM 2 or 3, CPU or GPU. incflo.advect_tracer = 1, incflo.use_temperature = 1, an EB
geometry, eb_flow.tracer = <ntrac values> and no eb_flow.temperature (or the reverse), with
eb_flow.velocity/vel_mag absent so that ReadParameters does not demand both lists. Any
incflo.diffusion_type. E.g. test_3d/benchmark.channel_sphere with incflo.use_temperature = 1,
incflo.mu_T = 1e-3 and eb_flow.tracer = 1.0 added.
Suggested fix
Give the temperature its own pair of MLEBABecLap ops and select them with a flag passed by
incflo::diffuse_temperature / compute_laps_T; the tracer keeps the existing ops. This is a pure
plumbing change (the per-component logic is untouched). Diff checked with git apply --check
(audit/notes/T6-scratch/019.diff):
--- a/src/diffusion/DiffusionScalarOp.H
+++ b/src/diffusion/DiffusionScalarOp.H
@@ -21,7 +21,8 @@
amrex::Vector<amrex::MultiFab*> const& eb_dirichlet,
amrex::Vector<int> const& use_chi,
amrex::Vector<amrex::BCRec> bcrec,
- amrex::Real dt);
+ amrex::Real dt,
+ bool is_temperature = false);
void diffuse_vel_components (
amrex::Vector<amrex::MultiFab*> const& vel,
@@ -33,7 +34,8 @@
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);
+ amrex::Vector<amrex::BCRec> bcrec,
+ bool is_temperature = false);
void compute_divtau (amrex::Vector<amrex::MultiFab*> const& a_divtau,
amrex::Vector<amrex::MultiFab const*> const& a_vel,
@@ -49,6 +51,11 @@
#ifdef AMREX_USE_EB
std::unique_ptr<amrex::MLEBABecLap> m_eb_scal_solve_op;
std::unique_ptr<amrex::MLEBABecLap> m_eb_scal_apply_op;
+ // The temperature gets its own EB ops: MLEBABecLap cannot revert an EB
+ // Dirichlet BC to Neumann, so the tracer and temperature passes must not
+ // share an operator when only one of them has an EB Dirichlet value.
+ std::unique_ptr<amrex::MLEBABecLap> m_eb_tem_solve_op;
+ std::unique_ptr<amrex::MLEBABecLap> m_eb_tem_apply_op;
std::unique_ptr<amrex::MLEBABecLap> m_eb_vel_solve_op;
std::unique_ptr<amrex::MLEBABecLap> m_eb_vel_apply_op;
#endif
--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -33,6 +33,15 @@
info_solve, ebfact);
m_eb_scal_solve_op->setMaxOrder(m_mg_maxorder);
+ if (m_incflo->m_use_temperature)
+ {
+ m_eb_tem_solve_op = std::make_unique<MLEBABecLap>(m_incflo->Geom(0,finest_level),
+ m_incflo->boxArray(0,finest_level),
+ m_incflo->DistributionMap(0,finest_level),
+ info_solve, ebfact);
+ m_eb_tem_solve_op->setMaxOrder(m_mg_maxorder);
+ }
+
if (!m_incflo->useTensorSolve())
{
m_eb_vel_solve_op = std::make_unique<MLEBABecLap>(m_incflo->Geom(0,finest_level),
@@ -51,6 +60,15 @@
m_incflo->DistributionMap(0,finest_level),
info_apply, ebfact);
m_eb_scal_apply_op->setMaxOrder(m_mg_maxorder);
+
+ if (m_incflo->m_use_temperature)
+ {
+ m_eb_tem_apply_op = std::make_unique<MLEBABecLap>(m_incflo->Geom(0,finest_level),
+ m_incflo->boxArray(0,finest_level),
+ m_incflo->DistributionMap(0,finest_level),
+ info_apply, ebfact);
+ m_eb_tem_apply_op->setMaxOrder(m_mg_maxorder);
+ }
}
if ( (m_incflo->need_divtau() && !m_incflo->useTensorSolve()) ||
@@ -135,8 +153,10 @@
Vector<MultiFab*> const& eb_dirichlet,
amrex::Vector<int> const& use_chi,
amrex::Vector<amrex::BCRec> bcrec,
- Real dt)
+ Real dt,
+ bool is_temperature)
{
+ amrex::ignore_unused(is_temperature);
//
// Solves
// [alpha a - beta div ( b grad )] sca = RHS
@@ -168,6 +188,12 @@
const int finest_level = m_incflo->finestLevel();
+#ifdef AMREX_USE_EB
+ // Tracer and temperature must not share an EB op (see the header)
+ MLEBABecLap* eb_op = is_temperature ? m_eb_tem_solve_op.get()
+ : m_eb_scal_solve_op.get();
+#endif
+
Vector<MultiFab> rhs_c(finest_level+1);
// Note only conservative uses this rhs_c container
for (int lev = 0; lev <= finest_level; ++lev) {
@@ -175,14 +201,14 @@
}
#ifdef AMREX_USE_EB
- if (m_eb_scal_solve_op)
+ if (eb_op)
{
- m_eb_scal_solve_op->setScalars(1.0, dt);
+ eb_op->setScalars(1.0, dt);
for (int lev = 0; lev <= finest_level; ++lev) {
if ( use_chi[0] ) {
- m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]);
+ eb_op->setACoeffs(lev, *chi[lev]);
} else {
- m_eb_scal_solve_op->setACoeffs(lev, 1.0);
+ eb_op->setACoeffs(lev, 1.0);
}
}
}
@@ -202,16 +228,16 @@
for (int comp = 0; comp < a_scalar[0]->nComp(); ++comp)
{
#ifdef AMREX_USE_EB
- if (m_eb_scal_solve_op)
+ if (eb_op)
{
- m_eb_scal_solve_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low,
+ eb_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low,
bcrec[0].lo()),
m_incflo->get_diffuse_scalar_bc(Orientation::high,
bcrec[0].hi()));
if ( m_incflo->m_has_mixedBC && comp>0 ) {
// Must reset scalars (and Acoef, done below) to reuse solver with Robin BC
- m_eb_scal_solve_op->setScalars(1.0, dt);
+ eb_op->setScalars(1.0, dt);
}
for (int lev = 0; lev <= finest_level; ++lev) {
@@ -219,9 +245,9 @@
if (comp > 0 &&
( m_incflo->m_has_mixedBC || (use_chi[comp] != use_chi[comp-1]) ) ) {
if ( use_chi[comp] ) {
- m_eb_scal_solve_op->setACoeffs(lev, *chi[lev]);
+ eb_op->setACoeffs(lev, *chi[lev]);
} else {
- m_eb_scal_solve_op->setACoeffs(lev, 1.0);
+ eb_op->setACoeffs(lev, 1.0);
}
}
@@ -230,12 +256,12 @@
// single component of eta that we are solving for
MultiFab phi (*eb_dirichlet[lev], amrex::make_alias, comp, 1);
MultiFab beta(*eta[lev] , amrex::make_alias, comp, 1);
- m_eb_scal_solve_op->setEBDirichlet(lev, phi, beta);
+ eb_op->setEBDirichlet(lev, phi, beta);
} // else use default homogeneous Neumann on EB
Array<MultiFab,AMREX_SPACEDIM> b = m_incflo->average_scalar_eta_to_faces(lev, comp, *eta[lev]);
- m_eb_scal_solve_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
+ eb_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
}
}
else
@@ -286,19 +312,19 @@
}
#ifdef AMREX_USE_EB
- if (m_eb_scal_solve_op) {
+ if (eb_op) {
if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]);
- m_eb_scal_solve_op->setLevelBC(lev, &phi[lev],
+ eb_op->setLevelBC(lev, &phi[lev],
&robin[0], &robin[1], &robin[2]);
}
else {
- m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]);
+ eb_op->setLevelBC(lev, &phi[lev]);
}
// For when we use the stencil for centroid values
- // m_eb_scal_solve_op->setPhiOnCentroid();
+ // eb_op->setPhiOnCentroid();
} else
#endif
{
@@ -307,7 +333,7 @@
}
#ifdef AMREX_USE_EB
- MLMG mlmg(m_eb_scal_solve_op ? static_cast<MLLinOp&>(*m_eb_scal_solve_op) : static_cast<MLLinOp&>(*m_reg_scal_solve_op));
+ MLMG mlmg(eb_op ? static_cast<MLLinOp&>(*eb_op) : static_cast<MLLinOp&>(*m_reg_scal_solve_op));
#else
MLMG mlmg(*m_reg_scal_solve_op);
#endif
@@ -492,15 +518,20 @@
Vector<MultiFab const*> const& a_scalar,
Vector<MultiFab const*> const& a_eta,
Vector<MultiFab*> const& eb_dirichlet,
- amrex::Vector<amrex::BCRec> bcrec)
+ amrex::Vector<amrex::BCRec> bcrec,
+ bool is_temperature)
{
BL_PROFILE("DiffusionScalarOp::compute_laps");
+ amrex::ignore_unused(is_temperature);
int finest_level = m_incflo->finestLevel();
int n_comp = a_laps[0]->nComp();
#ifdef AMREX_USE_EB
- if (m_eb_scal_apply_op)
+ // Tracer and temperature must not share an EB op (see the header)
+ MLEBABecLap* eb_op = is_temperature ? m_eb_tem_apply_op.get()
+ : m_eb_scal_apply_op.get();
+ if (eb_op)
{
Vector<MultiFab> laps_tmp(finest_level+1);
int tmp_comp = (m_incflo->m_redistribution_type == "StateRedist") ? 3 : 2;
@@ -514,13 +545,13 @@
// We want to return div (mu grad)) phi
- m_eb_scal_apply_op->setScalars(0.0, -1.0);
+ eb_op->setScalars(0.0, -1.0);
// For when we use the stencil for centroid values
- // m_eb_scal_apply_op->setPhiOnCentroid();
+ // eb_op->setPhiOnCentroid();
for (int comp = 0; comp < n_comp; ++comp) {
- m_eb_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low,
+ eb_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low,
bcrec[comp].lo()),
m_incflo->get_diffuse_scalar_bc(Orientation::high,
bcrec[comp].hi()));
@@ -529,7 +560,7 @@
if ( m_incflo->m_has_mixedBC && comp>0 ){
// Must reset scalars to reuse solver with Robin BC
- m_eb_scal_apply_op->setScalars(0.0, -1.0);
+ eb_op->setScalars(0.0, -1.0);
}
Vector<MultiFab> laps_comp;
@@ -540,35 +571,31 @@
// Use the same EB BC as the implicit solve in diffuse_scalar,
// otherwise the explicit and implicit diffusion terms are
- // inconsistent at the EB. NOTE that, exactly as for the solve
- // op, this op is shared by the tracer and temperature passes and
- // MLEBABecLap has no way to revert an EB Dirichlet BC back to
- // Neumann, so if only one of tracer_eb/temperature_eb is defined
- // the other pass will see the EB phi left by the first.
+ // inconsistent at the EB.
if (lev < static_cast<int>(eb_dirichlet.size()) && !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);
+ eb_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]);
- m_eb_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
+ eb_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_robinBC_MFs(lev, &scalar_comp[lev]);
- m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev],
+ eb_op->setLevelBC(lev, &scalar_comp[lev],
&robin[0], &robin[1], &robin[2]);
}
else {
- m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]);
+ eb_op->setLevelBC(lev, &scalar_comp[lev]);
}
}
- MLMG mlmg(*m_eb_scal_apply_op);
+ MLMG mlmg(*eb_op);
mlmg.apply(GetVecOfPtrs(laps_comp), GetVecOfPtrs(scalar_comp));
}
--- a/src/diffusion/incflo_diffusion.cpp
+++ b/src/diffusion/incflo_diffusion.cpp
@@ -80,7 +80,7 @@
Vector<MultiFab const*> const& eta)
{
get_diffusion_scalar_op()->compute_laps(laps, scalar, eta, get_temperature_eb(),
- get_temperature_bcrec());
+ get_temperature_bcrec(), true);
}
void
@@ -103,7 +103,7 @@
get_diffusion_scalar_op()->diffuse_scalar(temperature, rhocp, eta,
get_temperature_eb(),
{1} /* use rhocp */,
- get_temperature_bcrec(), dt_diff);
+ get_temperature_bcrec(), dt_diff, true);
}
void
Activity
Fix issues AMReX-Fluids#240-AMReX-Fluids#249 (AMReX-Fluids#266)