Skip to content

Commit

Permalink
progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Jul 15, 2024
1 parent d14417a commit 6033fc7
Show file tree
Hide file tree
Showing 8 changed files with 86 additions and 43 deletions.
17 changes: 12 additions & 5 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
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,0,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());

Expand Down Expand Up @@ -118,7 +119,9 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
}
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,0,0);
// z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,IntVect(nghost[0],nghost[1],0),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());

Expand Down Expand Up @@ -195,7 +198,9 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
}
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,0,0);
// z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,0,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());

Expand Down Expand Up @@ -273,8 +278,10 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
}
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,IntVect(nghost[0],nghost[1],0),
IntVect(nghost[0],nghost[1],0),m_geom.periodicity());
//z_nd_mf_loc.ParallelCopy(*m_z_phys_nd,0,0,1,IntVect(nghost[0],nghost[1],0),
// IntVect(nghost[0],nghost[1],0),m_geom.periodicity());
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());

Expand Down
1 change: 0 additions & 1 deletion Source/Diffusion/PBLModels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,6 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
const auto& ws10av_arr = most->get_mac_avg(level,4)->const_array(mfi);
const auto& t10av_arr = most->get_mac_avg(level,2)->const_array(mfi);
const auto& t_surf_arr = most->get_t_surf(level)->const_array(mfi);
const bool use_terrain = (z_phys_nd != nullptr);
const Array4<Real const> z_nd_arr = use_terrain ? z_phys_nd->array(mfi) : Array4<Real>{};
const Real most_zref = most->get_zref();

Expand Down
51 changes: 24 additions & 27 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,15 +211,15 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
// ********************************************************************************************
// Map factors
// ********************************************************************************************
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
BoxList bl2d_mf = ba.boxList();
for (auto& b : bl2d_mf) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));
BoxArray ba2d_mf(std::move(bl2d_mf));

mapfac_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,3);
mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(1,0,0)),dm,1,3);
mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d,IntVect(0,1,0)),dm,1,3);
mapfac_m[lev] = std::make_unique<MultiFab>(ba2d_mf,dm,1,3);
mapfac_u[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(1,0,0)),dm,1,3);
mapfac_v[lev] = std::make_unique<MultiFab>(convert(ba2d_mf,IntVect(0,1,0)),dm,1,3);
if (solverChoice.test_mapfactor) {
mapfac_m[lev]->setVal(0.5);
mapfac_u[lev]->setVal(0.5);
Expand Down Expand Up @@ -323,12 +323,12 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
{
lmask_lev[lev].resize(1);
auto ngv = lev_new[Vars::cons].nGrowVect(); ngv[2] = 0;
BoxList bl2d = ba.boxList();
for (auto& b : bl2d) {
BoxList bl2d_mask = ba.boxList();
for (auto& b : bl2d_mask) {
b.setRange(2,0);
}
BoxArray ba2d(std::move(bl2d));
lmask_lev[lev][0] = std::make_unique<iMultiFab>(ba2d,dm,1,ngv);
BoxArray ba2d_mask(std::move(bl2d_mask));
lmask_lev[lev][0] = std::make_unique<iMultiFab>(ba2d_mask,dm,1,ngv);
lmask_lev[lev][0]->setVal(1);
lmask_lev[lev][0]->FillBoundary(geom[lev].periodicity());
}
Expand Down Expand Up @@ -408,11 +408,9 @@ ERF::update_terrain_arrays (int lev, Real time)
{
if (solverChoice.use_terrain) {

if (lev == 0 && (init_type != "real" && init_type != "metgrid") ) {
prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time);
init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag);
}

//
// First interpolate from coarser level if there is one
//
if (lev > 0) {
Vector<MultiFab*> fmf = {z_phys_nd[lev].get(), z_phys_nd[lev].get()};
Vector<Real> ftime = {t_old[lev], t_new[lev]};
Expand All @@ -429,24 +427,23 @@ ERF::update_terrain_arrays (int lev, Real time)
geom[lev-1], geom[lev],
null_bc, 0, null_bc, 0, refRatio(lev-1),
mapper, domain_bcs_type, 0);
}

//
// Then if not using real/metgrid data, we
// 1) redefine the terrain at k=0 for every fine box which includes k=0
// 2) recreate z_phys_nd at every fine node using
// the data at the bottom of each fine grid
// which has been either been interpolated from the coarse grid (k>0)
// or set in init_custom_terrain (k=0)
//
//
// Then, if not using real/metgrid data,
// 1) redefine the terrain at k=0 for every fine box which includes k=0
// 2) recreate z_phys_nd at every fine node using
// the data at the bottom of each fine grid
// which has been either been interpolated from the coarse grid (k>0)
// or set in init_custom_terrain (k=0)
//
if (init_type != "real" && init_type != "metgrid") {
prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time);
if (init_type != "real" && init_type != "metgrid") {
init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag);
}
init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag,phys_bc_type);
}

make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);

make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
}
}
Expand Down
2 changes: 1 addition & 1 deletion Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ ERF::init_from_metgrid (int lev)
for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
// This defines only the z(i,j,0) values given the FAB filled from the NetCDF input
FArrayBox& z_phys_nd_fab = (*z_phys)[mfi];
init_terrain_from_metgrid(z_phys_nd_fab, NC_hgt_fab);
init_terrain_from_metgrid(z_phys_nd_fab, NC_hgt_fab,phys_bc_type);
} // mf

// This defines all the z(i,j,k) values given z(i,j,0) from above.
Expand Down
4 changes: 2 additions & 2 deletions Source/TimeIntegration/TI_fast_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,13 +47,13 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
// Make "old" fast geom -- store in z_phys_nd for convenience
if (verbose) Print() << "Making geometry at start of substep time: " << old_substep_time << std::endl;
prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_substep_time);
init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag);
init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag, phys_bc_type);
make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]);

// Make "new" fast geom
if (verbose) Print() << "Making geometry for end of substep time :" << new_substep_time << std::endl;
prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_substep_time);
init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag);
init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag, phys_bc_type);
make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]);
make_areas (fine_geom,*z_phys_nd_new[level], *ax_new[level], *ay_new[level], *az_new[level]);

Expand Down
6 changes: 3 additions & 3 deletions Source/TimeIntegration/TI_slow_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,19 +59,19 @@

if (verbose) Print() << "Re-making old geometry at old time : " << old_step_time << std::endl;
prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_step_time);
init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag);
init_terrain_grid (level,fine_geom,*z_phys_nd[level], zlevels_stag, phys_bc_type);
make_J (fine_geom,*z_phys_nd[level], *detJ_cc[level]);
make_areas (fine_geom,*z_phys_nd[level], *ax[level], *ay[level], *az[level]);

if (verbose) Print() << "Making src geometry at old_stage_time: " << old_stage_time << std::endl;
prob->init_custom_terrain(fine_geom,*z_phys_nd_src[level],old_stage_time);
init_terrain_grid (level,fine_geom,*z_phys_nd_src[level], zlevels_stag);
init_terrain_grid (level,fine_geom,*z_phys_nd_src[level], zlevels_stag, phys_bc_type);
make_J (fine_geom,*z_phys_nd_src[level], *detJ_cc_src[level]);
make_areas (fine_geom,*z_phys_nd_src[level], *ax_src[level], *ay_src[level], *az_src[level]);

if (verbose) Print() << "Making new geometry at new_stage_time: " << new_stage_time << std::endl;
prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_stage_time);
init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag);
init_terrain_grid (level,fine_geom,*z_phys_nd_new[level], zlevels_stag, phys_bc_type);
make_J (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]);
make_areas (fine_geom,*z_phys_nd_new[level], *ax_new[level], *ay_new[level], *az_new[level]);

