From f750b5b8714efefa994b331b5afe486573f7f517 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 15 Sep 2024 16:55:40 -0700 Subject: [PATCH] fix bcs with terrain (#1809) * fix bcs with terrain * fix physbc file also --- .../ERF_BoundaryConditions_cons.cpp | 12 ++++++------ .../ERF_BoundaryConditions_xvel.cpp | 16 +++++++++++----- .../ERF_BoundaryConditions_yvel.cpp | 11 +++++++++-- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 19 +++---------------- 4 files changed, 29 insertions(+), 29 deletions(-) diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp index d17b962d4..2f2531e4f 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp @@ -426,11 +426,11 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr // Loop over ghost cells in bottom XY-plane (valid box) Box xybx = bx; - int k0 = 0; - if (xybx.smallEnd(2) < 0) { - + int k0 = dom_lo.z; + if (xybx.smallEnd(2) < 0) + { xybx.setBig(2,dom_lo.z-1); - xybx.setSmall(2,bx.smallEnd()[2]); + xybx.setSmall(2,bx.smallEnd(2)); // Get the dz cell size Real dz = geomdata.CellSize(2); @@ -450,7 +450,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr // GradX at IJK location inside domain -- this relies on the assumption that we have // used foextrap for cell-centered quantities outside the domain to define the gradient as zero Real GradVarx, GradVary; - if (i < dom_lo.x-1 || i > dom_hi.x+1) { + if (i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) { GradVarx = 0.0; } else if (i+1 > bx_hi.x) { GradVarx = dxInv[0] * (dest_arr(i ,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp)); @@ -462,7 +462,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr // GradY at IJK location inside domain -- this relies on the assumption that we have // used foextrap for cell-centered quantities outside the domain to define the gradient as zero - if (j < dom_lo.y-1 || j > dom_hi.y+1) { + if (j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) { GradVary = 0.0; } else if (j+1 > bx_hi.y) { GradVary = dxInv[1] * (dest_arr(i,j ,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp)); diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_xvel.cpp index 1b82f7dc5..38f91e1c0 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_xvel.cpp @@ -248,13 +248,17 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, if (m_z_phys_nd) { const auto& bx_lo = lbound(bx); const auto& bx_hi = ubound(bx); + + const auto& zphys_lo = lbound(Box(z_phys_nd)); + const auto& zphys_hi = ubound(Box(z_phys_nd)); + // Neumann conditions (d/dn = 0) must be aware of the surface normal with terrain. // An additional source term arises from d/dx & d/dy & met_h_xi/eta/zeta. //===================================================================================== // Only modify scalars, U, or V // Loop over each component // Hit for Neumann condition at kmin - if(bcrs[0].lo(2) == ERFBCType::foextrap) { + if (bcrs[0].lo(2) == ERFBCType::foextrap) { // Loop over ghost cells in bottom XY-plane (valid box) Box xybx = bx; xybx.setBig(2,-1); @@ -269,7 +273,9 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, { // Clip indices for ghost-cells int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); + ii = amrex::min(amrex::max(ii,zphys_lo.x),zphys_hi.x); int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y); + jj = amrex::min(amrex::max(jj,zphys_lo.y),zphys_hi.y); // Get metrics Real met_h_xi = Compute_h_xi_AtIface (ii, jj, k0, dxInv, z_phys_nd); @@ -277,9 +283,9 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, Real met_h_zeta = Compute_h_zeta_AtIface(ii, jj, k0, dxInv, z_phys_nd); // GradX at IJK location inside domain -- this relies on the assumption that we have - // used foextrap for cell-centered quantities outside the domain to define the gradient as zero + // used foextrap for quantities outside the domain to define the gradient as zero Real GradVarx, GradVary; - if (i < dom_lo.x-1 || i > dom_hi.x+1) { + if ( i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) { GradVarx = 0.0; } else if (i+1 > bx_hi.x) { GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); @@ -290,8 +296,8 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, } // GradY at IJK location inside domain -- this relies on the assumption that we have - // used foextrap for cell-centered quantities outside the domain to define the gradient as zero - if (j < dom_lo.y-1 || j > dom_hi.y+1) { + // used foextrap for quantities outside the domain to define the gradient as zero + if ( j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) { GradVary = 0.0; } else if (j+1 > bx_hi.y) { GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); diff --git a/Source/BoundaryConditions/ERF_BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/ERF_BoundaryConditions_yvel.cpp index c1252214e..df98a9c4f 100644 --- a/Source/BoundaryConditions/ERF_BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/ERF_BoundaryConditions_yvel.cpp @@ -243,8 +243,13 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, } if (m_z_phys_nd) { + const auto& bx_lo = lbound(bx); const auto& bx_hi = ubound(bx); + + const auto& zphys_lo = lbound(Box(z_phys_nd)); + const auto& zphys_hi = ubound(Box(z_phys_nd)); + // Neumann conditions (d/dn = 0) must be aware of the surface normal with terrain. // An additional source term arises from d/dx & d/dy & met_h_xi/eta/zeta. //===================================================================================== @@ -266,7 +271,9 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, { // Clip indices for ghost-cells int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); + ii = amrex::min(amrex::max(ii,zphys_lo.x),zphys_hi.x); int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y); + jj = amrex::min(amrex::max(jj,zphys_lo.y),zphys_hi.y); // Get metrics Real met_h_xi = Compute_h_xi_AtJface (ii, jj, k0, dxInv, z_phys_nd); @@ -276,7 +283,7 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, // GradX at IJK location inside domain -- this relies on the assumption that we have // used foextrap for cell-centered quantities outside the domain to define the gradient as zero Real GradVarx, GradVary; - if (i < dom_lo.x-1 || i > dom_hi.x+1) { + if ( i < dom_lo.x-1 || i > dom_hi.x+1 || (i+1 > bx_hi.x && i-1 < bx_lo.x) ) { GradVarx = 0.0; } else if (i+1 > bx_hi.x) { GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); @@ -288,7 +295,7 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, // GradY at IJK location inside domain -- this relies on the assumption that we have // used foextrap for cell-centered quantities outside the domain to define the gradient as zero - if (j < dom_lo.y-1 || j > dom_hi.y+1) { + if ( j < dom_lo.y-1 || j > dom_hi.y+1 || (j+1 > bx_hi.y && j-1 < bx_lo.y) ) { GradVary = 0.0; } else if (j+1 > bx_hi.y) { GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index c791b4ae9..5a771e53b 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -31,13 +31,11 @@ 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) { + m_z_phys_nd->FillBoundary(m_geom.periodicity()); 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); } @@ -47,7 +45,6 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, 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()) @@ -107,23 +104,19 @@ 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) { + m_z_phys_nd->FillBoundary(m_geom.periodicity()); 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()) @@ -186,23 +179,19 @@ 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) { + m_z_phys_nd->FillBoundary(m_geom.periodicity()); 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()) @@ -278,8 +267,6 @@ 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,m_z_phys_nd->nGrowVect(), z_nd_mf_loc.nGrowVect()); }