diff --git a/EBGodunov/hydro_ebgodunov.H b/EBGodunov/hydro_ebgodunov.H index 3ab42107c..4aadbfff8 100644 --- a/EBGodunov/hydro_ebgodunov.H +++ b/EBGodunov/hydro_ebgodunov.H @@ -23,7 +23,9 @@ namespace EBGodunov { amrex::Vector const& h_bcrec, amrex::BCRec const* d_bcrec, const amrex::Geometry& geom, - amrex::Real dt, amrex::MultiFab const* velocity_on_eb_inflow = nullptr); + amrex::Real dt, + amrex::MultiFab const* velocity_on_eb_inflow = nullptr, + amrex::iMultiFab const* BC_MF = nullptr); void ComputeAdvectiveVel (AMREX_D_DECL(amrex::Box const& xbx, @@ -41,7 +43,8 @@ namespace EBGodunov { amrex::Array4 const& vel, amrex::Array4 const& flag, const amrex::Box& domain, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr = {}); void ExtrapVelToFacesOnBox ( amrex::Box const& bx, int ncomp, AMREX_D_DECL(amrex::Box const& xbx, @@ -78,7 +81,8 @@ namespace EBGodunov { AMREX_D_DECL(amrex::Array4 const& fcx, amrex::Array4 const& fcy, amrex::Array4 const& fcz), - amrex::Real* p, amrex::Array4 const& velocity_on_eb_inflow); + amrex::Real* p, amrex::Array4 const& velocity_on_eb_inflow, + amrex::Array4 const& bc_arr = {}); void ComputeEdgeState ( amrex::Box const& bx, int ncomp, @@ -107,7 +111,8 @@ namespace EBGodunov { amrex::Array4 const& fcz), amrex::Array4 const& ccent_arr, bool is_velocity, - amrex::Array4 const& values_on_eb_inflow); + amrex::Array4 const& values_on_eb_inflow, + amrex::Array4 const& bc_arr = {}); } // namespace ebgodunov diff --git a/EBGodunov/hydro_ebgodunov_bcs_K.H b/EBGodunov/hydro_ebgodunov_bcs_K.H index 33669cc5a..64747eafb 100644 --- a/EBGodunov/hydro_ebgodunov_bcs_K.H +++ b/EBGodunov/hydro_ebgodunov_bcs_K.H @@ -63,6 +63,8 @@ void SetXBCs (const int i, const int j, const int k, const int n, // upwinded intermediate edge state in Imx (which supplies hi) // With GPU, I think it's undertermined which Imx(i+..) you'll get here, // could be from PLM, could be the upwinded edgestate. + // Question is whether every thread in the stream has to finish a kernel + // before the next kernel is launched. // lo = hi_arr(2*domlo-i ,j,k,n); // hi = lo_arr(2*domlo-i-1,j,k,n); } diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp index 44fd10db3..77e3c8d7d 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp @@ -36,7 +36,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Array4 const& fcy, Array4 const& ccent_arr, bool is_velocity, - Array4 const& values_on_eb_inflow) + Array4 const& values_on_eb_inflow, + Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); Box const& ybx = amrex::surroundingNodes(bx,1); @@ -106,7 +107,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real lo = Ipx(i-1,j,k,n); Real hi = Imx(i ,j,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); @@ -118,7 +119,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real lo = Ipy(i,j-1,k,n); Real hi = Imy(i,j ,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); @@ -136,7 +137,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, { if (apy(i,j,k) > 0.) { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; l_yzlo = ylo(i,j,k,n); @@ -218,7 +219,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, sth += (fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.; } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -251,7 +252,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, { if (apx(i,j,k) > 0.) { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; l_xzlo = xlo(i,j,k,n); @@ -334,7 +335,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, sth += (fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp index 9cf1f23f5..321f5371d 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp @@ -39,7 +39,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Array4 const& fcz, Array4 const& ccent_arr, bool is_velocity, - Array4 const& values_on_eb_inflow) + Array4 const& values_on_eb_inflow, + Array4 const& bc_arr) { // bx is the cell-centered box on which we want to compute the advective update @@ -135,7 +136,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real uad = u_mac(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); @@ -153,7 +154,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real vad = v_mac(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); @@ -171,7 +172,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real wad = w_mac(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); @@ -199,7 +200,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(zylo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zylo, l_zyhi; EBGodunovCornerCouple::EBGodunov_corner_couple_zy(l_zylo, l_zyhi, i, j, k, n, l_dt, dy, iconserv[n], @@ -216,7 +217,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(yzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; EBGodunovCornerCouple::EBGodunov_corner_couple_yz(l_yzlo, l_yzhi, i, j, k, n, l_dt, dz, iconserv[n], @@ -302,7 +303,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, sth += (fq) ? Real(0.5)*l_dt*fq(i ,j,k,n) : Real(0.0); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -336,7 +337,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(xzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; EBGodunovCornerCouple::EBGodunov_corner_couple_xz(l_xzlo, l_xzhi, i, j, k, n, l_dt, dz, iconserv[n], @@ -353,7 +354,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(zxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zxlo, l_zxhi; EBGodunovCornerCouple::EBGodunov_corner_couple_zx(l_zxlo, l_zxhi, i, j, k, n, l_dt, dx, iconserv[n], @@ -439,7 +440,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, sth += (fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0.0); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -471,7 +472,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(xylo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xylo, l_xyhi; EBGodunovCornerCouple::EBGodunov_corner_couple_xy(l_xylo, l_xyhi, i, j, k, n, l_dt, dy, iconserv[n], @@ -488,7 +489,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Box(yxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yxlo, l_yxhi; EBGodunovCornerCouple::EBGodunov_corner_couple_yx(l_yxlo, l_yxhi, i, j, k, n, l_dt, dx, iconserv[n], @@ -573,7 +574,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, sth += (fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0.0); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp index 562c22e4c..2c0172407 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp @@ -23,7 +23,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, Vector const& h_bcrec, BCRec const* d_bcrec, const Geometry& geom, - Real l_dt, MultiFab const* velocity_on_eb_inflow) + Real l_dt, MultiFab const* velocity_on_eb_inflow, + iMultiFab const* BC_MF) { BL_PROFILE("EBGodunov::ExtrapVelToFaces()"); AMREX_ALWAYS_ASSERT(vel.hasEBFabFactory()); @@ -64,6 +65,9 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, Array4 const& a_vel = vel.const_array(mfi); Array4 const& a_f = vel_forces.const_array(mfi); + Array4 const& bc_arr = BC_MF ? BC_MF->const_array(mfi) + : Array4 {}; + // In 2-d: // 8*ncomp are: Imx, Ipx, Imy, Ipy, xlo/xhi, ylo/yhi // 2 are: u_ad, v_ad @@ -142,14 +146,14 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, #endif PLM::PredictVelOnXFace( xebx_g1, AMREX_SPACEDIM, Imx, Ipx, a_vel, a_vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); PLM::PredictVelOnYFace( yebx_g1, AMREX_SPACEDIM, Imy, Ipy, a_vel, a_vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); #if ( AMREX_SPACEDIM == 3 ) PLM::PredictVelOnZFace( zebx_g1, AMREX_SPACEDIM, Imz, Ipz, a_vel, a_vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); #endif @@ -159,7 +163,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, AMREX_D_DECL(Imx, Imy, Imz), AMREX_D_DECL(Ipx, Ipy, Ipz), a_vel, a_f, domain, l_dt, d_bcrec, - local_use_forces_in_trans); + local_use_forces_in_trans, + bc_arr); Godunov::ExtrapVelToFacesOnBox( bx, ncomp, AMREX_D_DECL(xbx, ybx, zbx), @@ -169,7 +174,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, AMREX_D_DECL(Imx, Imy, Imz), AMREX_D_DECL(Ipx, Ipy, Ipz), a_f, domain, dx, l_dt, d_bcrec, - local_use_forces_in_trans, p); + local_use_forces_in_trans, p, + bc_arr); } else { @@ -184,25 +190,25 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, EBPLM::PredictVelOnXFace( xebx_g2, Imx, Ipx, a_vel, a_vel, flagarr, vfrac_arr, AMREX_D_DECL(fcx,fcy,fcz),ccent_arr, - geom, l_dt, h_bcrec, d_bcrec ); + geom, l_dt, h_bcrec, d_bcrec, bc_arr ); EBPLM::PredictVelOnYFace( yebx_g2, Imy, Ipy, a_vel, a_vel, flagarr, vfrac_arr, AMREX_D_DECL(fcx,fcy,fcz),ccent_arr, - geom, l_dt, h_bcrec, d_bcrec ); + geom, l_dt, h_bcrec, d_bcrec, bc_arr ); #if (AMREX_SPACEDIM == 3) EBPLM::PredictVelOnZFace( zebx_g2, Imz, Ipz, a_vel, a_vel, flagarr, vfrac_arr, AMREX_D_DECL(fcx,fcy,fcz),ccent_arr, - geom, l_dt, h_bcrec, d_bcrec ); + geom, l_dt, h_bcrec, d_bcrec, bc_arr ); #endif EBGodunov::ComputeAdvectiveVel( AMREX_D_DECL(xebx_g2, yebx_g2, zebx_g2), AMREX_D_DECL(u_ad, v_ad, w_ad), AMREX_D_DECL(Imx, Imy, Imz), AMREX_D_DECL(Ipx, Ipy, Ipz), - a_vel, flagarr, domain, d_bcrec); + a_vel, flagarr, domain, d_bcrec, bc_arr); AMREX_D_TERM(Array4 const& apx = areafrac[0]->const_array(mfi);, Array4 const& apy = areafrac[1]->const_array(mfi);, @@ -226,7 +232,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, AMREX_D_DECL(fcx, fcy, fcz), p, velocity_on_eb_inflow ? - velocity_on_eb_inflow->const_array(mfi) : Array4{}); + velocity_on_eb_inflow->const_array(mfi) : Array4{}, + bc_arr); } Gpu::streamSynchronize(); // otherwise we might be using too much memory @@ -250,7 +257,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Array4 const& vel, Array4 const& flag, const Box& domain, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); const Dim3 dhi = amrex::ubound(domain); @@ -266,7 +274,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real lo = Ipx(i-1,j,k,n); Real hi = Imx(i ,j,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetXBCs(i, j, k, n, vel, lo, hi, bc.lo(0), bc.hi(0), dlo.x, dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -286,7 +294,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real lo = Ipy(i,j-1,k,n); Real hi = Imy(i,j ,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetYBCs(i, j, k, n, vel, lo, hi, bc.lo(1), bc.hi(1), dlo.y, dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -307,7 +315,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real lo = Ipz(i,j,k-1,n); Real hi = Imz(i,j,k ,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetZBCs(i, j, k, n, vel, lo, hi, bc.lo(2), bc.hi(2), dlo.z, dhi.z, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp index 9363adda3..2d97a4a54 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp @@ -34,7 +34,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, Array4 const& fcx, Array4 const& fcy, Real* p, - Array4 const& velocity_on_eb_inflow) + Array4 const& velocity_on_eb_inflow, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); const Dim3 dhi = amrex::ubound(domain); @@ -60,7 +61,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, Real lo = Ipx(i-1,j,k,n); Real hi = Imx(i ,j,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); @@ -72,7 +73,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, Real lo = Ipy(i,j-1,k,n); Real hi = Imy(i,j ,k,n); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); @@ -96,7 +97,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, if (flag(i,j,k).isConnected(0,-1,0)) { constexpr int n = 0; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; l_yzlo = ylo(i,j,k,n); @@ -121,7 +122,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, if (flag(i,j,k).isConnected(-1,0,0)) { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); // stl is on the left side of the lo-x side of cell (i,j) // sth is on the right side of the lo-x side of cell (i,j) @@ -221,7 +222,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, if (flag(i,j,k).isConnected(-1,0,0)) { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; l_xzlo = xlo(i,j,k,n); @@ -246,7 +247,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, if (flag(i,j,k).isConnected(0,-1,0)) { constexpr int n = 1; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); // stl is on the low side of the lo-y side of cell (i,j) // sth is on the high side of the lo-y side of cell (i,j) diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp index afd013314..ae03fd511 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp @@ -47,7 +47,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Array4 const& fcy, Array4 const& fcz, Real* p, - Array4 const& velocity_on_eb_inflow) + Array4 const& velocity_on_eb_inflow, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); const Dim3 dhi = amrex::ubound(domain); @@ -74,7 +75,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real hi = Imx(i ,j,k,n); Real uad = u_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetXBCs(i, j, k, n, q, lo, hi, bc.lo(0), bc.hi(0), dlo.x, dhi.x, true); @@ -91,7 +92,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real hi = Imy(i,j ,k,n); Real vad = v_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetYBCs(i, j, k, n, q, lo, hi, bc.lo(1), bc.hi(1), dlo.y, dhi.y, true); @@ -108,7 +109,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real hi = Imz(i,j,k ,n); Real wad = w_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); EBGodunovBC::SetZBCs(i, j, k, n, q, lo, hi, bc.lo(2), bc.hi(2), dlo.z, dhi.z, true); @@ -155,7 +156,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zylo, l_zyhi; EBGodunovCornerCouple::EBGodunov_corner_couple_zy(l_zylo, l_zyhi, i, j, k, n, l_dt, dy, false, @@ -173,7 +174,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; EBGodunovCornerCouple::EBGodunov_corner_couple_yz(l_yzlo, l_yzhi, i, j, k, n, l_dt, dz, false, @@ -193,7 +194,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, if (flag(i,j,k).isConnected(-1,0,0)) { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); // stl is on the lo side of the lo-x side of cell (i,j,k) // sth is on the hi side of the lo-x side of cell (i,j,k) @@ -306,7 +307,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; EBGodunovCornerCouple::EBGodunov_corner_couple_xz(l_xzlo, l_xzhi, i, j, k, n, l_dt, dz, false, @@ -324,7 +325,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zxlo, l_zxhi; EBGodunovCornerCouple::EBGodunov_corner_couple_zx(l_zxlo, l_zxhi, i, j, k, n, l_dt, dx, false, @@ -345,7 +346,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, if (flag(i,j,k).isConnected(0,-1,0)) { constexpr int n = 1; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); // stl is on the lo side of the lo-y side of cell (i,j,k) // sth is on the hi side of the lo-y side of cell (i,j,k) @@ -458,7 +459,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xylo, l_xyhi; EBGodunovCornerCouple::EBGodunov_corner_couple_xy(l_xylo, l_xyhi, i, j, k, n, l_dt, dy, false, @@ -480,7 +481,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yxlo, l_yxhi; EBGodunovCornerCouple::EBGodunov_corner_couple_yx(l_yxlo, l_yxhi, i, j, k, n, l_dt, dx, false, @@ -501,7 +502,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, if (flag(i,j,k).isConnected(0,0,-1)) { constexpr int n = 2; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); // stl is on the lo side of the lo-z side of cell (i,j,k) // sth is on the hi side of the lo-z side of cell (i,j,k) diff --git a/EBGodunov/hydro_ebgodunov_plm.H b/EBGodunov/hydro_ebgodunov_plm.H index 4b65a03f5..1b0192cde 100644 --- a/EBGodunov/hydro_ebgodunov_plm.H +++ b/EBGodunov/hydro_ebgodunov_plm.H @@ -27,7 +27,8 @@ void PredictVelOnXFace ( amrex::Box const& xebox, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr = {}); void PredictVelOnYFace ( amrex::Box const& yebox, amrex::Array4 const& Imy, @@ -43,7 +44,8 @@ void PredictVelOnYFace ( amrex::Box const& yebox, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr = {}); #if (AMREX_SPACEDIM==3) void PredictVelOnZFace ( amrex::Box const& zebox, @@ -60,7 +62,8 @@ void PredictVelOnZFace ( amrex::Box const& zebox, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr = {}); #endif @@ -78,7 +81,8 @@ void PredictStateOnXFace ( amrex::Box const& xebox, int ncomp, amrex::Geometry const& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc, bool is_velocity); + amrex::BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr = {}); void PredictStateOnYFace ( amrex::Box const& yebox, int ncomp, amrex::Array4 const& Imy, amrex::Array4 const& Ipy, @@ -93,7 +97,8 @@ void PredictStateOnYFace ( amrex::Box const& yebox, int ncomp, amrex::Geometry const& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc, bool is_velocity); + amrex::BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr = {}); #if (AMREX_SPACEDIM == 3) void PredictStateOnZFace ( amrex::Box const& zebox, int ncomp, @@ -109,7 +114,8 @@ void PredictStateOnZFace ( amrex::Box const& zebox, int ncomp, amrex::Geometry const& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc, bool is_velocity); + amrex::BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr = {}); #endif } diff --git a/EBGodunov/hydro_ebgodunov_plm.cpp b/EBGodunov/hydro_ebgodunov_plm.cpp index 0b37a5a3e..51a877ebf 100644 --- a/EBGodunov/hydro_ebgodunov_plm.cpp +++ b/EBGodunov/hydro_ebgodunov_plm.cpp @@ -7,6 +7,7 @@ #include #include +#include using namespace amrex; @@ -40,7 +41,8 @@ EBPLM::PredictVelOnXFace (Box const& xebox, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dx = geom.CellSize(0); const Real dtdx = dt/dx; @@ -57,17 +59,19 @@ EBPLM::PredictVelOnXFace (Box const& xebox, auto extdir_lohi_x = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::x)); auto extdir_lohi_y = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::y)); - bool has_extdir_or_ho_lo_x = extdir_lohi_x.first; - bool has_extdir_or_ho_hi_x = extdir_lohi_x.second; - bool has_extdir_or_ho_lo_y = extdir_lohi_y.first; - bool has_extdir_or_ho_hi_y = extdir_lohi_y.second; +// If we have a BC Array4, then we can't make these blanket statements about the domain face +// Best to assume we could have face-based BC somewhere + bool has_extdir_or_ho_lo_x = bc_arr ? true : extdir_lohi_x.first; + bool has_extdir_or_ho_hi_x = bc_arr ? true : extdir_lohi_x.second; + bool has_extdir_or_ho_lo_y = bc_arr ? true : extdir_lohi_y.first; + bool has_extdir_or_ho_hi_y = bc_arr ? true : extdir_lohi_y.second; #if (AMREX_SPACEDIM == 3) const int domain_klo = domain_box.smallEnd(2); const int domain_khi = domain_box.bigEnd(2); auto extdir_lohi_z = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::z)); - bool has_extdir_or_ho_lo_z = extdir_lohi_z.first; - bool has_extdir_or_ho_hi_z = extdir_lohi_z.second; + bool has_extdir_or_ho_lo_z = bc_arr ? true : extdir_lohi_z.first; + bool has_extdir_or_ho_hi_z = bc_arr ? true : extdir_lohi_z.second; #endif if ( (has_extdir_or_ho_lo_x && domain_ilo >= xebox.smallEnd(0)-1) || @@ -79,9 +83,7 @@ EBPLM::PredictVelOnXFace (Box const& xebox, (has_extdir_or_ho_lo_y && domain_jlo >= xebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= xebox.bigEnd(1) ) ) { - amrex::ParallelFor(xebox, ncomp, [q,ccvel,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imx,Ipx,dtdx,pbc,flag,ccc,vfrac,AMREX_D_DECL(fcx,fcy,fcz)] + amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -90,7 +92,9 @@ EBPLM::PredictVelOnXFace (Box const& xebox, // This means apx(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(-1,0,0)) { - const auto& bc = pbc[n]; + //const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); + bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || @@ -366,7 +370,8 @@ EBPLM::PredictVelOnYFace (Box const& yebox, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dy = geom.CellSize(1); const Real dtdy = dt/dy; @@ -382,17 +387,17 @@ EBPLM::PredictVelOnYFace (Box const& yebox, auto extdir_lohi_x = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::x)); auto extdir_lohi_y = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::y)); - bool has_extdir_or_ho_lo_x = extdir_lohi_x.first; - bool has_extdir_or_ho_hi_x = extdir_lohi_x.second; - bool has_extdir_or_ho_lo_y = extdir_lohi_y.first; - bool has_extdir_or_ho_hi_y = extdir_lohi_y.second; + bool has_extdir_or_ho_lo_x = bc_arr ? true : extdir_lohi_x.first; + bool has_extdir_or_ho_hi_x = bc_arr ? true : extdir_lohi_x.second; + bool has_extdir_or_ho_lo_y = bc_arr ? true : extdir_lohi_y.first; + bool has_extdir_or_ho_hi_y = bc_arr ? true : extdir_lohi_y.second; #if (AMREX_SPACEDIM == 3) const int domain_klo = domain_box.smallEnd(2); const int domain_khi = domain_box.bigEnd(2); auto extdir_lohi_z = has_extdir_or_ho(h_bcrec.data(), AMREX_SPACEDIM, static_cast(Direction::z)); - bool has_extdir_or_ho_lo_z = extdir_lohi_z.first; - bool has_extdir_or_ho_hi_z = extdir_lohi_z.second; + bool has_extdir_or_ho_lo_z = bc_arr ? true : extdir_lohi_z.first; + bool has_extdir_or_ho_hi_z = bc_arr ? true : extdir_lohi_z.second; #endif if ( (has_extdir_or_ho_lo_x && domain_ilo >= yebox.smallEnd(0)-1) || @@ -404,9 +409,7 @@ EBPLM::PredictVelOnYFace (Box const& yebox, (has_extdir_or_ho_lo_y && domain_jlo >= yebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= yebox.bigEnd(1) ) ) { - amrex::ParallelFor(yebox, ncomp, [q,ccvel,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imy,Ipy,dtdy,pbc,flag,vfrac,ccc,AMREX_D_DECL(fcx,fcy,fcz)] + amrex::ParallelFor(yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -415,7 +418,9 @@ EBPLM::PredictVelOnYFace (Box const& yebox, // This means apy(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(0,-1,0)) { - const auto& bc = pbc[n]; + //const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); + bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || @@ -692,7 +697,8 @@ EBPLM::PredictVelOnZFace (Box const& zebox, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dz = geom.CellSize(2); const Real dtdz = dt/dz; @@ -726,9 +732,7 @@ EBPLM::PredictVelOnZFace (Box const& zebox, (has_extdir_or_ho_lo_y && domain_jlo >= zebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= zebox.bigEnd(1) ) ) { - amrex::ParallelFor(zebox, ncomp, [q,ccvel,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imz,Ipz,dtdz,pbc,flag,vfrac,ccc,AMREX_D_DECL(fcx,fcy,fcz)] + amrex::ParallelFor(zebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -737,7 +741,9 @@ EBPLM::PredictVelOnZFace (Box const& zebox, // This means apz(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(0,0,-1)) { - const auto& bc = pbc[n]; + //const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); + bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || diff --git a/EBGodunov/hydro_ebgodunov_plm_fpu.cpp b/EBGodunov/hydro_ebgodunov_plm_fpu.cpp index 0d2c827d9..bffd5493c 100644 --- a/EBGodunov/hydro_ebgodunov_plm_fpu.cpp +++ b/EBGodunov/hydro_ebgodunov_plm_fpu.cpp @@ -7,6 +7,7 @@ #include #include +#include using namespace amrex; @@ -41,7 +42,8 @@ EBPLM::PredictStateOnXFace (Box const& xebox, int ncomp, Geometry const& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc, bool is_velocity) + BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr) { const Real dx = geom.CellSize(0); const Real dtdx = dt/dx; @@ -78,10 +80,7 @@ EBPLM::PredictStateOnXFace (Box const& xebox, int ncomp, (has_extdir_or_ho_lo_y && domain_jlo >= xebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= xebox.bigEnd(1) ) ) { - amrex::ParallelFor(xebox, ncomp, [q,umac,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imx,Ipx,dtdx,pbc,flag,vfrac,ccc,AMREX_D_DECL(fcx,fcy,fcz), - is_velocity] + amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -90,7 +89,7 @@ EBPLM::PredictStateOnXFace (Box const& xebox, int ncomp, // This means apx(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(-1,0,0)) { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || @@ -376,7 +375,8 @@ EBPLM::PredictStateOnYFace ( Box const& yebox, int ncomp, Geometry const& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc, bool is_velocity) + BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr) { const Real dy = geom.CellSize(1); const Real dtdy = dt/dy; @@ -413,10 +413,7 @@ EBPLM::PredictStateOnYFace ( Box const& yebox, int ncomp, (has_extdir_or_ho_lo_y && domain_jlo >= yebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= yebox.bigEnd(1) ) ) { - amrex::ParallelFor(yebox, ncomp, [q,vmac,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imy,Ipy,dtdy,pbc,flag,vfrac,ccc,AMREX_D_DECL(fcx,fcy,fcz), - is_velocity] + amrex::ParallelFor(yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -425,7 +422,7 @@ EBPLM::PredictStateOnYFace ( Box const& yebox, int ncomp, // This means apy(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(0,-1,0)) { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || @@ -715,7 +712,8 @@ EBPLM::PredictStateOnZFace ( Box const& zebox, int ncomp, Geometry const& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc, bool is_velocity) + BCRec const* pbc, bool is_velocity, + amrex::Array4 const& bc_arr) { const Real dz = geom.CellSize(1); const Real dtdz = dt/dz; @@ -747,10 +745,7 @@ EBPLM::PredictStateOnZFace ( Box const& zebox, int ncomp, (has_extdir_or_ho_lo_y && domain_jlo >= zebox.smallEnd(1)-1) || (has_extdir_or_ho_hi_y && domain_jhi <= zebox.bigEnd(1) ) ) { - amrex::ParallelFor(zebox, ncomp, [q,wmac,AMREX_D_DECL(domain_ilo,domain_jlo,domain_klo), - AMREX_D_DECL(domain_ihi,domain_jhi,domain_khi), - Imz,Ipz,dtdz,pbc,flag,vfrac,ccc,AMREX_D_DECL(fcx,fcy,fcz), - is_velocity] + amrex::ParallelFor(zebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real qpls(0.); @@ -759,7 +754,7 @@ EBPLM::PredictStateOnZFace ( Box const& zebox, int ncomp, // This means apz(i,j,k) > 0 and we have un-covered cells on both sides if (flag(i,j,k).isConnected(0,0,-1)) { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || diff --git a/Godunov/hydro_godunov.H b/Godunov/hydro_godunov.H index 3a22bf830..d9e88c687 100644 --- a/Godunov/hydro_godunov.H +++ b/Godunov/hydro_godunov.H @@ -24,7 +24,8 @@ void ExtrapVelToFaces ( amrex::MultiFab const& a_vel, const amrex::BCRec * d_bcrec, const amrex::Geometry& geom, amrex::Real l_dt, bool use_ppm, bool use_forces_in_trans, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::VanLeer, + amrex::iMultiFab const* BC_MF = nullptr); void ComputeAdvectiveVel (AMREX_D_DECL(amrex::Box const& xbx, amrex::Box const& ybx, @@ -43,7 +44,8 @@ void ComputeAdvectiveVel (AMREX_D_DECL(amrex::Box const& xbx, const amrex::Box& domain, amrex::Real dt, amrex::BCRec const* pbc, - bool use_forces_in_trans); + bool use_forces_in_trans, + amrex::Array4 const& bc_arr = {}); void ExtrapVelToFacesOnBox (amrex::Box const& bx, int ncomp, AMREX_D_DECL(amrex::Box const& xbx, @@ -68,7 +70,8 @@ void ExtrapVelToFacesOnBox (amrex::Box const& bx, int ncomp, amrex::Real dt, amrex::BCRec const* pbc, bool use_forces_in_trans, - amrex::Real* p); + amrex::Real* p, + amrex::Array4 const& bc_arr = {}); void ComputeEdgeState ( amrex::Box const& bx, int ncomp, amrex::Array4 const& q, @@ -86,7 +89,8 @@ void ComputeEdgeState ( amrex::Box const& bx, int ncomp, int const* iconserv, bool use_ppm, bool use_forces_in_trans, bool is_velocity, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::VanLeer, + amrex::Array4 const& bc_arr = {}); } diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp index 172170268..9c90dc77a 100644 --- a/Godunov/hydro_godunov_edge_state_2D.cpp +++ b/Godunov/hydro_godunov_edge_state_2D.cpp @@ -28,7 +28,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, const bool use_ppm, const bool use_forces_in_trans, const bool is_velocity, - const int limiter_type) + const int limiter_type, + amrex::Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); Box const& ybx = amrex::surroundingNodes(bx,1); @@ -100,15 +101,17 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); PLM::PredictStateOnXFace(i, j, k, n, l_dt, dx, Imx(i,j,k,n), Ipx(i-1,j,k,n), - q, umac(i,j,k), pbc[n], dlo.x, dhi.x, is_velocity); + q, umac(i,j,k), bc, dlo.x, dhi.x, is_velocity); }); amrex::ParallelFor(yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); PLM::PredictStateOnYFace(i, j, k, n, l_dt, dy, Imy(i,j,k,n), Ipy(i,j-1,k,n), - q, vmac(i,j,k), pbc[n], dlo.y, dhi.y, is_velocity); + q, vmac(i,j,k), bc, dlo.y, dhi.y, is_velocity); }); } @@ -128,7 +131,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, hi += 0.5*l_dt*fq(i ,j,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; @@ -145,7 +148,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, hi += 0.5*l_dt*fq(i,j ,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); @@ -163,7 +166,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(yzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; l_yzlo = ylo(i,j,k,n); @@ -217,7 +220,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, sth += (is_rz) ? -0.25 * l_dt * q(i,j,k,n)*( umac(i,j,k) + umac(i+1,j,k) ) / ( dx*(amrex::Math::abs(Real(i)+0.5)) ) : 0.; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -245,7 +248,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(xzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; l_xzlo = xlo(i,j,k,n); @@ -301,7 +304,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, sth += (is_rz) ? -0.25 * l_dt * q(i,j,k,n)*( umac(i,j ,k) + umac(i+1,j ,k) ) / ( dx*(amrex::Math::abs(Real(i)+0.5)) ) : 0.; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp index ea6076f46..78f08a1ef 100644 --- a/Godunov/hydro_godunov_edge_state_3D.cpp +++ b/Godunov/hydro_godunov_edge_state_3D.cpp @@ -30,7 +30,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, const bool use_ppm, const bool use_forces_in_trans, const bool is_velocity, - const int limiter_type) + const int limiter_type, + amrex::Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); Box const& ybx = amrex::surroundingNodes(bx,1); @@ -115,21 +116,24 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); PLM::PredictStateOnXFace(i, j, k, n, l_dt, dx, Imx(i,j,k,n), Ipx(i-1,j,k,n), - q, umac(i,j,k), pbc[n], dlo.x, dhi.x, is_velocity); + q, umac(i,j,k), bc, dlo.x, dhi.x, is_velocity); }); amrex::ParallelFor(yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); PLM::PredictStateOnYFace(i, j, k, n, l_dt, dy, Imy(i,j,k,n), Ipy(i,j-1,k,n), - q, vmac(i,j,k), pbc[n], dlo.y, dhi.y, is_velocity); + q, vmac(i,j,k), bc, dlo.y, dhi.y, is_velocity); }); amrex::ParallelFor(zebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); PLM::PredictStateOnZFace(i, j, k, n, l_dt, dz, Imz(i,j,k,n), Ipz(i,j,k-1,n), - q, wmac(i,j,k), pbc[n], dlo.z, dhi.z, is_velocity); + q, wmac(i,j,k), bc, dlo.z, dhi.z, is_velocity); }); } @@ -149,7 +153,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, hi += 0.5*l_dt*fq(i ,j,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; @@ -172,7 +176,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, hi += 0.5*l_dt*fq(i,j ,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); @@ -195,7 +199,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, hi += 0.5*l_dt*fq(i,j,k ,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); @@ -218,7 +222,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(zylo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zylo, l_zyhi; GodunovCornerCouple::AddCornerCoupleTermZY(l_zylo, l_zyhi, i, j, k, n, l_dt, dy, iconserv[n], @@ -235,7 +239,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(yzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; GodunovCornerCouple::AddCornerCoupleTermYZ(l_yzlo, l_yzhi, i, j, k, n, l_dt, dz, iconserv[n], @@ -255,43 +259,43 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, amrex::ParallelFor(xbx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - Real stl = xlo(i,j,k,n); - Real sth = xhi(i,j,k,n); + Real stl = xlo(i,j,k,n); + Real sth = xhi(i,j,k,n); // To match EBGodunov // Here we add dt/2 (-q u_x - (v q)_y - (w q)_z) to the term that is already - // q + dx/2 q_x + dt/2 (-u q_x) to get - // q + dx/2 q_x - dt/2 (u q_x + q u_x + (v q)_y + (w q)_z) which is equivalent to - // --> q + dx/2 q_x - dt/2 ( div (uvec q) ) - Real quxl = (umac(i,j,k) - umac(i-1,j,k)) * q(i-1,j,k,n); - stl += ( - (0.5*dtdx) * quxl - - (0.5*dtdy)*(yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) - - yzlo(i-1,j ,k ,n)*vmac(i-1,j ,k )) - - (0.5*dtdz)*(zylo(i-1,j ,k+1,n)*wmac(i-1,j ,k+1) - - zylo(i-1,j ,k ,n)*wmac(i-1,j ,k )) ); + // q + dx/2 q_x + dt/2 (-u q_x) to get + // q + dx/2 q_x - dt/2 (u q_x + q u_x + (v q)_y + (w q)_z) which is equivalent to + // --> q + dx/2 q_x - dt/2 ( div (uvec q) ) + Real quxl = (umac(i,j,k) - umac(i-1,j,k)) * q(i-1,j,k,n); + stl += ( - (0.5*dtdx) * quxl + - (0.5*dtdy)*(yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) + - yzlo(i-1,j ,k ,n)*vmac(i-1,j ,k )) + - (0.5*dtdz)*(zylo(i-1,j ,k+1,n)*wmac(i-1,j ,k+1) + - zylo(i-1,j ,k ,n)*wmac(i-1,j ,k )) ); - // Here we adjust for non-conservative by removing the q divu contribution to get - // q + dx/2 q_x - dt/2 ( div (uvec q) - q divu ) which is equivalent to - // --> q + dx/2 q_x - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : 0.; + // Here we adjust for non-conservative by removing the q divu contribution to get + // q + dx/2 q_x - dt/2 ( div (uvec q) - q divu ) which is equivalent to + // --> q + dx/2 q_x - dt/2 ( uvec dot grad q) + stl += (!iconserv[n]) ? 0.5*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : 0.; - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i-1,j,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i-1,j,k,n) : 0.; - // High side - Real quxh = (umac(i+1,j,k) - umac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdx) * quxh - - (0.5*dtdy)*(yzlo(i,j+1,k ,n)*vmac(i,j+1,k ) - - yzlo(i,j ,k ,n)*vmac(i,j ,k )) - - (0.5*dtdz)*(zylo(i,j ,k+1,n)*wmac(i,j ,k+1) - - zylo(i,j ,k ,n)*wmac(i,j ,k )) ); + // High side + Real quxh = (umac(i+1,j,k) - umac(i,j,k)) * q(i,j,k,n); + sth += ( - (0.5*dtdx) * quxh + - (0.5*dtdy)*(yzlo(i,j+1,k ,n)*vmac(i,j+1,k ) + - yzlo(i,j ,k ,n)*vmac(i,j ,k )) + - (0.5*dtdz)*(zylo(i,j ,k+1,n)*wmac(i,j ,k+1) + - zylo(i,j ,k ,n)*wmac(i,j ,k )) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i ,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? 0.5*l_dt* q(i ,j,k,n)*divu(i,j,k) : 0.; - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.; - auto bc = pbc[n]; - HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) { @@ -319,7 +323,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(xzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; GodunovCornerCouple::AddCornerCoupleTermXZ(l_xzlo, l_xzhi, i, j, k, n, l_dt, dz, iconserv[n], @@ -336,7 +340,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(zxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zxlo, l_zxhi; GodunovCornerCouple::AddCornerCoupleTermZX(l_zxlo, l_zxhi, i, j, k, n, l_dt, dx, iconserv[n], @@ -355,41 +359,41 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, amrex::ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - Real stl = ylo(i,j,k,n); - Real sth = yhi(i,j,k,n); + Real stl = ylo(i,j,k,n); + Real sth = yhi(i,j,k,n); - // Here we add dt/2 (-q v_y - (u q)_x - (w q)_z) to the term that is already - // q + dy/2 q_y + dt/2 (-v q_y) to get - // q + dy/2 q_y - dt/2 (v q_y + q v_y + (u q)_x + (w q)_z) which is equivalent to - // --> q + dy/2 q_y - dt/2 ( div (uvec q) ) - Real qvyl = (vmac(i,j,k) - vmac(i,j-1,k)) * q(i,j-1,k,n); - stl += ( - (0.5*dtdy) * qvyl - - (0.5*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) - - xzlo(i ,j-1,k ,n)*umac(i ,j-1,k )) - - (0.5*dtdz)*(zxlo(i ,j-1,k+1,n)*wmac(i ,j-1,k+1) - - zxlo(i ,j-1,k ,n)*wmac(i ,j-1,k )) ); + // Here we add dt/2 (-q v_y - (u q)_x - (w q)_z) to the term that is already + // q + dy/2 q_y + dt/2 (-v q_y) to get + // q + dy/2 q_y - dt/2 (v q_y + q v_y + (u q)_x + (w q)_z) which is equivalent to + // --> q + dy/2 q_y - dt/2 ( div (uvec q) ) + Real qvyl = (vmac(i,j,k) - vmac(i,j-1,k)) * q(i,j-1,k,n); + stl += ( - (0.5*dtdy) * qvyl + - (0.5*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) + - xzlo(i ,j-1,k ,n)*umac(i ,j-1,k )) + - (0.5*dtdz)*(zxlo(i ,j-1,k+1,n)*wmac(i ,j-1,k+1) + - zxlo(i ,j-1,k ,n)*wmac(i ,j-1,k )) ); - // Here we adjust for non-conservative by removing the q divu contribution to get - // q + dy/2 q_y - dt/2 ( div (uvec q) - q divu ) which is equivalent to - // --> q + dy/2 q_y - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : 0.; + // Here we adjust for non-conservative by removing the q divu contribution to get + // q + dy/2 q_y - dt/2 ( div (uvec q) - q divu ) which is equivalent to + // --> q + dy/2 q_y - dt/2 ( uvec dot grad q) + stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : 0.; - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j-1,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j-1,k,n) : 0.; - // High side - Real qvyh = (vmac(i,j+1,k) - vmac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdy) * qvyh - - (0.5*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) - - xzlo(i ,j,k ,n)*umac(i ,j,k )) - - (0.5*dtdz)*(zxlo(i ,j,k+1,n)*wmac(i ,j,k+1) - - zxlo(i ,j,k ,n)*wmac(i ,j,k )) ); + // High side + Real qvyh = (vmac(i,j+1,k) - vmac(i,j,k)) * q(i,j,k,n); + sth += ( - (0.5*dtdy) * qvyh + - (0.5*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) + - xzlo(i ,j,k ,n)*umac(i ,j,k )) + - (0.5*dtdz)*(zxlo(i ,j,k+1,n)*wmac(i ,j,k+1) + - zxlo(i ,j,k ,n)*wmac(i ,j,k )) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -418,7 +422,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(xylo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xylo, l_xyhi; GodunovCornerCouple::AddCornerCoupleTermXY(l_xylo, l_xyhi, i, j, k, n, l_dt, dy, iconserv[n], @@ -435,7 +439,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Box(yxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yxlo, l_yxhi; GodunovCornerCouple::AddCornerCoupleTermYX(l_yxlo, l_yxhi, i, j, k, n, l_dt, dx, iconserv[n], @@ -455,41 +459,41 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real stl = zlo(i,j,k,n); - Real sth = zhi(i,j,k,n); + Real sth = zhi(i,j,k,n); - // Here we add dt/2 (-q w_z - (u q)_x - (v q)_y) to the term that is already - // q + dz/2 q_z + dt/2 (-w q_z) to get - // q + dz/2 q_z - dt/2 (w q_z + q w_z + (u q)_x + (v q)_y) which is equivalent to - // --> q + dz/2 q_z - dt/2 ( div (uvec q) ) - Real qwzl = (wmac(i,j,k) - wmac(i,j,k-1)) * q(i,j,k-1,n); - stl += ( - (0.5*dtdz) * qwzl - - (0.5*dtdx)*(xylo(i+1,j ,k-1,n)*umac(i+1,j ,k-1) - -xylo(i ,j ,k-1,n)*umac(i ,j ,k-1)) - - (0.5*dtdy)*(yxlo(i ,j+1,k-1,n)*vmac(i ,j+1,k-1) - -yxlo(i ,j ,k-1,n)*vmac(i ,j ,k-1)) ); + // Here we add dt/2 (-q w_z - (u q)_x - (v q)_y) to the term that is already + // q + dz/2 q_z + dt/2 (-w q_z) to get + // q + dz/2 q_z - dt/2 (w q_z + q w_z + (u q)_x + (v q)_y) which is equivalent to + // --> q + dz/2 q_z - dt/2 ( div (uvec q) ) + Real qwzl = (wmac(i,j,k) - wmac(i,j,k-1)) * q(i,j,k-1,n); + stl += ( - (0.5*dtdz) * qwzl + - (0.5*dtdx)*(xylo(i+1,j ,k-1,n)*umac(i+1,j ,k-1) + -xylo(i ,j ,k-1,n)*umac(i ,j ,k-1)) + - (0.5*dtdy)*(yxlo(i ,j+1,k-1,n)*vmac(i ,j+1,k-1) + -yxlo(i ,j ,k-1,n)*vmac(i ,j ,k-1)) ); - // Here we adjust for non-conservative by removing the q divu contribution to get - // q + dz/2 q_z - dt/2 ( div (uvec q) - q divu ) which is equivalent to - // --> q + dz/2 q_z - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k-1,n)*divu(i,j,k-1) : 0.; + // Here we adjust for non-conservative by removing the q divu contribution to get + // q + dz/2 q_z - dt/2 ( div (uvec q) - q divu ) which is equivalent to + // --> q + dz/2 q_z - dt/2 ( uvec dot grad q) + stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k-1,n)*divu(i,j,k-1) : 0.; - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k-1,n) : 0.; + stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k-1,n) : 0.; - // High side - Real qwzh = (wmac(i,j,k+1) - wmac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdz) * qwzh - - (0.5*dtdx)*(xylo(i+1,j ,k,n)*umac(i+1,j ,k) - -xylo(i ,j ,k,n)*umac(i ,j ,k)) - - (0.5*dtdy)*(yxlo(i ,j+1,k,n)*vmac(i ,j+1,k) - -yxlo(i ,j ,k,n)*vmac(i ,j ,k)) ); + // High side + Real qwzh = (wmac(i,j,k+1) - wmac(i,j,k)) * q(i,j,k,n); + sth += ( - (0.5*dtdz) * qwzh + - (0.5*dtdx)*(xylo(i+1,j ,k,n)*umac(i+1,j ,k) + -xylo(i ,j ,k,n)*umac(i ,j ,k)) + - (0.5*dtdy)*(yxlo(i ,j+1,k,n)*vmac(i ,j+1,k) + -yxlo(i ,j ,k,n)*vmac(i ,j ,k)) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp index ebb4fbb57..eb300ca79 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp @@ -21,7 +21,8 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, const BCRec * d_bcrec, const Geometry& geom, Real l_dt, const bool use_ppm, const bool use_forces_in_trans, - const int limiter_type) + const int limiter_type, + iMultiFab const* BC_MF) { Box const& domain = geom.Domain(); const Real* dx = geom.CellSize(); @@ -42,9 +43,10 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, Array4 const& umac = a_umac.array(mfi); Array4 const& vmac = a_vmac.array(mfi); - Array4 const& vel = a_vel.const_array(mfi); Array4 const& f = a_forces.const_array(mfi); + Array4 const& bc_arr = BC_MF ? BC_MF->const_array(mfi) + : Array4 {}; scratch.resize(bxg1, (ncomp*4 + 1)*AMREX_SPACEDIM); Real* p = scratch.dataPtr(); @@ -87,7 +89,7 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, else { PLM::PredictVelOnXFace( Box(u_ad), AMREX_SPACEDIM, Imx, Ipx, vel, vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec); PLM::PredictVelOnYFace( Box(v_ad), AMREX_SPACEDIM, Imy, Ipy, vel, vel, geom, l_dt, h_bcrec, d_bcrec); } @@ -95,13 +97,15 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, ComputeAdvectiveVel( Box(u_ad), Box(v_ad), u_ad, v_ad, Imx, Imy, Ipx, Ipy, - vel, f, domain, l_dt, d_bcrec, use_forces_in_trans); + vel, f, domain, l_dt, d_bcrec, use_forces_in_trans, + bc_arr); ExtrapVelToFacesOnBox( bx, ncomp, xbx, ybx, umac, vmac, vel, u_ad, v_ad, Imx, Imy, Ipx, Ipy, - f, domain, dx, l_dt, d_bcrec, use_forces_in_trans, p); + f, domain, dx, l_dt, d_bcrec, use_forces_in_trans, p, + bc_arr); Gpu::streamSynchronize(); // otherwise we might be using too much memory } @@ -122,7 +126,8 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, const Box& domain, Real l_dt, BCRec const* pbc, - bool l_use_forces_in_trans ) + bool l_use_forces_in_trans, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); const Dim3 dhi = amrex::ubound(domain); @@ -142,7 +147,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, hi += 0.5*l_dt*f(i ,j,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -163,7 +168,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, hi += 0.5*l_dt*f(i,j ,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -192,7 +197,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, Real l_dt, BCRec const* pbc, bool l_use_forces_in_trans, - Real* p) + Real* p, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); @@ -228,7 +234,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, hi += 0.5*l_dt*f(i ,j,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); xlo(i,j,k,n) = lo; @@ -245,7 +251,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, hi += 0.5*l_dt*f(i,j ,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); ylo(i,j,k,n) = lo; @@ -263,7 +269,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; l_yzlo = ylo(i,j,k,n); @@ -280,7 +286,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real stl = xlo(i,j,k,n) - (0.25*l_dt/dy)*(v_ad(i-1,j+1,k )+v_ad(i-1,j,k))* (yzlo(i-1,j+1,k )-yzlo(i-1,j,k)); Real sth = xhi(i,j,k,n) - (0.25*l_dt/dy)*(v_ad(i ,j+1,k )+v_ad(i ,j,k))* @@ -319,7 +325,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; l_xzlo = xlo(i,j,k,n); @@ -338,7 +344,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real stl = ylo(i,j,k,n) - (0.25*l_dt/dx)*(u_ad(i+1,j-1,k )+u_ad(i,j-1,k))* (xzlo(i+1,j-1,k )-xzlo(i,j-1,k)); diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp index 62df05012..3835b3e91 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp @@ -20,10 +20,11 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, MultiFab& a_vmac, MultiFab& a_wmac, const Vector & h_bcrec, - const BCRec * d_bcrec, + const BCRec * d_bcrec, const Geometry& geom, Real l_dt, bool use_ppm, bool use_forces_in_trans, - const int limiter_type) + const int limiter_type, + amrex::iMultiFab const* BC_MF) { Box const& domain = geom.Domain(); const Real* dx = geom.CellSize(); @@ -49,6 +50,8 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, Array4 const& vel = a_vel.const_array(mfi); Array4 const& f = a_forces.const_array(mfi); + Array4 const& bc_arr = BC_MF ? BC_MF->const_array(mfi) + : Array4 {}; scratch.resize(bxg1, (ncomp*4 + 1)*AMREX_SPACEDIM); Real* p = scratch.dataPtr(); @@ -97,24 +100,26 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, else { PLM::PredictVelOnXFace( Box(u_ad), AMREX_SPACEDIM, Imx, Ipx, vel, vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); PLM::PredictVelOnYFace( Box(v_ad), AMREX_SPACEDIM, Imy, Ipy, vel, vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); PLM::PredictVelOnZFace( Box(w_ad), AMREX_SPACEDIM, Imz, Ipz, vel, vel, - geom, l_dt, h_bcrec, d_bcrec); + geom, l_dt, h_bcrec, d_bcrec, bc_arr); } ComputeAdvectiveVel( Box(u_ad), Box(v_ad), Box(w_ad), u_ad, v_ad, w_ad, Imx, Imy, Imz, Ipx, Ipy, Ipz, - vel, f, domain, l_dt, d_bcrec, use_forces_in_trans); + vel, f, domain, l_dt, d_bcrec, use_forces_in_trans, + bc_arr); ExtrapVelToFacesOnBox( bx, ncomp, xbx, ybx, zbx, umac, vmac, wmac, vel, u_ad, v_ad, w_ad, Imx, Imy, Imz, Ipx, Ipy, Ipz, - f, domain, dx, l_dt, d_bcrec, use_forces_in_trans, p); + f, domain, dx, l_dt, d_bcrec, use_forces_in_trans, p, + bc_arr); Gpu::streamSynchronize(); // otherwise we might be using too much memory } @@ -139,7 +144,8 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, const Box& domain, Real l_dt, BCRec const* pbc, - bool l_use_forces_in_trans) + bool l_use_forces_in_trans, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); const Dim3 dhi = amrex::ubound(domain); @@ -159,7 +165,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, hi += Real(0.5)*l_dt*f(i ,j,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -180,7 +186,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, hi += Real(0.5)*l_dt*f(i,j ,k,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -201,7 +207,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, hi += Real(0.5)*l_dt*f(i,j,k ,n); } - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; @@ -235,7 +241,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real l_dt, BCRec const* pbc, bool l_use_forces_in_trans, - Real* p) + Real* p, + Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); @@ -274,7 +281,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } Real uad = u_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); @@ -297,7 +304,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } Real vad = v_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); @@ -321,7 +328,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } Real wad = w_ad(i,j,k); - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); @@ -363,7 +370,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zylo, l_zyhi; GodunovCornerCouple::AddCornerCoupleTermZY(l_zylo, l_zyhi, i, j, k, n, l_dt, dy, false, @@ -381,7 +388,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yzlo, l_yzhi; GodunovCornerCouple::AddCornerCoupleTermYZ(l_yzlo, l_yzhi, i, j, k, n, l_dt, dz, false, @@ -400,7 +407,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real stl = xlo(i,j,k,n) - (Real(0.25)*l_dt/dy)*(v_ad(i-1,j+1,k )+v_ad(i-1,j,k))* (yzlo(i-1,j+1,k )-yzlo(i-1,j,k)) @@ -450,7 +457,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xzlo, l_xzhi; GodunovCornerCouple::AddCornerCoupleTermXZ(l_xzlo, l_xzhi, i, j, k, n, l_dt, dz, false, @@ -468,7 +475,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_zxlo, l_zxhi; GodunovCornerCouple::AddCornerCoupleTermZX(l_zxlo, l_zxhi, i, j, k, n, l_dt, dx, false, @@ -487,7 +494,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real stl = ylo(i,j,k,n) - (Real(0.25)*l_dt/dx)*(u_ad(i+1,j-1,k )+u_ad(i,j-1,k))* (xzlo(i+1,j-1,k )-xzlo(i,j-1,k)) @@ -539,7 +546,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_xylo, l_xyhi; GodunovCornerCouple::AddCornerCoupleTermXY(l_xylo, l_xyhi, i, j, k, n, l_dt, dy, false, @@ -561,7 +568,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; - const auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real l_yxlo, l_yxhi; GodunovCornerCouple::AddCornerCoupleTermYX(l_yxlo, l_yxhi, i, j, k, n, l_dt, dx, false, @@ -580,7 +587,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; - auto bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); Real stl = zlo(i,j,k,n) - (Real(0.25)*l_dt/dx)*(u_ad(i+1,j ,k-1)+u_ad(i,j,k-1))* (xylo(i+1,j ,k-1)-xylo(i,j,k-1)) - (Real(0.25)*l_dt/dy)*(v_ad(i ,j+1,k-1)+v_ad(i,j,k-1))* diff --git a/Godunov/hydro_godunov_plm.H b/Godunov/hydro_godunov_plm.H index 82a9fdd70..81f68c726 100644 --- a/Godunov/hydro_godunov_plm.H +++ b/Godunov/hydro_godunov_plm.H @@ -26,7 +26,8 @@ void PredictVelOnXFace ( amrex::Box const& xebox, int ncomp, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr = {}); void PredictVelOnYFace ( amrex::Box const& yebox, int ncomp, amrex::Array4 const& Imy, amrex::Array4 const& Ipy, @@ -35,7 +36,8 @@ void PredictVelOnYFace ( amrex::Box const& yebox, int ncomp, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr= {}); #if (AMREX_SPACEDIM==3) void PredictVelOnZFace ( amrex::Box const& zebox, int ncomp, @@ -45,7 +47,8 @@ void PredictVelOnZFace ( amrex::Box const& zebox, int ncomp, const amrex::Geometry& geom, amrex::Real dt, amrex::Vector const& h_bcrec, - amrex::BCRec const* pbc); + amrex::BCRec const* pbc, + amrex::Array4 const& bc_arr= {}); #endif diff --git a/Godunov/hydro_godunov_plm.cpp b/Godunov/hydro_godunov_plm.cpp index 7e57a6f0f..b33c67a02 100644 --- a/Godunov/hydro_godunov_plm.cpp +++ b/Godunov/hydro_godunov_plm.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -36,7 +37,8 @@ PLM::PredictVelOnXFace ( Box const& xebox, int ncomp, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dx = geom.CellSize(0); const Real dtdx = dt/dx; @@ -47,16 +49,18 @@ PLM::PredictVelOnXFace ( Box const& xebox, int ncomp, // At an ext_dir boundary, the boundary value is on the face, not cell center. auto extdir_lohi = has_extdir_or_ho(h_bcrec.data(), ncomp, static_cast(Direction::x)); - bool has_extdir_or_ho_lo = extdir_lohi.first; - bool has_extdir_or_ho_hi = extdir_lohi.second; + // If we have a BC Array4, then we can't make blanket statements about the domain face + // Best to assume we could have face-based BC somewhere + bool has_extdir_or_ho_lo = bc_arr ? true : extdir_lohi.first; + bool has_extdir_or_ho_hi = bc_arr ? true : extdir_lohi.second; if ((has_extdir_or_ho_lo && domain_ilo >= xebox.smallEnd(0)-1) || (has_extdir_or_ho_hi && domain_ihi <= xebox.bigEnd(0))) { - amrex::ParallelFor(xebox, ncomp, [q,vcc,domain_ilo,domain_ihi,Imx,Ipx,dtdx,pbc] + amrex::ParallelFor(xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_ilo = (bc.lo(0) == BCType::ext_dir) || (bc.lo(0) == BCType::hoextrap); bool extdir_or_ho_ihi = (bc.hi(0) == BCType::ext_dir) || @@ -99,7 +103,8 @@ PLM::PredictVelOnYFace (Box const& yebox, int ncomp, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dy = geom.CellSize(1); const Real dtdy = dt/dy; @@ -110,16 +115,16 @@ PLM::PredictVelOnYFace (Box const& yebox, int ncomp, // At an ext_dir boundary, the boundary value is on the face, not cell center. auto extdir_lohi = has_extdir_or_ho(h_bcrec.data(), ncomp, static_cast(Direction::y)); - bool has_extdir_or_ho_lo = extdir_lohi.first; - bool has_extdir_or_ho_hi = extdir_lohi.second; + bool has_extdir_or_ho_lo = bc_arr ? true : extdir_lohi.first; + bool has_extdir_or_ho_hi = bc_arr ? true : extdir_lohi.second; if ((has_extdir_or_ho_lo && domain_jlo >= yebox.smallEnd(1)-1) || (has_extdir_or_ho_hi && domain_jhi <= yebox.bigEnd(1))) { - amrex::ParallelFor(yebox, ncomp, [q,vcc,domain_jlo,domain_jhi,Imy,Ipy,dtdy,pbc] + amrex::ParallelFor(yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_jlo = (bc.lo(1) == BCType::ext_dir) || (bc.lo(1) == BCType::hoextrap); bool extdir_or_ho_jhi = (bc.hi(1) == BCType::ext_dir) || @@ -163,7 +168,8 @@ PLM::PredictVelOnZFace ( Box const& zebox, int ncomp, const Geometry& geom, Real dt, Vector const& h_bcrec, - BCRec const* pbc) + BCRec const* pbc, + Array4 const& bc_arr) { const Real dz = geom.CellSize(2); const Real dtdz = dt/dz; @@ -174,16 +180,16 @@ PLM::PredictVelOnZFace ( Box const& zebox, int ncomp, // At an ext_dir boundary, the boundary value is on the face, not cell center. auto extdir_lohi = has_extdir_or_ho(h_bcrec.data(), ncomp, static_cast(Direction::z)); - bool has_extdir_or_ho_lo = extdir_lohi.first; - bool has_extdir_or_ho_hi = extdir_lohi.second; + bool has_extdir_or_ho_lo = bc_arr ? true : extdir_lohi.first; + bool has_extdir_or_ho_hi = bc_arr ? true : extdir_lohi.second; if ((has_extdir_or_ho_lo && domain_klo >= zebox.smallEnd(2)-1) || (has_extdir_or_ho_hi && domain_khi <= zebox.bigEnd(2))) { - amrex::ParallelFor(zebox, ncomp, [q,vcc,domain_klo,domain_khi,Ipz,Imz,dtdz,pbc] + amrex::ParallelFor(zebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - const auto& bc = pbc[n]; + const auto bc = HydroBC::getBC(i, j, k, n, domain_box, pbc, bc_arr); bool extdir_or_ho_klo = (bc.lo(2) == BCType::ext_dir) || (bc.lo(2) == BCType::hoextrap); bool extdir_or_ho_khi = (bc.hi(2) == BCType::ext_dir) || diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 4787ec012..1d30e8cb4 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -16,6 +16,8 @@ namespace PPM { enum limiters {VanLeer, WENOZ, NoLimiter}; +static constexpr int default_limiter = VanLeer; + struct nolimiter { explicit nolimiter() = default; diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H index 67412dbbf..bde1d94a7 100644 --- a/Projections/hydro_MacProjector.H +++ b/Projections/hydro_MacProjector.H @@ -120,7 +120,10 @@ public: void setDomainBC (const amrex::Array& lobc, const amrex::Array& hibc); - void setLevelBC (int amrlev, const amrex::MultiFab* levelbcdata); + void setLevelBC (int amrlev, const amrex::MultiFab* levelbcdata, + const amrex::MultiFab* robin_a = nullptr, + const amrex::MultiFab* robin_b = nullptr, + const amrex::MultiFab* robin_f = nullptr); void setCoarseFineBC (const amrex::MultiFab* crse, int crse_ratio) { m_linop->setCoarseFineBC(crse, crse_ratio);} diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index be3e6478a..adfcabf94 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -205,11 +205,12 @@ MacProjector::setDomainBC (const Array& lobc, void -MacProjector::setLevelBC (int amrlev, const MultiFab* levelbcdata) +MacProjector::setLevelBC (int amrlev, const MultiFab* levelbcdata, const MultiFab* robin_a, + const MultiFab* robin_b, const MultiFab* robin_f) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_needs_domain_bcs, "setDomainBC must be called before setLevelBC"); - m_linop->setLevelBC(amrlev, levelbcdata); + m_linop->setLevelBC(amrlev, levelbcdata, robin_a, robin_b, robin_f); m_needs_level_bcs[amrlev] = false; } diff --git a/Utils/hydro_bcs_K.H b/Utils/hydro_bcs_K.H index 93cc7ea42..11b0d3e64 100644 --- a/Utils/hydro_bcs_K.H +++ b/Utils/hydro_bcs_K.H @@ -27,6 +27,38 @@ * */ namespace HydroBC{ +// +// Choose between single BC per domain face or position dependent BC array +// +[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +const amrex::BCRec getBC (const int i, const int j, const int k, const int n, + const amrex::Box& m_domain, const amrex::BCRec* bcr, + amrex::Array4 const& bca) +{ + if ( !bca ) { return bcr[n]; } + + int lo[AMREX_SPACEDIM]; + int hi[AMREX_SPACEDIM]; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) { + int index[] = {i,j,k}; + + for ( int dim = 0; dim < AMREX_SPACEDIM; dim++){ + if (index[dim] < m_domain.smallEnd(dim)) { index[dim] = m_domain.smallEnd(dim); } + if (index[dim] > m_domain.bigEnd(dim)) { index[dim] = m_domain.bigEnd(dim); } + if (dim == dir) { index[dim] = m_domain.smallEnd(dim)-1; } + } + lo[dir] = bca.contains(index[0], index[1], index[2]) ? + bca(index[0], index[1], index[2], n) : 0; +// FIXME?? if we don't contain (i,j,k) then it doesn't matter what the bc is there, because we don't touch it + + index[dir] = m_domain.bigEnd(dir)+1; + hi[dir] = bca.contains(index[0], index[1], index[2]) ? + bca(index[0], index[1], index[2], n) : 0; + } + + amrex::BCRec bc(lo, hi); + return bc; +} /** * diff --git a/Utils/hydro_compute_edgestate_and_flux.cpp b/Utils/hydro_compute_edgestate_and_flux.cpp index 9885af5f5..726696e7d 100644 --- a/Utils/hydro_compute_edgestate_and_flux.cpp +++ b/Utils/hydro_compute_edgestate_and_flux.cpp @@ -44,7 +44,8 @@ namespace { bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, std::string const& advection_type, - int limiter_type = PPM::VanLeer) + int limiter_type, + amrex::Array4 const& bc_arr) { #ifdef AMREX_USE_EB @@ -90,7 +91,7 @@ namespace { AMREX_D_DECL(apx,apy,apz), vfrac, AMREX_D_DECL(fcx,fcy,fcz), ccc, is_velocity, - values_on_eb_inflow); + values_on_eb_inflow, bc_arr); } else if (advection_type == "BDS") { @@ -122,7 +123,7 @@ namespace { geom, l_dt, d_bcrec, iconserv, godunov_use_ppm, godunov_use_forces_in_trans, - is_velocity, limiter_type); + is_velocity, limiter_type, bc_arr); } else if (advection_type == "BDS") { @@ -165,7 +166,8 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type) + int limiter_type, + amrex::Array4 const& bc_arr) { ComputeFluxesOnBoxFromState(bx, ncomp, mfi, q, @@ -177,7 +179,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, ebfact, /*values_on_eb_inflow*/ Array4{}, godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, fluxes_are_area_weighted, advection_type, - limiter_type); + limiter_type, bc_arr); } #endif @@ -209,7 +211,8 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type) + int limiter_type, + amrex::Array4 const& bc_arr) { ComputeFluxesOnBoxFromState(bx, ncomp, mfi, q, @@ -224,7 +227,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, #endif godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, fluxes_are_area_weighted, advection_type, - limiter_type); + limiter_type, bc_arr); } @@ -257,7 +260,8 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type) + int limiter_type, + amrex::Array4 const& bc_arr) { #ifdef AMREX_USE_EB @@ -284,7 +288,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, ebfact, values_on_eb_inflow, regular, #endif godunov_use_ppm, godunov_use_forces_in_trans, - is_velocity, advection_type, limiter_type); + is_velocity, advection_type, limiter_type, bc_arr); } // Compute fluxes. diff --git a/Utils/hydro_extrap_vel_to_faces.cpp b/Utils/hydro_extrap_vel_to_faces.cpp index 36d36be4c..6cdc7a9e1 100644 --- a/Utils/hydro_extrap_vel_to_faces.cpp +++ b/Utils/hydro_extrap_vel_to_faces.cpp @@ -53,15 +53,20 @@ HydroUtils::ExtrapVelToFaces ( amrex::MultiFab const& vel, #endif bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, - int limiter_type) + int limiter_type, + iMultiFab* BC_MF) { + amrex::ignore_unused(BC_MF); + if (advection_type == "Godunov") { #ifdef AMREX_USE_EB if (!ebfact.isAllRegular()) EBGodunov::ExtrapVelToFaces(vel, vel_forces, AMREX_D_DECL(u_mac, v_mac, w_mac), h_bcrec, d_bcrec, geom, dt, - velocity_on_eb_inflow); // Note that PPM not supported for EB + velocity_on_eb_inflow, + // Note that PPM is not supported for EB + BC_MF); else #endif Godunov::ExtrapVelToFaces(vel, vel_forces, diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index abdbc26f3..f0e1a6652 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -20,10 +20,9 @@ */ namespace HydroUtils { - /** * \brief Compute edge state and flux. Most general version for use with multilevel synchonization. - * + * All other versions ultimately call this one. */ void ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mfi, @@ -55,7 +54,8 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::default_limiter, + amrex::Array4 const& bc_arr = {}); /** * \brief Compute edge state and flux. For typical advection, and also allows for inflow on EB. * @@ -87,7 +87,8 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::default_limiter, + amrex::Array4 const& bc_arr = {}); /** * \brief Compute edge state and flux. For typical advection, but no inflow through EB. @@ -118,7 +119,8 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool godunov_use_ppm, bool godunov_use_forces_in_trans, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::default_limiter, + amrex::Array4 const& bc_arr = {}); #endif #ifdef AMREX_USE_EB @@ -135,7 +137,7 @@ ExtrapVelToFaces ( amrex::MultiFab const& vel, const amrex::EBFArrayBoxFactory& ebfact, bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::default_limiter); #endif void @@ -154,7 +156,8 @@ ExtrapVelToFaces ( amrex::MultiFab const& vel, #endif bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, - int limiter_type = PPM::VanLeer); + int limiter_type = PPM::default_limiter, + amrex::iMultiFab* BC_MF = nullptr); /** * \brief If convective, compute convTerm = u dot grad q = div (u q) - q div(u).