diff --git a/EBGodunov/hydro_ebgodunov.H b/EBGodunov/hydro_ebgodunov.H index 4aadbfff8..62efdac53 100644 --- a/EBGodunov/hydro_ebgodunov.H +++ b/EBGodunov/hydro_ebgodunov.H @@ -25,6 +25,7 @@ namespace EBGodunov { const amrex::Geometry& geom, amrex::Real dt, amrex::MultiFab const* velocity_on_eb_inflow = nullptr, + bool allow_inflow_on_outflow = false, amrex::iMultiFab const* BC_MF = nullptr); @@ -82,6 +83,7 @@ namespace EBGodunov { amrex::Array4 const& fcy, amrex::Array4 const& fcz), amrex::Real* p, amrex::Array4 const& velocity_on_eb_inflow, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); @@ -112,6 +114,7 @@ namespace EBGodunov { amrex::Array4 const& ccent_arr, bool is_velocity, amrex::Array4 const& values_on_eb_inflow, + bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr = {}); } // namespace ebgodunov diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp index 97d890970..c5b50050f 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp @@ -37,6 +37,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Array4 const& ccent_arr, bool is_velocity, Array4 const& values_on_eb_inflow, + bool allow_inflow_on_outflow, Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); @@ -222,15 +223,17 @@ 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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - if ( u_mac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - if ( u_mac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + if ( u_mac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + if ( u_mac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (u_mac(i,j,k) >= 0.) ? stl : sth; @@ -338,16 +341,19 @@ 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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - if ( v_mac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - if ( v_mac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + if ( v_mac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + if ( v_mac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } + Real temp = (v_mac(i,j,k) >= 0.) ? stl : sth; temp = (amrex::Math::abs(v_mac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; yedge(i,j,k,n) = temp; diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp index 57c7d13c7..ccb97d446 100644 --- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp @@ -40,6 +40,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp, Array4 const& ccent_arr, bool is_velocity, Array4 const& values_on_eb_inflow, + bool allow_inflow_on_outflow, Array4 const& bc_arr) { @@ -306,15 +307,17 @@ 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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - if ( u_mac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - if ( u_mac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + if ( u_mac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + if ( u_mac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (u_mac(i,j,k) >= 0.) ? stl : sth; @@ -443,16 +446,19 @@ 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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - if ( v_mac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - if ( v_mac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + if ( v_mac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + if ( v_mac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } + Real temp = (v_mac(i,j,k) >= 0.) ? stl : sth; temp = (amrex::Math::abs(v_mac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; yedge(i,j,k,n) = temp; @@ -577,16 +583,19 @@ 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); - if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) - { - if ( w_mac(i,j,k) >= 0. && n==ZVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) - { - if ( w_mac(i,j,k) <= 0. && n==ZVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) + { + if ( w_mac(i,j,k) >= 0. && n==ZVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) + { + if ( w_mac(i,j,k) <= 0. && n==ZVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } + Real temp = (w_mac(i,j,k) >= 0.) ? stl : sth; temp = (amrex::Math::abs(w_mac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; zedge(i,j,k,n) = temp; diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp index 7cb78f695..64c8668c8 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp @@ -24,6 +24,7 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, BCRec const* d_bcrec, const Geometry& geom, Real l_dt, MultiFab const* velocity_on_eb_inflow, + bool allow_inflow_on_outflow, iMultiFab const* BC_MF) { BL_PROFILE("EBGodunov::ExtrapVelToFaces()"); @@ -175,7 +176,7 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, AMREX_D_DECL(Ipx, Ipy, Ipz), a_f, domain, dx, l_dt, d_bcrec, local_use_forces_in_trans, p, - bc_arr); + allow_inflow_on_outflow,bc_arr); } else { @@ -233,6 +234,7 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel, p, velocity_on_eb_inflow ? velocity_on_eb_inflow->const_array(mfi) : Array4{}, + allow_inflow_on_outflow, bc_arr); } diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp index 2d97a4a54..a38f82e3d 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp @@ -35,6 +35,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp, Array4 const& fcy, Real* p, Array4 const& velocity_on_eb_inflow, + bool allow_inflow_on_outflow, Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); @@ -190,16 +191,17 @@ 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); - // Prevent backflow - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; bool ltm = ( (stl <= 0. && sth >= 0.) || (amrex::Math::abs(stl+sth) < small_vel) ); @@ -312,16 +314,17 @@ 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); - // Prevent backflow - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp index ae03fd511..d5145961a 100644 --- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp +++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp @@ -48,6 +48,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, Array4 const& fcz, Real* p, Array4 const& velocity_on_eb_inflow, + bool allow_inflow_on_outflow, Array4 const& bc_arr) { const Dim3 dlo = amrex::lbound(domain); @@ -267,16 +268,17 @@ 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); - // Prevent backflow - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; @@ -419,16 +421,17 @@ 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); - // Prevent backflow - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; @@ -574,17 +577,17 @@ 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); - - // Prevent backflow - if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; diff --git a/Godunov/hydro_godunov.H b/Godunov/hydro_godunov.H index d9e88c687..127722c5c 100644 --- a/Godunov/hydro_godunov.H +++ b/Godunov/hydro_godunov.H @@ -25,6 +25,7 @@ void ExtrapVelToFaces ( amrex::MultiFab const& a_vel, const amrex::Geometry& geom, amrex::Real l_dt, bool use_ppm, bool use_forces_in_trans, int limiter_type = PPM::VanLeer, + bool allow_inflow_on_outflow = false, amrex::iMultiFab const* BC_MF = nullptr); void ComputeAdvectiveVel (AMREX_D_DECL(amrex::Box const& xbx, @@ -71,6 +72,7 @@ void ExtrapVelToFacesOnBox (amrex::Box const& bx, int ncomp, amrex::BCRec const* pbc, bool use_forces_in_trans, amrex::Real* p, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); void ComputeEdgeState ( amrex::Box const& bx, int ncomp, @@ -90,8 +92,8 @@ void ComputeEdgeState ( amrex::Box const& bx, int ncomp, bool use_ppm, bool use_forces_in_trans, bool is_velocity, int limiter_type = PPM::VanLeer, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); - } #endif diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp index 2c77a4e44..68748c315 100644 --- a/Godunov/hydro_godunov_edge_state_2D.cpp +++ b/Godunov/hydro_godunov_edge_state_2D.cpp @@ -29,6 +29,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, const bool use_forces_in_trans, const bool is_velocity, const int limiter_type, + const bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); @@ -229,15 +230,17 @@ 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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (umac(i,j,k) >= 0.) ? stl : sth; @@ -313,15 +316,17 @@ 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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (vmac(i,j,k) >= 0.) ? stl : sth; diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp index af7d3c094..6a573c599 100644 --- a/Godunov/hydro_godunov_edge_state_3D.cpp +++ b/Godunov/hydro_godunov_edge_state_3D.cpp @@ -31,6 +31,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, const bool use_forces_in_trans, const bool is_velocity, const int limiter_type, + const bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { Box const& xbx = amrex::surroundingNodes(bx,0); @@ -303,15 +304,17 @@ 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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (umac(i,j,k) >= 0.) ? stl : sth; @@ -402,15 +405,17 @@ 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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (vmac(i,j,k) >= 0.) ? stl : sth; @@ -473,10 +478,10 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // --> 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)) ); + -(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 @@ -497,20 +502,20 @@ 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); - if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) - { - if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) - { - if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) + { + if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) + { + if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real temp = (wmac(i,j,k) >= 0.) ? stl : sth; diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp index d4d6ef22f..0f60235ff 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp @@ -22,6 +22,7 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, const Geometry& geom, Real l_dt, const bool use_ppm, const bool use_forces_in_trans, const int limiter_type, + const bool allow_inflow_on_outflow, iMultiFab const* BC_MF) { Box const& domain = geom.Domain(); @@ -111,7 +112,7 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, u_ad, v_ad, Imx, Imy, Ipx, Ipy, f, domain, dx, l_dt, d_bcrec, use_forces_in_trans, p, - bc_arr); + allow_inflow_on_outflow, bc_arr); Gpu::streamSynchronize(); // otherwise we might be using too much memory } @@ -204,6 +205,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, BCRec const* pbc, bool l_use_forces_in_trans, Real* p, + bool allow_inflow_on_outflow, Array4 const& bc_arr) { @@ -306,15 +308,17 @@ Godunov::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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; @@ -365,15 +369,17 @@ Godunov::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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp index 920f19995..6de8c3bf2 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp @@ -24,6 +24,7 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_vel, const Geometry& geom, Real l_dt, bool use_ppm, bool use_forces_in_trans, const int limiter_type, + const bool allow_inflow_on_outflow, amrex::iMultiFab const* BC_MF) { Box const& domain = geom.Domain(); @@ -125,7 +126,7 @@ Godunov::ExtrapVelToFaces ( MultiFab const& a_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, - bc_arr); + allow_inflow_on_outflow, bc_arr); Gpu::streamSynchronize(); // otherwise we might be using too much memory } @@ -248,6 +249,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp, BCRec const* pbc, bool l_use_forces_in_trans, Real* p, + const bool allow_inflow_on_outflow, Array4 const& bc_arr) { @@ -432,15 +434,17 @@ Godunov::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); - if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) - { - sth = amrex::min(sth,0.0_rt); - stl = sth; - } - if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) - { - stl = amrex::max(stl,0.0_rt); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) + { + sth = amrex::min(sth,0.0_rt); + stl = sth; + } + if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) + { + stl = amrex::max(stl,0.0_rt); + sth = stl; + } } Real st = ( (stl+sth) >= Real(0.0)) ? stl : sth; @@ -520,15 +524,17 @@ Godunov::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); - if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) - { - sth = amrex::min(sth,Real(0.0)); - stl = sth; - } - if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) - { - stl = amrex::max(stl,Real(0.0)); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) + { + sth = amrex::min(sth,Real(0.0)); + stl = sth; + } + if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) + { + stl = amrex::max(stl,Real(0.0)); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; @@ -611,15 +617,17 @@ Godunov::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); - if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) - { - sth = amrex::min(sth,Real(0.0)); - stl = sth; - } - if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) - { - stl = amrex::max(stl,Real(0.0)); - sth = stl; + if (!allow_inflow_on_outflow) { + if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) + { + sth = amrex::min(sth,Real(0.0)); + stl = sth; + } + if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) + { + stl = amrex::max(stl,Real(0.0)); + sth = stl; + } } Real st = ( (stl+sth) >= 0.) ? stl : sth; diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index 71ac154a8..83e4a0c45 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -737,15 +737,17 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx; - if (vel_edge(i+1,j,k) > small_vel) + if (vel_edge(i+1,j,k) > small_vel) { Ip = sp - (half*sigmap)*((sp - sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - else + } else { Ip = S(i,j,k,n); + } - if(vel_edge(i,j,k) < -small_vel) + if(vel_edge(i,j,k) < -small_vel) { Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); - else + } else { Im = S(i,j,k,n); + } } template @@ -781,15 +783,17 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx; - if (vel_edge(i,j+1,k) > small_vel) + if (vel_edge(i,j+1,k) > small_vel) { Ip = sp - (half*sigmap)*((sp - sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - else + } else { Ip = S(i,j,k,n); + } - if (vel_edge(i,j,k) < -small_vel) + if (vel_edge(i,j,k) < -small_vel) { Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); - else + } else { Im = S(i,j,k,n); + } } @@ -827,15 +831,17 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx; amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx; - if(vel_edge(i,j,k+1) > small_vel) + if(vel_edge(i,j,k+1) > small_vel) { Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); - else + } else { Ip = S(i,j,k,n); + } - if(vel_edge(i,j,k) < -small_vel) + if(vel_edge(i,j,k) < -small_vel) { Im = sm + (half*sigmam)*((sp-sm) + (amrex::Real(1.0) - two3rds*sigmam)*s6); - else + } else { Im = S(i,j,k,n); + } } #endif diff --git a/MOL/hydro_mol.H b/MOL/hydro_mol.H index 603109427..17504b05c 100644 --- a/MOL/hydro_mol.H +++ b/MOL/hydro_mol.H @@ -97,7 +97,8 @@ void ExtrapVelToFaces ( const amrex::MultiFab& vel, amrex::MultiFab& wmac ), const amrex::Geometry& a_geom, amrex::Vector const& h_bcrec, - const amrex::BCRec* d_bcrec ); + const amrex::BCRec* d_bcrec, + bool allow_inflow_on_outflow = false); /** * \brief For Computing the pre-MAC edge states to be MAC-projected. @@ -115,7 +116,8 @@ void ExtrapVelToFacesBox ( AMREX_D_DECL( amrex::Box const& ubx, amrex::Array4 const& vcc, const amrex::Geometry& geom, amrex::Vector const& h_bcrec, - const amrex::BCRec* d_bcrec ); + const amrex::BCRec* d_bcrec, + bool allow_inflow_on_outflow = false); } diff --git a/MOL/hydro_mol_extrap_vel_to_faces.cpp b/MOL/hydro_mol_extrap_vel_to_faces.cpp index ec1289040..10040eb3f 100644 --- a/MOL/hydro_mol_extrap_vel_to_faces.cpp +++ b/MOL/hydro_mol_extrap_vel_to_faces.cpp @@ -21,7 +21,8 @@ MOL::ExtrapVelToFaces ( const MultiFab& a_vel, MultiFab& a_wmac ), const Geometry& a_geom, const Vector& h_bcrec, - BCRec const* d_bcrec) + BCRec const* d_bcrec, + bool allow_inflow_on_outflow) { BL_PROFILE("MOL::ExtrapVelToFaces"); @@ -42,7 +43,8 @@ MOL::ExtrapVelToFaces ( const MultiFab& a_vel, Array4 const& vcc = a_vel.const_array(mfi); ExtrapVelToFacesBox( AMREX_D_DECL(ubx,vbx,wbx), AMREX_D_DECL(u,v,w), - vcc,a_geom,h_bcrec, d_bcrec); + vcc,a_geom,h_bcrec, d_bcrec, + allow_inflow_on_outflow); } } diff --git a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp index 1355b2cd3..66ae09a15 100644 --- a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp +++ b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp @@ -27,8 +27,6 @@ namespace { } } - - void MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Box const& vbx, @@ -39,7 +37,8 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, Array4 const& vcc, const Geometry& geom, Vector const& h_bcrec, - BCRec const* d_bcrec) + BCRec const* d_bcrec, + bool allow_inflow_on_outflow) { BL_PROFILE("MOL::ExtrapVelToFacesBox"); @@ -69,7 +68,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, if ((has_extdir_or_ho_lo && domain_ilo >= ubx.smallEnd(0)-1) || (has_extdir_or_ho_hi && domain_ihi <= ubx.bigEnd(0))) { - amrex::ParallelFor(ubx, [vcc,domain_ilo,domain_ihi,u,d_bcrec] + amrex::ParallelFor(ubx, [vcc,domain_ilo,domain_ihi,u,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool extdir_or_ho_ilo = (d_bcrec[0].lo(0) == BCType::ext_dir) || @@ -90,15 +89,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetXEdgeBCs(i, j, k, n, vcc, 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) ) - { - upls = amrex::min(upls,Real(0.0)); - umns = upls; - } - if ( (i==domain_ihi+1) && (d_bcrec[0].hi(0) == BCType::foextrap || d_bcrec[0].hi(0) == BCType::hoextrap) ) - { - umns = amrex::max(umns,Real(0.0)); - upls = umns; + if (!allow_inflow_on_outflow) { + if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) + { + upls = amrex::min(upls,Real(0.0)); + umns = upls; + } + if ( (i==domain_ihi+1) && (d_bcrec[0].hi(0) == BCType::foextrap || d_bcrec[0].hi(0) == BCType::hoextrap) ) + { + umns = amrex::max(umns,Real(0.0)); + upls = umns; + } } Real u_val(0); @@ -127,7 +128,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, } else { - amrex::ParallelFor(ubx, [vcc,domain_ilo,domain_ihi,u,d_bcrec] + amrex::ParallelFor(ubx, [vcc,domain_ilo,domain_ihi,u,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; @@ -137,15 +138,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetXEdgeBCs(i, j, k, n, vcc, 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) ) - { - upls = amrex::min(upls,Real(0.0)); - umns = upls; - } - if ( (i==domain_ihi+1) && (d_bcrec[0].hi(0) == BCType::foextrap || d_bcrec[0].hi(0) == BCType::hoextrap) ) - { - umns = amrex::max(umns,Real(0.0)); - upls = umns; + if (!allow_inflow_on_outflow) { + if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) ) + { + upls = amrex::min(upls,Real(0.0)); + umns = upls; + } + if ( (i==domain_ihi+1) && (d_bcrec[0].hi(0) == BCType::foextrap || d_bcrec[0].hi(0) == BCType::hoextrap) ) + { + umns = amrex::max(umns,Real(0.0)); + upls = umns; + } } Real u_val(0); @@ -177,7 +180,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, if ((has_extdir_or_ho_lo && domain_jlo >= vbx.smallEnd(1)-1) || (has_extdir_or_ho_hi && domain_jhi <= vbx.bigEnd(1))) { - amrex::ParallelFor(vbx, [vcc,domain_jlo,domain_jhi,v,d_bcrec] + amrex::ParallelFor(vbx, [vcc,domain_jlo,domain_jhi,v,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool extdir_or_ho_jlo = (d_bcrec[1].lo(1) == BCType::ext_dir) || @@ -197,15 +200,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetYEdgeBCs(i, j, k, n, vcc, 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) ) - { - vpls = amrex::min(vpls,Real(0.0)); - vmns = vpls; - } - if ( (j==domain_jhi+1) && (d_bcrec[1].hi(1) == BCType::foextrap || d_bcrec[1].hi(1) == BCType::hoextrap) ) - { - vmns = amrex::max(vmns,Real(0.0)); - vpls = vmns; + if (!allow_inflow_on_outflow) { + if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) + { + vpls = amrex::min(vpls,Real(0.0)); + vmns = vpls; + } + if ( (j==domain_jhi+1) && (d_bcrec[1].hi(1) == BCType::foextrap || d_bcrec[1].hi(1) == BCType::hoextrap) ) + { + vmns = amrex::max(vmns,Real(0.0)); + vpls = vmns; + } } Real v_val(0); @@ -234,7 +239,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, } else { - amrex::ParallelFor(vbx, [vcc,domain_jlo,domain_jhi,v,d_bcrec] + amrex::ParallelFor(vbx, [vcc,domain_jlo,domain_jhi,v,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 1; @@ -244,15 +249,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetYEdgeBCs(i, j, k, n, vcc, 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) ) - { - vpls = amrex::min(vpls,Real(0.0)); - vmns = vpls; - } - if ( (j==domain_jhi+1) && (d_bcrec[1].hi(1) == BCType::foextrap || d_bcrec[1].hi(1) == BCType::hoextrap) ) - { - vmns = amrex::max(vmns,Real(0.0)); - vpls = vmns; + if (!allow_inflow_on_outflow) { + if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) ) + { + vpls = amrex::min(vpls,Real(0.0)); + vmns = vpls; + } + if ( (j==domain_jhi+1) && (d_bcrec[1].hi(1) == BCType::foextrap || d_bcrec[1].hi(1) == BCType::hoextrap) ) + { + vmns = amrex::max(vmns,Real(0.0)); + vpls = vmns; + } } Real v_val(0); @@ -284,7 +291,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, if ((has_extdir_or_ho_lo && domain_klo >= wbx.smallEnd(2)-1) || (has_extdir_or_ho_hi && domain_khi <= wbx.bigEnd(2))) { - amrex::ParallelFor(wbx, [vcc,domain_klo,domain_khi,w,d_bcrec] + amrex::ParallelFor(wbx, [vcc,domain_klo,domain_khi,w,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool extdir_or_ho_klo = (d_bcrec[2].lo(2) == BCType::ext_dir) || @@ -304,15 +311,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetZEdgeBCs(i, j, k, n, vcc, 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) ) - { - wpls = amrex::min(wpls,Real(0.0)); - wmns = wpls; - } - if ( (k==domain_khi+1) && (d_bcrec[2].hi(2) == BCType::foextrap || d_bcrec[2].hi(2) == BCType::hoextrap) ) - { - wmns = amrex::max(wmns,Real(0.0)); - wpls = wmns; + if (!allow_inflow_on_outflow) { + if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) + { + wpls = amrex::min(wpls,Real(0.0)); + wmns = wpls; + } + if ( (k==domain_khi+1) && (d_bcrec[2].hi(2) == BCType::foextrap || d_bcrec[2].hi(2) == BCType::hoextrap) ) + { + wmns = amrex::max(wmns,Real(0.0)); + wpls = wmns; + } } Real w_val(0); @@ -341,7 +350,7 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, } else { - amrex::ParallelFor(wbx, [vcc,domain_klo,domain_khi,w,d_bcrec] + amrex::ParallelFor(wbx, [vcc,domain_klo,domain_khi,w,d_bcrec,allow_inflow_on_outflow] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 2; @@ -351,15 +360,17 @@ MOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx, HydroBC::SetZEdgeBCs(i, j, k, n, vcc, 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) ) - { - wpls = amrex::min(wpls,Real(0.0)); - wmns = wpls; - } - if ( (k==domain_khi+1) && (d_bcrec[2].hi(2) == BCType::foextrap || d_bcrec[2].hi(2) == BCType::hoextrap) ) - { - wmns = amrex::max(wmns,Real(0.0)); - wpls = wmns; + if (!allow_inflow_on_outflow) { + if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) ) + { + wpls = amrex::min(wpls,Real(0.0)); + wmns = wpls; + } + if ( (k==domain_khi+1) && (d_bcrec[2].hi(2) == BCType::foextrap || d_bcrec[2].hi(2) == BCType::hoextrap) ) + { + wmns = amrex::max(wmns,Real(0.0)); + wpls = wmns; + } } Real w_val(0); diff --git a/Utils/hydro_compute_edgestate_and_flux.cpp b/Utils/hydro_compute_edgestate_and_flux.cpp index 726696e7d..3d49e092a 100644 --- a/Utils/hydro_compute_edgestate_and_flux.cpp +++ b/Utils/hydro_compute_edgestate_and_flux.cpp @@ -45,6 +45,7 @@ namespace { bool is_velocity, std::string const& advection_type, int limiter_type, + bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { @@ -91,7 +92,9 @@ namespace { AMREX_D_DECL(apx,apy,apz), vfrac, AMREX_D_DECL(fcx,fcy,fcz), ccc, is_velocity, - values_on_eb_inflow, bc_arr); + values_on_eb_inflow, + allow_inflow_on_outflow, + bc_arr); } else if (advection_type == "BDS") { @@ -123,7 +126,7 @@ namespace { geom, l_dt, d_bcrec, iconserv, godunov_use_ppm, godunov_use_forces_in_trans, - is_velocity, limiter_type, bc_arr); + is_velocity, limiter_type, allow_inflow_on_outflow, bc_arr); } else if (advection_type == "BDS") { @@ -167,6 +170,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type, + bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { @@ -179,7 +183,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, bc_arr); + limiter_type, allow_inflow_on_outflow, bc_arr); } #endif @@ -212,6 +216,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type, + bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { @@ -227,7 +232,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, bc_arr); + limiter_type, allow_inflow_on_outflow, bc_arr); } @@ -261,6 +266,7 @@ HydroUtils::ComputeFluxesOnBoxFromState (Box const& bx, int ncomp, MFIter& mfi, bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type, + bool allow_inflow_on_outflow, amrex::Array4 const& bc_arr) { @@ -288,7 +294,8 @@ 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, bc_arr); + is_velocity, advection_type, limiter_type, + allow_inflow_on_outflow, bc_arr); } // Compute fluxes. diff --git a/Utils/hydro_extrap_vel_to_faces.cpp b/Utils/hydro_extrap_vel_to_faces.cpp index 6cdc7a9e1..f6cac52e0 100644 --- a/Utils/hydro_extrap_vel_to_faces.cpp +++ b/Utils/hydro_extrap_vel_to_faces.cpp @@ -27,13 +27,14 @@ HydroUtils::ExtrapVelToFaces ( amrex::MultiFab const& vel, const EBFArrayBoxFactory& ebfact, bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, - int limiter_type) + int limiter_type, + bool allow_inflow_on_outflow) { ExtrapVelToFaces(vel, vel_forces, AMREX_D_DECL(u_mac,v_mac,w_mac), h_bcrec, d_bcrec, geom, dt, ebfact, /*velocity_on_eb_inflow*/ nullptr, godunov_ppm, godunov_use_forces_in_trans, - advection_type, limiter_type); + advection_type, limiter_type, allow_inflow_on_outflow); } #endif @@ -54,6 +55,7 @@ HydroUtils::ExtrapVelToFaces ( amrex::MultiFab const& vel, bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, int limiter_type, + bool allow_inflow_on_outflow, iMultiFab* BC_MF) { amrex::ignore_unused(BC_MF); @@ -72,7 +74,8 @@ HydroUtils::ExtrapVelToFaces ( amrex::MultiFab const& vel, Godunov::ExtrapVelToFaces(vel, vel_forces, AMREX_D_DECL(u_mac, v_mac, w_mac), h_bcrec, d_bcrec, - geom, dt, godunov_ppm, godunov_use_forces_in_trans, limiter_type); + geom, dt, godunov_ppm, godunov_use_forces_in_trans, + limiter_type, allow_inflow_on_outflow); } else if (advection_type == "MOL") { @@ -81,7 +84,7 @@ HydroUtils::ExtrapVelToFaces ( amrex::MultiFab const& vel, EBMOL::ExtrapVelToFaces(vel, AMREX_D_DECL(u_mac, v_mac, w_mac), geom, h_bcrec, d_bcrec); else #endif - MOL::ExtrapVelToFaces(vel, AMREX_D_DECL(u_mac, v_mac, w_mac), geom, h_bcrec, d_bcrec); + MOL::ExtrapVelToFaces(vel, AMREX_D_DECL(u_mac, v_mac, w_mac), geom, h_bcrec, d_bcrec, allow_inflow_on_outflow); } else { amrex::Abort("Dont know this advection_type in HydroUtils::ExtrapVelToFaces"); } diff --git a/Utils/hydro_utils.H b/Utils/hydro_utils.H index af966dfb7..7111e109e 100644 --- a/Utils/hydro_utils.H +++ b/Utils/hydro_utils.H @@ -66,6 +66,7 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type = PPM::default_limiter, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); /** * \brief Compute edge state and flux. For typical advection, and also allows for inflow on EB. @@ -99,6 +100,7 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type = PPM::default_limiter, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); /** @@ -131,6 +133,7 @@ ComputeFluxesOnBoxFromState ( amrex::Box const& bx, int ncomp, amrex::MFIter& mf bool is_velocity, bool fluxes_are_area_weighted, std::string const& advection_type, int limiter_type = PPM::default_limiter, + bool allow_inflow_on_outflow = false, amrex::Array4 const& bc_arr = {}); #endif @@ -148,7 +151,8 @@ 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::default_limiter); + int limiter_type = PPM::default_limiter, + bool allow_inflow_on_outflow = false); #endif void @@ -168,6 +172,7 @@ ExtrapVelToFaces ( amrex::MultiFab const& vel, bool godunov_ppm, bool godunov_use_forces_in_trans, std::string const& advection_type, int limiter_type = PPM::default_limiter, + bool allow_inflow_on_outflow = false, amrex::iMultiFab* BC_MF = nullptr); /**