Skip to content

Fix issues #179, #181-#185, #188-#190 #204

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鈥檒l occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
asalmgren merged 1 commit into from
Sep 27, 2026
Merged

Fix issues #179, #181-#185, #188-#190 #204

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
35 changes: 31 additions & 4 deletions Source/MacProj.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -561,6 +561,14 @@ MacProj::mac_sync_compute (int level,
forcing_term = std::make_unique<MultiFab>(grids, dmap, num_state_comps, NavierStokesBase::nghost_force());
divu_fp.reset(ns_level.getDivCond(NavierStokesBase::nghost_force(),prev_time));

// Get divu to time n+1/2, exactly as velocity_advection and
// scalar_advection do before predicting the edge states, so that the
// sync rebuilds the edge states the way the advance made them.
{
std::unique_ptr<MultiFab> dsdt(ns_level.getDsdt(NavierStokesBase::nghost_force(),prev_time));
MultiFab::Saxpy(*divu_fp, 0.5*dt, *dsdt, 0, 0, 1, NavierStokesBase::nghost_force());
}

MultiFab& Gp = ns_level.get_old_data(Gradp_Type);

visc_terms.setVal(0.0); // Initialize to make calls below safe
Expand Down Expand Up @@ -770,12 +778,20 @@ MacProj::mac_sync_compute (int level,
bool do_crse_add = false;
bool do_fine_add = update_fluxreg;

//
// ComputeAofs always builds an Array4 from the state MultiFab, even though
// it does not use the values when the edge states are known. So it must be
// handed a defined MultiFab; an alias of Sync has the right BoxArray and
// DistributionMap.
//
MultiFab Smf(Sync, amrex::make_alias, Sync_indx, ncomp);

//
// Compute the mac sync correction.
//
ns_level.ComputeAofs(Sync, /*Ssync_comp*/ Sync_indx, state_comp,
ncomp,
/*State*/ MultiFab(), /*S_comp*/ int(),//not used when known_edgestates
/*State*/ Smf, /*S_comp*/ 0, //not used when known_edgestates
/*forcing*/ nullptr, /*f_comp*/ int(), //not used when known_edgestates
/*constraint divU*/ nullptr, //not used when known_edgestates
fluxes, /*flux_comp*/ 0,
Expand Down Expand Up @@ -988,14 +1004,25 @@ MacProj::test_umac_periodic (int level,
Vector<IntVect> pshifts(27);
std::vector< std::pair<int,Box> > isects;

//
// MultiFabCopyDescriptor and the FArrayBox comparison below both operate
// on the host, so give them host-accessible copies of u_mac.
//
Array<MultiFab,AMREX_SPACEDIM> h_umac;

for (int dim = 0; dim < AMREX_SPACEDIM; dim++)
{
if (geom.isPeriodic(dim))
{
Box eDomain = amrex::surroundingNodes(geom.Domain(),dim);

mfid[dim] = mfcd.RegisterMultiFab(&u_mac[dim]);
h_umac[dim].define(u_mac[dim].boxArray(), u_mac[dim].DistributionMap(),
u_mac[dim].nComp(), u_mac[dim].nGrowVect(),
MFInfo().SetArena(The_Pinned_Arena()));
amrex::dtoh_memcpy(h_umac[dim], u_mac[dim]);
Gpu::streamSynchronize();

mfid[dim] = mfcd.RegisterMultiFab(&h_umac[dim]);

// How to combine pirm into one global pirm?
// don't think std::vector::push_back() is thread safe
Expand Down Expand Up @@ -1061,11 +1088,11 @@ MacProj::test_umac_periodic (int level,
AMREX_ASSERT(pirm_i.m_srcBox.sameSize(pirm_i.m_dstBox));
AMREX_ASSERT(u_mac[dim].DistributionMap()[pirm_i.m_idx] == ParallelDescriptor::MyProc());

diff.resize(pirm_i.m_srcBox, 1);
diff.resize(pirm_i.m_srcBox, 1, The_Pinned_Arena());

mfcd.FillFab(mfid[dim], pirm_i.m_fbid, diff);

diff.minus<RunOn::Host>(u_mac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);
diff.minus<RunOn::Host>(h_umac[dim][pirm_i.m_idx],pirm_i.m_dstBox,diff.box(),0,0,1);

const Real max_norm = diff.norm<RunOn::Host>(0);

Expand Down
167 changes: 119 additions & 48 deletions Source/NS_derive.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,25 +131,51 @@ namespace derive_functions
// Need to check if there are covered cells in neighbours --
// -- if so, use one-sided difference computation (but still quadratic)
if (!flag_fab(i,j,k).isConnected( 1,0,0)) {
vx = - (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i-1,j,k,1)
+ c2 * dat_arr(i-2,j,k,1)) * idx;
if (flag_fab(i,j,k).isConnected(-1,0,0)) {
if (flag_fab(i-1,j,k).isConnected(-1,0,0)) {
vx = - (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i-1,j,k,1)
+ c2 * dat_arr(i-2,j,k,1)) * idx;
} else {
// Only one uncovered cell that way, drop to linear
vx = (dat_arr(i ,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
}
}
// Covered on both sides, leave the derivative at zero
} else if (!flag_fab(i,j,k).isConnected(-1,0,0)) {
vx = (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i+1,j,k,1)
+ c2 * dat_arr(i+2,j,k,1)) * idx;
if (flag_fab(i+1,j,k).isConnected( 1,0,0)) {
vx = (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i+1,j,k,1)
+ c2 * dat_arr(i+2,j,k,1)) * idx;
} else {
// Only one uncovered cell that way, drop to linear
vx = (dat_arr(i+1,j,k,1) - dat_arr(i ,j,k,1)) * idx;
}
} else {
vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
}
// Do the same in y-direction
if (!flag_fab(i,j,k).isConnected( 0,1,0)) {
uy = - (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j-1,k,0)
+ c2 * dat_arr(i,j-2,k,0)) * idy;
if (flag_fab(i,j,k).isConnected(0,-1,0)) {
if (flag_fab(i,j-1,k).isConnected(0,-1,0)) {
uy = - (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j-1,k,0)
+ c2 * dat_arr(i,j-2,k,0)) * idy;
} else {
// Only one uncovered cell that way, drop to linear
uy = (dat_arr(i,j ,k,0) - dat_arr(i,j-1,k,0)) * idy;
}
}
// Covered on both sides, leave the derivative at zero
} else if (!flag_fab(i,j,k).isConnected(0,-1,0)) {
uy = (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j+1,k,0)
+ c2 * dat_arr(i,j+2,k,0)) * idy;
if (flag_fab(i,j+1,k).isConnected( 0,1,0)) {
uy = (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j+1,k,0)
+ c2 * dat_arr(i,j+2,k,0)) * idy;
} else {
// Only one uncovered cell that way, drop to linear
uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j ,k,0)) * idy;
}
} else {
uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy;
}
Expand All @@ -167,59 +193,104 @@ namespace derive_functions
// -- if so, use one-sided difference computation (but still quadratic)
if (!flag_fab(i,j,k).isConnected( 1,0,0)) {
// Covered cell to the right, go fish left
vx = - (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i-1,j,k,1)
+ c2 * dat_arr(i-2,j,k,1)) * idx;
wx = - (c0 * dat_arr(i ,j,k,2)
+ c1 * dat_arr(i-1,j,k,2)
+ c2 * dat_arr(i-2,j,k,2)) * idx;
if (flag_fab(i,j,k).isConnected(-1,0,0)) {
if (flag_fab(i-1,j,k).isConnected(-1,0,0)) {
vx = - (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i-1,j,k,1)
+ c2 * dat_arr(i-2,j,k,1)) * idx;
wx = - (c0 * dat_arr(i ,j,k,2)
+ c1 * dat_arr(i-1,j,k,2)
+ c2 * dat_arr(i-2,j,k,2)) * idx;
} else {
// Only one uncovered cell that way, drop to linear
vx = (dat_arr(i ,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
wx = (dat_arr(i ,j,k,2) - dat_arr(i-1,j,k,2)) * idx;
}
}
// Covered on both sides, leave the derivative at zero
} else if (!flag_fab(i,j,k).isConnected(-1,0,0)) {
// Covered cell to the left, go fish right
vx = (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i+1,j,k,1)
+ c2 * dat_arr(i+2,j,k,1)) * idx;
wx = (c0 * dat_arr(i ,j,k,2)
+ c1 * dat_arr(i+1,j,k,2)
+ c2 * dat_arr(i+2,j,k,2)) * idx;
if (flag_fab(i+1,j,k).isConnected( 1,0,0)) {
vx = (c0 * dat_arr(i ,j,k,1)
+ c1 * dat_arr(i+1,j,k,1)
+ c2 * dat_arr(i+2,j,k,1)) * idx;
wx = (c0 * dat_arr(i ,j,k,2)
+ c1 * dat_arr(i+1,j,k,2)
+ c2 * dat_arr(i+2,j,k,2)) * idx;
} else {
// Only one uncovered cell that way, drop to linear
vx = (dat_arr(i+1,j,k,1) - dat_arr(i ,j,k,1)) * idx;
wx = (dat_arr(i+1,j,k,2) - dat_arr(i ,j,k,2)) * idx;
}
} else {
// No covered cells right or left, use standard stencil
vx = 0.5 * (dat_arr(i+1,j,k,1) - dat_arr(i-1,j,k,1)) * idx;
wx = 0.5 * (dat_arr(i+1,j,k,2) - dat_arr(i-1,j,k,2)) * idx;
}
// Do the same in y-direction
if (!flag_fab(i,j,k).isConnected(0, 1,0)) {
uy = - (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j-1,k,0)
+ c2 * dat_arr(i,j-2,k,0)) * idy;
wy = - (c0 * dat_arr(i,j ,k,2)
+ c1 * dat_arr(i,j-1,k,2)
+ c2 * dat_arr(i,j-2,k,2)) * idy;
if (flag_fab(i,j,k).isConnected(0,-1,0)) {
if (flag_fab(i,j-1,k).isConnected(0,-1,0)) {
uy = - (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j-1,k,0)
+ c2 * dat_arr(i,j-2,k,0)) * idy;
wy = - (c0 * dat_arr(i,j ,k,2)
+ c1 * dat_arr(i,j-1,k,2)
+ c2 * dat_arr(i,j-2,k,2)) * idy;
} else {
// Only one uncovered cell that way, drop to linear
uy = (dat_arr(i,j ,k,0) - dat_arr(i,j-1,k,0)) * idy;
wy = (dat_arr(i,j ,k,2) - dat_arr(i,j-1,k,2)) * idy;
}
}
// Covered on both sides, leave the derivative at zero
} else if (!flag_fab(i,j,k).isConnected(0,-1,0)) {
uy = (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j+1,k,0)
+ c2 * dat_arr(i,j+2,k,0)) * idy;
wy = (c0 * dat_arr(i,j ,k,2)
+ c1 * dat_arr(i,j+1,k,2)
+ c2 * dat_arr(i,j+2,k,2)) * idy;
if (flag_fab(i,j+1,k).isConnected(0, 1,0)) {
uy = (c0 * dat_arr(i,j ,k,0)
+ c1 * dat_arr(i,j+1,k,0)
+ c2 * dat_arr(i,j+2,k,0)) * idy;
wy = (c0 * dat_arr(i,j ,k,2)
+ c1 * dat_arr(i,j+1,k,2)
+ c2 * dat_arr(i,j+2,k,2)) * idy;
} else {
// Only one uncovered cell that way, drop to linear
uy = (dat_arr(i,j+1,k,0) - dat_arr(i,j ,k,0)) * idy;
wy = (dat_arr(i,j+1,k,2) - dat_arr(i,j ,k,2)) * idy;
}
} else {
uy = 0.5 * (dat_arr(i,j+1,k,0) - dat_arr(i,j-1,k,0)) * idy;
wy = 0.5 * (dat_arr(i,j+1,k,2) - dat_arr(i,j-1,k,2)) * idy;
}
// Do the same in z-direction
if (!flag_fab(i,j,k).isConnected(0,0, 1)) {
uz = - (c0 * dat_arr(i,j,k ,0)
+ c1 * dat_arr(i,j,k-1,0)
+ c2 * dat_arr(i,j,k-2,0)) * idz;
vz = - (c0 * dat_arr(i,j,k ,1)
+ c1 * dat_arr(i,j,k-1,1)
+ c2 * dat_arr(i,j,k-2,1)) * idz;
if (flag_fab(i,j,k).isConnected(0,0,-1)) {
if (flag_fab(i,j,k-1).isConnected(0,0,-1)) {
uz = - (c0 * dat_arr(i,j,k ,0)
+ c1 * dat_arr(i,j,k-1,0)
+ c2 * dat_arr(i,j,k-2,0)) * idz;
vz = - (c0 * dat_arr(i,j,k ,1)
+ c1 * dat_arr(i,j,k-1,1)
+ c2 * dat_arr(i,j,k-2,1)) * idz;
} else {
// Only one uncovered cell that way, drop to linear
uz = (dat_arr(i,j,k ,0) - dat_arr(i,j,k-1,0)) * idz;
vz = (dat_arr(i,j,k ,1) - dat_arr(i,j,k-1,1)) * idz;
}
}
// Covered on both sides, leave the derivative at zero
} else if (!flag_fab(i,j,k).isConnected(0,0,-1)) {
uz = (c0 * dat_arr(i,j,k ,0)
+ c1 * dat_arr(i,j,k+1,0)
+ c2 * dat_arr(i,j,k+2,0)) * idz;
vz = (c0 * dat_arr(i,j,k ,1)
+ c1 * dat_arr(i,j,k+1,1)
+ c2 * dat_arr(i,j,k+2,1)) * idz;
if (flag_fab(i,j,k+1).isConnected(0,0, 1)) {
uz = (c0 * dat_arr(i,j,k ,0)
+ c1 * dat_arr(i,j,k+1,0)
+ c2 * dat_arr(i,j,k+2,0)) * idz;
vz = (c0 * dat_arr(i,j,k ,1)
+ c1 * dat_arr(i,j,k+1,1)
+ c2 * dat_arr(i,j,k+2,1)) * idz;
} else {
// Only one uncovered cell that way, drop to linear
uz = (dat_arr(i,j,k+1,0) - dat_arr(i,j,k ,0)) * idz;
vz = (dat_arr(i,j,k+1,1) - dat_arr(i,j,k ,1)) * idz;
}
} else {
uz = 0.5 * (dat_arr(i,j,k+1,0) - dat_arr(i,j,k-1,0)) * idz;
vz = 0.5 * (dat_arr(i,j,k+1,1) - dat_arr(i,j,k-1,1)) * idz;
Expand Down
19 changes: 13 additions & 6 deletions Source/NS_getForce.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,10 +92,13 @@ NavierStokesBase::getForce (FArrayBox& force,
#endif

// Compute min/max
for (int n=0; n<ncomp; n++) {
amrex::Print() << "State comp " << scomp+n << " min/max "
<< State.min<RunOn::Gpu>(scomp+n) << " / "
<< State.max<RunOn::Gpu>(scomp+n) << '\n';
// State may not hold all of the state components (scalar callers pass
// a velocity-only FAB), so index it by its own layout rather than by
// the absolute state component scomp+n.
for (int n=0; n<State.nComp(); n++) {
amrex::Print() << "State comp " << n << " min/max "
<< State.min<RunOn::Gpu>(n) << " / "
<< State.max<RunOn::Gpu>(n) << '\n';
}
for (int n=auxScomp; n<Aux.nComp(); n++) {
amrex::Print() << "aux comp " << n << " min/max "
Expand Down Expand Up @@ -178,10 +181,14 @@ NavierStokesBase::getForce (FArrayBox& force,

if (ParallelDescriptor::IOProcessor() && getForceVerbose) {
// Compute min/max
// For scalar-only calls force holds ncomp components starting at
// component 0 (see scomp_scal above), so index force by its own layout
// while still reporting the absolute state component in the label.
for (int n=0; n<ncomp; n++) {
const int fcomp = ( scomp<AMREX_SPACEDIM ) ? scomp+n : n;
amrex::Print() << "Force comp " << scomp+n << " min/max "
<< force.min<RunOn::Gpu>(scomp+n) << " / "
<< force.max<RunOn::Gpu>(scomp+n) << '\n';
<< force.min<RunOn::Gpu>(fcomp) << " / "
<< force.max<RunOn::Gpu>(fcomp) << '\n';
}

amrex::Print() << "NavierStokesBase::getForce(): Leaving..."
Expand Down
28 changes: 18 additions & 10 deletions Source/NavierStokesBase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2774,17 +2774,15 @@ NavierStokesBase::scalar_advection_update (Real dt,
// Must do FillPatch here instead of MF iterator because we need the
// boundary values in the old data (especially at inflow)
//
const int index_new_s = Density;
const int index_new_rho = Density;
const int index_old_s = index_new_s - Density;
const int index_old_rho = index_new_rho - Density;

FillPatchIterator S_fpi(*this,S_old,1,prev_time,State_Type,Density,1);
MultiFab& Smf=S_fpi.get_mf();

ConservativeScalMinMax(S_new, index_new_s, index_new_rho,
Smf, index_old_s, index_old_rho);

//
// Clamp the new density to the min/max of the old-time density in
// the neighborhood. (The conservative s/rho form with s = rho is
// an identity, i.e. a no-op.)
//
ConvectiveScalMinMax(S_new, Density, Smf, 0);
}

++sComp;
Expand Down Expand Up @@ -4277,8 +4275,13 @@ NavierStokesBase::ConservativeScalMinMax ( amrex::MultiFab& Snew, const in
amrex::ParallelFor(bx, [=]
AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
#ifdef AMREX_USE_EB
// Leave covered cells alone; their stencil holds no valid data.
if ( vfrac(i,j,k) == 0. ) { return; }
#endif

Real smn = std::numeric_limits<Real>::max();
Real smx = std::numeric_limits<Real>::min();
Real smx = std::numeric_limits<Real>::lowest();

#if (AMREX_SPACEDIM==3)
int ks = -1;
Expand Down Expand Up @@ -4335,8 +4338,13 @@ NavierStokesBase::ConvectiveScalMinMax ( amrex::MultiFab& Snew, const int
amrex::ParallelFor(bx, [=]
AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
#ifdef AMREX_USE_EB
// Leave covered cells alone; their stencil holds no valid data.
if ( vfrac(i,j,k) == 0. ) { return; }
#endif

Real smn = std::numeric_limits<Real>::max();
Real smx = std::numeric_limits<Real>::min();
Real smx = std::numeric_limits<Real>::lowest();

#if (AMREX_SPACEDIM==3)
int ks = -1;
Expand Down
3 changes: 2 additions & 1 deletion Source/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,8 @@ main (int argc,
//
if (Amr::RegridOnRestart())
{
if ( (amrptr->levelSteps(0) >= max_step ) ||
if ( ( (max_step >= 0) &&
(amrptr->levelSteps(0) >= max_step) ) ||
( (stop_time >= 0.0) &&
(amrptr->cumTime() >= stop_time) ) )
{
Expand Down
Loading
Loading