Skip to content

Commit

Permalink
fix bcs with terrain (erf-model#1809)
Browse files Browse the repository at this point in the history 
* fix bcs with terrain

* fix physbc file also
  • Loading branch information
@asalmgren
asalmgren authored Sep 15, 2024
1 parent b203466 commit f750b5b
Show file tree
Hide file tree
Showing 4 changed files with 29 additions and 29 deletions.
12 changes: 6 additions & 6 deletions Source/BoundaryConditions/ERF_BoundaryConditions_cons.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -426,11 +426,11 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& 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);
Expand All @@ -450,7 +450,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& 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));
Expand All @@ -462,7 +462,7 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4<Real>& 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));
Expand Down
16 changes: 11 additions & 5 deletions Source/BoundaryConditions/ERF_BoundaryConditions_xvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,13 +248,17 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4<Real>& 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<var>/dn = 0) must be aware of the surface normal with terrain.
// An additional source term arises from d<var>/dx & d<var>/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);
Expand All @@ -269,17 +273,19 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4<Real>& 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);
Real met_h_eta = Compute_h_eta_AtIface (ii, jj, k0, dxInv, z_phys_nd);
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));
Expand All @@ -290,8 +296,8 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4<Real>& 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));
Expand Down
11 changes: 9 additions & 2 deletions Source/BoundaryConditions/ERF_BoundaryConditions_yvel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,8 +243,13 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& 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<var>/dn = 0) must be aware of the surface normal with terrain.
// An additional source term arises from d<var>/dx & d<var>/dy & met_h_xi/eta/zeta.
//=====================================================================================
Expand All @@ -266,7 +271,9 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& 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);
Expand All @@ -276,7 +283,7 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& 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));
Expand All @@ -288,7 +295,7 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4<Real>& 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));
Expand Down
19 changes: 3 additions & 16 deletions Source/BoundaryConditions/ERF_PhysBCFunct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
Expand All @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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());
}
Expand Down

0 comments on commit f750b5b

Please sign in to comment.