diff --git a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp index a49db9391..6abb2221d 100644 --- a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp @@ -12,7 +12,11 @@ using namespace amrex; */ void -ERF::fill_from_wrfbdy (const Vector& mfs, const Real time) +ERF::fill_from_wrfbdy (const Vector& mfs, + const Real time, + bool cons_only, + int icomp_cons, + int ncomp_cons) { int lev = 0; @@ -46,10 +50,13 @@ ERF::fill_from_wrfbdy (const Vector& mfs, const Real time) ind_map.push_back( {0} ); // zvel // Nvars to loop over - Vector comp_var = {NVAR, 1, 1, 1}; + Vector 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]; @@ -61,16 +68,18 @@ ERF::fill_from_wrfbdy (const Vector& 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; @@ -138,7 +147,6 @@ ERF::fill_from_wrfbdy (const Vector& 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 @@ -158,14 +166,14 @@ ERF::fill_from_wrfbdy (const Vector& 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); }); diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 9890fe906..ba9cc3c7e 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -92,7 +92,7 @@ ERF::FillPatch (int lev, Real time, const Vector& 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); @@ -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); diff --git a/Source/ERF.H b/Source/ERF.H index 7eb43bb7f..02beaddb7 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -297,7 +297,10 @@ public: #ifdef ERF_USE_NETCDF void fill_from_wrfbdy (const amrex::Vector& mfs, - amrex::Real time); + amrex::Real time, + bool cons_only = false, + int icomp_cons = 0, + int ncomp_cons = NVAR); #endif #ifdef ERF_USE_PARTICLES diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index 81218520f..0996568d3 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -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")); } @@ -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")); } @@ -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]); } @@ -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