Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update AMReX-Hydro routines to take allow_inflow_on_outflow flag; thi… #137

Merged
merged 2 commits into from
Aug 1, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions EBGodunov/hydro_ebgodunov.H
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
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
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
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
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
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
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,
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