Skip to content

Commit

Permalink
This version runs both compressible and anelastic with factor 2 refin…
Browse files Browse the repository at this point in the history 
…ement of the DC problem to 900s
  • Loading branch information
@asalmgren
asalmgren committed Oct 15, 2024
1 parent 4ce47d3 commit 88d2efd
Show file tree
Hide file tree
Showing 13 changed files with 141 additions and 153 deletions.
20 changes: 11 additions & 9 deletions Source/BoundaryConditions/ERF_BoundaryConditions_cons.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::impose_lateral_cons_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
int icomp, int ncomp, int ngz)
int icomp, int ncomp, IntVect ng)
{
BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs);
const auto& dom_lo = lbound(domain);
Expand Down Expand Up @@ -70,6 +70,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
{
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);

ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
Expand Down Expand Up @@ -108,6 +109,7 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
{
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);

ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
Expand Down Expand Up @@ -147,10 +149,10 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
// Populate ghost cells on lo-x and hi-x domain boundaries
Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1);
if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,ngz);
if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ngz);
if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ngz);
if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ngz);
if (bx_xlo.smallEnd(2) != domain.smallEnd(2)) bx_xlo.growLo(2,ng[2]);
if (bx_xlo.bigEnd(2) != domain.bigEnd(2)) bx_xlo.growHi(2,ng[2]);
if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ng[2]);
if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ng[2]);
ParallelFor(
bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
Expand Down Expand Up @@ -196,10 +198,10 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4<Real>& dest_arr,
// Populate ghost cells on lo-y and hi-y domain boundaries
Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,ngz);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ngz);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ngz);
if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ngz);
if (bx_ylo.smallEnd(2) != domain.smallEnd(2)) bx_ylo.growLo(2,ng[2]);
if (bx_ylo.bigEnd(2) != domain.bigEnd(2)) bx_ylo.growHi(2,ng[2]);
if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ng[2]);
if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ng[2]);
ParallelFor(
bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
{
Expand Down
78 changes: 9 additions & 69 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,24 +84,18 @@ ERF::FillPatch (int lev, Real time,
FillPatchSingleLevel(*mfs_vel[Vars::cons], ngvect_cons, time, fmf, IntVect(0,0,0), ftime,
0, 0, ncomp, geom[lev]);

(*physbcs_cons[lev])(*mfs_vel[Vars::cons],0,ncomp,ngvect_cons,time,BCVars::cons_bc);

if (!cons_only) {
fmf = {&vars_old[lev][Vars::xvel], &vars_new[lev][Vars::xvel]};
FillPatchSingleLevel(*mfs_vel[Vars::xvel], ngvect_vels, time, fmf,
IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]);
(*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc);

fmf = {&vars_old[lev][Vars::yvel], &vars_new[lev][Vars::yvel]};
FillPatchSingleLevel(*mfs_vel[Vars::yvel], ngvect_vels, time, fmf,
IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]);
(*physbcs_v[lev])(*mfs_vel[Vars::yvel],0,1,ngvect_vels,time,BCVars::xvel_bc);

fmf = {&vars_old[lev][Vars::zvel], &vars_new[lev][Vars::zvel]};
FillPatchSingleLevel(*mfs_vel[Vars::zvel], ngvect_vels, time, fmf,
IntVect(0,0,0), ftime, 0, 0, 1, geom[lev]);
(*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel],
ngvect_vels,time,BCVars::zvel_bc);
} // !cons_only

} else {
Expand All @@ -115,23 +109,17 @@ ERF::FillPatch (int lev, Real time,
mapper = &cell_cons_interp;

// Impose physical bc's on coarse data (note time and 0 are not used)
// Note that we call FillBoundary inside the physbcs call
(*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc);
(*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc);

// Make sure internal ghost cells are filled as well
vars_old[lev-1][Vars::cons].FillBoundary(geom[lev-1].periodicity());
vars_new[lev-1][Vars::cons].FillBoundary(geom[lev-1].periodicity());

// Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled
FillPatchTwoLevels(mf_c, ngvect_cons, IntVect(0,0,0),
time, cmf, ctime, fmf, ftime,
0, 0, mf_c.nComp(), geom[lev-1], geom[lev],
refRatio(lev-1), mapper, domain_bcs_type,
BCVars::cons_bc);

// Impose physical bc's on fine data
(*physbcs_cons[lev])(mf_c,0,mf_c.nComp(),ngvect_cons,time,BCVars::cons_bc);

if (!cons_only)
{
mapper = &face_cons_linear_interp;
Expand All @@ -143,13 +131,10 @@ ERF::FillPatch (int lev, Real time,
// **********************************************************************

// Impose physical bc's on coarse data (note time and 0 are not used)
// Note that we call FillBoundary inside the physbcs call
(*physbcs_u[lev-1])(vars_old[lev-1][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc);
(*physbcs_u[lev-1])(vars_new[lev-1][Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc);

// Make sure internal ghost cells are filled as well
vars_old[lev-1][Vars::xvel].FillBoundary(geom[lev-1].periodicity());
vars_new[lev-1][Vars::xvel].FillBoundary(geom[lev-1].periodicity());

fmf = {&vars_old[lev ][Vars::xvel], &vars_new[lev ][Vars::xvel]};
cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]};

Expand All @@ -160,19 +145,13 @@ ERF::FillPatch (int lev, Real time,
refRatio(lev-1), mapper, domain_bcs_type,
BCVars::xvel_bc);

// Impose physical bc's on fine data
(*physbcs_u[lev])(vars_new[lev][Vars::xvel],0,mf_u.nComp(),ngvect_vels,time,BCVars::xvel_bc);

// **********************************************************************

// Impose physical bc's on coarse data (note time and 0 are not used)
// Note that we call FillBoundary inside the physbcs call
(*physbcs_v[lev-1])(vars_old[lev-1][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc);
(*physbcs_v[lev-1])(vars_new[lev-1][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc);

// Make sure internal ghost cells are filled as well
vars_old[lev-1][Vars::yvel].FillBoundary(geom[lev-1].periodicity());
vars_new[lev-1][Vars::yvel].FillBoundary(geom[lev-1].periodicity());

fmf = {&vars_old[lev ][Vars::yvel], &vars_new[lev ][Vars::yvel]};
cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]};

Expand All @@ -183,12 +162,10 @@ ERF::FillPatch (int lev, Real time,
refRatio(lev-1), mapper, domain_bcs_type,
BCVars::yvel_bc);

// Impose physical bc's on fine data
(*physbcs_v[lev])(vars_new[lev][Vars::yvel],0,1,ngvect_vels,time,BCVars::yvel_bc);

// **********************************************************************

// Impose physical bc's on coarse data (note time and 0 are not used)
// Note that we call FillBoundary inside the physbcs call
(*physbcs_w[lev-1])(vars_old[lev-1][Vars::zvel],
vars_old[lev-1][Vars::xvel],
vars_old[lev-1][Vars::yvel],
Expand All @@ -198,10 +175,6 @@ ERF::FillPatch (int lev, Real time,
vars_new[lev-1][Vars::yvel],
ngvect_vels,time,BCVars::zvel_bc);

// Make sure internal ghost cells are filled as well
vars_old[lev-1][Vars::zvel].FillBoundary(geom[lev-1].periodicity());
vars_new[lev-1][Vars::zvel].FillBoundary(geom[lev-1].periodicity());

fmf = {&vars_old[lev ][Vars::zvel], &vars_new[lev ][Vars::zvel]};
cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]};

Expand All @@ -211,12 +184,6 @@ ERF::FillPatch (int lev, Real time,
0, 0, 1, geom[lev-1], geom[lev],
refRatio(lev-1), mapper, domain_bcs_type,
BCVars::zvel_bc);


// Impose physical bc's on fine data -- note the u and v have been filled above
(*physbcs_w[lev])(*mfs_vel[Vars::zvel],*mfs_vel[Vars::xvel],*mfs_vel[Vars::yvel],
ngvect_vels,time,BCVars::zvel_bc);

} // !cons_only
} // lev > 0

Expand All @@ -236,6 +203,7 @@ ERF::FillPatch (int lev, Real time,
if (m_r2d) fill_from_bndryregs(mfs_vel,time);

// We call these even if init_type == real because these will fill the vertical bcs
// Note that we call FillBoundary inside the physbcs call
(*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,ngvect_cons,time,BCVars::cons_bc);
if (!cons_only) {
(*physbcs_u[lev])(*mfs_vel[Vars::xvel],0,1,ngvect_vels,time,BCVars::xvel_bc);
Expand All @@ -245,37 +213,6 @@ ERF::FillPatch (int lev, Real time,
}
}

/*
* Fill ghost cells of qmoist
*
* @param[in] lev level of refinement at which to fill the data
* @param[in] time time at which the data should be filled
* @param[out] mf MultiFab to be filled (qmoist[lev])
*/
void
ERF::FillPatchMoistVars (int lev, MultiFab& mf)
{
BL_PROFILE_VAR("ERF::FillPatchMoistVars()",ERF_FillPatchMoistVars);
// ***************************************************************************
// Physical bc's at domain boundary
// ***************************************************************************
int icomp_cons = 0;
int ncomp_cons = 1; // We only fill qv, the first component

// Note that we are filling qv, stored in qmoist[lev], with the input data (if there is any), stored
// in RhoQ1_comp.

if (!use_real_bcs) {
Real time = Real(0.0);
IntVect ngvect_cons = mf.nGrowVect();
int bccomp_cons = BCVars::RhoQ1_bc_comp;

(*physbcs_cons[lev])(mf,icomp_cons,ncomp_cons,ngvect_cons,time,bccomp_cons);
}

mf.FillBoundary(geom[lev].periodicity());
}

/*
* Fill valid and ghost data
* This version fills mfs in valid regions with the values in "mfs" when it is passed in;
Expand Down Expand Up @@ -549,7 +486,10 @@ ERF::FillIntermediatePatch (int lev, Real time,
domain_bcs_type);
}

mfs_mom[Vars::cons]->FillBoundary(geom[lev].periodicity());
mfs_mom[IntVars::cons]->FillBoundary(geom[lev].periodicity());
mfs_mom[IntVars::xmom]->FillBoundary(geom[lev].periodicity());
mfs_mom[IntVars::ymom]->FillBoundary(geom[lev].periodicity());
mfs_mom[IntVars::zmom]->FillBoundary(geom[lev].periodicity());
}

/*
Expand Down
2 changes: 1 addition & 1 deletion Source/BoundaryConditions/ERF_PhysBCFunct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ public:

void impose_lateral_cons_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain,
int icomp, int ncomp, int ngz);
int icomp, int ncomp, amrex::IntVect ng);
void impose_vertical_cons_bcs (const amrex::Array4<amrex::Real>& dest_arr,
const amrex::Box& bx, const amrex::Box& domain,
const amrex::Array4<amrex::Real const>& z_nd,
Expand Down
73 changes: 57 additions & 16 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,12 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
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.
//
mf.FillBoundary(m_geom.periodicity());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand Down Expand Up @@ -77,9 +83,11 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,

if (!m_use_real_bcs)
{
impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,nghost[2]);
// We send a box with ghost cells in the lateral directions only
impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,nghost);
}

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

Expand Down Expand Up @@ -118,6 +126,12 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
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.
//
mf.FillBoundary(m_geom.periodicity());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -133,7 +147,10 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
//
// These are the boxes we use to test on relative to the domain
//
Box xbx1 = surroundingNodes(bx,0); xbx1.grow(IntVect(nghost[0],nghost[1],0));
Box xbx1 = surroundingNodes(bx,0); xbx1.grow(nghost);
if(xbx1.smallEnd(2) < domain.smallEnd(2)) xbx1.setSmall(2,domain.smallEnd(2));
if(xbx1.bigEnd(2) > domain.bigEnd(2)) xbx1.setBig(2,domain.bigEnd(2));

Box xbx2 = surroundingNodes(bx,0); xbx2.grow(nghost);

Array4<const Real> z_nd_arr;
Expand All @@ -157,7 +174,6 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,

impose_vertical_xvel_bcs(velx_arr,xbx2,domain,z_nd_arr,dxInv,bccomp,time);
}

} // MFIter
} // OpenMP
} // operator()
Expand Down Expand Up @@ -193,6 +209,12 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
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.
//
mf.FillBoundary(m_geom.periodicity());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -208,7 +230,10 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
//
// These are the boxes we use to test on relative to the domain
//
Box ybx1 = surroundingNodes(bx,1); ybx1.grow(IntVect(nghost[0],nghost[1],0));
Box ybx1 = surroundingNodes(bx,1); ybx1.grow(nghost);
if (ybx1.smallEnd(2) < domain.smallEnd(2)) ybx1.setSmall(2,domain.smallEnd(2));
if (ybx1.bigEnd(2) > domain.bigEnd(2)) ybx1.setBig(2,domain.bigEnd(2));

Box ybx2 = surroundingNodes(bx,1); ybx2.grow(nghost);

Array4<const Real> z_nd_arr;
Expand Down Expand Up @@ -272,6 +297,12 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
}
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.
//
mf.FillBoundary(m_geom.periodicity());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -287,29 +318,39 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
//
// These are the boxes we use to test on relative to the domain
//
Box zbx = surroundingNodes(bx,2); zbx.grow(0,nghost[0]);
zbx.grow(1,nghost[1]);
Box zbx1 = surroundingNodes(bx,2); zbx1.grow(nghost);
if (zbx1.smallEnd(2) < domain.smallEnd(2)) zbx1.setSmall(2,domain.smallEnd(2));
if (zbx1.bigEnd(2) > domain.bigEnd(2)) zbx1.setBig(2,domain.bigEnd(2)+1);

Box zbx2 = surroundingNodes(bx,2); zbx2.grow(nghost);

//
// These are the boxes we use to test on relative to the domain
//
// Box zbx = surroundingNodes(bx,2); zbx.grow(0,nghost[0]);
// zbx.grow(1,nghost[1]);
Array4<const Real> z_nd_arr;

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

Array4<const Real> const& velx_arr = xvel.const_array(mfi);
Array4<const Real> const& vely_arr = yvel.const_array(mfi);
Array4< Real> const& velz_arr = mf.array(mfi);

if (!m_use_real_bcs)
if (!gdomainz.contains(zbx2))
{
if (!gdomainz.contains(zbx))
Array4<const Real> const& velx_arr = xvel.const_array(mfi);
Array4<const Real> const& vely_arr = yvel.const_array(mfi);
Array4< Real> const& velz_arr = mf.array(mfi);

if (!m_use_real_bcs)
{
impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv,bccomp_w);
if (!gdomainz.contains(zbx1))
{
impose_lateral_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx1,domain,z_nd_arr,dxInv,bccomp_w);
}
}
} // m_use_real_bcs

if (!gdomainz.contains(zbx)) {
impose_vertical_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx,domain,z_nd_arr,dxInv,
impose_vertical_zvel_bcs(velz_arr,velx_arr,vely_arr,zbx2,domain,z_nd_arr,dxInv,
bccomp_u, bccomp_v, bccomp_w, m_terrain_type);
}
} // MFIter
Expand Down
Loading

0 comments on commit 88d2efd

Please sign in to comment.