Skip to content

Fix issues #221, #223, #225, #228, #231, #232, #234, #235, #237, #238… #257

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 1 commit into from
Oct 8, 2026
Merged

Fix issues #221, #223, #225, #228, #231, #232, #234, #235, #237, #238… #257

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion Exec/eb_run2d/regtest.2d.bubble
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,4 +133,4 @@ prob.blob_radius = 0.2
prob.density_ic = 2.0

#mag_vort is stepping out of bounds
#amr.derive_plot_vars = mag_vort diveru avg_pressure
#amr.derive_plot_vars = mag_vort avg_pressure
10 changes: 0 additions & 10 deletions Source/MacProj.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,6 @@ public:
amrex::Real dt,
amrex::Real prev_time,
const amrex::MultiFab& divu,
int have_divu,
const amrex::BCRec& density_math_bc,
bool increment_vel_register = true );

Expand Down Expand Up @@ -104,14 +103,6 @@ public:
void check_div_cond (int level,
amrex::MultiFab U_edge[]) const;
//
// Boundary conditions.
//
void set_outflow_bcs (int level,
amrex::MultiFab* mac_phi,
const amrex::MultiFab* u_mac,
const amrex::MultiFab& S,
const amrex::MultiFab& divu);
//
// Pointers to amr,amrlevel.
//
amrex::Amr* parent;
Expand All @@ -132,7 +123,6 @@ public:
int finest_level_allocated;

static int verbose;
static int do_outflow_bcs;
static amrex::Real mac_tol;
static amrex::Real mac_abs_tol;
static amrex::Real mac_sync_tol;
Expand Down
144 changes: 8 additions & 136 deletions Source/MacProj.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
#include <AMReX_ParmParse.H>
#include <MacProj.H>
#include <NavierStokesBase.H>
#include <OutFlowBC.H>
#include <hydro_MacProjector.H>

#ifdef AMREX_USE_EB
Expand All @@ -25,7 +24,6 @@ int MacProj::verbose;
Real MacProj::mac_tol;
Real MacProj::mac_abs_tol;
Real MacProj::mac_sync_tol;
int MacProj::do_outflow_bcs;
int MacProj::check_umac_periodicity;
int MacProj::max_order = 4;
int MacProj::agglomeration = 1;
Expand All @@ -49,7 +47,6 @@ MacProj::Initialize ()
MacProj::mac_tol = 1.0e-12;
MacProj::mac_abs_tol = 1.0e-16;
MacProj::mac_sync_tol = 1.0e-10;
MacProj::do_outflow_bcs = 1;
//
// Only check umac periodicity when debugging. Can be overridden on input.
//
Expand All @@ -65,7 +62,6 @@ MacProj::Initialize ()
pp.query("mac_tol", mac_tol);
pp.query("mac_abs_tol", mac_abs_tol);
pp.query("mac_sync_tol", mac_sync_tol);
pp.query("do_outflow_bcs", do_outflow_bcs);
pp.query("check_umac_periodicity", check_umac_periodicity);
pp.query("umac_periodic_test_Tol", umac_periodic_test_Tol);

Expand Down Expand Up @@ -184,38 +180,6 @@ MacProj::cleanup (int level)
//
// Projection functions follow ...
//
static
bool
grids_on_side_of_domain (const BoxArray& grids,
const Box& domain,
const Orientation& outFace)
{
const int idir = outFace.coordDir();

if (outFace.isLow())
{
for (int igrid = 0; igrid < grids.size(); igrid++)
{
if (grids[igrid].smallEnd(idir) == domain.smallEnd(idir))
{
return true;
}
}
}

if (outFace.isHigh())
{
for (int igrid = 0; igrid < grids.size(); igrid++)
{
if (grids[igrid].bigEnd(idir) == domain.bigEnd(idir))
{
return true;
}
}
}

return false;
}

