Skip to content

Commit

Permalink
add option for direction-dependent bc's (#136)
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren authored Aug 8, 2024
1 parent 555ba57 commit 186df80
Show file tree
Hide file tree
Showing 17 changed files with 448 additions and 235 deletions.
118 changes: 83 additions & 35 deletions EBGodunov/hydro_ebgodunov_bcs_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<const amrex::Real> &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 )
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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);
Expand All @@ -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;
}
Expand All @@ -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);
Expand All @@ -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<const amrex::Real> &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 )
Expand All @@ -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;
}
Expand All @@ -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);
Expand All @@ -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;
}
Expand All @@ -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);
Expand All @@ -201,26 +233,34 @@ 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<const amrex::Real> &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<const amrex::Real> &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;


// 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;
}
Expand All @@ -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);
Expand All @@ -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;
}
Expand All @@ -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);
Expand Down
22 changes: 14 additions & 8 deletions EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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);
Expand Down Expand Up @@ -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) )
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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) )
Expand Down
Loading

0 comments on commit 186df80

Please sign in to comment.