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

Rework phys_bc to use BCRec #469

Merged
merged 14 commits into from
Jul 18, 2024
15 changes: 8 additions & 7 deletions Source/Maestro.H
Original file line number Diff line number Diff line change
Expand Up @@ -1734,7 +1734,7 @@ class Maestro : public amrex::AmrCore {
amrex::Array4<amrex::Real> const ury,
amrex::Array4<amrex::Real> const uimhy,
const amrex::Box& domainBox,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;

zingale marked this conversation as resolved.
Show resolved Hide resolved
void VelPredVelocities(const amrex::MFIter& mfi,
amrex::Array4<const amrex::Real> const utilde,
Expand All @@ -1755,7 +1755,8 @@ class Maestro : public amrex::AmrCore {
amrex::Array4<const amrex::Real> const force,
amrex::Array4<const amrex::Real> const w0_cart_in,
const amrex::Box& domainBox,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;

#else
void VelPredInterface(const amrex::MFIter& mfi,
amrex::Array4<const amrex::Real> const utilde,
Expand All @@ -1779,7 +1780,7 @@ class Maestro : public amrex::AmrCore {
amrex::Array4<amrex::Real> const urz,
amrex::Array4<amrex::Real> const uimhz,
const amrex::Box& domainBox,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;

void VelPredTransverse(const amrex::MFIter& mfi,
amrex::Array4<const amrex::Real> const utilde,
Expand All @@ -1802,7 +1803,7 @@ class Maestro : public amrex::AmrCore {
amrex::Array4<amrex::Real> const wimhxy,
amrex::Array4<amrex::Real> const wimhyx,
const amrex::Box& domainBox,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;

void VelPredVelocities(const amrex::MFIter& mfi,
amrex::Array4<const amrex::Real> const utilde,
Expand Down Expand Up @@ -1836,7 +1837,8 @@ class Maestro : public amrex::AmrCore {
amrex::Array4<const amrex::Real> const force,
amrex::Array4<const amrex::Real> const w0_cart_in,
const amrex::Box& domainBox,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;

#endif
////////////////////////

Expand Down Expand Up @@ -1912,8 +1914,7 @@ class Maestro : public amrex::AmrCore {
amrex::Vector<amrex::iMultiFab> cell_cc_to_r;

/// stores domain boundary conditions.
/// These muse be vectors (rather than arrays) so we can ParmParse them
IntVector phys_bc;
amrex::BCRec phys_bc;

/// Boundary condition objects needed for FillPatch routines.
/// This is essentially an array (over components)
Expand Down
22 changes: 14 additions & 8 deletions Source/MaestroFillData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -194,8 +194,6 @@ void Maestro::FillUmacGhost(
const auto domlo = domainBox.loVect3d();
const auto domhi = domainBox.hiVect3d();

const int* AMREX_RESTRICT physbc_p = phys_bc.dataPtr();

// get references to the MultiFabs at level lev
MultiFab& sold_mf =
sold[lev]; // need a cell-centered MF for the MFIter
Expand All @@ -216,11 +214,13 @@ void Maestro::FillUmacGhost(
#if (AMREX_SPACEDIM == 3)
const Array4<Real> wmac = umac_in[lev][2].array(mfi);
#endif
int bclo = phys_bc.lo(0);
int bchi = phys_bc.hi(0);

ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
// lo x-faces
if (i == domlo[0] - 1) {
switch (physbc_p[0]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bclo) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j, k) = umac(i + 1, j, k);
vmac(i, j, k) = 0.0;
Expand Down Expand Up @@ -258,7 +258,7 @@ void Maestro::FillUmacGhost(

// hi x-faces
if (i == domhi[0] + 2) {
switch (physbc_p[AMREX_SPACEDIM]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bchi) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j, k) = umac(i - 1, j, k);
vmac(i - 1, j, k) = 0.0;
Expand Down Expand Up @@ -297,10 +297,13 @@ void Maestro::FillUmacGhost(

Gpu::synchronize();

bclo = phys_bc.lo(1);
bchi = phys_bc.hi(1);

ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
// lo y-faces
if (j == domlo[1] - 1) {
switch (physbc_p[1]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bclo) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j, k) = 0.0;
vmac(i, j, k) = vmac(i, j + 1, k);
Expand Down Expand Up @@ -338,7 +341,7 @@ void Maestro::FillUmacGhost(

// hi y-faces
if (j == domhi[1] + 2) {
switch (physbc_p[AMREX_SPACEDIM + 1]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bchi) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j - 1, k) = 0.0;
vmac(i, j, k) = vmac(i, j - 1, k);
Expand Down Expand Up @@ -379,10 +382,13 @@ void Maestro::FillUmacGhost(

Gpu::synchronize();

bclo = phys_bc.lo(2);
bchi = phys_bc.hi(2);

ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
// lo z-faces
if (k == domlo[2] - 1) {
switch (physbc_p[2]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bclo) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j, k) = 0.0;
vmac(i, j, k) = 0.0;
Expand Down Expand Up @@ -412,7 +418,7 @@ void Maestro::FillUmacGhost(

// hi z-faces
if (k == domhi[2] + 2) {
switch (physbc_p[2 + AMREX_SPACEDIM]) { // NOLINT(bugprone-switch-missing-default-case)
switch (bchi) { // NOLINT(bugprone-switch-missing-default-case)
case amrex::PhysBCType::inflow:
umac(i, j, k - 1) = 0.0;
vmac(i, j, k - 1) = 0.0;
Expand Down
4 changes: 2 additions & 2 deletions Source/MaestroMacProj.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -434,14 +434,14 @@ void Maestro::SetMacSolverBCs(MLABecLaplacian& mlabec) {
mlmg_lobc[idim] = mlmg_hibc[idim] = LinOpBCType::Periodic;
} else {
// lo-side BCs
if (phys_bc[idim] == amrex::PhysBCType::outflow) {
if (phys_bc.lo(idim) == amrex::PhysBCType::outflow) {
mlmg_lobc[idim] = LinOpBCType::Dirichlet;
} else {
mlmg_lobc[idim] = LinOpBCType::Neumann;
}

// hi-side BCs
if (phys_bc[AMREX_SPACEDIM + idim] == amrex::PhysBCType::outflow) {
if (phys_bc.hi(idim) == amrex::PhysBCType::outflow) {
mlmg_hibc[idim] = LinOpBCType::Dirichlet;
} else {
mlmg_hibc[idim] = LinOpBCType::Neumann;
Expand Down
20 changes: 10 additions & 10 deletions Source/MaestroMakeUtrans.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,8 @@ void Maestro::MakeUtrans(
}

// create utrans
int bclo = phys_bc[0];
int bchi = phys_bc[AMREX_SPACEDIM];
int bclo = phys_bc.lo(0);
int bchi = phys_bc.hi(0);

ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real ulx = 0.0;
Expand Down Expand Up @@ -182,8 +182,8 @@ void Maestro::MakeUtrans(
}

// create vtrans
int bclo = phys_bc[1];
int bchi = phys_bc[AMREX_SPACEDIM + 1];
int bclo = phys_bc.lo(1);
int bchi = phys_bc.hi(1);

ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real vly = 0.0;
Expand Down Expand Up @@ -298,8 +298,8 @@ void Maestro::MakeUtrans(
}

// create utrans
int bclo = phys_bc[0];
int bchi = phys_bc[AMREX_SPACEDIM];
int bclo = phys_bc.lo(0);
int bchi = phys_bc.hi(0);

ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real ulx = 0.0;
Expand Down Expand Up @@ -413,8 +413,8 @@ void Maestro::MakeUtrans(
}

// create vtrans
int bclo = phys_bc[1];
int bchi = phys_bc[AMREX_SPACEDIM + 1];
int bclo = phys_bc.lo(1);
int bchi = phys_bc.hi(1);

ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real vly = 0.0;
Expand Down Expand Up @@ -529,8 +529,8 @@ void Maestro::MakeUtrans(
}

// create wtrans
int bclo = phys_bc[2];
int bchi = phys_bc[AMREX_SPACEDIM + 2];
int bclo = phys_bc.lo(2);
int bchi = phys_bc.hi(2);

ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real wlz = 0.0;
Expand Down
12 changes: 6 additions & 6 deletions Source/MaestroNodalProj.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -155,13 +155,13 @@ void Maestro::NodalProj(int proj_type, Vector<MultiFab>& rhcc,
if (Geom(0).isPeriodic(idim)) {
mlmg_lobc[idim] = mlmg_hibc[idim] = LinOpBCType::Periodic;
} else {
if (phys_bc[idim] == amrex::PhysBCType::outflow) {
if (phys_bc.lo(idim) == amrex::PhysBCType::outflow) {
mlmg_lobc[idim] = LinOpBCType::Dirichlet;
} else {
mlmg_lobc[idim] = LinOpBCType::Neumann;
}

if (phys_bc[AMREX_SPACEDIM + idim] == amrex::PhysBCType::outflow) {
if (phys_bc.hi(idim) == amrex::PhysBCType::outflow) {
mlmg_hibc[idim] = LinOpBCType::Dirichlet;
} else {
mlmg_hibc[idim] = LinOpBCType::Neumann;
Expand Down Expand Up @@ -505,8 +505,8 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {
const Box& domainBox = geom[lev].Domain();

for (int idir = 0; idir < BL_SPACEDIM; idir++) {
if (phys_bc[idir] != amrex::PhysBCType::inflow &&
phys_bc[AMREX_SPACEDIM + idir] != amrex::PhysBCType::inflow) {
if (phys_bc.lo(idir) != amrex::PhysBCType::inflow &&
phys_bc.hi(idir) != amrex::PhysBCType::inflow) {
vel[lev].setBndry(0.0, idir, 1);
} else {
#ifdef _OPENMP
Expand All @@ -521,7 +521,7 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {

BoxList bxlist(reg);

if (phys_bc[idir] == amrex::PhysBCType::inflow &&
if (phys_bc.lo(idir) == amrex::PhysBCType::inflow &&
reg.smallEnd(idir) == domainBox.smallEnd(idir)) {
Box bx; // bx is the region we *protect* from zero'ing

Expand All @@ -545,7 +545,7 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {
bxlist.push_back(bx);
}

if (phys_bc[AMREX_SPACEDIM + idir] == amrex::PhysBCType::inflow &&
if (phys_bc.hi(idir) == amrex::PhysBCType::inflow &&
reg.bigEnd(idir) == domainBox.bigEnd(idir)) {
Box bx; // bx is the region we *protect* from zero'ing

Expand Down
61 changes: 29 additions & 32 deletions Source/MaestroPlot.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1681,45 +1681,42 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,

Array4<const Real> const u = vel[lev].array(mfi);
Array4<Real> const vort = vorticity[lev].array(mfi);
GpuArray<int, AMREX_SPACEDIM * 2> physbc;
for (int n = 0; n < AMREX_SPACEDIM * 2; ++n) {
physbc[n] = phys_bc[n];
}


#if (AMREX_SPACEDIM == 2)
bool fix_lo_x = (phys_bc.lo(0) == amrex::PhysBCType::inflow || phys_bc.lo(0) == amrex::PhysBCType::slipwall ||
phys_bc.lo(0) == amrex::PhysBCType::noslipwall);
bool fix_hi_x = (phys_bc.hi(0) == amrex::PhysBCType::inflow || phys_bc.hi(0) == amrex::PhysBCType::slipwall ||
phys_bc.hi(0) == amrex::PhysBCType::noslipwall);
bool fix_lo_y = (phys_bc.lo(1) == amrex::PhysBCType::inflow || phys_bc.lo(1) == amrex::PhysBCType::slipwall ||
phys_bc.lo(1) == amrex::PhysBCType::noslipwall);
bool fix_hi_y = (phys_bc.hi(1) == amrex::PhysBCType::inflow || phys_bc.hi(1) == amrex::PhysBCType::slipwall ||
phys_bc.hi(1) == amrex::PhysBCType::noslipwall);

ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 1)) / hx;
Real uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;

if (i == ilo && (physbc[0] == amrex::PhysBCType::inflow || physbc[0] == amrex::PhysBCType::slipwall ||
physbc[0] == amrex::PhysBCType::noslipwall)) {
if (i == ilo && fix_lo_x) {
vx = (u(i + 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
4.0 * u(i - 1, j, k, 1)) /
hx;
uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;

} else if (i == ihi + 1 &&
(physbc[AMREX_SPACEDIM] == amrex::PhysBCType::inflow ||
physbc[AMREX_SPACEDIM] == amrex::PhysBCType::slipwall ||
physbc[AMREX_SPACEDIM] == amrex::PhysBCType::noslipwall)) {
} else if (i == ihi + 1 && fix_hi_x) {
vx = -(u(i - 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
4.0 * u(i + 1, j, k, 1)) /
hx;
uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
}

if (j == jlo && (physbc[1] == amrex::PhysBCType::inflow || physbc[1] == amrex::PhysBCType::slipwall ||
physbc[1] == amrex::PhysBCType::noslipwall)) {
if (j == jlo && fix_lo_y) {
vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 0)) / hx;
uy = (u(i, j + 1, k, 0) + 3.0 * u(i, j, k, 0) -
4.0 * u(i, j - 1, k, 0)) /
hy;

} else if (j == jhi + 1 &&
(physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::inflow ||
physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::slipwall ||
physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::noslipwall)) {
} else if (j == jhi + 1 && fix_hi_y) {
vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 1)) / hx;
uy = -(u(i, j - 1, k, 0) + 3.0 * u(i, j, k, 0) -
4.0 * u(i, j + 1, k, 0)) /
Expand All @@ -1729,7 +1726,22 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,
vort(i, j, k) = vx - uy;
});

#else
# else
bool fix_lo_x =
(phys_bc.lo(0) == amrex::PhysBCType::inflow || phys_bc.lo(0) == amrex::PhysBCType::noslipwall);
bool fix_hi_x = (phys_bc.hi(0) == amrex::PhysBCType::inflow ||
phys_bc.hi(0) == amrex::PhysBCType::noslipwall);

bool fix_lo_y =
(phys_bc.lo(1) == amrex::PhysBCType::inflow || phys_bc.lo(1) == amrex::PhysBCType::noslipwall);
bool fix_hi_y = (phys_bc.hi(1) == amrex::PhysBCType::inflow ||
phys_bc.hi(1) == amrex::PhysBCType::noslipwall);

bool fix_lo_z =
(phys_bc.lo(2) == amrex::PhysBCType::inflow || phys_bc.lo(2) == amrex::PhysBCType::noslipwall);
bool fix_hi_z = (phys_bc.hi(2) == amrex::PhysBCType::inflow ||
phys_bc.hi(2) == amrex::PhysBCType::noslipwall);

ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
Real uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
Real uz = 0.5 * (u(i, j, k + 1, 0) - u(i, j, k - 1, 0)) / hz;
Expand All @@ -1738,21 +1750,6 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,
Real wx = 0.5 * (u(i + 1, j, k, 2) - u(i - 1, j, k, 2)) / hx;
Real wy = 0.5 * (u(i, j + 1, k, 2) - u(i, j - 1, k, 2)) / hy;

bool fix_lo_x =
(physbc[0] == amrex::PhysBCType::inflow || physbc[0] == amrex::PhysBCType::noslipwall);
bool fix_hi_x = (physbc[AMREX_SPACEDIM] == amrex::PhysBCType::inflow ||
physbc[AMREX_SPACEDIM] == amrex::PhysBCType::noslipwall);

bool fix_lo_y =
(physbc[1] == amrex::PhysBCType::inflow || physbc[1] == amrex::PhysBCType::noslipwall);
bool fix_hi_y = (physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::inflow ||
physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::noslipwall);

bool fix_lo_z =
(physbc[2] == amrex::PhysBCType::inflow || physbc[2] == amrex::PhysBCType::noslipwall);
bool fix_hi_z = (physbc[AMREX_SPACEDIM + 2] == amrex::PhysBCType::inflow ||
physbc[AMREX_SPACEDIM + 2] == amrex::PhysBCType::noslipwall);

// First do all the faces
if (fix_lo_x && i == ilo) {
vx = (u(i + 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
Expand Down
Loading