Skip to content

Commit

Permalink
Add support for mixed bcs in (EB)godunov PLM and MacProjector (#123)
Browse files Browse the repository at this point in the history 
For the MacProjector, allow specifying Robin BCs (since mixed BCs can be
viewed as a special case of this).

For (EB)Godunov, this has been implemented by allowing the option to
pass a position dependent BC MF/Array4. The BC MF is cell-centered and
carries the amrex::BCType type in the first ghost cell (precedent set by
MLMG's treatment of Robin BCs) and must fully specify the BC on all
faces. If a position dependent BCs are passed in, they take precedent
and single bc per face BCRecs are ignored.
  • Loading branch information
@cgilet
cgilet authored Feb 28, 2024
1 parent 2b89ad7 commit 1ec1247
Show file tree
Hide file tree
Showing 24 changed files with 403 additions and 294 deletions.
13 changes: 9 additions & 4 deletions EBGodunov/hydro_ebgodunov.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,9 @@ namespace EBGodunov {
amrex::Vector<amrex::BCRec> const& h_bcrec,
amrex::BCRec const* d_bcrec,
const amrex::Geometry& geom,
amrex::Real dt, amrex::MultiFab const* velocity_on_eb_inflow = nullptr);
amrex::Real dt,
amrex::MultiFab const* velocity_on_eb_inflow = nullptr,
amrex::iMultiFab const* BC_MF = nullptr);


void ComputeAdvectiveVel (AMREX_D_DECL(amrex::Box const& xbx,
Expand All @@ -41,7 +43,8 @@ namespace EBGodunov {
amrex::Array4<amrex::Real const> const& vel,
amrex::Array4<amrex::EBCellFlag const> const& flag,
const amrex::Box& domain,
amrex::BCRec const* pbc);
amrex::BCRec const* pbc,
amrex::Array4<int const> const& bc_arr = {});

void ExtrapVelToFacesOnBox ( amrex::Box const& bx, int ncomp,
AMREX_D_DECL(amrex::Box const& xbx,
Expand Down Expand Up @@ -78,7 +81,8 @@ namespace EBGodunov {
AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
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);
amrex::Real* p, amrex::Array4<amrex::Real const> const& velocity_on_eb_inflow,
amrex::Array4<int const> const& bc_arr = {});


void ComputeEdgeState ( amrex::Box const& bx, int ncomp,
Expand Down Expand Up @@ -107,7 +111,8 @@ namespace EBGodunov {
amrex::Array4<amrex::Real const> const& fcz),
amrex::Array4<amrex::Real const> const& ccent_arr,
bool is_velocity,
amrex::Array4<amrex::Real const> const& values_on_eb_inflow);
amrex::Array4<amrex::Real const> const& values_on_eb_inflow,
amrex::Array4<int const> const& bc_arr = {});

} // namespace ebgodunov

Expand Down
2 changes: 2 additions & 0 deletions EBGodunov/hydro_ebgodunov_bcs_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,8 @@ void SetXBCs (const int i, const int j, const int k, const int n,
// upwinded intermediate edge state in Imx (which supplies hi)
// With GPU, I think it's undertermined which Imx(i+..) you'll get here,
// could be from PLM, could be the upwinded edgestate.
// Question is whether every thread in the stream has to finish a kernel
// before the next kernel is launched.
// lo = hi_arr(2*domlo-i ,j,k,n);
// hi = lo_arr(2*domlo-i-1,j,k,n);
}
Expand Down
15 changes: 8 additions & 7 deletions EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Array4<Real const> const& fcy,
Array4<Real const> const& ccent_arr,
bool is_velocity,
Array4<Real const> const& values_on_eb_inflow)
Array4<Real const> const& values_on_eb_inflow,
Array4<int const> const& bc_arr)
{
Box const& xbx = amrex::surroundingNodes(bx,0);
Box const& ybx = amrex::surroundingNodes(bx,1);
Expand Down Expand Up @@ -106,7 +107,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Real lo = Ipx(i-1,j,k,n);
Real hi = Imx(i ,j,k,n);

auto bc = pbc[n];
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);

Expand All @@ -118,7 +119,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Real lo = Ipy(i,j-1,k,n);
Real hi = Imy(i,j ,k,n);

auto bc = pbc[n];
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);

