Skip to content

Commit

Permalink
fix bcs with terrain and clean up time integrator
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Nov 1, 2024
1 parent 7262816 commit a93bf1b
Show file tree
Hide file tree
Showing 13 changed files with 89 additions and 146 deletions.
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_BoundaryConditions_basestate.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,7 +259,7 @@ void ERFPhysBCFunct_base::impose_vertical_basestate_bcs (const Array4<Real>& des
{
dest_arr(i,j,k,BaseState::r0_comp) = dest_arr(i,j,dom_lo.z,BaseState::r0_comp);
dest_arr(i,j,k,BaseState::p0_comp) = p_0 -
dest_arr(i,j,k,BaseState::r0_comp) * hz * 9.81;
dest_arr(i,j,k,BaseState::r0_comp) * hz * CONST_GRAV;
dest_arr(i,j,k,BaseState::pi0_comp) = dest_arr(i,j,dom_lo.z,BaseState::pi0_comp);
dest_arr(i,j,k,BaseState::th0_comp) = dest_arr(i,j,dom_lo.z,BaseState::th0_comp);
}
Expand Down
4 changes: 2 additions & 2 deletions Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
const Array4<Real const>& z_phys_nd,
const GpuArray<Real,AMREX_SPACEDIM> dxInv,
int icomp, int ncomp)
int icomp, int ncomp, bool do_terrain_adjustment)
{
BL_PROFILE_VAR("impose_vertical_cons_bcs()",impose_vertical_cons_bcs);
const auto& dom_lo = lbound(domain);
Expand Down Expand Up @@ -416,7 +416,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
);
}

if (m_z_phys_nd) {
if (do_terrain_adjustment && m_z_phys_nd) {
const auto& bx_lo = lbound(bx);
const auto& bx_hi = ubound(bx);
const BCRec* bc_ptr_h = bcrs.data();
Expand Down
14 changes: 9 additions & 5 deletions Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,11 +61,13 @@ ERF::FillIntermediatePatch (int lev, Real time,
}
}

// amrex::Print() << "CONS ONLY " << cons_only << std::endl;
// amrex::Print() << "ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl;
// amrex::Print() << "LEVEL " << lev << " CONS ONLY " << cons_only <<
// " ICOMP NCOMP " << icomp_cons << " " << ncomp_cons << " NGHOST " << ng_cons << std::endl;

AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes);
AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes);
if (!cons_only) {
AMREX_ALWAYS_ASSERT(mfs_mom.size() == IntVars::NumTypes);
AMREX_ALWAYS_ASSERT(mfs_vel.size() == Vars::NumTypes);
}

// Enforce no penetration for thin immersed body
if (!cons_only) {
Expand Down Expand Up @@ -122,7 +124,9 @@ ERF::FillIntermediatePatch (int lev, Real time,
Vector<Real> ftime = {time,time};

// Impose physical bc's on fine data (note time and 0 are not used)
(*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true);
bool do_fb = true; bool do_terrain_adjustment = false;
(*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,
do_fb, do_terrain_adjustment);

if ( (icomp_cons+ncomp_cons > 1) && (interpolation_type == StateInterpType::Perturbational) )
{
Expand Down
4 changes: 2 additions & 2 deletions Source/BoundaryConditions/ERF_PhysBCFunct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ public:
*/
void operator() (amrex::MultiFab& mf, int icomp, int ncomp,
amrex::IntVect const& nghost, const amrex::Real time, int bccomp_cons,
bool do_fb);
bool do_fb = true, bool do_terrain_adjustment = true);

void impose_lateral_cons_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain,
Expand All @@ -62,7 +62,7 @@ public:
const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<amrex::Real const>& z_nd,
const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dxInv,
int icomp, int ncomp);
int icomp, int ncomp, bool do_terrain_adjustment = true);

private:
int m_lev;
Expand Down
68 changes: 6 additions & 62 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using namespace amrex;

