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.

Sign up for GitHub

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

Add support for mixed bcs in (EB)godunov PLM and MacProjector #123

Merged
merged 12 commits into from
Feb 28, 2024
Merged

Add support for mixed bcs in (EB)godunov PLM and MacProjector #123

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
a731ab3
Add option to pass Robin BCs
cgilet Jan 24, 2024
1612b96
Fix typo
cgilet Jan 24, 2024
9e60faa
WIP - passing through position dep BCs, use exclusively if exist
cgilet Feb 8, 2024
f1d5904
Define and use a PPM::default_limiter, was defacto VanLeer
cgilet Feb 8, 2024
fb5331e
WIP - compiles now
cgilet Feb 8, 2024
cb1542e
WIP - extrapolated velocity handles mixed BCs. Edgestates broken
cgilet Feb 12, 2024
cbf5380
WIP - EBGod appears to work with passed position dependent bc array
cgilet Feb 12, 2024
e39e4f5
Add poisition dependent BC option for 3D (EB)Godunov. Needs cleanup
cgilet Feb 12, 2024
fdbb31e
Make regular 2D Godunov accept position de[endent BC array. Regular
cgilet Feb 14, 2024
04f6a5f
Clean up
cgilet Feb 21, 2024
18b9649
Fix spot in regular Godunov PLM that missed the bc_arr
cgilet Feb 21, 2024
1db329e
Remove extra semicolon HIP compiler complained about and consistent
cgilet Feb 21, 2024
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
13 changes: 9 additions & 4 deletions EBGodunov/hydro_ebgodunov.H
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading