Skip to content

Commit

Permalink
Update AMReX-Hydro routines to take allow_inflow_on_outflow flag; thi…
Browse files Browse the repository at this point in the history 
…s defaults to false but allows the user to override
  • Loading branch information
@asalmgren
asalmgren committed Aug 1, 2024
1 parent aca1032 commit f541dce
Show file tree
Hide file tree
Showing 18 changed files with 373 additions and 285 deletions.
3 changes: 3 additions & 0 deletions EBGodunov/hydro_ebgodunov.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);


Expand Down Expand Up @@ -82,6 +83,7 @@ namespace EBGodunov {
amrex::Array4<amrex::Real const> const& fcy,
amrex::Array4<amrex::Real const> const& fcz),
amrex::Real* p, amrex::Array4<amrex::Real const> const& velocity_on_eb_inflow,
bool allow_inflow_on_outflow = false,
amrex::Array4<int const> const& bc_arr = {});


Expand Down Expand Up @@ -112,6 +114,7 @@ namespace EBGodunov {
amrex::Array4<amrex::Real const> const& ccent_arr,
bool is_velocity,
amrex::Array4<amrex::Real const> const& values_on_eb_inflow,
bool allow_inflow_on_outflow,
amrex::Array4<int const> const& bc_arr = {});

} // namespace ebgodunov
Expand Down
42 changes: 24 additions & 18 deletions EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Array4<Real const> const& ccent_arr,
bool is_velocity,
Array4<Real const> const& values_on_eb_inflow,
bool allow_inflow_on_outflow,
Array4<int const> const& bc_arr)
{
Box const& xbx = amrex::surroundingNodes(bx,0);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
63 changes: 36 additions & 27 deletions EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Array4<Real const> const& ccent_arr,
bool is_velocity,
Array4<Real const> const& values_on_eb_inflow,
bool allow_inflow_on_outflow,
Array4<int const> const& bc_arr)
{

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 3 additions & 1 deletion EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()");
Expand Down Expand Up @@ -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
{
Expand Down Expand Up @@ -233,6 +234,7 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
p,
velocity_on_eb_inflow ?
velocity_on_eb_inflow->const_array(mfi) : Array4<Real const>{},
allow_inflow_on_outflow,
bc_arr);
}

Expand Down
43 changes: 23 additions & 20 deletions EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
Array4<Real const> const& fcy,
Real* p,
Array4<Real const> const& velocity_on_eb_inflow,
bool allow_inflow_on_outflow,
Array4<int const> const& bc_arr)
{
const Dim3 dlo = amrex::lbound(domain);
Expand Down Expand Up @@ -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) );
Expand Down Expand Up @@ -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;
Expand Down
65 changes: 34 additions & 31 deletions EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
Array4<Real const> const& fcz,
Real* p,
Array4<Real const> const& velocity_on_eb_inflow,
bool allow_inflow_on_outflow,
Array4<int const> const& bc_arr)
{
const Dim3 dlo = amrex::lbound(domain);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
4 changes: 3 additions & 1 deletion Godunov/hydro_godunov.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -71,6 +72,7 @@ void ExtrapVelToFacesOnBox (amrex::Box const& bx, int ncomp,
amrex::BCRec const* pbc,
bool use_forces_in_trans,
amrex::Real* p,
const bool allow_inflow_on_outflow = false,
amrex::Array4<int const> const& bc_arr = {});

void ComputeEdgeState ( amrex::Box const& bx, int ncomp,
Expand All @@ -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<int const> const& bc_arr = {});

}

#endif
Expand Down
Loading

0 comments on commit f541dce

Please sign in to comment.