Skip to content

Commit

Permalink
fix the bc fill for when we have fine patches and terrain (#1684)
Browse files Browse the repository at this point in the history 
* fix the bc fill for when we have fine patches and terrain

* more fixes for multilevel

* enlarge diffusion arrays in vertical -- needed for multilevel

* fix spelling

* fix whitespace

* progress
  • Loading branch information
@asalmgren
asalmgren authored Jul 15, 2024
1 parent 0277aff commit 0ea0a01
Show file tree
Hide file tree
Showing 12 changed files with 195 additions and 61 deletions.
11 changes: 9 additions & 2 deletions Source/BoundaryConditions/BoundaryConditions_cons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,13 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);

Box per_grown_domain(domain);
int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
per_grown_domain.grow(IntVect(growx,growy,0));
const auto& perdom_lo = lbound(per_grown_domain);
const auto& perdom_hi = ubound(per_grown_domain);

GeometryData const& geomdata = m_geom.data();

// xlo: ori = 0
Expand Down Expand Up @@ -317,8 +324,8 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& dest_arr
ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
// Clip indices for ghost-cells
int ii = amrex::min(amrex::max(i,dom_lo.x),dom_hi.x);
int jj = amrex::min(amrex::max(j,dom_lo.y),dom_hi.y);
int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);

// Get metrics
Real met_h_xi = Compute_h_xi_AtCellCenter (ii,jj,k0,dxInv,z_phys_nd);
Expand Down
11 changes: 9 additions & 2 deletions Source/BoundaryConditions/BoundaryConditions_xvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,13 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4<Real>& dest_arr,
const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);

Box per_grown_domain(domain);
int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
per_grown_domain.grow(IntVect(growx,growy,0));
const auto& perdom_lo = lbound(per_grown_domain);
const auto& perdom_hi = ubound(per_grown_domain);

GeometryData const& geomdata = m_geom.data();

// Based on BCRec for the domain, we need to make BCRec for this Box
Expand Down Expand Up @@ -262,8 +269,8 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4<Real>& dest_arr,
ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
// Clip indices for ghost-cells
int ii = amrex::min(amrex::max(i,dom_lo.x),dom_hi.x);
int jj = amrex::min(amrex::max(j,dom_lo.y),dom_hi.y);
int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);

// Get metrics
Real met_h_xi = Compute_h_xi_AtIface (ii, jj, k0, dxInv, z_phys_nd);
Expand Down
11 changes: 9 additions & 2 deletions Source/BoundaryConditions/BoundaryConditions_yvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,13 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& dest_arr,
const auto& dom_lo = lbound(domain);
const auto& dom_hi = ubound(domain);

Box per_grown_domain(domain);
int growx = (m_geom.isPeriodic(0)) ? 1 : 0;
int growy = (m_geom.isPeriodic(1)) ? 1 : 0;
per_grown_domain.grow(IntVect(growx,growy,0));
const auto& perdom_lo = lbound(per_grown_domain);
const auto& perdom_hi = ubound(per_grown_domain);

// Based on BCRec for the domain, we need to make BCRec for this Box
// bccomp is used as starting index for m_domain_bcs_type
// 0 is used as starting index for bcrs
Expand Down Expand Up @@ -260,8 +267,8 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& dest_arr,
ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
// Clip indices for ghost-cells
int ii = amrex::min(amrex::max(i,dom_lo.x),dom_hi.x);
int jj = amrex::min(amrex::max(j,dom_lo.y),dom_hi.y);
int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x);
int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y);

// Get metrics
Real met_h_xi = Compute_h_xi_AtJface (ii, jj, k0, dxInv, z_phys_nd);
Expand Down
79 changes: 75 additions & 4 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,24 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,
}
}

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;
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());
}
z_nd_mf_loc.FillBoundary(m_geom.periodicity());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -53,7 +71,7 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp,

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

if (!gdomain.contains(cbx2))
Expand Down Expand Up @@ -89,6 +107,24 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,
}
}

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;
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,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());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -111,7 +147,7 @@ void ERFPhysBCFunct_u::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,

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

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

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;
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,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());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -172,7 +226,7 @@ void ERFPhysBCFunct_v::operator() (MultiFab& mf, int /*icomp*/, int /*ncomp*/,

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

if (!gdomainy.contains(ybx2))
Expand Down Expand Up @@ -214,6 +268,23 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,
// We want to make sure we impose the z-vels at k=0 if the box includes k=0
if (gdomainz.smallEnd(2) == 0) gdomainz.setSmall(2,1);

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,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());

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand All @@ -235,7 +306,7 @@ void ERFPhysBCFunct_w::operator() (MultiFab& mf, MultiFab& xvel, MultiFab& yvel,

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

Array4<const Real> const& velx_arr = xvel.const_array(mfi);;
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
79 changes: 40 additions & 39 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 All @@ -353,25 +353,29 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
BoxArray ba23 = convert(ba, IntVect(0,1,1));

if (l_use_diff) {
Tau11_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
Tau22_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
Tau33_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
Tau12_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,0) );
Tau13_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,0) );
Tau23_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,0) );
//
// NOTE: We require ghost cells in the vertical when allowing grids that don't
// cover the entire vertical extent of the domain at this level
//
Tau11_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
Tau22_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
Tau33_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
Tau12_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
Tau13_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
Tau23_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
if (l_use_terrain) {
Tau21_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,0) );
Tau31_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,0) );
Tau32_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,0) );
Tau21_lev[lev] = std::make_unique<MultiFab>( ba12, dm, 1, IntVect(1,1,1) );
Tau31_lev[lev] = std::make_unique<MultiFab>( ba13, dm, 1, IntVect(1,1,1) );
Tau32_lev[lev] = std::make_unique<MultiFab>( ba23, dm, 1, IntVect(1,1,1) );
} else {
Tau21_lev[lev] = nullptr;
Tau31_lev[lev] = nullptr;
Tau32_lev[lev] = nullptr;
}
SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
SFS_hfx1_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx2_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx3_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_diss_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,0) );
SFS_diss_lev[lev] = std::make_unique<MultiFab>( ba , dm, 1, IntVect(1,1,1) );
SFS_hfx1_lev[lev]->setVal(0.);
SFS_hfx2_lev[lev]->setVal(0.);
SFS_hfx3_lev[lev]->setVal(0.);
Expand Down Expand Up @@ -404,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 @@ -425,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: 3 additions & 1 deletion Source/TimeIntegration/ERF_ComputeTimestep.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,9 @@ ERF::ComputeDt (int step)
dt_0 = amrex::min(dt_0, n_factor*dt_tmp[lev]);
if (step == 0){
dt_0 *= init_shrink;
Print() << "Timestep 0: shrink initial dt by " << init_shrink << std::endl;
if (verbose) {
Print() << "Timestep 0: shrink initial dt at level " << lev << " by " << init_shrink << std::endl;
}
}
}

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
Loading

0 comments on commit 0ea0a01

Please sign in to comment.