Skip to content

Commit

Permalink
more changes to support direction_dependent bc
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Aug 6, 2024
1 parent 836cba0 commit ecd625f
Show file tree
Hide file tree
Showing 9 changed files with 174 additions and 75 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
2 changes: 2 additions & 0 deletions EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
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 @@ -267,6 +268,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
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
1 change: 1 addition & 0 deletions EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -354,6 +354,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
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);
Expand Down
9 changes: 6 additions & 3 deletions EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) );
Expand All @@ -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) );
Expand All @@ -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) );
Expand Down
13 changes: 8 additions & 5 deletions EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -111,7 +113,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,

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;
Expand Down
1 change: 1 addition & 0 deletions Godunov/hydro_godunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,6 +270,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
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);
Expand Down
Loading

0 comments on commit ecd625f

Please sign in to comment.