//
// Compute the level advance mac projection.
Expand All @@ -228,7 +192,6 @@ MacProj::mac_project (int level,
Real dt,
Real time,
const MultiFab& divu,
int have_divu,
const BCRec& density_math_bc,
bool increment_vel_register )
{
Expand Down Expand Up @@ -262,10 +225,6 @@ MacProj::mac_project (int level,
const MultiFab& rhotime = ns.get_rho(time);
MultiFab::Copy(S, rhotime, 0, Density, 1, 1);

if (OutFlowBC::HasOutFlowBC(phys_bc) && have_divu && do_outflow_bcs) {
set_outflow_bcs(level, mac_phi, u_mac, S, divu);
}

//
// Set up the mac projection
//
Expand Down Expand Up @@ -861,100 +820,6 @@ MacProj::check_div_cond (int level,
}
}

void
MacProj::set_outflow_bcs (int level,
MultiFab* mac_phi,
const MultiFab* /*u_mac*/,
const MultiFab& /*S*/,
const MultiFab& /*divu*/)
{
//
// This code is very similar to the outflow BC stuff in the Projection
// class except that here the the phi to be solved for lives on the
// out-directed faces. The projection equation to satisfy is
//
// (1/r)(d/dr)[r/rho dphi/dr] = dv/dr - S
//
bool hasOutFlow;
Orientation outFaces[2*AMREX_SPACEDIM];
int numOutFlowFaces;

OutFlowBC::GetOutFlowFaces(hasOutFlow,outFaces,phys_bc,numOutFlowFaces);

const BoxArray& grids = LevelData[level]->boxArray();
const Geometry& geom = parent->Geom(level);
const Box& domain = parent->Geom(level).Domain();
//
// Create 1-wide cc box just outside boundary to hold phi.
//
BoxList ccBoxList, phiBoxList;
// numOutFlowFaces gives the number of outflow faces on the entire
// problem domain
// nOutFlowTouched gives the number of outflow faces a level touches, so
// nOutFlowTouched = numOutFlowFaces for level 0, but
// nOutFlowTouched <= numOutFlowFaces for levels > 0, since
// finer levels may not span the entire problem domain
int nOutFlowTouched = 0;
for (int iface = 0; iface < numOutFlowFaces; iface++)
{
if (grids_on_side_of_domain(grids,geom.Domain(),outFaces[iface]))
{
nOutFlowTouched++;
const int outDir = outFaces[iface].coordDir();

Box ccBndBox;
if (outFaces[iface].faceDir() == Orientation::high)
{
ccBndBox = amrex::adjCellHi(domain,outDir,2);
ccBndBox.shift(outDir,-2);
}
else
{
ccBndBox = amrex::adjCellLo(domain,outDir,2);
ccBndBox.shift(outDir,2);
}
ccBoxList.push_back(ccBndBox);

Box phiBox = amrex::adjCell(domain,outFaces[iface],1);
phiBoxList.push_back(phiBox);

const Box& valid_ccBndBox = ccBndBox & domain;
const BoxArray uncovered_outflow_ba = amrex::complementIn(valid_ccBndBox,grids);

if ((!uncovered_outflow_ba.empty()) &&
grids.intersects(valid_ccBndBox))
{
amrex::Error("MacProj: Cannot yet handle partially refined outflow");
}
}
}

if ( !ccBoxList.isEmpty() )
{
BoxArray phiBoxArray(phiBoxList);
phiBoxList.clear();

//
// Must do this kind of copy instead of mac_phi->copy(phidat);
// because we're copying onto the ghost cells of the FABs,
// not the valid regions.
//
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for ( int iface = 0; iface < nOutFlowTouched; ++iface )
{
for (MFIter mfi(*mac_phi); mfi.isValid(); ++mfi)
{
Box ovlp = (*mac_phi)[mfi].box() & phiBoxArray[iface];
if (ovlp.ok()) {
(*mac_phi)[mfi].setVal<RunOn::Gpu>(0,ovlp,0,1);
}
}
}
}
}

//
// Structure used by test_umac_periodic().
//
Expand Down Expand Up @@ -1188,7 +1053,14 @@ MacProj::mlmg_mac_solve (Amr* a_parent, const MultiFab* cphi, const BCRec& a_phy
set_mac_solve_bc(mlmg_lobc, mlmg_hibc, a_phys_bc, geom);

macproj.setDomainBC(mlmg_lobc, mlmg_hibc);
if (level > 0 && cphi)
//
// cphi is null for the sync solve, which wants a homogeneous Dirichlet
// coarse/fine BC. MLMG still needs to be told the true refinement
// ratio, however, or it assumes 2 and mis-positions the coarse ghost
// value. So register the ratio on every level above 0, with or without
// coarse data.
//
if (level > 0)
{
macproj.setCoarseFineBC(cphi, a_parent->refRatio(level-1)[0]);
}
Expand Down
10 changes: 10 additions & 0 deletions Source/NS_LES.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,16 @@ NavierStokesBase::calc_mut_LES(MultiFab* mu_LES[AMREX_SPACEDIM], const Real time
MultiFab** tensorflux = fb.get();
std::array<MultiFab*,AMREX_SPACEDIM> grad_Uvel{AMREX_D_DECL(tensorflux[0], tensorflux[1], tensorflux[2])};

//
// FluxBoxes does not initialize the face MultiFabs, and compVelGrad skips
// covered FABs (covered tiles when tiling). The loop below evaluates the
// LES formula everywhere, so zero the gradients first to get mu_t = 0 on
// covered faces rather than whatever the arena held.
//
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
grad_Uvel[idim]->setVal(0.0);
}

tensorop.compVelGrad(0,{grad_Uvel},{Uvel},MLLinOp::Location::FaceCenter);

//
Expand Down
20 changes: 20 additions & 0 deletions Source/NavierStokes.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,14 @@ NavierStokes::Initialize ()

NavierStokesBase::Initialize();

{
//
// Documented in RunningProblems.rst; only used with EB.
//
ParmParse pp("ns");
pp.query("set_plot_coveredCell_val", set_plot_coveredCell_val);
}

//
// Set number of state variables.
//
Expand Down Expand Up @@ -95,6 +103,18 @@ NavierStokes::Initialize_bcs ()
pp.getarr("hi_bc",hi_bc,0,AMREX_SPACEDIM);
for (int i = 0; i < AMREX_SPACEDIM; i++)
{
//
// phys_bc indexes the six-entry tables in NS_BC.H, so anything
// outside [interior,noslipwall] would read past them. In particular
// PhysBCType also defines inflowoutflow=6, which IAMR does not
// support, and the integer input path is otherwise unchecked.
//
if ( lo_bc[i] < PhysBCType::interior || lo_bc[i] > PhysBCType::noslipwall ||
hi_bc[i] < PhysBCType::interior || hi_bc[i] > PhysBCType::noslipwall )
{
amrex::Abort("NavierStokes::Initialize_bcs: ns.lo_bc/ns.hi_bc must be in [0,5]; see the list of physical BC types in the documentation.");
}

phys_bc.setLo(i,lo_bc[i]);
phys_bc.setHi(i,hi_bc[i]);
}
Expand Down
50 changes: 17 additions & 33 deletions Source/NavierStokesBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,37 +55,6 @@ struct HomExtDirFill
}
};

//
// A dummy function because FillPatch requires something to exist for filling dirichlet boundary conditions,
// even if we know we cannot have an ext_dir BC.
// u_mac BCs are only either periodic (BCType::int_dir) or first order extrapolation (FOEXTRAP).
//
struct umacFill
{
AMREX_GPU_DEVICE
void operator()(
const amrex::IntVect& /*iv*/,
amrex::Array4<amrex::Real> const& /*dummy*/,
const int /*dcomp*/,
const int numcomp,
amrex::GeometryData const& /*geom*/,
const amrex::Real /*time*/,
const amrex::BCRec* bcr,
const int bcomp,
const int /*orig_comp*/) const
{
// Abort if this function is expected to fill an ext_dir BC.
for (int n = bcomp; n < bcomp+numcomp; ++n) {
const amrex::BCRec& bc = bcr[n];
if ( AMREX_D_TERM( bc.lo(0) == amrex::BCType::ext_dir || bc.hi(0) == amrex::BCType::ext_dir,
|| bc.lo(1) == amrex::BCType::ext_dir || bc.hi(1) == amrex::BCType::ext_dir,
|| bc.lo(2) == amrex::BCType::ext_dir || bc.hi(2) == amrex::BCType::ext_dir ) ) {
amrex::Abort("NavierStokesBase::umacFill: umac should not have BCType::ext_dir");
}
}
}
};


BCRec NavierStokesBase::phys_bc;
Projection* NavierStokesBase::projector = nullptr;
Expand Down Expand Up @@ -2146,7 +2115,7 @@ NavierStokesBase::mac_project (Real time,

Vector<BCRec> density_math_bc = fetchBCArray(State_Type,Density,1);

mac_projector->mac_project(level,u_mac,S_old,dt,time,*divu,have_divu,
mac_projector->mac_project(level,u_mac,S_old,dt,time,*divu,
density_math_bc[0], increment_vel_register);

create_umac_grown(ngrow, divu);
Expand Down Expand Up @@ -2758,7 +2727,7 @@ NavierStokesBase::set_state_in_checkpoint (Vector<int>& state_in_checkpoint)
// Abort if any of the NSB::*_in_checkpoint variables haven't been set by user.
//
if ( gradp_in_checkpoint<0 || average_in_checkpoint<0 )
Abort("\n\n Checkpoint file is missing one or more state types. Set both\n ns.gradp_in_checkpoint and ns.avg_in_checkpoint to identify missing\n data. Set to 1 if present in checkpoint, 0 if not present. If unsure,\n try setting both to 0.\n\n If you just activated Time Averaging, you should add \n ns.avg_in_checkpoint=0 ns.gradp_in_checkpoint=1 \n\n");
Abort("\n\n Checkpoint file is missing one or more state types. Set both\n ns.gradp_in_checkpoint and ns.avg_in_checkpoint to identify missing\n data. Set to 1 if present in checkpoint, 0 if not present. These must\n flag exactly the missing types: the checkpoint is read as a sequential\n stream, so marking a type that is actually present shifts every later\n read and silently restarts from mis-assigned data.\n\n If you just activated Time Averaging, you should add \n ns.avg_in_checkpoint=0 ns.gradp_in_checkpoint=1 \n\n");

//
// Tell AmrLevel which types are in the checkpoint, so it knows what to copy.
Expand Down Expand Up @@ -3987,6 +3956,21 @@ NavierStokesBase::post_timestep_particle (int crse_iteration)
n = timestamp_indices.size();
nextras = timestamp_num_extras();

//
// These index the FillPatched state below. NUM_STATE is not
// yet known when read_particle_params() reads the list, so
// validate here: BaseFab::copy only range-checks the
// component with an AMREX_ASSERT, so a bad entry would read
// outside the FAB in a release build.
//
for (int i = 0; i < n; ++i)
{
if (timestamp_indices[i] < 0 || timestamp_indices[i] >= NUM_STATE)
{
amrex::Abort("NavierStokesBase::post_timestep_particle: particles.timestamp_indices entries must be in [0,NUM_STATE)");
}
}

int sz = n + nextras;
tindices.reserve(sz);

Expand Down
Loading
Loading