void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
IntVect const& nghost, const Real /*time*/, int /*bccomp*/,
bool do_fb)
bool do_fb, bool do_terrain_adjustment)
{
BL_PROFILE("ERFPhysBCFunct_cons::()");

Expand All @@ -32,21 +32,6 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
}
}

MultiFab z_nd_mf_loc;
if (m_z_phys_nd) {
m_z_phys_nd->FillBoundary(m_geom.periodicity());
BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
for (auto& b : bl_z_phys) {
b.setSmall(2,0);
b.setBig(2,1);
}
BoxArray ba_z(std::move(bl_z_phys));

z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0));
z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
z_nd_mf_loc.nGrowVect());
}

//
// We fill all of the interior and periodic ghost cells first, so we can fill
// those directly inside the lateral and vertical calls.
Expand Down Expand Up @@ -77,7 +62,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,

if (m_z_phys_nd)
{
z_nd_arr = z_nd_mf_loc.const_array(mfi);
z_nd_arr = m_z_phys_nd->const_array(mfi);
}

if (!gdomain.contains(cbx2))
Expand All @@ -91,7 +76,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
}

// We send the full FAB box with ghost cells
impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp);
impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,do_terrain_adjustment);
}

} // MFIter
Expand All @@ -116,20 +101,6 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
}
}

MultiFab z_nd_mf_loc;
if (m_z_phys_nd) {
m_z_phys_nd->FillBoundary(m_geom.periodicity());
BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
for (auto& b : bl_z_phys) {
b.setSmall(2,0);
b.setBig(2,1);
}
BoxArray ba_z(std::move(bl_z_phys));
z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0]+1,nghost[1],0));
z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
z_nd_mf_loc.nGrowVect());
}

//
// We fill all of the interior and periodic ghost cells first, so we can fill
// those directly inside the lateral and vertical calls.
Expand Down Expand Up @@ -163,7 +134,7 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,

if (m_z_phys_nd)
{
z_nd_arr = z_nd_mf_loc.const_array(mfi);
z_nd_arr = m_z_phys_nd->const_array(mfi);
}

if (!gdomainx.contains(xbx2))
Expand Down Expand Up @@ -202,20 +173,6 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
}
}

MultiFab z_nd_mf_loc;
if (m_z_phys_nd) {
m_z_phys_nd->FillBoundary(m_geom.periodicity());
BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
for (auto& b : bl_z_phys) {
b.setSmall(2,0);
b.setBig(2,1);
}
BoxArray ba_z(std::move(bl_z_phys));
z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1]+1,0));
z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
z_nd_mf_loc.nGrowVect());
}

//
// We fill all of the interior and periodic ghost cells first, so we can fill
// those directly inside the lateral and vertical calls.
Expand Down Expand Up @@ -249,7 +206,7 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,

if (m_z_phys_nd)
{
z_nd_arr = z_nd_mf_loc.const_array(mfi);
z_nd_arr = m_z_phys_nd->const_array(mfi);
}

if (!gdomainy.contains(ybx2))
Expand Down Expand Up @@ -295,19 +252,6 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,

Box ndomain = convert(domain,IntVect(1,1,1));

MultiFab z_nd_mf_loc;
if (m_z_phys_nd) {
BoxList bl_z_phys = convert(mf.boxArray(),IntVect(1,1,1)).boxList();
for (auto& b : bl_z_phys) {
b &= ndomain;
}
BoxArray ba_z(std::move(bl_z_phys));
z_nd_mf_loc.define(ba_z,mf.DistributionMap(),1,IntVect(nghost[0],nghost[1],0));
z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,m_z_phys_nd->nGrowVect(),
z_nd_mf_loc.nGrowVect());
}
z_nd_mf_loc.FillBoundary(m_geom.periodicity());

//
// We fill all of the interior and periodic ghost cells first, so we can fill
// those directly inside the lateral and vertical calls.
Expand Down Expand Up @@ -339,7 +283,7 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,

