Skip to content

Fix issues #240-#249 #266

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
asalmgren merged 8 commits into from
Sep 21, 2026
Merged

Fix issues #240-#249 #266

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
304eb2a
more precision changes
asalmgren May 18, 2026
1f80483
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren May 18, 2026
ff13b1d
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren Jul 15, 2026
f5fa182
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren Aug 21, 2026
62c0b00
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren Aug 21, 2026
ab031be
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren Sep 21, 2026
8dc0ae5
Merge branch 'development' of github.com:AMReX-Codes/incflo into deve…
asalmgren Sep 21, 2026
a876fe8
Fix issues #240-#249
asalmgren Sep 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 9 additions & 2 deletions src/diffusion/DiffusionScalarOp.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ public:
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,
Expand All @@ -33,7 +34,8 @@ public:
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,
Expand All @@ -49,6 +51,11 @@ private:
#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
Expand Down
109 changes: 68 additions & 41 deletions src/diffusion/DiffusionScalarOp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,15 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo)
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),
Expand All @@ -51,6 +60,15 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo)
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()) ||
Expand Down Expand Up @@ -135,8 +153,10 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
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
Expand Down Expand Up @@ -168,21 +188,27 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,

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) {
rhs_c[lev].define(a_scalar[lev]->boxArray(), a_scalar[lev]->DistributionMap(), 1, 0);
}

#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);
}
}
}
Expand All @@ -202,26 +228,26 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
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,
bcrec[0].lo()),
m_incflo->get_diffuse_scalar_bc(Orientation::high,
bcrec[0].hi()));
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) {
// Only reset Acoeff if necessary.
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);
}
}

Expand All @@ -230,12 +256,12 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
// 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
Expand Down Expand Up @@ -286,19 +312,19 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
}

#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],
&robin[0], &robin[1], &robin[2]);
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
{
Expand All @@ -307,7 +333,7 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& a_scalar,
}

#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
Expand Down Expand Up @@ -492,15 +518,20 @@ 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)
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;
Expand All @@ -514,22 +545,22 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,


// 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,
bcrec[comp].lo()),
m_incflo->get_diffuse_scalar_bc(Orientation::high,
bcrec[comp].hi()));
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()));

int eta_comp = comp;

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;
Expand All @@ -540,35 +571,31 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,

// 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],
&robin[0], &robin[1], &robin[2]);
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));
}

Expand Down
4 changes: 2 additions & 2 deletions src/diffusion/incflo_diffusion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ incflo::compute_laps_T(Vector<MultiFab *> const& laps,
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
Expand All @@ -103,7 +103,7 @@ incflo::diffuse_temperature(Vector<MultiFab *> const& temperature,
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
Expand Down
11 changes: 8 additions & 3 deletions src/incflo.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ void incflo::InitData ()

// xxxxx TODO averagedown ???

if (m_check_int > 0) { WriteCheckPointFile(); }
if (m_check_int > 0) { WriteCheckPointFile(); m_last_chk = 0; }

// Plot initial distribution
if (m_plot_int > 0 || m_plot_per_exact > 0 || m_plot_per_approx > 0)
Expand All @@ -98,19 +98,24 @@ void incflo::InitData ()
// Read starting configuration from chk file.
ReadCheckpointFile();

// The checkpoint just read is the checkpoint for this step; without this
// the final-output block in Evolve() rewrites it (renaming the original to
// *.old.*) when the restart does not take any steps.
m_last_chk = m_nstep;

#ifdef INCFLO_USE_PARTICLES
particleData.Redistribute();
#endif

if (m_plotfile_on_restart)
{
WritePlotFile();
m_last_plt = 0;
m_last_plt = m_nstep;
}
if (m_smallplotfile_on_restart)
{
WriteSmallPlotFile();
m_last_smallplt = 0;
m_last_smallplt = m_nstep;
}
}

Expand Down
4 changes: 3 additions & 1 deletion src/incflo_update_density.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,9 @@ void incflo::update_density (StepType step_type)
{
BL_PROFILE("incflo::update_density");

int ng = (step_type == StepType::Corrector) ? 0 : 1;
// One ghost cell in both passes: ApplyCCProjection averages density_nph to
// faces and reads its ghost cells, so the corrector must refresh them too.
const int ng = 1;

Real l_dt = m_dt;

Expand Down
2 changes: 1 addition & 1 deletion src/incflo_update_temperature.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ void incflo::update_temperature (StepType step_type, Vector<MultiFab>& tem_eta,
{
compute_temperature_diff_coeff(new_time, GetVecOfPtrs(tem_eta));
if (m_diff_type == DiffusionType::Explicit) {
compute_laps_T(get_laps_new(), get_temperature_new_const(), GetVecOfConstPtrs(tem_eta));
compute_laps_T(get_laps_tem_new(), get_temperature_new_const(), GetVecOfConstPtrs(tem_eta));
}
}

Expand Down
14 changes: 9 additions & 5 deletions src/particles/incflo_PCInit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,11 +238,12 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par
Real z_ctr = cyl_center[2];
#endif

// The cylinder axis matters here exactly as it does in AdvectWithFlow;
// assuming a z-parallel axis culls against the wrong axis for direction 0/1.
// The cylinder axis matters here exactly as it does in AdvectWithFlow,
// which accepts 0, 1 and 2 in both 2D and 3D; the culling below must
// use the same distance for each direction.
int cyl_direction;
pp.get("direction",cyl_direction);
AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction < AMREX_SPACEDIM);
AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= 2);

// Remove particles that are outside of the cylinder
for (ParIterType pti(*this, lev); pti.isValid(); ++pti)
Expand All @@ -264,8 +265,11 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par
: (cyl_direction == 1) ? std::sqrt(x*x + z*z)
: std::sqrt(x*x + y*y);
#else
// In 2D the only cylinder axis is the out-of-plane one
Real r = std::sqrt(x*x + y*y);
// Same convention as AdvectWithFlow and amrex::EB2::CylinderIF:
// direction 2 is a disk, directions 0 and 1 are slabs.
Real r = (cyl_direction == 2) ? std::sqrt(x*x + y*y)
: (cyl_direction == 1) ? std::abs(x)
: std::abs(y);
#endif

if (r > cyl_radius) {
Expand Down
Loading
Loading