diff --git a/EBGodunov/hydro_ebgodunov_bcs_K.H b/EBGodunov/hydro_ebgodunov_bcs_K.H index 64747eafb..d2fe7e67c 100644 --- a/EBGodunov/hydro_ebgodunov_bcs_K.H +++ b/EBGodunov/hydro_ebgodunov_bcs_K.H @@ -25,8 +25,8 @@ namespace EBGodunovBC { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetXBCs (const int i, const int j, const int k, const int n, const amrex::Array4 &s, - amrex::Real &lo, - amrex::Real &hi, + amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi, const bool is_velocity ) @@ -37,7 +37,14 @@ void SetXBCs (const int i, const int j, const int k, const int n, // Low X if (i <= domlo) { - if (bclo==BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) ) + { + lo = s(domlo-1, j, k, n); + hi = s(domlo-1, j, k, n); + } + else if ( bclo == BCType::ext_dir || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo = s(domlo-1,j,k,n); // For turbulent inflow, there are times when the inflow face @@ -46,7 +53,8 @@ void SetXBCs (const int i, const int j, const int k, const int n, // tangential components to transport values from the interior. if( n == XVEL && is_velocity ) hi=lo; } - else if (bclo == BCType::foextrap || bclo == BCType::hoextrap) + else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || + (bclo == BCType::direction_dependent && velp < 0.0) ) { lo = hi; } else if (bclo == BCType::reflect_even) @@ -71,8 +79,8 @@ void SetXBCs (const int i, const int j, const int k, const int n, else if (bclo == BCType::reflect_odd) { if ( i==domlo ) { - hi = 0.; - lo = 0.; + hi = Real(0.); + lo = Real(0.); } else { Abort("EBGodunovBC::SetXBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // lo = -hi_arr(2*domlo-i ,j,k,n); @@ -83,12 +91,20 @@ void SetXBCs (const int i, const int j, const int k, const int n, // High X else if (i > domhi) { - if (bchi==BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == XVEL) && (s(domhi+1,j,k,n) <= 0.0) ) + { + hi = s(domhi+1, j, k, n); + lo = s(domhi+1, j, k, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { hi = s(domhi+1,j,k,n) ; if( n ==XVEL && is_velocity ) lo=hi; } - else if (bchi == BCType::foextrap || bchi == BCType::hoextrap ) + else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || + (bchi== BCType::direction_dependent && velm > 0.0) ) { hi = lo; } @@ -104,8 +120,8 @@ void SetXBCs (const int i, const int j, const int k, const int n, else if (bchi == BCType::reflect_odd) { if ( i==domhi+1 ) { - hi = 0.; - lo = 0.; + hi = Real(0.); + lo = Real(0.); } else { Abort("EBGodunovBC::SetXBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // hi = -lo_arr(2*(domhi+1)-i-1,j,k,n); @@ -120,8 +136,8 @@ void SetXBCs (const int i, const int j, const int k, const int n, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetYBCs (const int i, const int j, const int k, const int n, const amrex::Array4 &s, - amrex::Real &lo, - amrex::Real &hi, + amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi, const bool is_velocity ) @@ -132,12 +148,20 @@ void SetYBCs (const int i, const int j, const int k, const int n, // Low Y if (j <= domlo) { - if (bclo==BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == YVEL) && (s(i,domlo-1,k,n) >= 0.0) ) + { + lo = s(i, domlo-1, k, n); + hi = s(i, domlo-1, k, n); + } + else if ( bclo == BCType::ext_dir || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo = s(i,domlo-1,k,n); if ( n == YVEL && is_velocity ) hi = lo; } - else if (bclo == BCType::foextrap || bclo == BCType::hoextrap) + else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || + (bclo == BCType::direction_dependent && velp < 0.0) ) { lo = hi; } @@ -153,8 +177,8 @@ void SetYBCs (const int i, const int j, const int k, const int n, else if(bclo == BCType::reflect_odd) { if ( j==domlo ) { - hi = 0.; - lo = 0.; + hi = Real(0.); + lo = Real(0.); } else { Abort("EBGodunovBC::SetYBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // lo = -hi_arr(i,2*domlo-j ,k,n); @@ -165,12 +189,20 @@ void SetYBCs (const int i, const int j, const int k, const int n, // High Y else if (j > domhi) { - if (bchi==BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == YVEL) && (s(i,domhi+1,k,n) <= 0.0) ) + { + hi = s(i, domhi+1, k, n); + lo = s(i, domhi+1, k, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { hi = s(i,domhi+1,k,n); if( n == YVEL && is_velocity ) lo = hi ; } - else if (bchi == BCType::foextrap || bchi == BCType::hoextrap) + else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || + (bchi== BCType::direction_dependent && velm > 0.0) ) { hi = lo; } @@ -186,8 +218,8 @@ void SetYBCs (const int i, const int j, const int k, const int n, else if (bchi == BCType::reflect_odd) { if ( j==domhi+1 ) { - hi = 0.; - lo = 0.; + hi = Real(0.); + lo = Real(0.); } else { Abort("EBGodunovBC::SetYBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // hi = -lo_arr(i,2*(domhi+1)-j-1,k,n); @@ -201,13 +233,13 @@ void SetYBCs (const int i, const int j, const int k, const int n, #if (AMREX_SPACEDIM==3) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void SetZBCs(const int i, const int j, const int k, const int n, - const amrex::Array4 &s, - amrex::Real &lo, - amrex::Real &hi, - const int bclo, const int bchi, - const int domlo, const int domhi, - const bool is_velocity) +void SetZBCs (const int i, const int j, const int k, const int n, + const amrex::Array4 &s, + amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, + const int bclo, const int bchi, + const int domlo, const int domhi, + const bool is_velocity) { using namespace amrex; @@ -215,12 +247,20 @@ void SetZBCs(const int i, const int j, const int k, const int n, // Low Z if (k <= domlo) { - if (bclo==BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == ZVEL) && (s(i,j,domlo-1,n) >= 0.0) ) + { + lo = s(i, j, domlo-1, n); + hi = s(i, j, domlo-1, n); + } + else if ( bclo == BCType::ext_dir || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo =s(i,j,domlo-1,n); if ( n == ZVEL && is_velocity ) hi = lo; } - else if (bclo == BCType::foextrap || bclo == BCType::hoextrap) + else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || + (bclo == BCType::direction_dependent && velp < 0.0) ) { lo = hi; } @@ -234,8 +274,8 @@ void SetZBCs(const int i, const int j, const int k, const int n, { Abort("EBGodunovBC::SetZBCs not yet implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // if ( k==domlo ) { - // hi = 0.; - // lo = 0.; + // hi = Real(0.); + // lo = Real(0.); // } else { // lo = -hi_arr(i,j,2*domlo-k ,n); // hi = -lo_arr(i,j,2*domlo-k-1,n); @@ -245,12 +285,20 @@ void SetZBCs(const int i, const int j, const int k, const int n, // High Z else if (k > domhi) { - if (bchi==BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == ZVEL) && (s(i,j,domhi+1,n) <= 0.0) ) + { + hi = s(i, j, domhi+1, n); + lo = s(i, j, domhi+1, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { hi = s(i,j,domhi+1,n); if ( n == ZVEL && is_velocity ) lo = hi ; } - else if (bchi == BCType::foextrap || bchi == BCType::hoextrap) + else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || + (bchi== BCType::direction_dependent && velm > 0.0) ) { hi = lo; } @@ -264,8 +312,8 @@ void SetZBCs(const int i, const int j, const int k, const int n, { Abort("EBGodunovBC::SetZBCs not yet implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs."); // if ( k==domhi+1 ) { - // hi = 0.; - // lo = 0.; + // hi = Real(0.); + // lo = Real(0.); // } else { // hi = -lo_arr(i,j,2*(domhi+1)-k-1,n); // lo = -hi_arr(i,j,2*(domhi+1)-k ,n); diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp index c5b50050f..c813017b3 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp @@ -110,7 +110,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; @@ -122,7 +123,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -143,9 +145,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, l_yzlo = ylo(i,j,k,n); l_yzhi = yhi(i,j,k,n); - Real vad = v_mac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + Real vad = v_mac(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); @@ -221,7 +224,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, } 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -261,9 +265,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, l_xzlo = xlo(i,j,k,n); l_xzhi = xhi(i,j,k,n); - Real uad = u_mac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = u_mac(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); @@ -339,7 +344,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, } 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if (!allow_inflow_on_outflow) { 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 ccb97d446..226f2b681 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp @@ -138,8 +138,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real uad = u_mac(i,j,k); 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; @@ -157,7 +157,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -171,15 +172,15 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Real lo = Ipz(i,j,k-1,n); Real hi = Imz(i,j,k ,n); - Real wad = w_mac(i,j,k); - 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); + HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, w_mac(i,j,k), w_mac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); zlo(i,j,k,n) = lo; zhi(i,j,k,n) = hi; + Real wad = w_mac(i,j,k); Real st = (wad >= 0.) ? lo : hi; Real fuz = (amrex::Math::abs(wad) < small_vel) ? 0. : 1.; Imz(i,j,k,n) = fuz*st + (Real(1.0) - fuz)*Real(0.5)*(hi + lo); @@ -208,9 +209,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, v_mac, yed); - Real wad = w_mac(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, w_mac(i,j,k), w_mac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + Real wad = w_mac(i,j,k); Real st = (wad >= 0.) ? l_zylo : l_zyhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0); zylo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zyhi + l_zylo); @@ -225,9 +227,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, w_mac, zed); - Real vad = v_mac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + Real vad = v_mac(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0); yzlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yzhi + l_yzlo); @@ -305,7 +308,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, } 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -347,9 +351,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, xlo(i,j,k,n), xhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, w_mac, zed); - Real uad = u_mac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = u_mac(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xzlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xzhi + l_xzlo); @@ -364,9 +369,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, u_mac, xed); - Real wad = w_mac(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, w_mac(i,j,k), w_mac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + Real wad = w_mac(i,j,k); Real st = (wad >= 0.) ? l_zxlo : l_zxhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; zxlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zxhi + l_zxlo); @@ -444,7 +450,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, } 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -485,9 +492,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, xlo(i,j,k,n), xhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, v_mac, yed); - Real uad = u_mac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, u_mac(i,j,k), u_mac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = u_mac(i,j,k); Real st = (uad >= 0.) ? l_xylo : l_xyhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xylo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xyhi + l_xylo); @@ -502,9 +510,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, u_mac, xed); - Real vad = v_mac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, v_mac(i,j,k), v_mac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + Real vad = v_mac(i,j,k); Real st = (vad >= 0.) ? l_yxlo : l_yxhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yxlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yxhi + l_yxlo); @@ -581,7 +590,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, } 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); + HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, w_mac(i,j,k), w_mac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); if (!allow_inflow_on_outflow) { 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 64c8668c8..ffbd52207 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp @@ -277,7 +277,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real hi = Imx(i ,j,k,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); + EBGodunovBC::SetXBCs(i, j, k, n, vel, lo, hi, lo, hi, + bc.lo(0), bc.hi(0), dlo.x, dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -297,7 +298,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real hi = Imy(i,j ,k,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); + EBGodunovBC::SetYBCs(i, j, k, n, vel, lo, hi, lo, hi, + bc.lo(1), bc.hi(1), dlo.y, dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -318,7 +320,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx, Real hi = Imz(i,j,k ,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); + EBGodunovBC::SetZBCs(i, j, k, n, vel, lo, hi, lo, hi, + bc.lo(2), bc.hi(2), dlo.z, dhi.z, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp index a38f82e3d..f82400e83 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp @@ -64,7 +64,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; @@ -76,7 +77,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -103,9 +105,10 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, l_yzlo = ylo(i,j,k,n); l_yzhi = yhi(i,j,k,n); - Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yhat(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); @@ -189,7 +192,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, } } - HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -230,9 +234,10 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, l_xzlo = xlo(i,j,k,n); l_xzhi = xhi(i,j,k,n); - Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + Real uad = u_ad(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xhat(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); @@ -312,7 +317,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, } } - HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp index d5145961a..f7ecdb56f 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp @@ -75,14 +75,15 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real lo = Ipx(i-1,j,k,n); Real hi = Imx(i ,j,k,n); - Real uad = u_ad(i,j,k); 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); + EBGodunovBC::SetXBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), bc.hi(0), dlo.x, dhi.x, true); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; + Real uad = u_ad(i,j,k); Real st = (uad >= 0.) ? lo : hi; Real fu = (amrex::Math::abs(uad) < small_vel) ? Real(0.0) : Real(1.0); Imx(i, j, k, n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (hi + lo); // store xedge @@ -92,14 +93,15 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real lo = Ipy(i,j-1,k,n); Real hi = Imy(i,j ,k,n); - Real vad = v_ad(i,j,k); 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); + EBGodunovBC::SetYBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), bc.hi(1), dlo.y, dhi.y, true); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; + Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? lo : hi; Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0); Imy(i, j, k, n) = fu*st + (Real(1.0) - fu)*Real(0.5)*(hi + lo); // store yedge @@ -109,14 +111,15 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real lo = Ipz(i,j,k-1,n); Real hi = Imz(i,j,k ,n); - Real wad = w_ad(i,j,k); 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); + EBGodunovBC::SetZBCs(i, j, k, n, q, lo, hi, w_ad(i,j,k), w_ad(i,j,k), + bc.lo(2), bc.hi(2), dlo.z, dhi.z, true); zlo(i,j,k,n) = lo; zhi(i,j,k,n) = hi; + Real wad = w_ad(i,j,k); Real st = (wad >= 0.) ? lo : hi; Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0); Imz(i, j, k, n) = fu*st + (Real(1.0) - fu)*Real(0.5)*(hi + lo); // store zedge @@ -164,10 +167,11 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge); - Real wad = w_ad(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, w_ad(i,j,k), w_ad(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + Real wad = w_ad(i,j,k); Real st = (wad >= 0.) ? l_zylo : l_zyhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0); zylo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zyhi + l_zylo); @@ -182,9 +186,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge); - Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0); yzlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yzhi + l_yzlo); @@ -266,7 +271,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } } - HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -316,10 +322,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, xlo(i,j,k,n), xhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge); - Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); - + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + Real uad = u_ad(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : Real(1.0); xzlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xzhi + l_xzlo); @@ -334,10 +340,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge); - Real wad = w_ad(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); - + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, w_ad(i,j,k), w_ad(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + Real wad = w_ad(i,j,k); Real st = (wad >= 0.) ? l_zxlo : l_zxhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : Real(1.0); zxlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zxhi + l_zxlo); @@ -419,7 +425,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } } - HydroBC::SetYEdgeBCs( i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -470,7 +477,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge); Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = (uad >= 0.) ? l_xylo : l_xyhi; @@ -491,10 +499,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge); - Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); - + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? l_yxlo : l_yxhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : Real(1.0); yxlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yxhi + l_yxlo); @@ -575,7 +583,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, } } - HydroBC::SetZEdgeBCs( i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, w_ad(i,j,k), w_ad(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); if (!allow_inflow_on_outflow) { if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) diff --git a/EBMOL/hydro_ebmol_edge_state_K.H b/EBMOL/hydro_ebmol_edge_state_K.H index dd7f438fd..19a9f33a1 100644 --- a/EBMOL/hydro_ebmol_edge_state_K.H +++ b/EBMOL/hydro_ebmol_edge_state_K.H @@ -141,7 +141,8 @@ amrex::Real hydro_ebmol_xedge_state_extdir ( AMREX_D_DECL(int i, int j, int k), qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin); - HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, umac(i,j,k), umac(i,j,k), + d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity); if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { @@ -249,7 +250,8 @@ amrex::Real hydro_ebmol_xedge_state ( AMREX_D_DECL(int i, int j, int k), int n, qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin); - HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, umac(i,j,k), umac(i,j,k), + d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity); if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { @@ -403,7 +405,8 @@ amrex::Real hydro_ebmol_yedge_state_extdir ( AMREX_D_DECL(int i, int j, int k), + delta_y * slopes_eb_lo[1]; #endif - HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, vmac(i,j,k), vmac(i,j,k), + d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity); if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { @@ -515,7 +518,8 @@ amrex::Real hydro_ebmol_yedge_state ( AMREX_D_DECL(int i, int j, int k), int n, qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin); - HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, vmac(i,j,k), vmac(i,j,k), + d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity); if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { @@ -658,7 +662,8 @@ amrex::Real hydro_ebmol_zedge_state_extdir ( int i, int j, int k, int n, qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin); - HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, wmac(i,j,k), wmac(i,j,k), + d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity); if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { @@ -752,7 +757,8 @@ amrex::Real hydro_ebmol_zedge_state ( int i, int j, int k, int n, qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin); - HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, wmac(i,j,k), wmac(i,j,k), + d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity); if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { diff --git a/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp b/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp index 284702435..d1b475d38 100644 --- a/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp +++ b/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp @@ -178,7 +178,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, #endif umns = amrex::max(amrex::min(umns, cc_umax), cc_umin); - HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); + HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, umns, upls, + d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) { @@ -267,7 +268,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, #endif umns = amrex::max(amrex::min(umns, cc_umax), cc_umin); - HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); + HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, umns, upls, + d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) { @@ -413,7 +415,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, vmns = amrex::max(amrex::min(vmns, cc_vmax), cc_vmin); - HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); + HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, vmns, vpls, + d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) { @@ -505,7 +508,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, vmns = amrex::max(amrex::min(vmns, cc_vmax), cc_vmin); - HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); + HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, vmns, vpls, + d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) { @@ -633,7 +637,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, wmns = amrex::max(amrex::min(wmns, cc_wmax), cc_wmin); - HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); + HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, wmns, wpls, + d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) { @@ -713,7 +718,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, wmns = amrex::max(amrex::min(wmns, cc_wmax), cc_wmin); - HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); + HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, wmns, wpls, + d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) { diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp index 68748c315..033483475 100644 --- a/Godunov/hydro_godunov_edge_state_2D.cpp +++ b/Godunov/hydro_godunov_edge_state_2D.cpp @@ -140,7 +140,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; }, @@ -157,7 +158,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -179,7 +181,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, l_yzlo = ylo(i,j,k,n); l_yzhi = yhi(i,j,k,n); Real vad = vmac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; @@ -228,7 +231,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -263,9 +267,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, l_xzlo = xlo(i,j,k,n); l_xzhi = xhi(i,j,k,n); - Real uad = umac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = umac(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); @@ -314,7 +319,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if (!allow_inflow_on_outflow) { 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 6a573c599..6648d8af4 100644 --- a/Godunov/hydro_godunov_edge_state_3D.cpp +++ b/Godunov/hydro_godunov_edge_state_3D.cpp @@ -162,7 +162,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; Real st = (uval) ? lo : hi; @@ -185,7 +186,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -208,7 +210,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wmac(i,j,k), wmac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); zlo(i,j,k,n) = lo; zhi(i,j,k,n) = hi; @@ -236,9 +239,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, vmac, Imy); - Real wad = wmac(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wmac(i,j,k), wmac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + Real wad = wmac(i,j,k); Real st = (wad >= 0.) ? l_zylo : l_zyhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; zylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zyhi + l_zylo); @@ -253,9 +257,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, wmac, Imz); - Real vad = vmac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + Real vad = vmac(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); @@ -302,7 +307,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -339,9 +345,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, xlo(i,j,k,n), xhi(i,j,k,n), q, divu, wmac, Imz); - Real uad = umac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = umac(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); @@ -356,9 +363,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, zlo(i,j,k,n), zhi(i,j,k,n), q, divu, umac, Imx); - Real wad = wmac(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wmac(i,j,k), wmac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); + Real wad = wmac(i,j,k); Real st = (wad >= 0.) ? l_zxlo : l_zxhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; zxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zxhi + l_zxlo); @@ -403,7 +411,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -440,9 +449,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, xlo(i,j,k,n), xhi(i,j,k,n), q, divu, vmac, Imy); - Real uad = umac(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, umac(i,j,k), umac(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity); + Real uad = umac(i,j,k); Real st = (uad >= 0.) ? l_xylo : l_xyhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xyhi + l_xylo); @@ -457,9 +467,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, ylo(i,j,k,n), yhi(i,j,k,n), q, divu, umac, Imx); - Real vad = vmac(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vmac(i,j,k), vmac(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity); + Real vad = vmac(i,j,k); Real st = (vad >= 0.) ? l_yxlo : l_yxhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yxhi + l_yxlo); @@ -503,7 +514,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; 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); + HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wmac(i,j,k), wmac(i,j,k), + bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity); if (!allow_inflow_on_outflow) { 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 0f60235ff..85735cdbc 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp @@ -155,7 +155,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, } 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); + HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -176,7 +176,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, } 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); + HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -243,7 +243,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, } 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; @@ -260,7 +261,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, } 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -283,9 +285,10 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, l_yzlo = ylo(i,j,k,n); l_yzhi = yhi(i,j,k,n); - Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; yzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); @@ -306,7 +309,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, sth += 0.5 * l_dt * f(i ,j,k,n); } - HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -341,16 +345,15 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, l_xzlo = xlo(i,j,k,n); l_xzhi = xhi(i,j,k,n); - Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); - + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k), + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + Real uad = u_ad(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; xzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); }); - amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; @@ -367,7 +370,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, sth += 0.5 * l_dt * f(i,j ,k,n); } - HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k), + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp index 6de8c3bf2..31c136493 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp @@ -173,7 +173,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, } 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); + HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -194,7 +194,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, } 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); + HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -215,7 +215,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, } 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); + HydroBC::SetZEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); Real st = ( (lo+hi) >= 0.) ? lo : hi; bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); @@ -291,7 +291,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real uad = u_ad(i,j,k); 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); + HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; @@ -314,7 +314,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real vad = v_ad(i,j,k); 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); + HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; @@ -338,7 +338,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Real wad = w_ad(i,j,k); 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); + HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); zlo(i,j,k,n) = lo; zhi(i,j,k,n) = hi; @@ -386,7 +386,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, v_ad, yedge); Real wad = w_ad(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad, + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); Real st = (wad >= 0.) ? l_zylo : l_zyhi; @@ -404,7 +405,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, w_ad, zedge); Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; @@ -432,7 +434,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, sth += Real(0.5) * l_dt * f(i ,j,k,n); } - HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + Real uad = u_ad(i,j,k); + HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) @@ -475,7 +478,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, w_ad, zedge); Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; @@ -493,7 +497,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, u_ad, xedge); Real wad = w_ad(i,j,k); - HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad, + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); Real st = (wad >= 0.) ? l_zxlo : l_zxhi; Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; @@ -522,7 +527,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, sth += Real(0.5) * l_dt * f(i,j ,k,n); } - HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + Real vad = v_ad(i,j,k); + HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) @@ -566,8 +572,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, v_ad, yedge); Real uad = u_ad(i,j,k); - HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); - + HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, + bc.lo(0), dlo.x, bc.hi(0), dhi.x, true); Real st = (uad >= 0.) ? l_xylo : l_xyhi; Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; @@ -588,7 +594,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, q, divu, u_ad, xedge); Real vad = v_ad(i,j,k); - HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); + HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, + bc.lo(1), dlo.y, bc.hi(1), dhi.y, true); Real st = (vad >= 0.) ? l_yxlo : l_yxhi; @@ -615,7 +622,9 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, sth += Real(0.5) * l_dt * f(i,j,k ,n); } - HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); + Real wad = w_ad(i,j,k); + HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, + bc.lo(2), dlo.z, bc.hi(2), dhi.z, true); if (!allow_inflow_on_outflow) { if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) diff --git a/Godunov/hydro_godunov_plm.H b/Godunov/hydro_godunov_plm.H index 81f68c726..fbe5a68af 100644 --- a/Godunov/hydro_godunov_plm.H +++ b/Godunov/hydro_godunov_plm.H @@ -73,7 +73,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, int order = 4; - if (i == domain_ilo && (bc.lo(0) == BCType::ext_dir)) + if (i == domain_ilo && ((bc.lo(0) == BCType::ext_dir) || + (bc.lo(0) == BCType::direction_dependent && umac >= 0.0)) ) { umns = S(i-1,j,k,n); @@ -88,7 +89,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, } } - else if (i == domain_ihi+1 && (bc.hi(0) == BCType::ext_dir)) + else if (i == domain_ihi+1 && ((bc.hi(0) == BCType::ext_dir) || + (bc.hi(0) == BCType::direction_dependent && umac <= 0.0)) ) { upls = S(i ,j,k,n); @@ -140,7 +142,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, int order = 4; - if (j == domain_jlo && (bc.lo(1) == BCType::ext_dir)) + if (j == domain_jlo && ((bc.lo(1) == BCType::ext_dir) || + (bc.lo(1) == BCType::direction_dependent && vmac >= 0.0)) ) { vmns = S(i,j-1,k,n); if ( n==YVEL && is_velocity ) @@ -153,7 +156,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, amrex_calc_yslope_extdir(i,j ,k,n,order,S, extdir_or_ho_jlo, extdir_or_ho_jhi, domain_jlo, domain_jhi); } } - else if (j == domain_jhi+1 && (bc.hi(1) == BCType::ext_dir)) + else if (j == domain_jhi+1 && ((bc.hi(1) == BCType::ext_dir) || + (bc.hi(1) == BCType::direction_dependent && vmac <= 0.0)) ) { vpls = S(i,j ,k,n); if ( n==YVEL && is_velocity ) @@ -205,7 +209,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, int order = 4; - if (k == domain_klo && (bc.lo(2) == BCType::ext_dir)) + if (k == domain_klo && ((bc.lo(2) == BCType::ext_dir) || + (bc.lo(2) == BCType::direction_dependent && wmac >= 0.0)) ) { wmns = S(i,j,k-1,n); if ( n == ZVEL && is_velocity ) @@ -218,7 +223,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, amrex_calc_zslope_extdir(i,j,k ,n,order,S, extdir_or_ho_klo, extdir_or_ho_khi, domain_klo, domain_khi); } } - else if (k == domain_khi+1 && (bc.hi(2) == BCType::ext_dir)) + else if (k == domain_khi+1 && ((bc.hi(2) == BCType::ext_dir) || + (bc.hi(2) == BCType::direction_dependent && wmac <= 0.0)) ) { wpls = S(i,j,k ,n); if ( n == ZVEL && is_velocity ) diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 83e4a0c45..739270189 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -10,6 +10,7 @@ #include #include +#include #include namespace PPM { @@ -281,12 +282,14 @@ void SetXBCs ( const int i, const int j, const int k, const int n, amrex::Real &sm, amrex::Real &sp, amrex::Real &sedge1, amrex::Real &sedge2, const amrex::Array4 &s, + const amrex::Real velm, const amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi) { using namespace amrex; - if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) + if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { if (i == domlo) { @@ -326,7 +329,8 @@ void SetXBCs ( const int i, const int j, const int k, const int n, } } - if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) + if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { if (i == domhi) { @@ -373,12 +377,14 @@ void SetYBCs ( const int i, const int j, const int k, const int n, amrex::Real &sm, amrex::Real &sp, amrex::Real &sedge1, amrex::Real &sedge2, const amrex::Array4 &s, + const amrex::Real velm, const amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi) { using namespace amrex; - if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) + if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || + (bclo == BCType::direction_dependent && velm >= 0.0) ) { if (j == domlo) { @@ -418,7 +424,8 @@ void SetYBCs ( const int i, const int j, const int k, const int n, } } - if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) + if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || + (bchi == BCType::direction_dependent && velp <= 0.0) ) { if (j == domhi) { @@ -465,12 +472,14 @@ void SetZBCs ( const int i, const int j, const int k, const int n, amrex::Real &sm, amrex::Real &sp, amrex::Real &sedge1, amrex::Real &sedge2, const amrex::Array4 &s, + const amrex::Real velm, const amrex::Real velp, const int bclo, const int bchi, const int domlo, const int domhi) { using namespace amrex; - if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) + if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || + (bclo == BCType::direction_dependent && velm >= 0.0) ) { if (k == domlo) { @@ -510,7 +519,8 @@ void SetZBCs ( const int i, const int j, const int k, const int n, } } - if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) + if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || + (bchi == BCType::direction_dependent && velp <= 0.0) ) { if (k == domhi) { @@ -582,7 +592,8 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(0), bc.hi(0), domlo, domhi); + SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, + bc.lo(0), bc.hi(0), domlo, domhi); amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); @@ -629,7 +640,8 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(1), bc.hi(1), domlo, domhi); + SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, + bc.lo(1), bc.hi(1), domlo, domhi); amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); @@ -677,7 +689,8 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(2), bc.hi(2), domlo, domhi); + SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, + bc.lo(2), bc.hi(2), domlo, domhi); amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); @@ -730,7 +743,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(0), bc.hi(0), domlo, domhi); + SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(0), bc.hi(0), domlo, domhi); amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp); @@ -776,7 +790,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(1), bc.hi(1), domlo, domhi); + SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(1), bc.hi(1), domlo, domhi); amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); @@ -825,7 +840,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, amrex::Real sm, sp; amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2); - SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(2), bc.hi(2), domlo, domhi); + SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), + bc.lo(2), bc.hi(2), domlo, domhi); amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx; @@ -919,9 +935,7 @@ void PredictStateOnFaces (amrex::Box const& bx, #endif }); } - - -} +} // namespace #endif /** @} */ diff --git a/MOL/hydro_mol_edge_state_K.H b/MOL/hydro_mol_edge_state_K.H index 6ff324dde..7adf93888 100644 --- a/MOL/hydro_mol_edge_state_K.H +++ b/MOL/hydro_mol_edge_state_K.H @@ -45,7 +45,8 @@ amrex::Real hydro_mol_xedge_state_extdir ( int i, int j, int k, int n, amrex::Real qmns = q(i-1,j,k,n) + 0.5 * amrex_calc_xslope_extdir( i-1, j, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity); + HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,umac(i,j,k),umac(i,j,k), + d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity); if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { @@ -89,7 +90,8 @@ amrex::Real hydro_mol_xedge_state ( int i, int j, int k, int n, amrex::Real qpls = q(i ,j,k,n) - 0.5 * amrex_calc_xslope( i , j, k, n, order, q ); amrex::Real qmns = q(i-1,j,k,n) + 0.5 * amrex_calc_xslope( i-1, j, k, n, order, q ); - HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity); + HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,umac(i,j,k),umac(i,j,k), + d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity); if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { @@ -152,7 +154,8 @@ amrex::Real hydro_mol_yedge_state_extdir ( int i, int j, int k, int n, amrex::Real qmns = q(i,j-1,k,n) + 0.5 * amrex_calc_yslope_extdir( i, j-1, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity); + HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,vmac(i,j,k),vmac(i,j,k), + d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity); if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { @@ -198,7 +201,8 @@ amrex::Real hydro_mol_yedge_state ( int i, int j, int k, int n, amrex::Real qpls = q(i,j ,k,n) - 0.5 * amrex_calc_yslope( i, j , k, n, order, q ); amrex::Real qmns = q(i,j-1,k,n) + 0.5 * amrex_calc_yslope( i, j-1, k, n, order, q ); - HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity); + HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,vmac(i,j,k),vmac(i,j,k), + d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity); if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { @@ -263,7 +267,8 @@ amrex::Real hydro_mol_zedge_state_extdir ( int i, int j, int k, int n, amrex::Real qmns = q(i,j,k-1,n) + 0.5 * amrex_calc_zslope_extdir( i, j, k-1, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity); + HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,wmac(i,j,k),wmac(i,j,k), + d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity); if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { @@ -309,7 +314,8 @@ amrex::Real hydro_mol_zedge_state ( int i, int j, int k, int n, amrex::Real qpls = q(i,j,k ,n) - 0.5 * amrex_calc_zslope( i, j, k , n, order, q ); amrex::Real qmns = q(i,j,k-1,n) + 0.5 * amrex_calc_zslope( i, j, k-1, n, order, q ); - HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity); + HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,wmac(i,j,k),wmac(i,j,k), + d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity); if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { diff --git a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp index 66ae09a15..2204f7b89 100644 --- a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp +++ b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp @@ -87,7 +87,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real umns = vcc_mns + Real(0.5) * amrex_calc_xslope_extdir( i-1,j,k,0,order,vcc,extdir_or_ho_ilo, extdir_or_ho_ihi, domain_ilo, domain_ihi); - HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); + HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, umns, upls, + d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); if (!allow_inflow_on_outflow) { if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) @@ -136,7 +137,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real upls = vcc(i ,j,k,0) - Real(0.5) * amrex_calc_xslope(i ,j,k,0,order,vcc); Real umns = vcc(i-1,j,k,0) + Real(0.5) * amrex_calc_xslope(i-1,j,k,0,order,vcc); - HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); + HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, umns, upls, + d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true); if (!allow_inflow_on_outflow) { if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) @@ -198,7 +200,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real vmns = vcc_mns + Real(0.5) * amrex_calc_yslope_extdir( i,j-1,k,1,order,vcc,extdir_or_ho_jlo,extdir_or_ho_jhi,domain_jlo,domain_jhi); - HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); + HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, vmns, vpls, + d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); if (!allow_inflow_on_outflow) { if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) @@ -247,7 +250,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real vpls = vcc(i,j ,k,1) - Real(0.5) * amrex_calc_yslope(i,j ,k,1,order,vcc); Real vmns = vcc(i,j-1,k,1) + Real(0.5) * amrex_calc_yslope(i,j-1,k,1,order,vcc); - HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); + HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, vmns, vpls, + d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true); if (!allow_inflow_on_outflow) { if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) @@ -309,7 +313,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real wmns = vcc_mns + Real(0.5) * amrex_calc_zslope_extdir( i,j,k-1,2,order,vcc,extdir_or_ho_klo,extdir_or_ho_khi,domain_klo,domain_khi); - HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); + HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, wmns, wpls, + d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); if (!allow_inflow_on_outflow) { if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) @@ -358,7 +363,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Real wpls = vcc(i,j,k ,2) - Real(0.5) * amrex_calc_zslope(i,j,k ,2,order,vcc); Real wmns = vcc(i,j,k-1,2) + Real(0.5) * amrex_calc_zslope(i,j,k-1,2,order,vcc); - HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); + HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, wmns, wpls, + d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true); if (!allow_inflow_on_outflow) { if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) diff --git a/Utils/hydro_bcs_K.H b/Utils/hydro_bcs_K.H index 050d6f38e..6fd578e79 100644 --- a/Utils/hydro_bcs_K.H +++ b/Utils/hydro_bcs_K.H @@ -77,19 +77,32 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetXEdgeBCs( int i, int j, int k, int n, const amrex::Array4 &s, amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, int bclo, int domlo, int bchi, int domhi, bool is_velocity ) { - using namespace amrex; // This function is for setting BCs ON domain faces only. + // NOTE: this is called in two ways: + // 1) in the extrapolation of velocity components only -- + // in this case velm and velp may be predicted + // 2) for all components in the construction of edge states + // in this case velm and velp are identical and equal to + // the MAC velocity on the current edge AMREX_ASSERT( domlo <= i && i <= domhi+1 ); if (i == domlo) { - if (bclo == BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) ) + { + lo = s(domlo-1, j, k, n); + hi = s(domlo-1, j, k, n); + } + else if ( (bclo == BCType::ext_dir) || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo = s(domlo-1, j, k, n); // For turbulent inflow, there are times when the inflow face @@ -99,32 +112,41 @@ void SetXEdgeBCs( int i, int j, int k, int n, if ( n==XVEL && is_velocity ) hi=lo; } else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even ) + bclo == BCType::reflect_even || + (bclo == BCType::direction_dependent && velp < 0.0)) { lo = hi; } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; + hi = Real(0.); + lo = Real(0.); } } else if (i == domhi+1) { - if (bchi == BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == XVEL) && (s(domhi+1,j,k,n) <= 0.0) ) + { + hi = s(domhi+1, j, k, n); + lo = s(domhi+1, j, k, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { hi = s(domhi+1, j, k, n); if ( n==XVEL && is_velocity ) lo = hi; } else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even ) + bchi == BCType::reflect_even || + (bchi == BCType::direction_dependent && velm > 0.0)) { hi = lo; } else if (bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; + lo = Real(0.); + hi = Real(0.); } } else @@ -141,12 +163,13 @@ void SetXEdgeBCs( int i, int j, int k, int n, */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void SetYEdgeBCs ( int i, int j, int k, int n, - const amrex::Array4 &s, - amrex::Real &lo, amrex::Real &hi, - int bclo, int domlo, - int bchi, int domhi, - bool is_velocity ) +void SetYEdgeBCs (int i, int j, int k, int n, + const amrex::Array4 &s, + amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, + int bclo, int domlo, + int bchi, int domhi, + bool is_velocity ) { using namespace amrex; @@ -155,38 +178,54 @@ void SetYEdgeBCs ( int i, int j, int k, int n, if (j == domlo) { - if (bclo == BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == YVEL) && (s(i,domlo-1,k,n) >= 0.0) ) + { + lo = s(i, domlo-1, k, n); + hi = s(i, domlo-1, k, n); + } + else if ( (bclo == BCType::ext_dir) || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo = s(i, domlo-1, k, n); if ( n==YVEL && is_velocity ) hi=lo; } else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even ) + bclo == BCType::reflect_even || + (bclo == BCType::direction_dependent && velp < 0.0)) { lo = hi; } else if (bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; + hi = Real(0.); + lo = Real(0.); } } else if (j == domhi+1) { - if (bchi == BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == YVEL) && (s(i,domhi+1,k,n) <= 0.0) ) + { + hi = s(i, domhi+1, k, n); + lo = s(i, domhi+1, k, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { hi = s(i, domhi+1, k, n); if ( n==YVEL && is_velocity ) lo=hi; } else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even ) + bchi == BCType::reflect_even || + (bchi == BCType::direction_dependent && velm > 0.0)) { hi = lo; } - else if(bchi == BCType::reflect_odd) + else if (bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; + lo = Real(0.); + hi = Real(0.); } } else @@ -208,12 +247,13 @@ void SetYEdgeBCs ( int i, int j, int k, int n, */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void SetZEdgeBCs ( int i, int j, int k, int n, - const amrex::Array4 &s, - amrex::Real &lo, amrex::Real &hi, - int bclo, int domlo, - int bchi, int domhi, - bool is_velocity ) +void SetZEdgeBCs (int i, int j, int k, int n, + const amrex::Array4 &s, + amrex::Real &lo, amrex::Real &hi, + amrex::Real velm, amrex::Real velp, + int bclo, int domlo, + int bchi, int domhi, + bool is_velocity ) { using namespace amrex; @@ -222,38 +262,54 @@ void SetZEdgeBCs ( int i, int j, int k, int n, if (k == domlo) { - if (bclo == BCType::ext_dir) + if ( (bclo == BCType::direction_dependent) && is_velocity && + (n == ZVEL) && (s(i,j,domlo-1,n) >= 0.0) ) + { + lo = s(i, j, domlo-1, n); + hi = s(i, j, domlo-1, n); + } + else if ( (bclo == BCType::ext_dir) || + (bclo == BCType::direction_dependent && velp >= 0.0) ) { lo = s(i, j, domlo-1, n); if ( n==ZVEL && is_velocity ) hi=lo; } else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap || - bclo == BCType::reflect_even ) + bclo == BCType::reflect_even || + (bclo == BCType::direction_dependent && velp < 0.0)) { lo = hi; } else if(bclo == BCType::reflect_odd) { - hi = 0.; - lo = hi; + hi = Real(0.); + lo = Real(0.); } } else if (k == domhi+1) { - if (bchi == BCType::ext_dir) + if ( (bchi == BCType::direction_dependent) && is_velocity && + (n == ZVEL) && (s(i,j,domhi+1,n) <= 0.0) ) + { + hi = s(i, j, domhi+1, n); + lo = s(i, j, domhi+1, n); + } + else if ( (bchi == BCType::ext_dir) || + (bchi == BCType::direction_dependent && velm <= 0.0) ) { - hi = s(i,j,domhi+1, n); + hi = s(i, j, domhi+1, n); if ( n==ZVEL && is_velocity ) lo=hi; } else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap || - bchi == BCType::reflect_even ) + bchi == BCType::reflect_even || + (bchi == BCType::direction_dependent && velm > 0.0)) { hi = lo; } else if(bchi == BCType::reflect_odd) { - lo = 0.; - hi = lo; + lo = Real(0.); + hi = Real(0.); } } else