Skip to content

Commit

Permalink
fix bounds of terrain_mf
Browse files Browse the repository at this point in the history
  • Loading branch information
asalmgren committed Jan 9, 2025
1 parent 017a80b commit e0d49c7
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 14 deletions.
2 changes: 1 addition & 1 deletion Source/ERF_MakeNewArrays.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -492,7 +492,7 @@ ERF::init_zphys (int lev, Real time)
Box bx(surroundingNodes(geom[lev].Domain())); bx.grow(2);
BoxArray ba(makeSlab(bx,2,0));
DistributionMapping dm(ba);
MultiFab terrain_mf(ba,dm,1,0);
MultiFab terrain_mf(ba,dm,1,1);
terrain_mf.setVal(-1.e23);

//
Expand Down
25 changes: 12 additions & 13 deletions Source/ERF_ProbCommon.H
Original file line number Diff line number Diff line change
Expand Up @@ -296,7 +296,7 @@ public:
read_custom_terrain(const std::string& fname,
const bool is_usgs,
const amrex::Geometry& geom,
amrex::MultiFab& z_phys_nd,
amrex::MultiFab& terrain_mf,
const amrex::Real& /*time*/)
{
// Read terrain file
Expand Down Expand Up @@ -341,9 +341,9 @@ public:
m_zterrain.push_back(value3);
counter += 1;
}
AMREX_ASSERT(m_xterrain.size() == static_cast<long unsigned int>(nx*ny));
AMREX_ASSERT(m_yterrain.size() == static_cast<long unsigned int>(ny));
AMREX_ASSERT(m_zterrain.size() == static_cast<long unsigned int>(nx*ny));
AMREX_ASSERT(m_xterrain.size() == static_cast<long int>(nx*ny));
AMREX_ASSERT(m_yterrain.size() == static_cast<long int>(ny));
AMREX_ASSERT(m_zterrain.size() == static_cast<long int>(nx*ny));

} else {
file >> nx; file >> ny;
Expand Down Expand Up @@ -382,8 +382,7 @@ public:
amrex::Real* d_yt = d_yterrain.data();
amrex::Real* d_zt = d_zterrain.data();

// Populate z_phys data
int ngrow = z_phys_nd.nGrow();
int ngrow = terrain_mf.nGrow();
auto dx = geom.CellSizeArray();
auto ProbLoArr = geom.ProbLoArray();

Expand All @@ -393,13 +392,13 @@ public:
int ihi = geom.Domain().bigEnd(0) + 1;
int jhi = geom.Domain().bigEnd(1) + 1;

for (amrex::MFIter mfi(z_phys_nd,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
for (amrex::MFIter mfi(terrain_mf,amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
// Grown box with no z range
amrex::Box xybx = mfi.growntilebox(ngrow);
xybx.setRange(2,0);

amrex::Array4<amrex::Real> const& z_arr = z_phys_nd.array(mfi);
amrex::Array4<amrex::Real> const& z_arr = terrain_mf.array(mfi);

amrex::ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/)
{
Expand Down Expand Up @@ -484,32 +483,32 @@ public:
* Note: Terrain functionality can also be used to provide grid stretching.
*
* @param[in] geom container for geometric information
* @param[out] z_phys_nd height coordinate at nodes
* @param[out] terrain_mf height coordinate at nodes
* @param[in] time current time
*/
virtual void
init_custom_terrain (const amrex::Geometry& /*geom*/,
amrex::MultiFab& z_phys_nd,
amrex::MultiFab& terrain_mf,
const amrex::Real& /*time*/)
{
// Note that this only sets the terrain value at the ground IF k=0 is in the box
amrex::Print() << "Initializing flat terrain at z=0" << std::endl;

// Number of ghost cells
int ngrow = z_phys_nd.nGrow();
int ngrow = terrain_mf.nGrow();

// Bottom of domain
int k0 = 0;

for ( amrex::MFIter mfi(z_phys_nd, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
for ( amrex::MFIter mfi(terrain_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
amrex::Box zbx = mfi.nodaltilebox(2);
if (zbx.smallEnd(2) <= 0) {
// Grown box with no z range
amrex::Box xybx = mfi.growntilebox(ngrow);
xybx.setRange(2,0);

amrex::Array4<amrex::Real> const& z_arr = z_phys_nd.array(mfi);
amrex::Array4<amrex::Real> const& z_arr = terrain_mf.array(mfi);

ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) {
z_arr(i,j,k0) = 0.0;
Expand Down

0 comments on commit e0d49c7

Please sign in to comment.