Skip to content

Commit

Permalink
Fix for failing WPS test. (#1261)
Browse files Browse the repository at this point in the history 
* Fix for failing WPS test.

* Streamline indice selection.

* We need to write and read ghost cells for qmoist since we don't directly impose boundary conditions on those vars.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Oct 6, 2023
1 parent f0d8d23 commit d31ad0a
Show file tree
Hide file tree
Showing 4 changed files with 68 additions and 49 deletions.
30 changes: 19 additions & 11 deletions Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,11 @@ using namespace amrex;
*/

void
ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs, const Real time)
ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs,
const Real time,
bool cons_only,
int icomp_cons,
int ncomp_cons)
{
int lev = 0;

Expand Down Expand Up @@ -46,10 +50,13 @@ ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs, const Real time)
ind_map.push_back( {0} ); // zvel

// Nvars to loop over
Vector<int> comp_var = {NVAR, 1, 1, 1};
Vector<int> comp_var = {ncomp_cons, 1, 1, 1};

// End of vars loop
int var_idx_end = (cons_only) ? Vars::cons + 1 : Vars::NumTypes;

// Loop over all variable types
for (int var_idx = Vars::cons; var_idx < Vars::NumTypes; ++var_idx)
for (int var_idx = Vars::cons; var_idx < var_idx_end; ++var_idx)
{
MultiFab& mf = *mfs[var_idx];

Expand All @@ -61,16 +68,18 @@ ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs, const Real time)
const auto& dom_lo = amrex::lbound(domain);
const auto& dom_hi = amrex::ubound(domain);

// Offset only applys to cons (we may fill a subset of these vars)
int offset = (var_idx == Vars::cons) ? icomp_cons : 0;

// Loop over each component
for (int comp_idx(0); comp_idx < comp_var[var_idx]; ++comp_idx)
for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx)
{
int width;
int width = wrfbdy_set_width;;

// Variable can be read from wrf bdy
//------------------------------------
if (is_read[var_idx][comp_idx])
{
width = wrfbdy_set_width;
int ivar = ind_map[var_idx][comp_idx];
IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0;

Expand Down Expand Up @@ -138,7 +147,6 @@ ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs, const Real time)
// Variable not read from wrf bdy
//------------------------------------
} else {
width = wrfbdy_width - 1;
IntVect ng_vect = mf.nGrowVect(); ng_vect[2] = 0;

#ifdef AMREX_USE_OMP
Expand All @@ -158,14 +166,14 @@ ERF::fill_from_wrfbdy (const Vector<MultiFab*>& mfs, const Real time)
ParallelFor(bx_xlo, bx_xhi,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
int jj = std::max(j , dom_lo.y);
jj = std::min(jj, dom_hi.y);
int jj = std::max(j , dom_lo.y+width);
jj = std::min(jj, dom_hi.y-width);
dest_arr(i,j,k,comp_idx) = dest_arr(dom_lo.x+width,jj,k,comp_idx);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
int jj = std::max(j , dom_lo.y);
jj = std::min(jj, dom_hi.y);
int jj = std::max(j , dom_lo.y+width);
jj = std::min(jj, dom_hi.y-width);
dest_arr(i,j,k,comp_idx) = dest_arr(dom_hi.x-width,jj,k,comp_idx);
});

Expand Down
7 changes: 5 additions & 2 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)

#ifdef ERF_USE_NETCDF
// We call this here because it is an ERF routine
if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs,time);
if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs, time);
#endif

if (m_r2d) fill_from_bndryregs(mfs,time);
Expand Down Expand Up @@ -245,8 +245,11 @@ ERF::FillIntermediatePatch (int lev, Real time,
IntVect ngvect_vels = IntVect(ng_vel ,ng_vel ,ng_vel);

#ifdef ERF_USE_NETCDF
// NOTE: This routine needs to be aware of what FillIntermediatePatch is operating on
// --- i.e., cons_only and which cons indices (icomp_cons & ncomp_cons)

// We call this here because it is an ERF routine
if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs,time);
if (init_type=="real" && lev==0) fill_from_wrfbdy(mfs, time, cons_only, icomp_cons, ncomp_cons);
#endif

if (m_r2d) fill_from_bndryregs(mfs,time);
Expand Down
5 changes: 4 additions & 1 deletion Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,10 @@ public:

#ifdef ERF_USE_NETCDF
void fill_from_wrfbdy (const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::Real time);
amrex::Real time,
bool cons_only = false,
int icomp_cons = 0,
int ncomp_cons = NVAR);
#endif

#ifdef ERF_USE_PARTICLES
Expand Down
75 changes: 40 additions & 35 deletions Source/IO/Checkpoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,23 +126,26 @@ ERF::WriteCheckpointFile () const
MultiFab::Copy(zvel,vars_new[lev][Vars::zvel],0,0,1,0);
VisMF::Write(zvel, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "ZFace"));

IntVect ng;
#ifdef ERF_USE_MOISTURE
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),0);
MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),0);
// We must read and write qmoist with ghost cells because we don't directly impose BCs on these vars
ng = qmoist[lev].nGrowVect();
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
MultiFab::Copy(moist_vars,qmoist[lev],0,0,qmoist[lev].nComp(),ng);
VisMF::Write(moist_vars, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MoistVars"));
#endif

// Note that we write the ghost cells of the base state (unlike above)
IntVect ngvect_base = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ngvect_base);
MultiFab::Copy(base,base_state[lev],0,0,base.nComp(),ngvect_base);
ng = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng);
MultiFab::Copy(base,base_state[lev],0,0,base.nComp(),ng);
VisMF::Write(base, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState"));

if (solverChoice.use_terrain) {
// Note that we also write the ghost cells of z_phys_nd
IntVect ngvect = z_phys_nd[lev]->nGrowVect();
MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ngvect);
MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ngvect);
ng = z_phys_nd[lev]->nGrowVect();
MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng);
MultiFab::Copy(z_height,*z_phys_nd[lev],0,0,1,ng);
VisMF::Write(z_height, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z_Phys_nd"));
}

Expand All @@ -153,19 +156,19 @@ ERF::WriteCheckpointFile () const
}
BoxArray ba2d(std::move(bl2d));

IntVect ngvect_mf = mapfac_m[lev]->nGrowVect();
MultiFab mf_m(ba2d,dmap[lev],1,ngvect_mf);
MultiFab::Copy(mf_m,*mapfac_m[lev],0,0,1,ngvect_mf);
ng = mapfac_m[lev]->nGrowVect();
MultiFab mf_m(ba2d,dmap[lev],1,ng);
MultiFab::Copy(mf_m,*mapfac_m[lev],0,0,1,ng);
VisMF::Write(mf_m, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_m"));

ngvect_mf = mapfac_u[lev]->nGrowVect();
MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ngvect_mf);
MultiFab::Copy(mf_u,*mapfac_u[lev],0,0,1,ngvect_mf);
ng = mapfac_u[lev]->nGrowVect();
MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ng);
MultiFab::Copy(mf_u,*mapfac_u[lev],0,0,1,ng);
VisMF::Write(mf_u, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_u"));

ngvect_mf = mapfac_v[lev]->nGrowVect();
MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ngvect_mf);
MultiFab::Copy(mf_v,*mapfac_v[lev],0,0,1,ngvect_mf);
ng = mapfac_v[lev]->nGrowVect();
MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng);
MultiFab::Copy(mf_v,*mapfac_v[lev],0,0,1,ng);
VisMF::Write(mf_v, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "MapFactor_v"));
}

Expand Down Expand Up @@ -327,25 +330,27 @@ ERF::ReadCheckpointFile ()
VisMF::Read(zvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "ZFace"));
MultiFab::Copy(vars_new[lev][Vars::zvel],zvel,0,0,1,0);

IntVect ng;
#ifdef ERF_USE_MOISTURE
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),0);
VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),0);
ng = qmoist[lev].nGrowVect();
MultiFab moist_vars(grids[lev],dmap[lev],qmoist[lev].nComp(),ng);
VisMF::Read(moist_vars, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MoistVars"));
MultiFab::Copy(qmoist[lev],moist_vars,0,0,qmoist[lev].nComp(),ng);
#endif

// Note that we read the ghost cells of the base state (unlike above)
IntVect ngvect_base = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ngvect_base);
ng = base_state[lev].nGrowVect();
MultiFab base(grids[lev],dmap[lev],base_state[lev].nComp(),ng);
VisMF::Read(base, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState"));
MultiFab::Copy(base_state[lev],base,0,0,base.nComp(),ngvect_base);
MultiFab::Copy(base_state[lev],base,0,0,base.nComp(),ng);
base_state[lev].FillBoundary(geom[lev].periodicity());

if (solverChoice.use_terrain) {
// Note that we also read the ghost cells of z_phys_nd
IntVect ngvect = z_phys_nd[lev]->nGrowVect();
MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ngvect);
ng = z_phys_nd[lev]->nGrowVect();
MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng);
VisMF::Read(z_height, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Z_Phys_nd"));
MultiFab::Copy(*z_phys_nd[lev],z_height,0,0,1,ngvect);
MultiFab::Copy(*z_phys_nd[lev],z_height,0,0,1,ng);
update_terrain_arrays(lev, t_new[lev]);
}

Expand All @@ -356,20 +361,20 @@ ERF::ReadCheckpointFile ()
}
BoxArray ba2d(std::move(bl2d));

IntVect ngvect_mf = mapfac_m[lev]->nGrowVect();
MultiFab mf_m(ba2d,dmap[lev],1,ngvect_mf);
ng = mapfac_m[lev]->nGrowVect();
MultiFab mf_m(ba2d,dmap[lev],1,ng);
VisMF::Read(mf_m, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_m"));
MultiFab::Copy(*mapfac_m[lev],mf_m,0,0,1,ngvect_mf);
MultiFab::Copy(*mapfac_m[lev],mf_m,0,0,1,ng);

ngvect_mf = mapfac_u[lev]->nGrowVect();
MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ngvect_mf);
ng = mapfac_u[lev]->nGrowVect();
MultiFab mf_u(convert(ba2d,IntVect(1,0,0)),dmap[lev],1,ng);
VisMF::Read(mf_u, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_u"));
MultiFab::Copy(*mapfac_u[lev],mf_u,0,0,1,ngvect_mf);
MultiFab::Copy(*mapfac_u[lev],mf_u,0,0,1,ng);

ngvect_mf = mapfac_v[lev]->nGrowVect();
MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ngvect_mf);
ng = mapfac_v[lev]->nGrowVect();
MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng);
VisMF::Read(mf_v, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_v"));
MultiFab::Copy(*mapfac_v[lev],mf_v,0,0,1,ngvect_mf);
MultiFab::Copy(*mapfac_v[lev],mf_v,0,0,1,ng);
}

#ifdef ERF_USE_PARTICLES
Expand Down

0 comments on commit d31ad0a

Please sign in to comment.