if (m_z_phys_nd)
{
z_nd_arr = z_nd_mf_loc.const_array(mfi);
z_nd_arr = m_z_phys_nd->const_array(mfi);
}

//
Expand Down
16 changes: 12 additions & 4 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
// Define dmap[lev] to be dm
SetDistributionMap(lev, dm);

// amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl;
if (verbose) {
amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl;
}

if (lev == 0) init_bcs();

Expand Down Expand Up @@ -179,7 +181,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
{
AMREX_ALWAYS_ASSERT(lev > 0);

// amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl;
if (verbose) {
amrex::Print() <<" NEW BA FROM COARSE AT LEVEL " << lev << " " << ba << std::endl;
}

//********************************************************************************************
// This allocates all kinds of things, including but not limited to: solution arrays,
Expand Down Expand Up @@ -283,15 +287,19 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
void
ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm)
{
// amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl;
if (verbose) {
amrex::Print() <<" REMAKING WITH NEW BA AT LEVEL " << lev << " " << ba << std::endl;
}

AMREX_ALWAYS_ASSERT(lev > 0);
AMREX_ALWAYS_ASSERT(solverChoice.terrain_type != TerrainType::Moving);

BoxArray ba_old(vars_new[lev][Vars::cons].boxArray());
DistributionMapping dm_old(vars_new[lev][Vars::cons].DistributionMap());

// amrex::Print() <<" OLD BA AT LEVEL " << lev << " " << ba_old << std::endl;
if (verbose) {
amrex::Print() <<" OLD BA AT LEVEL " << lev << " " << ba_old << std::endl;
}

int ncomp_cons = vars_new[lev][Vars::cons].nComp();
IntVect ngrow_state = vars_new[lev][Vars::cons].nGrowVect();
Expand Down
17 changes: 11 additions & 6 deletions Source/IO/ERF_Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1410,10 +1410,15 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p

int lev0 = 0;
int desired_ratio = std::max(std::max(ref_ratio[lev0][0],ref_ratio[lev0][1]),ref_ratio[lev0][2]);
bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 0) || (ref_ratio[lev0][1] == 0) ) ||
(ref_ratio[lev0][2] == 0) );
bool any_ratio_one = ( ( (ref_ratio[lev0][0] == 1) || (ref_ratio[lev0][1] == 1) ) ||
(ref_ratio[lev0][2] == 1) );
for (int lev = 1; lev < finest_level; lev++) {
any_ratio_one = any_ratio_one ||
( ( (ref_ratio[lev][0] == 1) || (ref_ratio[lev][1] == 1) ) ||
(ref_ratio[lev][2] == 1) );
}

if (any_ratio_one == 1 && m_expand_plotvars_to_unif_rr)
if (any_ratio_one && m_expand_plotvars_to_unif_rr)
{
Vector<IntVect> r2(finest_level);
Vector<Geometry> g2(finest_level+1);
Expand All @@ -1435,9 +1440,9 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p

for (int lev = 1; lev <= finest_level; ++lev) {
if (lev > 1) {
r2[lev-1][0] *= desired_ratio / ref_ratio[lev-1][0];
r2[lev-1][1] *= desired_ratio / ref_ratio[lev-1][1];
r2[lev-1][2] *= desired_ratio / ref_ratio[lev-1][2];
r2[lev-1][0] = r2[lev-2][0] * desired_ratio / ref_ratio[lev-1][0];
r2[lev-1][1] = r2[lev-2][1] * desired_ratio / ref_ratio[lev-1][1];
r2[lev-1][2] = r2[lev-2][2] * desired_ratio / ref_ratio[lev-1][2];
}

mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0);
Expand Down
42 changes: 8 additions & 34 deletions Source/TimeIntegration/ERF_MRI.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,11 +56,8 @@ private:
int force_stage1_single_substep;

/**
* \brief The pre_update function is called by the integrator on stage data before using it to evaluate a right-hand side.
* \brief The post_update function is called by the integrator on stage data at the end of the stage
* \brief The no_substep function is called when we have no acoustic substepping
*/
std::function<void (T&, int)> pre_update;
std::function<void (T&, amrex::Real, int, int)> post_update;
std::function<void (T&, T&, T&, amrex::Real, amrex::Real, int)> no_substep;


Expand Down Expand Up @@ -164,16 +161,6 @@ public:
return slow_fast_timestep_ratio;
}

