From 83e7a36f4a35fae8efe35c75c1f93fdaca226753 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 2 Sep 2024 21:20:26 -0700 Subject: [PATCH] Update for many scalars (#1782) * update so passive scalars use same bcs * fix * missed a fix * fix bc access in diffusion operations * fix to BC_cons * more fixes for N scalars * fix typo --- .../BoundaryConditions_bndryreg.cpp | 16 +- .../BoundaryConditions_cons.cpp | 354 +++++++++++------- .../BoundaryConditions_xvel.cpp | 232 ++++++------ .../BoundaryConditions_yvel.cpp | 246 ++++++------ .../BoundaryConditions_zvel.cpp | 219 ++++++----- Source/BoundaryConditions/ERF_PhysBCFunct.H | 44 +-- Source/BoundaryConditions/ERF_PhysBCFunct.cpp | 6 +- Source/Diffusion/DiffusionSrcForState_N.cpp | 202 ++++++---- Source/Diffusion/DiffusionSrcForState_T.cpp | 145 ++++--- Source/ERF.H | 4 +- Source/IO/ERF_ReadBndryPlanes.H | 4 +- Source/IO/ERF_ReadBndryPlanes.cpp | 6 +- Source/IndexDefines.H | 15 +- Source/Initialization/ERF_init_bcs.cpp | 39 +- 14 files changed, 842 insertions(+), 690 deletions(-) diff --git a/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp b/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp index 8882d95e4..988ef47d1 100644 --- a/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_bndryreg.cpp @@ -74,14 +74,18 @@ ERF::fill_from_bndryregs (const Vector& mfs, const Real time) ParallelFor( bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].lo(0) == ERFBCType::ext_dir_ingested) { + int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+n; + if (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_ingested) { int jb = std::min(std::max(j,dom_lo.y),dom_hi.y); int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); dest_arr(i,j,k,icomp+n) = bdatxlo(dom_lo.x-1,jb,kb,bccomp+n); } }, bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].hi(0) == ERFBCType::ext_dir_ingested) { + int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+n; + if (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_ingested) { int jb = std::min(std::max(j,dom_lo.y),dom_hi.y); int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); dest_arr(i,j,k,icomp+n) = bdatxhi(dom_hi.x+1,jb,kb,bccomp+n); @@ -100,14 +104,18 @@ ERF::fill_from_bndryregs (const Vector& mfs, const Real time) ParallelFor( bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].lo(1) == ERFBCType::ext_dir_ingested) { + int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+n; + if (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_ingested) { int ib = std::min(std::max(i,dom_lo.x),dom_hi.x); int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); dest_arr(i,j,k,icomp+n) = bdatylo(ib,dom_lo.y-1,kb,bccomp+n); } }, bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].hi(1) == ERFBCType::ext_dir_ingested) { + int bc_comp = (icomp+n >= RhoScalar_comp && icomp+n < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+n; + if (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_ingested) { int ib = std::min(std::max(i,dom_lo.x),dom_hi.x); int kb = std::min(std::max(k,dom_lo.z),dom_hi.z); dest_arr(i,j,k,icomp+n) = bdatyhi(ib,dom_hi.y+1,kb,bccomp+n); diff --git a/Source/BoundaryConditions/BoundaryConditions_cons.cpp b/Source/BoundaryConditions/BoundaryConditions_cons.cpp index e1593a3da..acd86160d 100644 --- a/Source/BoundaryConditions/BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_cons.cpp @@ -12,11 +12,10 @@ using namespace amrex; * @param[in] icomp index into the MultiFab -- this can be any value from 0 to NVAR-1 * @param[in] ncomp the number of components -- this can be any value from 1 to NVAR * as long as icomp+ncomp <= NVAR-1. - * @param[in] bccomp index into m_domain_bcs_type */ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - int icomp, int ncomp, int bccomp, int ngz) + int icomp, int ncomp, int ngz) { BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs); const auto& dom_lo = lbound(domain); @@ -30,24 +29,36 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, // zhi: ori = 5 // Based on BCRec for the domain, we need to make BCRec for this Box - // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - Vector bcrs(icomp+ncomp); - setBC(bx, domain, bccomp, 0, icomp+ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(ncomp); - Gpu::DeviceVector bcrs_d(icomp+ncomp); - Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); - const BCRec* bc_ptr = bcrs_d.data(); + GpuArray,NBCVAR_max> l_bc_extdir_vals_d; + + const int* bxlo = bx.loVect(); + const int* bxhi = bx.hiVect(); + const int* dlo = domain.loVect(); + const int* dhi = domain.hiVect(); - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; - for (int i = 0; i < icomp+ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int nc = 0; nc < ncomp; nc++) + { + int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+nc; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir] + ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir )); + bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir] + ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir )); + } + + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[nc][ori] = m_bc_extdir_vals[bc_comp][ori]; + } + } - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_neumann_vals_d; - for (int i = 0; i < icomp+ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_neumann_vals_d[i][ori] = m_bc_neumann_vals[bccomp+i][ori]; + Gpu::DeviceVector bcrs_d(ncomp); + Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); + const BCRec* bc_ptr = bcrs_d.data(); GeometryData const& geomdata = m_geom.data(); bool is_periodic_in_x = geomdata.isPeriodic(0); @@ -59,20 +70,28 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].lo(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][0]; - } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::ext_dir_prim) { + bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(0); + + if (l_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][0]; + } else if (l_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(dom_lo.x,j,k,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][0]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][0]; } }, - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].hi(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][3]; - } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::ext_dir_prim) { + bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(0); + + if (h_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][3]; + } else if (h_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(dom_hi.x,j,k,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][3]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][3]; } } ); @@ -83,20 +102,26 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); ParallelFor( - bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].lo(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][1]; - } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::ext_dir_prim) { + bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(1); + if (l_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][1]; + } else if (l_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(i,dom_lo.y,k,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][1]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][1]; } }, - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].hi(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][4]; - } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::ext_dir_prim) { + bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(1); + if (h_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][4]; + } else if (h_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(i,dom_hi.y,k,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][4]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][4]; } } ); @@ -114,32 +139,38 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, if (bx_xhi.smallEnd(2) != domain.smallEnd(2)) bx_xhi.growLo(2,ngz); if (bx_xhi.bigEnd(2) != domain.bigEnd(2)) bx_xhi.growHi(2,ngz); ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(0); int iflip = dom_lo.x - 1 - i; - if (bc_ptr[icomp+n].lo(0) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_lo.x,j,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_lo.x,j,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(iflip,j,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(iflip,j,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(0) == ERFBCType::hoextrapcc) { - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(dom_lo.x,j,k,icomp+n) - dest_arr(dom_lo.x+1,j,k,icomp+n) ; + if (l_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_lo.x,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp); + } else if (l_bc_type == ERFBCType::hoextrapcc) { + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(dom_lo.x,j,k,dest_comp) - dest_arr(dom_lo.x+1,j,k,dest_comp) ; } }, - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(0); int iflip = 2*dom_hi.x + 1 - i; - if (bc_ptr[icomp+n].hi(0) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_hi.x,j,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(dom_hi.x,j,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(iflip,j,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(iflip,j,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(0) == ERFBCType::hoextrapcc) { - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(dom_hi.x,j,k,icomp+n) - dest_arr(dom_hi.x-1,j,k,icomp+n) ; + if (h_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(dom_hi.x,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(iflip,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(iflip,j,k,dest_comp); + } else if (h_bc_type == ERFBCType::hoextrapcc) { + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(dom_hi.x,j,k,dest_comp) - dest_arr(dom_hi.x-1,j,k,dest_comp) ; } } ); @@ -155,33 +186,39 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, if (bx_yhi.smallEnd(2) != domain.smallEnd(2)) bx_yhi.growLo(2,ngz); if (bx_yhi.bigEnd(2) != domain.bigEnd(2)) bx_yhi.growHi(2,ngz); ParallelFor( - bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(1); int jflip = dom_lo.y - 1 - j; - if (bc_ptr[icomp+n].lo(1) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_lo.y,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_lo.y,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,jflip,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(i,jflip,k,icomp+n); - } else if (bc_ptr[icomp+n].lo(1) == ERFBCType::hoextrapcc) { - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,dom_lo.y,k,icomp+n) - dest_arr(i,dom_lo.y+1,k,icomp+n) ; + if (l_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp); + } else if (l_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_lo.y,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp); + } else if (l_bc_type == ERFBCType::hoextrapcc) { + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(i,dom_lo.y,k,dest_comp) - dest_arr(i,dom_lo.y+1,k,dest_comp) ; } }, - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(1); int jflip = 2*dom_hi.y + 1 - j; - if (bc_ptr[icomp+n].hi(1) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_hi.y,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,dom_hi.y,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,jflip,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(i,jflip,k,icomp+n); - } else if (bc_ptr[icomp+n].hi(1) == ERFBCType::hoextrapcc) { - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,dom_hi.y,k,icomp+n) - dest_arr(i,dom_hi.y-1,k,icomp+n); + if (h_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp); + } else if (h_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,dom_hi.y,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,jflip,k,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,jflip,k,dest_comp); + } else if (h_bc_type == ERFBCType::hoextrapcc) { + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(i,dom_hi.y,k,dest_comp) - dest_arr(i,dom_hi.y-1,k,dest_comp); } } ); @@ -200,13 +237,12 @@ void ERFPhysBCFunct_cons::impose_lateral_cons_bcs (const Array4& dest_arr, * @param[in] icomp the index of the first component to be filled * @param[in] ncomp the number of components -- this can be any value from 1 to NVAR * as long as icomp+ncomp <= NVAR-1. - * @param[in] bccomp index into m_domain_bcs_type */ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, const Array4& z_phys_nd, const GpuArray dxInv, - int icomp, int ncomp, int bccomp) + int icomp, int ncomp) { BL_PROFILE_VAR("impose_lateral_cons_bcs()",impose_lateral_cons_bcs); const auto& dom_lo = lbound(domain); @@ -229,44 +265,64 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr // zhi: ori = 5 // Based on BCRec for the domain, we need to make BCRec for this Box - // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - Vector bcrs(icomp+ncomp); - setBC(bx, domain, bccomp, 0, icomp+ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(ncomp); + GpuArray,NBCVAR_max> l_bc_extdir_vals_d; + GpuArray,NBCVAR_max> l_bc_neumann_vals_d; + + const int* bxlo = bx.loVect(); + const int* bxhi = bx.hiVect(); + const int* dlo = domain.loVect(); + const int* dhi = domain.hiVect(); + + for (int nc = 0; nc < ncomp; nc++) + { + int bc_comp = (icomp+nc >= RhoScalar_comp && icomp+nc < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : icomp+nc; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + bcrs[nc].setLo(dir, ( bxlo[dir]<=dlo[dir] + ? m_domain_bcs_type[bc_comp].lo(dir) : BCType::int_dir )); + bcrs[nc].setHi(dir, ( bxhi[dir]>=dhi[dir] + ? m_domain_bcs_type[bc_comp].hi(dir) : BCType::int_dir )); + } + + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[nc][ori] = m_bc_extdir_vals[bc_comp][ori]; + l_bc_neumann_vals_d[nc][ori] = m_bc_neumann_vals[bc_comp][ori]; + } + } Gpu::DeviceVector bcrs_d(icomp+ncomp); Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); const BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; - for (int i = 0; i < icomp+ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; - - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_neumann_vals_d; - for (int i = 0; i < icomp+ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_neumann_vals_d[i][ori] = m_bc_neumann_vals[bccomp+i][ori]; - - { Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1); Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1); ParallelFor( - bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].lo(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][2]; - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::ext_dir_prim) { + bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(2); + + if (l_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][2]; + + } else if (l_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(i,j,dom_lo.z,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][2]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][2]; } }, - bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - if (bc_ptr[icomp+n].hi(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k,icomp+n) = l_bc_extdir_vals_d[icomp+n][5]; - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::ext_dir_prim) { + bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(2); + if (h_bc_type == ERFBCType::ext_dir) { + dest_arr(i,j,k,dest_comp) = l_bc_extdir_vals_d[n][5]; + } else if (h_bc_type == ERFBCType::ext_dir_prim) { Real rho = dest_arr(i,j,dom_hi.z,Rho_comp); - dest_arr(i,j,k,icomp+n) = rho * l_bc_extdir_vals_d[icomp+n][5]; + dest_arr(i,j,k,dest_comp) = rho * l_bc_extdir_vals_d[n][5]; } } @@ -278,45 +334,51 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1); // Populate ghost cells on lo-z and hi-z domain boundaries ParallelFor( - bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(2); int kflip = dom_lo.z - 1 - i; - if (bc_ptr[icomp+n].lo(2) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_lo.z,icomp+n); - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_lo.z,icomp+n); - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,kflip,icomp+n); - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(i,j,kflip,icomp+n); - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::neumann) { + if (l_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp); + } else if (l_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp); + } else if (l_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp); + } else if (l_bc_type == ERFBCType::neumann) { Real delta_z = (dom_lo.z - k) / dxInv[2]; - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_lo.z,icomp+n) - - delta_z*l_bc_neumann_vals_d[icomp+n][2]*dest_arr(i,j,dom_lo.z,Rho_comp); - } else if (bc_ptr[icomp+n].lo(2) == ERFBCType::hoextrapcc) { - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,j,dom_lo.z,icomp+n) - dest_arr(i,j,dom_lo.z+1,icomp+n); + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_lo.z,dest_comp) - + delta_z*l_bc_neumann_vals_d[n][2]*dest_arr(i,j,dom_lo.z,Rho_comp); + } else if (l_bc_type == ERFBCType::hoextrapcc) { + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(i,j,dom_lo.z,dest_comp) - dest_arr(i,j,dom_lo.z+1,dest_comp); } }, - bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + int dest_comp = icomp+n; + int h_bc_type = bc_ptr[n].hi(2); int kflip = 2*dom_hi.z + 1 - i; - if (bc_ptr[icomp+n].hi(2) == ERFBCType::foextrap) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_hi.z,icomp+n); - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::open) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_hi.z,icomp+n); - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::reflect_even) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,kflip,icomp+n); - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::reflect_odd) { - dest_arr(i,j,k,icomp+n) = -dest_arr(i,j,kflip,icomp+n); - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::neumann) { + if (h_bc_type == ERFBCType::foextrap) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp); + } else if (h_bc_type == ERFBCType::open) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_even) { + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,kflip,dest_comp); + } else if (h_bc_type == ERFBCType::reflect_odd) { + dest_arr(i,j,k,dest_comp) = -dest_arr(i,j,kflip,dest_comp); + } else if (h_bc_type == ERFBCType::neumann) { Real delta_z = (k - dom_hi.z) / dxInv[2]; if( (icomp+n) == Rho_comp ) { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_hi.z,icomp+n) + - delta_z*l_bc_neumann_vals_d[icomp+n][5]; + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) + + delta_z*l_bc_neumann_vals_d[n][5]; } else { - dest_arr(i,j,k,icomp+n) = dest_arr(i,j,dom_hi.z,icomp+n) + - delta_z*l_bc_neumann_vals_d[icomp+n][5]*dest_arr(i,j,dom_hi.z,Rho_comp); + dest_arr(i,j,k,dest_comp) = dest_arr(i,j,dom_hi.z,dest_comp) + + delta_z*l_bc_neumann_vals_d[n][5]*dest_arr(i,j,dom_hi.z,Rho_comp); } - } else if (bc_ptr[icomp+n].hi(2) == ERFBCType::hoextrapcc){ - dest_arr(i,j,k,icomp+n) = 2.0*dest_arr(i,j,dom_hi.z,icomp+n) - dest_arr(i,j,dom_hi.z-1,icomp+n); + } else if (h_bc_type == ERFBCType::hoextrapcc){ + dest_arr(i,j,k,dest_comp) = 2.0*dest_arr(i,j,dom_hi.z,dest_comp) - dest_arr(i,j,dom_hi.z-1,dest_comp); } } ); @@ -331,9 +393,11 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr //===================================================================================== // Only modify scalars, U, or V // Loop over each component - for (int n = icomp; n < icomp+ncomp; n++) { + for (int n = 0; n < ncomp; n++) { // Hit for Neumann condition at kmin - if( bcrs[n].lo(2) == ERFBCType::foextrap) + int dest_comp = icomp+n; + int l_bc_type = bc_ptr[n].lo(2); + if(l_bc_type == ERFBCType::foextrap) { // Loop over ghost cells in bottom XY-plane (valid box) Box xybx = bx; @@ -365,11 +429,11 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr if (i < dom_lo.x-1 || i > dom_hi.x+1) { GradVarx = 0.0; } else if (i+1 > bx_hi.x) { - GradVarx = dxInv[0] * (dest_arr(i ,j,k0,n) - dest_arr(i-1,j,k0,n)); + GradVarx = dxInv[0] * (dest_arr(i ,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp)); } else if (i-1 < bx_lo.x) { - GradVarx = dxInv[0] * (dest_arr(i+1,j,k0,n) - dest_arr(i ,j,k0,n)); + GradVarx = dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i ,j,k0,dest_comp)); } else { - GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0,n) - dest_arr(i-1,j,k0,n)); + GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0,dest_comp) - dest_arr(i-1,j,k0,dest_comp)); } // GradY at IJK location inside domain -- this relies on the assumption that we have @@ -377,18 +441,18 @@ void ERFPhysBCFunct_cons::impose_vertical_cons_bcs (const Array4& dest_arr if (j < dom_lo.y-1 || j > dom_hi.y+1) { GradVary = 0.0; } else if (j+1 > bx_hi.y) { - GradVary = dxInv[1] * (dest_arr(i,j ,k0,n) - dest_arr(i,j-1,k0,n)); + GradVary = dxInv[1] * (dest_arr(i,j ,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp)); } else if (j-1 < bx_lo.y) { - GradVary = dxInv[1] * (dest_arr(i,j+1,k0,n) - dest_arr(i,j ,k0,n)); + GradVary = dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j ,k0,dest_comp)); } else { - GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0,n) - dest_arr(i,j-1,k0,n)); + GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0,dest_comp) - dest_arr(i,j-1,k0,dest_comp)); } // Prefactor Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. ); // Accumulate in bottom ghost cell (EXTRAP already populated) - dest_arr(i,j,k,n) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); + dest_arr(i,j,k,dest_comp) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); }); } // box includes k0 } // foextrap diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index 2c3eb8edd..a42da9d39 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -31,19 +31,18 @@ void ERFPhysBCFunct_u::impose_lateral_xvel_bcs (const Array4& dest_arr, // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; - Vector bcrs(ncomp); - setBC(enclosedCells(bx), domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(1); + setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs); - Gpu::DeviceVector bcrs_d(ncomp); + Gpu::DeviceVector bcrs_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); const BCRec* bc_ptr = bcrs_d.data(); GpuArray,1> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } GeometryData const& geomdata = m_geom.data(); bool is_periodic_in_x = geomdata.isPeriodic(0); @@ -57,58 +56,58 @@ void ERFPhysBCFunct_u::impose_lateral_xvel_bcs (const Array4& dest_arr, Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+2); Box bx_xlo_face(bx); bx_xlo_face.setSmall(0,dom_lo.x ); bx_xlo_face.setBig(0,dom_lo.x ); Box bx_xhi_face(bx); bx_xhi_face.setSmall(0,dom_hi.x+1); bx_xhi_face.setBig(0,dom_hi.x+1); - ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + ParallelFor(bx_xlo, bx_xlo_face, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = dom_lo.x - i; - if (bc_ptr[n].lo(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][0]; - } else if (bc_ptr[n].lo(0) == ERFBCType::foextrap) { + if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0]; + } else if (bc_ptr[0].lo(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::open) { + } else if (bc_ptr[0].lo(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::neumann_int) { + } else if (bc_ptr[0].lo(0) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/3.0; } }, // We only set the values on the domain faces themselves if EXT_DIR - bx_xlo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (bc_ptr[n].lo(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][0]; - } else if (bc_ptr[n].lo(0) == ERFBCType::neumann_int) { + if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0]; + } else if (bc_ptr[0].lo(0) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k))/3.0; } } ); - ParallelFor( - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + ParallelFor(bx_xhi, bx_xhi_face, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = 2*(dom_hi.x + 1) - i; - if (bc_ptr[n].hi(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][3]; - } else if (bc_ptr[n].hi(0) == ERFBCType::foextrap) { + if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3]; + } else if (bc_ptr[0].hi(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x+1,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::open) { + } else if (bc_ptr[0].hi(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_hi.x+1,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::neumann_int) { + } else if (bc_ptr[0].hi(0) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/3.0; } }, // We only set the values on the domain faces themselves if EXT_DIR - bx_xhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (bc_ptr[n].hi(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][3]; - } else if (bc_ptr[n].hi(0) == ERFBCType::neumann_int) { + if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3]; + } else if (bc_ptr[0].hi(0) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k))/3.0; } } @@ -121,32 +120,32 @@ void ERFPhysBCFunct_u::impose_lateral_xvel_bcs (const Array4& dest_arr, Real* xvel_bc_ptr = m_u_bc_data; Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); - ParallelFor( - bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int jflip = dom_lo.y - 1 - j; - if (bc_ptr[n].lo(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][1]; - } else if (bc_ptr[n].lo(1) == ERFBCType::foextrap) { + if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1]; + } else if (bc_ptr[0].lo(1) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::open) { + } else if (bc_ptr[0].lo(1) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,jflip,k); } }, - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int jflip = 2*dom_hi.y + 1 - j; - if (bc_ptr[n].hi(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[n][4]; - } else if (bc_ptr[n].hi(1) == ERFBCType::foextrap) { + if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (xvel_bc_ptr) ? xvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4]; + } else if (bc_ptr[0].hi(1) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::open) { + } else if (bc_ptr[0].hi(1) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,jflip,k); } } @@ -193,54 +192,53 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; - Vector bcrs(ncomp); - setBC(enclosedCells(bx), domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(1); + setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs); - Gpu::DeviceVector bcrs_d(ncomp); + Gpu::DeviceVector bcrs_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); const BCRec* bc_ptr = bcrs_d.data(); GpuArray,1> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } { // Populate ghost cells on lo-z and hi-z domain boundaries Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1); Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1); - ParallelFor( - bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + ParallelFor(bx_zlo, bx_zhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int kflip = dom_lo.z - 1 - k; - if (bc_ptr[n].lo(2) == ERFBCType::ext_dir) { + if (bc_ptr[0].lo(2) == ERFBCType::ext_dir) { #ifdef ERF_USE_TERRAIN_VELOCITY dest_arr(i,j,k) = prob->compute_terrain_velocity(time); #else - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][2]; #endif - } else if (bc_ptr[n].lo(2) == ERFBCType::foextrap) { + } else if (bc_ptr[0].lo(2) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); - } else if (bc_ptr[n].lo(2) == ERFBCType::open) { + } else if (bc_ptr[0].lo(2) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); - } else if (bc_ptr[n].lo(2) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr[n].lo(2) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,j,kflip); } }, - bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int kflip = 2*dom_hi.z + 1 - k; - if (bc_ptr[n].hi(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; - } else if (bc_ptr[n].hi(2) == ERFBCType::foextrap) { + if (bc_ptr[0].hi(2) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5]; + } else if (bc_ptr[0].hi(2) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); - } else if (bc_ptr[n].hi(2) == ERFBCType::open) { + } else if (bc_ptr[0].hi(2) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); - } else if (bc_ptr[n].hi(2) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr[n].hi(2) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,j,kflip); } } @@ -255,61 +253,61 @@ void ERFPhysBCFunct_u::impose_vertical_xvel_bcs (const Array4& dest_arr, //===================================================================================== // Only modify scalars, U, or V // Loop over each component - for (int n = 0; n < ncomp; n++) { - // Hit for Neumann condition at kmin - if(bcrs[n].lo(2) == ERFBCType::foextrap) { - // Loop over ghost cells in bottom XY-plane (valid box) - Box xybx = bx; - xybx.setBig(2,-1); - xybx.setSmall(2,bx.smallEnd()[2]); - int k0 = 0; + // Hit for Neumann condition at kmin + if(bcrs[0].lo(2) == ERFBCType::foextrap) { + // Loop over ghost cells in bottom XY-plane (valid box) + Box xybx = bx; + xybx.setBig(2,-1); + xybx.setSmall(2,bx.smallEnd()[2]); + int k0 = 0; - // Get the dz cell size - Real dz = geomdata.CellSize(2); + // Get the dz cell size + Real dz = geomdata.CellSize(2); - // Fill all the Neumann srcs with terrain - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); - int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y); + // Fill all the Neumann srcs with terrain + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); + int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_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); + // 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 - Real GradVarx, GradVary; - if (i < dom_lo.x-1 || i > dom_hi.x+1) - GradVarx = 0.0; - else if (i+1 > bx_hi.x) - GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); - else if (i-1 < bx_lo.x) - GradVarx = dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i ,j,k0)); - else - GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0)); + // 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) { + GradVarx = 0.0; + } else if (i+1 > bx_hi.x) { + GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); + } else if (i-1 < bx_lo.x) { + GradVarx = dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i ,j,k0)); + } else { + GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0)); + } - // 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) - GradVary = 0.0; - else if (j+1 > bx_hi.y) - GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); - else if (j-1 < bx_lo.y) - GradVary = dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j ,k0)); - else - GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0)); + // 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) { + GradVary = 0.0; + } else if (j+1 > bx_hi.y) { + GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); + } else if (j-1 < bx_lo.y) { + GradVary = dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j ,k0)); + } else { + GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0)); + } - // Prefactor - Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. ); + // Prefactor + Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. ); - // Accumulate in bottom ghost cell (EXTRAP already populated) - dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); + // Accumulate in bottom ghost cell (EXTRAP already populated) + dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); }); - } // foextrap - } // ncomp + } // foextrap } //m_z_phys_nd Gpu::streamSynchronize(); } diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index 2fe5825ef..c1252214e 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -22,9 +22,8 @@ void ERFPhysBCFunct_v::impose_lateral_yvel_bcs (const Array4& dest_arr, // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; - Vector bcrs(ncomp); - setBC(enclosedCells(bx), domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(1); + setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs); // xlo: ori = 0 // ylo: ori = 1 @@ -33,15 +32,15 @@ void ERFPhysBCFunct_v::impose_lateral_yvel_bcs (const Array4& dest_arr, // yhi: ori = 4 // zhi: ori = 5 - Gpu::DeviceVector bcrs_d(ncomp); + Gpu::DeviceVector bcrs_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); const BCRec* bc_ptr = bcrs_d.data(); GpuArray, 1> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } GeometryData const& geomdata = m_geom.data(); bool is_periodic_in_x = geomdata.isPeriodic(0); @@ -54,32 +53,32 @@ void ERFPhysBCFunct_v::impose_lateral_yvel_bcs (const Array4& dest_arr, Real* yvel_bc_ptr = m_v_bc_data; Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); - ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = dom_lo.x - 1- i; - if (bc_ptr[n].lo(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][0]; - } else if (bc_ptr[n].lo(0) == ERFBCType::foextrap) { + if (bc_ptr[0].lo(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0]; + } else if (bc_ptr[0].lo(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::open) { + } else if (bc_ptr[0].lo(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr[n].lo(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); } }, - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = 2*dom_hi.x + 1 - i; - if (bc_ptr[n].hi(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][3]; - } else if (bc_ptr[n].hi(0) == ERFBCType::foextrap) { + if (bc_ptr[0].hi(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3]; + } else if (bc_ptr[0].hi(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::open) { + } else if (bc_ptr[0].hi(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr[n].hi(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].hi(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); } } @@ -94,58 +93,58 @@ void ERFPhysBCFunct_v::impose_lateral_yvel_bcs (const Array4& dest_arr, Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+2); Box bx_ylo_face(bx); bx_ylo_face.setSmall(1,dom_lo.y ); bx_ylo_face.setBig(1,dom_lo.y ); Box bx_yhi_face(bx); bx_yhi_face.setSmall(1,dom_hi.y+1); bx_yhi_face.setBig(1,dom_hi.y+1); - ParallelFor( - bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + ParallelFor(bx_ylo, bx_ylo_face, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int jflip = dom_lo.y-j; - if (bc_ptr[n].lo(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][1]; - } else if (bc_ptr[n].lo(1) == ERFBCType::foextrap) { + if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1]; + } else if (bc_ptr[0].lo(1) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::open) { + } else if (bc_ptr[0].lo(1) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(1) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,jflip,k); - } else if (bc_ptr[n].lo(1) == ERFBCType::neumann_int) { + } else if (bc_ptr[0].lo(1) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/3.0; } }, // We only set the values on the domain faces themselves if EXT_DIR - bx_ylo_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (bc_ptr[n].lo(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][1]; - } else if (bc_ptr[n].lo(1) == ERFBCType::neumann_int) { + if (bc_ptr[0].lo(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1]; + } else if (bc_ptr[0].lo(1) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(i,dom_lo.y+1,k) - dest_arr(i,dom_lo.y+2,k))/3.0; } } ); - ParallelFor( - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + ParallelFor(bx_yhi, bx_yhi_face, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { - int jflip = 2*(dom_hi.y + 1) - j; - if (bc_ptr[n].hi(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][4]; - } else if (bc_ptr[n].hi(1) == ERFBCType::foextrap) { - dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::open) { - dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::reflect_even) { - dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::reflect_odd) { - dest_arr(i,j,k) = -dest_arr(i,jflip,k); - } else if (bc_ptr[n].hi(1) == ERFBCType::neumann_int) { - dest_arr(i,j,k) = (4.0*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/3.0; + int jflip = 2*(dom_hi.y + 1) - j; + if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4]; + } else if (bc_ptr[0].hi(1) == ERFBCType::foextrap) { + dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k); + } else if (bc_ptr[0].hi(1) == ERFBCType::open) { + dest_arr(i,j,k) = dest_arr(i,dom_hi.y+1,k); + } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_even) { + dest_arr(i,j,k) = dest_arr(i,jflip,k); + } else if (bc_ptr[0].hi(1) == ERFBCType::reflect_odd) { + dest_arr(i,j,k) = -dest_arr(i,jflip,k); + } else if (bc_ptr[0].hi(1) == ERFBCType::neumann_int) { + dest_arr(i,j,k) = (4.0*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/3.0; } }, // We only set the values on the domain faces themselves if EXT_DIR - bx_yhi_face, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (bc_ptr[n].hi(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[n][4]; - } else if (bc_ptr[n].hi(1) == ERFBCType::neumann_int) { + if (bc_ptr[0].hi(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (yvel_bc_ptr) ? yvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4]; + } else if (bc_ptr[0].hi(1) == ERFBCType::neumann_int) { dest_arr(i,j,k) = (4.0*dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k))/3.0; } } @@ -185,9 +184,8 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; - Vector bcrs(ncomp); - setBC(enclosedCells(bx), domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs); + Vector bcrs(1); + setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs); // xlo: ori = 0 // ylo: ori = 1 @@ -196,15 +194,15 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, // yhi: ori = 4 // zhi: ori = 5 - Gpu::DeviceVector bcrs_d(ncomp); + Gpu::DeviceVector bcrs_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs.begin(), bcrs.end(), bcrs_d.begin()); const BCRec* bc_ptr = bcrs_d.data(); GpuArray, 1> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } GeometryData const& geomdata = m_geom.data(); @@ -212,32 +210,32 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, // Populate ghost cells on lo-z and hi-z domain boundaries Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1); Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1); - ParallelFor( - bx_zlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + ParallelFor(bx_zlo, bx_zhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int kflip = dom_lo.z - 1 - k; - if (bc_ptr[n].lo(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][2]; - } else if (bc_ptr[n].lo(2) == ERFBCType::foextrap) { + if (bc_ptr[0].lo(2) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][2]; + } else if (bc_ptr[0].lo(2) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); - } else if (bc_ptr[n].lo(2) == ERFBCType::open) { + } else if (bc_ptr[0].lo(2) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z); - } else if (bc_ptr[n].lo(2) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr[n].lo(2) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].lo(2) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,j,kflip); } }, - bx_zhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int kflip = 2*dom_hi.z + 1 - k; - if (bc_ptr[n].hi(2) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = l_bc_extdir_vals_d[n][5]; - } else if (bc_ptr[n].hi(2) == ERFBCType::foextrap) { + if (bc_ptr[0].hi(2) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = l_bc_extdir_vals_d[0][5]; + } else if (bc_ptr[0].hi(2) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); - } else if (bc_ptr[n].hi(2) == ERFBCType::open) { + } else if (bc_ptr[0].hi(2) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z); - } else if (bc_ptr[n].hi(2) == ERFBCType::reflect_even) { + } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(i,j,kflip); - } else if (bc_ptr[n].hi(2) == ERFBCType::reflect_odd) { + } else if (bc_ptr[0].hi(2) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(i,j,kflip); } } @@ -252,61 +250,61 @@ void ERFPhysBCFunct_v::impose_vertical_yvel_bcs (const Array4& dest_arr, //===================================================================================== // Only modify scalars, U, or V // Loop over each component - for (int n = 0; n < ncomp; n++) { - // Hit for Neumann condition at kmin - if(bcrs[n].lo(2) == ERFBCType::foextrap) { - // Loop over ghost cells in bottom XY-plane (valid box) - Box xybx = bx; - xybx.setBig(2,-1); - xybx.setSmall(2,bx.smallEnd()[2]); - int k0 = 0; + // Hit for Neumann condition at kmin + if(bcrs[0].lo(2) == ERFBCType::foextrap) { + // Loop over ghost cells in bottom XY-plane (valid box) + Box xybx = bx; + xybx.setBig(2,-1); + xybx.setSmall(2,bx.smallEnd()[2]); + int k0 = 0; - // Get the dz cell size - Real dz = geomdata.CellSize(2); + // Get the dz cell size + Real dz = geomdata.CellSize(2); - // Fill all the Neumann srcs with terrain - ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - // Clip indices for ghost-cells - int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); - int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y); + // Fill all the Neumann srcs with terrain + ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // Clip indices for ghost-cells + int ii = amrex::min(amrex::max(i,perdom_lo.x),perdom_hi.x); + int jj = amrex::min(amrex::max(j,perdom_lo.y),perdom_hi.y); - // Get metrics - Real met_h_xi = Compute_h_xi_AtJface (ii, jj, k0, dxInv, z_phys_nd); - Real met_h_eta = Compute_h_eta_AtJface (ii, jj, k0, dxInv, z_phys_nd); - Real met_h_zeta = Compute_h_zeta_AtJface(ii, jj, k0, dxInv, z_phys_nd); + // Get metrics + Real met_h_xi = Compute_h_xi_AtJface (ii, jj, k0, dxInv, z_phys_nd); + Real met_h_eta = Compute_h_eta_AtJface (ii, jj, k0, dxInv, z_phys_nd); + Real met_h_zeta = Compute_h_zeta_AtJface(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 - Real GradVarx, GradVary; - if (i < dom_lo.x-1 || i > dom_hi.x+1) - GradVarx = 0.0; - else if (i+1 > bx_hi.x) - GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); - else if (i-1 < bx_lo.x) - GradVarx = dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i ,j,k0)); - else - GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0)); + // 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) { + GradVarx = 0.0; + } else if (i+1 > bx_hi.x) { + GradVarx = dxInv[0] * (dest_arr(i ,j,k0) - dest_arr(i-1,j,k0)); + } else if (i-1 < bx_lo.x) { + GradVarx = dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i ,j,k0)); + } else { + GradVarx = 0.5 * dxInv[0] * (dest_arr(i+1,j,k0) - dest_arr(i-1,j,k0)); + } - // 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) - GradVary = 0.0; - else if (j+1 > bx_hi.y) - GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); - else if (j-1 < bx_lo.y) - GradVary = dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j ,k0)); - else - GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0)); + // 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) { + GradVary = 0.0; + } else if (j+1 > bx_hi.y) { + GradVary = dxInv[1] * (dest_arr(i,j ,k0) - dest_arr(i,j-1,k0)); + } else if (j-1 < bx_lo.y) { + GradVary = dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j ,k0)); + } else { + GradVary = 0.5 * dxInv[1] * (dest_arr(i,j+1,k0) - dest_arr(i,j-1,k0)); + } - // Prefactor - Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. ); + // Prefactor + Real met_fac = met_h_zeta / ( met_h_xi*met_h_xi + met_h_eta*met_h_eta + 1. ); - // Accumulate in bottom ghost cell (EXTRAP already populated) - dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); - }); - } // foextrap - } // ncomp + // Accumulate in bottom ghost cell (EXTRAP already populated) + dest_arr(i,j,k) -= dz * met_fac * ( met_h_xi * GradVarx + met_h_eta * GradVary ); + }); + } // foextrap } //m_z_phys_nd Gpu::streamSynchronize(); } diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index d50c9aecd..935177505 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -28,7 +28,6 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; Vector bcrs_w(1); setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs_w); @@ -39,7 +38,7 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a // yhi: ori = 4 // zhi: ori = 5 - Gpu::DeviceVector bcrs_w_d(ncomp); + Gpu::DeviceVector bcrs_w_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs_w.begin(), bcrs_w.end(), bcrs_w_d.begin()); const BCRec* bc_ptr_w = bcrs_w_d.data(); @@ -47,9 +46,9 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a bool l_use_terrain = (m_z_phys_nd != nullptr); - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } GeometryData const& geomdata = m_geom.data(); bool is_periodic_in_x = geomdata.isPeriodic(0); @@ -63,38 +62,38 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a Real* zvel_bc_ptr = m_w_bc_data; Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1); Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+1); - ParallelFor( - bx_xlo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + ParallelFor(bx_xlo, bx_xhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = dom_lo.x - 1 - i; - if (bc_ptr_w[n].lo(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[n][0]; + if (bc_ptr_w[0].lo(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[0][0]; if (l_use_terrain) { dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); } - } else if (bc_ptr_w[n].lo(0) == ERFBCType::foextrap) { + } else if (bc_ptr_w[0].lo(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr_w[n].lo(0) == ERFBCType::open) { + } else if (bc_ptr_w[0].lo(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_lo.x,j,k); - } else if (bc_ptr_w[n].lo(0) == ERFBCType::reflect_even) { + } else if (bc_ptr_w[0].lo(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr_w[n].lo(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr_w[0].lo(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); } }, - bx_xhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { + [=] AMREX_GPU_DEVICE (int i, int j, int k) { int iflip = 2*dom_hi.x + 1 - i; - if (bc_ptr_w[n].hi(0) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[n][3]; + if (bc_ptr_w[0].hi(0) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[0][3]; if (l_use_terrain) { dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); } - } else if (bc_ptr_w[n].hi(0) == ERFBCType::foextrap) { + } else if (bc_ptr_w[0].hi(0) == ERFBCType::foextrap) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); - } else if (bc_ptr_w[n].hi(0) == ERFBCType::open) { + } else if (bc_ptr_w[0].hi(0) == ERFBCType::open) { dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k); - } else if (bc_ptr_w[n].hi(0) == ERFBCType::reflect_even) { + } else if (bc_ptr_w[0].hi(0) == ERFBCType::reflect_even) { dest_arr(i,j,k) = dest_arr(iflip,j,k); - } else if (bc_ptr_w[n].hi(0) == ERFBCType::reflect_odd) { + } else if (bc_ptr_w[0].hi(0) == ERFBCType::reflect_odd) { dest_arr(i,j,k) = -dest_arr(iflip,j,k); } } @@ -107,42 +106,42 @@ void ERFPhysBCFunct_w::impose_lateral_zvel_bcs (const Array4& dest_a Real* zvel_bc_ptr = m_w_bc_data; Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1); Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1); - ParallelFor(bx_ylo, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int jflip = dom_lo.y - 1 - j; - if (bc_ptr_w[n].lo(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[n][1]; - if (l_use_terrain) { - dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + ParallelFor(bx_ylo, bx_yhi, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int jflip = dom_lo.y - 1 - j; + if (bc_ptr_w[0].lo(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[0][1]; + if (l_use_terrain) { + dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + } + } else if (bc_ptr_w[0].lo(1) == ERFBCType::foextrap) { + dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); + } else if (bc_ptr_w[0].lo(1) == ERFBCType::open) { + dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); + } else if (bc_ptr_w[0].lo(1) == ERFBCType::reflect_even) { + dest_arr(i,j,k) = dest_arr(i,jflip,k); + } else if (bc_ptr_w[0].lo(1) == ERFBCType::reflect_odd) { + dest_arr(i,j,k) = -dest_arr(i,jflip,k); } - } else if (bc_ptr_w[n].lo(1) == ERFBCType::foextrap) { - dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr_w[n].lo(1) == ERFBCType::open) { - dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k); - } else if (bc_ptr_w[n].lo(1) == ERFBCType::reflect_even) { - dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr_w[n].lo(1) == ERFBCType::reflect_odd) { - dest_arr(i,j,k) = -dest_arr(i,jflip,k); - } - }, - bx_yhi, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - int jflip = 2*dom_hi.y + 1 - j; - if (bc_ptr_w[n].hi(1) == ERFBCType::ext_dir) { - dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[n][4]; - if (l_use_terrain) { - dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int jflip = 2*dom_hi.y + 1 - j; + if (bc_ptr_w[0].hi(1) == ERFBCType::ext_dir) { + dest_arr(i,j,k) = (zvel_bc_ptr) ? zvel_bc_ptr[k] : l_bc_extdir_vals_d[0][4]; + if (l_use_terrain) { + dest_arr(i,j,k) = WFromOmega(i,j,k,dest_arr(i,j,k),xvel_arr,yvel_arr,z_phys_nd,dxInv); + } + } else if (bc_ptr_w[0].hi(1) == ERFBCType::foextrap) { + dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); + } else if (bc_ptr_w[0].hi(1) == ERFBCType::open) { + dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); + } else if (bc_ptr_w[0].hi(1) == ERFBCType::reflect_even) { + dest_arr(i,j,k) = dest_arr(i,jflip,k); + } else if (bc_ptr_w[0].hi(1) == ERFBCType::reflect_odd) { + dest_arr(i,j,k) = -dest_arr(i,jflip,k); } - } else if (bc_ptr_w[n].hi(1) == ERFBCType::foextrap) { - dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); - } else if (bc_ptr_w[n].hi(1) == ERFBCType::open) { - dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k); - } else if (bc_ptr_w[n].hi(1) == ERFBCType::reflect_even) { - dest_arr(i,j,k) = dest_arr(i,jflip,k); - } else if (bc_ptr_w[n].hi(1) == ERFBCType::reflect_odd) { - dest_arr(i,j,k) = -dest_arr(i,jflip,k); - } - }); + }); } - Gpu::streamSynchronize(); } @@ -183,7 +182,6 @@ void ERFPhysBCFunct_w::impose_vertical_zvel_bcs (const Array4& dest_arr, // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; Vector bcrs_u(1), bcrs_v(1), bcrs_w(1); setBC(enclosedCells(bx), domain, bccomp_u, 0, 1, m_domain_bcs_type, bcrs_u); setBC(enclosedCells(bx), domain, bccomp_v, 0, 1, m_domain_bcs_type, bcrs_v); @@ -199,10 +197,8 @@ void ERFPhysBCFunct_w::impose_vertical_zvel_bcs (const Array4& dest_arr, GpuArray,1> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) { - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp_w+i][ori]; - } + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp_w][ori]; } // ******************************************************* @@ -291,7 +287,6 @@ void ERFPhysBCFunct_w_no_terrain::impose_lateral_zvel_bcs (const Array4 bcrs_w(1); setBC(enclosedCells(bx), domain, bccomp, 0, 1, m_domain_bcs_type, bcrs_w); @@ -302,15 +297,15 @@ void ERFPhysBCFunct_w_no_terrain::impose_lateral_zvel_bcs (const Array4 bcrs_w_d(ncomp); + Gpu::DeviceVector bcrs_w_d(1); Gpu::copyAsync(Gpu::hostToDevice, bcrs_w.begin(), bcrs_w.end(), bcrs_w_d.begin()); const BCRec* bc_ptr_w = bcrs_w_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp+i][ori]; + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp][ori]; + } GeometryData const& geomdata = m_geom.data(); bool is_periodic_in_x = geomdata.isPeriodic(0); @@ -323,32 +318,32 @@ void ERFPhysBCFunct_w_no_terrain::impose_lateral_zvel_bcs (const Array4& // Based on BCRec for the domain, we need to make BCRec for this Box // bccomp is used as starting index for m_domain_bcs_type // 0 is used as starting index for bcrs - int ncomp = 1; Vector bcrs_w(1); setBC(enclosedCells(bx), domain, bccomp_w, 0, 1, m_domain_bcs_type, bcrs_w); // We use these for the asserts below const BCRec* bc_ptr_w_h = bcrs_w.data(); - GpuArray,AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM> l_bc_extdir_vals_d; - for (int i = 0; i < ncomp; i++) { - for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { - l_bc_extdir_vals_d[i][ori] = m_bc_extdir_vals[bccomp_w+i][ori]; - } + for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) { + l_bc_extdir_vals_d[0][ori] = m_bc_extdir_vals[bccomp_w][ori]; } // ******************************************************* diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.H b/Source/BoundaryConditions/ERF_PhysBCFunct.H index d42e2c493..b11ac0373 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.H +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.H @@ -26,8 +26,8 @@ public: ERFPhysBCFunct_cons (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_neumann_vals, std::unique_ptr& z_phys_nd, const bool use_real_bcs) : m_lev(lev), m_geom(geom), @@ -56,20 +56,20 @@ public: void impose_lateral_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, - int icomp, int ncomp, int bccomp, int ngz); + int icomp, int ncomp, int ngz); void impose_vertical_cons_bcs (const amrex::Array4& dest_arr, const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& z_nd, const amrex::GpuArray dxInv, - int icomp, int ncomp, int bccomp); + int icomp, int ncomp); private: int m_lev; amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; amrex::MultiFab* m_z_phys_nd; bool m_use_real_bcs; }; @@ -80,8 +80,8 @@ public: ERFPhysBCFunct_u (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_neumann_vals, std::unique_ptr& z_phys_nd, const bool use_real_bcs, amrex::Real* u_bc_data) @@ -123,8 +123,8 @@ private: amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; amrex::MultiFab* m_z_phys_nd; bool m_use_real_bcs; amrex::Real* m_u_bc_data; @@ -136,8 +136,8 @@ public: ERFPhysBCFunct_v (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_neumann_vals, std::unique_ptr& z_phys_nd, const bool use_real_bcs, amrex::Real* v_bc_data) @@ -180,8 +180,8 @@ private: amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; amrex::MultiFab* m_z_phys_nd; bool m_use_real_bcs; amrex::Real* m_v_bc_data; @@ -193,8 +193,8 @@ public: ERFPhysBCFunct_w (const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_neumann_vals, const TerrainType& terrain_type, std::unique_ptr& z_phys_nd, const bool use_real_bcs, amrex::Real* w_bc_data) @@ -245,8 +245,8 @@ private: amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; TerrainType m_terrain_type; amrex::MultiFab* m_z_phys_nd; bool m_use_real_bcs; @@ -260,8 +260,8 @@ public: const int lev, const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_extdir_vals, - amrex::Array,AMREX_SPACEDIM+NVAR_max> bc_neumann_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_extdir_vals, + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> bc_neumann_vals, const bool use_real_bcs, amrex::Real* w_bc_data) : m_lev(lev), @@ -300,8 +300,8 @@ private: amrex::Geometry m_geom; amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; - amrex::Array,AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; bool m_use_real_bcs; amrex::Real* m_w_bc_data; }; diff --git a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp index d4155577d..7e2509caa 100644 --- a/Source/BoundaryConditions/ERF_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ERF_PhysBCFunct.cpp @@ -14,7 +14,7 @@ using namespace amrex; */ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, - IntVect const& nghost, const Real /*time*/, int bccomp) + IntVect const& nghost, const Real /*time*/, int /*bccomp*/) { BL_PROFILE("ERFPhysBCFunct_cons::()"); @@ -80,10 +80,10 @@ void ERFPhysBCFunct_cons::operator() (MultiFab& mf, int icomp, int ncomp, if (!m_use_real_bcs) { - impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,bccomp,nghost[2]); + impose_lateral_cons_bcs(cons_arr,cbx1,domain,icomp,ncomp,nghost[2]); } - impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp,bccomp); + impose_vertical_cons_bcs(cons_arr,cbx2,domain,z_nd_arr,dxInv,icomp,ncomp); } } // MFIter diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 90fb1f990..aa13818f4 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -95,12 +95,11 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, const Box zbx = surroundingNodes(bx,2); const int end_comp = start_comp + num_comp - 1; - const int qty_offset = RhoTheta_comp; // Theta, KE, QKE, Scalar - Vector alpha_eff(NVAR_max, 0.0); + Vector alpha_eff(NPRIMVAR_max, 0.0); if (l_consA) { - for (int i = 0; i < NVAR_max; ++i) { + for (int i = 0; i < NPRIMVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.alpha_T; @@ -132,7 +131,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, } } } else { - for (int i = 0; i < NVAR_max; ++i) { + for (int i = 0; i < NPRIMVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T; @@ -198,20 +197,25 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, if (l_consA && l_turb) { ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; - - bool ext_dir_on_xlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir_prim)) + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; + const int prim_scal_index = prim_index; + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_xlo = ( ((bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim)) && i == dom_lo.x); - bool ext_dir_on_xhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_xhi = ( ((bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim)) && i == dom_hi.x+1); Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index]) - + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index]) + + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) ); if (ext_dir_on_xlo) { xflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i-1, j, k, prim_index) @@ -222,25 +226,31 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, - 3. * cell_prim(i-1, j, k, prim_index) + (1./3.) * cell_prim(i-2, j, k, prim_index) ) * dx_inv; } else { - xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0); + xflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - + cell_prim(i-1, j, k, prim_index)) * dx_inv * mf_u(i,j,0); } }); ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; - - bool ext_dir_on_ylo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir_prim)) + const int prim_index = qty_index - 1; + // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; + const int prim_scal_index = prim_index; + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_ylo = ( ((bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim)) && j == dom_lo.y); - bool ext_dir_on_yhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_yhi = ( ((bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim)) && j == dom_hi.y+1); Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index]) - + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index]) + + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) ); if (ext_dir_on_ylo) { yflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j-1, k, prim_index) @@ -257,18 +267,23 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; + // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; + const int prim_scal_index = prim_index; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && @@ -307,13 +322,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_xlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_xlo = ( ((bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim)) && i == dom_lo.x); - bool ext_dir_on_xhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_xhi = ( ((bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim)) && i == dom_hi.x+1); Real rhoAlpha = d_alpha_eff[prim_index]; @@ -335,13 +353,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_ylo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_ylo = ( ((bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim)) && j == dom_lo.y); - bool ext_dir_on_yhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_yhi = ( ((bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim)) && j == dom_hi.y+1); Real rhoAlpha = d_alpha_eff[prim_index]; @@ -363,20 +384,23 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == dom_lo.z); if (ext_dir_on_zlo) { @@ -412,13 +436,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_xlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_xlo = ( ((bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim)) && i == dom_lo.x); - bool ext_dir_on_xhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_xhi = ( ((bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim)) && i == dom_hi.x+1); Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ); @@ -439,13 +466,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_ylo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_ylo = ( ((bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim)) && j == dom_lo.y); - bool ext_dir_on_yhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_yhi = ( ((bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim)) && j == dom_hi.y+1); Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ); @@ -466,19 +496,22 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == dom_lo.z); if (ext_dir_on_zlo) { @@ -515,13 +548,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_xlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(0) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_xlo = ( ((bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(0) == ERFBCType::ext_dir_prim)) && i == dom_lo.x); - bool ext_dir_on_xhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(0) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_xhi = ( ((bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(0) == ERFBCType::ext_dir_prim)) && i == dom_hi.x+1); Real rhoAlpha = d_alpha_eff[prim_index]; @@ -541,13 +577,16 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; - bool ext_dir_on_ylo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(1) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_ylo = ( ((bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(1) == ERFBCType::ext_dir_prim)) && j == dom_lo.y); - bool ext_dir_on_yhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(1) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_yhi = ( ((bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(1) == ERFBCType::ext_dir_prim)) && j == dom_hi.y+1); Real rhoAlpha = d_alpha_eff[prim_index]; @@ -567,15 +606,18 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].hi(2) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 74cc3a7c5..9bb0e1eec 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -108,12 +108,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Box zbx3 = zbx; const int end_comp = start_comp + num_comp - 1; - const int qty_offset = RhoTheta_comp; // Theta, KE, QKE, Scalar - Vector alpha_eff(NVAR_max, 0.0); + Vector alpha_eff(NPRIMVAR_max, 0.0); if (l_consA) { - for (int i = 0; i < NVAR_max; ++i) { + for (int i = 0; i < NPRIMVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.alpha_T; @@ -145,7 +144,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, } } } else { - for (int i = 0; i < NVAR_max; ++i) { + for (int i = 0; i < NPRIMVAR_max; ++i) { switch (i) { case PrimTheta_comp: alpha_eff[PrimTheta_comp] = diffChoice.rhoAlpha_T; @@ -212,18 +211,23 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; + const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index]) - + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_index]) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_scal_index]) + + mu_turb(i-1, j, k, d_eddy_diff_idx[prim_scal_index]) ); Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ax(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -240,18 +244,22 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; + const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index]) - + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_index]) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_scal_index]) + + mu_turb(i, j-1, k, d_eddy_diff_idy[prim_scal_index]) ); Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ay(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -268,24 +276,28 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; + const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); Real met_h_zeta = az(i,j,k); Real GradCz; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir_prim) ) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim) ) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) @@ -313,7 +325,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; rhoAlpha += 0.5 * ( mu_turb(i , j, k, d_eddy_diff_idx[prim_index]) @@ -322,8 +334,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ax(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -340,7 +355,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; rhoAlpha += 0.5 * ( mu_turb(i, j , k, d_eddy_diff_idy[prim_index]) @@ -349,8 +364,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ay(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -367,7 +385,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; @@ -377,14 +395,17 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_zeta = az(i,j,k); Real GradCz; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) @@ -412,7 +433,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i-1, j, k, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; @@ -420,8 +441,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ax(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -438,7 +462,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j-1, k, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; @@ -446,8 +470,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); Real met_h_zeta = ay(i,j,k); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -464,7 +491,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; @@ -472,14 +499,17 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, Real met_h_zeta = az(i,j,k); Real GradCz; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) @@ -507,7 +537,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; @@ -515,8 +545,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, met_h_xi = Compute_h_xi_AtIface (i,j,k,cellSizeInv,z_nd); met_h_zeta = Compute_h_zeta_AtIface(i,j,k,cellSizeInv,z_nd); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i-1, j, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i-1, j, k-1, prim_index) ); @@ -533,7 +566,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(ybx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; @@ -541,8 +574,11 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, met_h_eta = Compute_h_eta_AtJface (i,j,k,cellSizeInv,z_nd); met_h_zeta = Compute_h_zeta_AtJface(i,j,k,cellSizeInv,z_nd); + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); bool most_on_zlo = ( use_most && rot_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); Real GradCz = 0.25 * dz_inv * ( cell_prim(i, j, k+1, prim_index) + cell_prim(i, j-1, k+1, prim_index) - cell_prim(i, j, k-1, prim_index) - cell_prim(i, j-1, k-1, prim_index) ); @@ -559,7 +595,7 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { const int qty_index = start_comp + n; - const int prim_index = qty_index - qty_offset; + const int prim_index = qty_index - 1; Real rhoAlpha = d_alpha_eff[prim_index]; @@ -567,14 +603,17 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); Real GradCz; - bool ext_dir_on_zlo = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::ext_dir_prim)) + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[BCVars::cons_bc+qty_index].lo(5) == ERFBCType::ext_dir_prim)) + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) && k == dom_hi.z+1); bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && k == 0); + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); if (ext_dir_on_zlo) { GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) diff --git a/Source/ERF.H b/Source/ERF.H index f21ec6e82..6ef021ef9 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -796,10 +796,10 @@ private: amrex::Array domain_bc_type; // These hold the Dirichlet values at walls which need them ... - amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals; + amrex::Array, AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals; // These hold the Neumann values at walls which need them ... - amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_neumann_vals; + amrex::Array, AMREX_SPACEDIM+NBCVAR_max> m_bc_neumann_vals; // These are the "physical" boundary condition types (e.g. "inflow") amrex::GpuArray phys_bc_type; diff --git a/Source/IO/ERF_ReadBndryPlanes.H b/Source/IO/ERF_ReadBndryPlanes.H index c7be2cbc7..a05f0c96d 100644 --- a/Source/IO/ERF_ReadBndryPlanes.H +++ b/Source/IO/ERF_ReadBndryPlanes.H @@ -27,11 +27,11 @@ public: void read_input_files (amrex::Real time, amrex::Real dt, - amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals); + amrex::Array, AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals); void read_file (int idx, amrex::Vector>& data_to_fill, - amrex::Array, AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals); + amrex::Array, AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals); // Return the pointer to PlaneVectors at time "time" amrex::Vector>& interp_in_time (const amrex::Real& time); diff --git a/Source/IO/ERF_ReadBndryPlanes.cpp b/Source/IO/ERF_ReadBndryPlanes.cpp index 0d87044da..0d73b7f1a 100644 --- a/Source/IO/ERF_ReadBndryPlanes.cpp +++ b/Source/IO/ERF_ReadBndryPlanes.cpp @@ -267,7 +267,7 @@ void ReadBndryPlanes::read_time_file () */ void ReadBndryPlanes::read_input_files (Real time, Real dt, - Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals) + Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals) { BL_PROFILE("ERF::ReadBndryPlanes::read_input_files"); @@ -341,7 +341,7 @@ void ReadBndryPlanes::read_input_files (Real time, */ void ReadBndryPlanes::read_file (const int idx, Vector>& data_to_fill, - Array,AMREX_SPACEDIM+NVAR_max> m_bc_extdir_vals) + Array,AMREX_SPACEDIM+NBCVAR_max> m_bc_extdir_vals) { const int t_step = m_in_timesteps[idx]; const std::string chkname1 = m_filename + Concatenate("/bndry_output", t_step); @@ -353,7 +353,7 @@ void ReadBndryPlanes::read_file (const int idx, BoxArray ba(domain); DistributionMapping dm{ba}; - GpuArray, AMREX_SPACEDIM+NVAR_max> l_bc_extdir_vals_d; + GpuArray, AMREX_SPACEDIM+NBCVAR_max> l_bc_extdir_vals_d; for (int i = 0; i < BCVars::NumTypes; i++) { diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 8f23788cc..0b972037a 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -23,6 +23,14 @@ // but not to allocate actual solution data #define NVAR_max (NDRY + NSCALARS + NMOIST_max) +// This is the number of components if we assume there is only one passive scalar +// We can use this when we assume all passive scalars have the same bc's +#define NBCVAR_max (NDRY + 1 + NMOIST_max) + +// This is the number of components if we assume there is only one passive scalar +// We can use this when we assume all passive scalars have the same diffusion coefficients +#define NPRIMVAR_max (NDRY + 1 + NMOIST_max) + // Cell-centered state variables #define Rho_comp 0 #define RhoTheta_comp (Rho_comp+1) @@ -51,6 +59,7 @@ #define PrimQ6_comp (RhoQ6_comp-1) // NOTE: we still use this indexing even if no moisture +// NOTE: We assume a single boundary condition for all passive scalars namespace BCVars { enum { cons_bc = 0, @@ -61,9 +70,9 @@ namespace BCVars { RhoScalar_bc_comp, RhoQ1_bc_comp, RhoQ2_bc_comp, - xvel_bc = NVAR_max, - yvel_bc = NVAR_max+1, - zvel_bc = NVAR_max+2, + xvel_bc = NBCVAR_max, + yvel_bc = NBCVAR_max+1, + zvel_bc = NBCVAR_max+2, NumTypes }; } diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index 0f72695ac..ad7cea2e4 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -275,8 +275,8 @@ void ERF::init_bcs () // // ***************************************************************************** { - domain_bcs_type.resize(AMREX_SPACEDIM+NVAR_max); - domain_bcs_type_d.resize(AMREX_SPACEDIM+NVAR_max); + domain_bcs_type.resize(AMREX_SPACEDIM+NBCVAR_max); + domain_bcs_type_d.resize(AMREX_SPACEDIM+NBCVAR_max); for (OrientationIter oit; oit; ++oit) { Orientation ori = oit(); @@ -394,6 +394,7 @@ void ERF::init_bcs () // // Here we translate the physical boundary conditions -- one type per face -- // into logical boundary conditions for each cell-centered variable + // NOTE: all "scalars" share the same type of boundary condition // // ***************************************************************************** { @@ -405,11 +406,11 @@ void ERF::init_bcs () if ( bct == ERF_BC::symmetry ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::reflect_even); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::reflect_even); } } @@ -417,11 +418,11 @@ void ERF::init_bcs () else if ( bct == ERF_BC::outflow ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); } } @@ -429,11 +430,11 @@ void ERF::init_bcs () else if ( bct == ERF_BC::ho_outflow ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::hoextrapcc); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::hoextrapcc); } } @@ -441,17 +442,17 @@ void ERF::init_bcs () else if ( bct == ERF_BC::open ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) + for (int i = 0; i < NBCVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::open); } else { - for (int i = 0; i < NVAR_max; i++) + for (int i = 0; i < NBCVAR_max; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::open); } } else if ( bct == ERF_BC::no_slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) { @@ -465,7 +466,7 @@ void ERF::init_bcs () domain_bcs_type[BCVars::RhoTheta_bc_comp].setLo(dir, ERFBCType::neumann); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); } if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) { @@ -483,7 +484,7 @@ void ERF::init_bcs () else if (bct == ERF_BC::slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) { @@ -496,7 +497,7 @@ void ERF::init_bcs () domain_bcs_type[BCVars::Rho_bc_comp].setLo(dir, ERFBCType::neumann); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::foextrap); } if (m_bc_extdir_vals[BCVars::RhoTheta_bc_comp][ori] > 0.) { @@ -513,7 +514,7 @@ void ERF::init_bcs () else if (bct == ERF_BC::inflow) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir); if (input_bndry_planes && dir < 2 && ( ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || @@ -527,7 +528,7 @@ void ERF::init_bcs () } } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir); if (input_bndry_planes && dir < 2 && ( ( (BCVars::cons_bc+i == BCVars::Rho_bc_comp) && m_r2d->ingested_density()) || @@ -546,11 +547,11 @@ void ERF::init_bcs () else if (bct == ERF_BC::periodic) { if (side == Orientation::low) { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::int_dir); } } else { - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::int_dir); } } @@ -558,7 +559,7 @@ void ERF::init_bcs () else if ( bct == ERF_BC::MOST ) { AMREX_ALWAYS_ASSERT(dir == 2 && side == Orientation::low); - for (int i = 0; i < NVAR_max; i++) { + for (int i = 0; i < NBCVAR_max; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::foextrap); } }