From d6127ecb3b2259c241313a976781b5c240d046b0 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Thu, 5 Oct 2023 14:58:21 -0700 Subject: [PATCH 1/3] Fix for failing WPS test. --- .../BoundaryConditions_wrfbdy.cpp | 32 ++++++++++++------- Source/BoundaryConditions/ERF_FillPatch.cpp | 7 ++-- Source/ERF.H | 5 ++- 3 files changed, 30 insertions(+), 14 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp index a49db9391..5d7cc9417 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,14 @@ 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 = Vars::NumTypes; + if (cons_only) var_idx_end = Vars::cons + 1; // 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 +69,19 @@ 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 = 0; + if (var_idx == Vars::cons) offset = icomp_cons; + // 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 +149,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 +168,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 From 37a67a8a9ae1ebcb740ec5136530179d813031f0 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 6 Oct 2023 10:10:17 -0700 Subject: [PATCH 2/3] Streamline indice selection. --- Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp index 5d7cc9417..6abb2221d 100644 --- a/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_wrfbdy.cpp @@ -53,8 +53,7 @@ ERF::fill_from_wrfbdy (const Vector& mfs, Vector comp_var = {ncomp_cons, 1, 1, 1}; // End of vars loop - int var_idx_end = Vars::NumTypes; - if (cons_only) var_idx_end = Vars::cons + 1; + int var_idx_end = (cons_only) ? Vars::cons + 1 : Vars::NumTypes; // Loop over all variable types for (int var_idx = Vars::cons; var_idx < var_idx_end; ++var_idx) @@ -70,8 +69,7 @@ ERF::fill_from_wrfbdy (const Vector& mfs, const auto& dom_hi = amrex::ubound(domain); // Offset only applys to cons (we may fill a subset of these vars) - int offset = 0; - if (var_idx == Vars::cons) offset = icomp_cons; + int offset = (var_idx == Vars::cons) ? icomp_cons : 0; // Loop over each component for (int comp_idx(offset); comp_idx < (comp_var[var_idx]+offset); ++comp_idx) From 585c7e0ebebd2dc18e2c0f8d4ff80cefd4e29171 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Fri, 6 Oct 2023 13:02:18 -0700 Subject: [PATCH 3/3] We need to write and read ghost cells for qmoist since we don't directly impose boundary conditions on those vars. --- Source/IO/Checkpoint.cpp | 75 +++++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 35 deletions(-) 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