Expand All @@ -136,7 +137,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
{
if (apy(i,j,k) > 0.)
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_yzlo, l_yzhi;

l_yzlo = ylo(i,j,k,n);
Expand Down Expand Up @@ -218,7 +219,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
sth += (fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.;
}

auto bc = pbc[n];
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) )
Expand Down Expand Up @@ -251,7 +252,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
{
if (apx(i,j,k) > 0.)
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_xzlo, l_xzhi;

l_xzlo = xlo(i,j,k,n);
Expand Down Expand Up @@ -334,7 +335,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
sth += (fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.;
}

auto bc = pbc[n];
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) )
Expand Down
27 changes: 14 additions & 13 deletions EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Array4<Real const> const& fcz,
Array4<Real const> const& ccent_arr,
bool is_velocity,
Array4<Real const> const& values_on_eb_inflow)
Array4<Real const> const& values_on_eb_inflow,
Array4<int const> const& bc_arr)
{

// bx is the cell-centered box on which we want to compute the advective update
Expand Down Expand Up @@ -135,7 +136,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,

Real uad = u_mac(i,j,k);

auto bc = pbc[n];
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);

Expand All @@ -153,7 +154,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,

Real vad = v_mac(i,j,k);

auto bc = pbc[n];
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);

Expand All @@ -171,7 +172,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,

Real wad = w_mac(i,j,k);

auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);

HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);

Expand Down Expand Up @@ -199,7 +200,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(zylo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_zylo, l_zyhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_zy(l_zylo, l_zyhi,
i, j, k, n, l_dt, dy, iconserv[n],
Expand All @@ -216,7 +217,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(yzlo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_yzlo, l_yzhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_yz(l_yzlo, l_yzhi,
i, j, k, n, l_dt, dz, iconserv[n],
Expand Down Expand Up @@ -302,7 +303,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
sth += (fq) ? Real(0.5)*l_dt*fq(i ,j,k,n) : Real(0.0);
}

auto bc = pbc[n];
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) )
Expand Down Expand Up @@ -336,7 +337,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(xzlo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_xzlo, l_xzhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_xz(l_xzlo, l_xzhi,
i, j, k, n, l_dt, dz, iconserv[n],
Expand All @@ -353,7 +354,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(zxlo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_zxlo, l_zxhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_zx(l_zxlo, l_zxhi,
i, j, k, n, l_dt, dx, iconserv[n],
Expand Down Expand Up @@ -439,7 +440,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
sth += (fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0.0);
}

auto bc = pbc[n];
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) )
Expand Down Expand Up @@ -471,7 +472,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(xylo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_xylo, l_xyhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_xy(l_xylo, l_xyhi,
i, j, k, n, l_dt, dy, iconserv[n],
Expand All @@ -488,7 +489,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
Box(yxlo), ncomp,
[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const auto bc = pbc[n];
const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
Real l_yxlo, l_yxhi;
EBGodunovCornerCouple::EBGodunov_corner_couple_yx(l_yxlo, l_yxhi,
i, j, k, n, l_dt, dx, iconserv[n],
Expand Down Expand Up @@ -573,7 +574,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
sth += (fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0.0);
}

auto bc = pbc[n];
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) )
Expand Down
38 changes: 23 additions & 15 deletions EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
Vector<BCRec> const& h_bcrec,
BCRec const* d_bcrec,
const Geometry& geom,
Real l_dt, MultiFab const* velocity_on_eb_inflow)
Real l_dt, MultiFab const* velocity_on_eb_inflow,
iMultiFab const* BC_MF)
{
BL_PROFILE("EBGodunov::ExtrapVelToFaces()");
AMREX_ALWAYS_ASSERT(vel.hasEBFabFactory());
Expand Down Expand Up @@ -64,6 +65,9 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
Array4<Real const> const& a_vel = vel.const_array(mfi);
Array4<Real const> const& a_f = vel_forces.const_array(mfi);

Array4<int const> const& bc_arr = BC_MF ? BC_MF->const_array(mfi)
: Array4<int const> {};

// In 2-d:
// 8*ncomp are: Imx, Ipx, Imy, Ipy, xlo/xhi, ylo/yhi
// 2 are: u_ad, v_ad
Expand Down Expand Up @@ -142,14 +146,14 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
#endif

PLM::PredictVelOnXFace( xebx_g1, AMREX_SPACEDIM, Imx, Ipx, a_vel, a_vel,
geom, l_dt, h_bcrec, d_bcrec);
geom, l_dt, h_bcrec, d_bcrec, bc_arr);

PLM::PredictVelOnYFace( yebx_g1, AMREX_SPACEDIM, Imy, Ipy, a_vel, a_vel,
geom, l_dt, h_bcrec, d_bcrec);
geom, l_dt, h_bcrec, d_bcrec, bc_arr);