Expand Down
3 changes: 2 additions & 1 deletion Source/Utils/TerrainMetrics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,8 @@ void init_zlevels (amrex::Vector<amrex::Real>& zlevels_stag,

void init_terrain_grid (int lev, const amrex::Geometry& geom,
amrex::MultiFab& z_phys_nd,
amrex::Vector<amrex::Real> const& z_levels_h);
amrex::Vector<amrex::Real> const& z_levels_h,
amrex::GpuArray<ERF_BC, AMREX_SPACEDIM*2>& phys_bc_type);

//*****************************************************************************************
// Compute terrain metric terms at cell-center
Expand Down
45 changes: 42 additions & 3 deletions Source/Utils/TerrainMetrics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,8 @@ init_zlevels (Vector<Real>& zlevels_stag,

void
init_terrain_grid (int lev, const Geometry& geom, MultiFab& z_phys_nd,
Vector<Real> const& z_levels_h)
Vector<Real> const& z_levels_h,
GpuArray<ERF_BC, AMREX_SPACEDIM*2>& phys_bc_type)
{
// z_nd is nodal in all directions
const Box& domain = geom.Domain();
Expand Down Expand Up @@ -397,10 +398,48 @@ init_terrain_grid (int lev, const Geometry& geom, MultiFab& z_phys_nd,
}
} //switch

// Fill ghost layers and corners (including periodic)
z_phys_nd.FillBoundary(geom.periodicity());
// Fill ghost layers and corners (including periodic)
z_phys_nd.FillBoundary(geom.periodicity());


if (phys_bc_type[Orientation(0,Orientation::low )] == ERF_BC::symmetry ||
phys_bc_type[Orientation(0,Orientation::high)] == ERF_BC::symmetry ||
phys_bc_type[Orientation(1,Orientation::low )] == ERF_BC::symmetry ||
phys_bc_type[Orientation(1,Orientation::high)] == ERF_BC::symmetry) {

const auto& dom_lo = lbound(convert(geom.Domain(),IntVect(1,1,1)));
const auto& dom_hi = ubound(convert(geom.Domain(),IntVect(1,1,1)));

for (MFIter mfi(z_phys_nd,true); mfi.isValid(); ++mfi) {
const Box& bx = mfi.growntilebox();
const Array4< Real> z_nd_arr = z_phys_nd.array(mfi);
if (phys_bc_type[Orientation(0,Orientation::low)] == ERF_BC::symmetry && bx.smallEnd(0) == dom_lo.x) {
ParallelFor(makeSlab(bx,0,1), [=] AMREX_GPU_DEVICE (int , int j, int k)
{
z_nd_arr(dom_lo.x-1,j,k) = z_nd_arr(dom_lo.x+1,j,k);
});
}
if (phys_bc_type[Orientation(0,Orientation::high)] == ERF_BC::symmetry && bx.bigEnd(0) == dom_hi.x) {
ParallelFor(makeSlab(bx,0,1), [=] AMREX_GPU_DEVICE (int , int j, int k)
{
z_nd_arr(dom_hi.x+1,j,k) = z_nd_arr(dom_hi.x-1,j,k);
});
}
if (phys_bc_type[Orientation(1,Orientation::low)] == ERF_BC::symmetry && bx.smallEnd(1) == dom_lo.y) {
ParallelFor(makeSlab(bx,1,1), [=] AMREX_GPU_DEVICE (int i, int , int k)
{
z_nd_arr(i,dom_lo.y-1,k) = z_nd_arr(i,dom_lo.y+1,k);
});
}
if (phys_bc_type[Orientation(1,Orientation::high)] == ERF_BC::symmetry && bx.bigEnd(1) == dom_hi.y) {
ParallelFor(makeSlab(bx,1,1), [=] AMREX_GPU_DEVICE (int i, int , int k)
{
z_nd_arr(i,dom_hi.y+1,k) = z_nd_arr(i,dom_hi.y-1,k);
});
}
}
}

//********************************************************************************
// Populate domain boundary cells in z-direction
//********************************************************************************
Expand Down

0 comments on commit 6033fc7

Please sign in to comment.