void set_pre_update (std::function<void (T&, int)> F)
{
pre_update = F;
}

void set_post_update (std::function<void (T&, amrex::Real, int, int)> F)
{
post_update = F;
}

void set_no_substep (std::function<void (T&, T&, T&, amrex::Real, amrex::Real, int)> F)
{
no_substep = F;
Expand Down Expand Up @@ -247,8 +234,6 @@ public:
// RK3 for compressible integrator
for (int nrk = 0; nrk < 3; nrk++)
{
// amrex::Print() << "Starting RK3: Step " << nrk+1 << std::endl;

// Capture the time we got to in the previous RK step
old_time_stage = time_stage;

Expand Down Expand Up @@ -285,12 +270,6 @@ public:
// step 2 starts with S_stage = S^* and we always start substepping at the old time
// step 3 starts with S_stage = S^** and we always start substepping at the old time

// All pre_update does is call cons_to_prim, and we have done this with the old
// data already before starting the RK steps
if (nrk > 0) {
pre_update(S_new, S_new[IntVars::cons].nGrow());
}

// S_scratch also holds the average momenta over the fast iterations --
// to be used to update the slow variables -- we will initialize with
// the momenta used in the first call to the slow_rhs, then update
Expand Down Expand Up @@ -327,10 +306,6 @@ public:
// (because we did update the fast variables in the substepping)
// ****************************************************
slow_rhs_post(*F_slow, S_old, S_new, *S_sum, *S_scratch, time, old_time_stage, time_stage, nrk);

// Call the post-update hook for S_new after all the fast steps completed
// This will update S_prim that is used in the slow RHS
post_update(S_new, time + nsubsteps*dtau, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow());
} // nrk

} else {
Expand All @@ -343,12 +318,6 @@ public:
if (nrk == 0) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; }
if (nrk == 1) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; }

// All pre_update does is call cons_to_prim, and we have done this with the old
// data already before starting the RK steps
if (nrk > 0) {
pre_update(S_new, S_new[IntVars::cons].nGrow());
}

// S_scratch also holds the average momenta over the fast iterations --
// to be used to update the slow variables -- we will initialize with
// the momenta used in the first call to the slow_rhs, then update
Expand All @@ -358,9 +327,14 @@ public:

no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk);

// ****************************************************
// Evaluate F_slow(S_stage) only for the slow variables
// Note that we are using the current stage versions (in S_new) of the slow variables
// (because we didn't update the slow variables in the substepping)
// but we are using the "new" versions (in S_sum) of the velocities
// (because we did update the fast variables in the substepping)
// ****************************************************
slow_rhs_post(*F_slow, S_old, S_new, *S_sum, *S_scratch, time, old_time_stage, time_stage, nrk);

post_update(S_new, time + nsubsteps*dtau, S_new[IntVars::cons].nGrow(), S_new[IntVars::xmom].nGrow());
} // nrk
}

Expand Down
4 changes: 3 additions & 1 deletion Source/TimeIntegration/ERF_TI_fast_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,9 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
const Real new_substep_time)
{
BL_PROFILE("fast_rhs_fun");
if (verbose) amrex::Print() << "Calling fast rhs at level " << level << " with dt = " << dtau << std::endl;
if (verbose) amrex::Print() << "Fast time integration at level " << level
<< " from " << old_substep_time << " to " << new_substep_time
<< " with dt = " << dtau << std::endl;

// Define beta_s here so that it is consistent between where we make the fast coefficients
// and where we use them
Expand Down
Loading

0 comments on commit a93bf1b

Please sign in to comment.