#if ( AMREX_SPACEDIM == 3 )
PLM::PredictVelOnZFace( zebx_g1, AMREX_SPACEDIM, Imz, Ipz, a_vel, a_vel,
geom, l_dt, h_bcrec, d_bcrec);
geom, l_dt, h_bcrec, d_bcrec, bc_arr);
#endif


Expand All @@ -159,7 +163,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
AMREX_D_DECL(Imx, Imy, Imz),
AMREX_D_DECL(Ipx, Ipy, Ipz),
a_vel, a_f, domain, l_dt, d_bcrec,
local_use_forces_in_trans);
local_use_forces_in_trans,
bc_arr);

Godunov::ExtrapVelToFacesOnBox( bx, ncomp,
AMREX_D_DECL(xbx, ybx, zbx),
Expand All @@ -169,7 +174,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
AMREX_D_DECL(Imx, Imy, Imz),
AMREX_D_DECL(Ipx, Ipy, Ipz),
a_f, domain, dx, l_dt, d_bcrec,
local_use_forces_in_trans, p);
local_use_forces_in_trans, p,
bc_arr);
}
else
{
Expand All @@ -184,25 +190,25 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
EBPLM::PredictVelOnXFace( xebx_g2, Imx, Ipx, a_vel, a_vel,
flagarr, vfrac_arr,
AMREX_D_DECL(fcx,fcy,fcz),ccent_arr,
geom, l_dt, h_bcrec, d_bcrec );
geom, l_dt, h_bcrec, d_bcrec, bc_arr );

EBPLM::PredictVelOnYFace( yebx_g2, Imy, Ipy, a_vel, a_vel,
flagarr, vfrac_arr,
AMREX_D_DECL(fcx,fcy,fcz),ccent_arr,
geom, l_dt, h_bcrec, d_bcrec );
geom, l_dt, h_bcrec, d_bcrec, bc_arr );

#if (AMREX_SPACEDIM == 3)
EBPLM::PredictVelOnZFace( zebx_g2, Imz, Ipz, a_vel, a_vel,
flagarr, vfrac_arr,
AMREX_D_DECL(fcx,fcy,fcz),ccent_arr,
geom, l_dt, h_bcrec, d_bcrec );
geom, l_dt, h_bcrec, d_bcrec, bc_arr );
#endif

EBGodunov::ComputeAdvectiveVel( AMREX_D_DECL(xebx_g2, yebx_g2, zebx_g2),
AMREX_D_DECL(u_ad, v_ad, w_ad),
AMREX_D_DECL(Imx, Imy, Imz),
AMREX_D_DECL(Ipx, Ipy, Ipz),
a_vel, flagarr, domain, d_bcrec);
a_vel, flagarr, domain, d_bcrec, bc_arr);

AMREX_D_TERM(Array4<Real const> const& apx = areafrac[0]->const_array(mfi);,
Array4<Real const> const& apy = areafrac[1]->const_array(mfi);,
Expand All @@ -226,7 +232,8 @@ EBGodunov::ExtrapVelToFaces ( MultiFab const& vel,
AMREX_D_DECL(fcx, fcy, fcz),
p,
velocity_on_eb_inflow ?
velocity_on_eb_inflow->const_array(mfi) : Array4<Real const>{});
velocity_on_eb_inflow->const_array(mfi) : Array4<Real const>{},
bc_arr);
}

Gpu::streamSynchronize(); // otherwise we might be using too much memory
Expand All @@ -250,7 +257,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
Array4<Real const> const& vel,
Array4<EBCellFlag const> const& flag,
const Box& domain,
BCRec const* pbc)
BCRec const* pbc,
Array4<int const> const& bc_arr)
{
const Dim3 dlo = amrex::lbound(domain);
const Dim3 dhi = amrex::ubound(domain);
Expand All @@ -266,7 +274,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
Real lo = Ipx(i-1,j,k,n);
Real hi = Imx(i ,j,k,n);

auto bc = pbc[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);

Real st = ( (lo+hi) >= 0.) ? lo : hi;
Expand All @@ -286,7 +294,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
Real lo = Ipy(i,j-1,k,n);
Real hi = Imy(i,j ,k,n);

auto bc = pbc[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);

Real st = ( (lo+hi) >= 0.) ? lo : hi;
Expand All @@ -307,7 +315,7 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
Real lo = Ipz(i,j,k-1,n);
Real hi = Imz(i,j,k ,n);

auto bc = pbc[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);

Real st = ( (lo+hi) >= 0.) ? lo : hi;
Expand Down
Loading

0 comments on commit 1ec1247

Please sign in to comment.