Skip to content

Commit

Permalink
remove manifold-specfic changes, will go in later PR
Browse files Browse the repository at this point in the history
  • Loading branch information
baperry2 committed Sep 17, 2024
1 parent af7813d commit 472c78d
Show file tree
Hide file tree
Showing 2 changed files with 110 additions and 124 deletions.
1 change: 0 additions & 1 deletion Source/PeleLMeX.H
Original file line number Diff line number Diff line change
Expand Up @@ -606,7 +606,6 @@ public:
* velocity, ensuring they sum up to zero. \param a_fluxes diffusion fluxes to
* be updated \param a_spec species rhoYs state data on all levels
*/
template <typename EOSType>
void adjustSpeciesFluxes(
const amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>>&
a_spfluxes,
Expand Down
233 changes: 110 additions & 123 deletions Source/PeleLMeX_Diffusion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -204,126 +204,6 @@ PeleLM::computeDifferentialDiffusionTerms(
#endif
}

template <typename EOSType>
void
PeleLM::adjustSpeciesFluxes(
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& a_spfluxes,
Vector<MultiFab const*> const& a_spec)
{

BL_PROFILE("PeleLMeX::adjustSpeciesFluxes()");

// Get the species BCRec
auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES);

for (int lev = 0; lev <= finest_level; ++lev) {

const Box& domain = geom[lev].Domain();

#ifdef AMREX_USE_EB
//------------------------------------------------------------------------
// Get the edge species state needed for EB
int nGrow = 1;
const auto& ba = a_spec[lev]->boxArray();
const auto& dm = a_spec[lev]->DistributionMap();
const auto& ebfact = EBFactory(lev);
Array<MultiFab, AMREX_SPACEDIM> edgstate{AMREX_D_DECL(
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(0)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact),
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(1)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact),
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(2)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact))};
EB_interp_CellCentroid_to_FaceCentroid(
*a_spec[lev], GetArrOfPtrs(edgstate), 0, 0, NUM_SPECIES, geom[lev],
bcRecSpec);
auto const& areafrac = ebfact.getAreaFrac();
#endif

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*a_spec[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
const Box& ebx = mfi.nodaltilebox(idim);
const Box& edomain = amrex::surroundingNodes(domain, idim);
auto const& rhoY = a_spec[lev]->const_array(mfi);
auto const& flux_dir = a_spfluxes[lev][idim]->array(mfi);

const auto bc_lo = bcRecSpec[0].lo(idim);
const auto bc_hi = bcRecSpec[0].hi(idim);

#ifdef AMREX_USE_EB
auto const& flagfab = ebfact.getMultiEBCellFlagFab()[mfi];

if (flagfab.getType(amrex::grow(ebx, 0)) != FabType::covered) {
// No cut cells in tile + nghost-cell width halo -> use non-eb routine
if (flagfab.getType(amrex::grow(ebx, nGrow)) == FabType::regular) {
amrex::ParallelFor(
ebx, [idim, rhoY, flux_dir, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux(i, j, k, idim, on_lo, on_hi, rhoY, flux_dir);
});
} else {
auto const& rhoYed_ar = edgstate[idim].const_array(mfi);
auto const& areafrac_ar = areafrac[idim]->const_array(mfi);
amrex::ParallelFor(
ebx,
[idim, rhoY, flux_dir, rhoYed_ar, areafrac_ar, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux_eb(
i, j, k, idim, on_lo, on_hi, rhoY, rhoYed_ar, areafrac_ar,
flux_dir);
});
}
}
#else
amrex::ParallelFor(
ebx, [idim, rhoY, flux_dir, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux(i, j, k, idim, on_lo, on_hi, rhoY, flux_dir);
});
#endif
}
}
}
}

template <>
void
PeleLM::adjustSpeciesFluxes<pele::physics::eos::Manifold>(
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& /*a_spfluxes*/,
Vector<MultiFab const*> const& /*a_spec*/)
{
// Manifold Model: "Species" don't sum to unity, so no need to adjust the
// fluxes
BL_PROFILE("PeleLM::adjustSpeciesFluxes()");
}

void
PeleLM::computeDifferentialDiffusionFluxes(
const TimeStamp& a_time,
Expand Down Expand Up @@ -410,8 +290,7 @@ PeleLM::computeDifferentialDiffusionFluxes(
}

// Adjust species diffusion fluxes to ensure their sum is zero
adjustSpeciesFluxes<pele::physics::PhysicsType::eos_type>(
a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time)));
adjustSpeciesFluxes(a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time)));
//----------------------------------------------------------------

//----------------------------------------------------------------
Expand Down Expand Up @@ -760,6 +639,114 @@ PeleLM::addSoretTerm(
}
}

void
PeleLM::adjustSpeciesFluxes(
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& a_spfluxes,
Vector<MultiFab const*> const& a_spec)
{

BL_PROFILE("PeleLMeX::adjustSpeciesFluxes()");

// Get the species BCRec
auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES);

for (int lev = 0; lev <= finest_level; ++lev) {

const Box& domain = geom[lev].Domain();

#ifdef AMREX_USE_EB
//------------------------------------------------------------------------
// Get the edge species state needed for EB
int nGrow = 1;
const auto& ba = a_spec[lev]->boxArray();
const auto& dm = a_spec[lev]->DistributionMap();
const auto& ebfact = EBFactory(lev);
Array<MultiFab, AMREX_SPACEDIM> edgstate{AMREX_D_DECL(
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(0)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact),
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(1)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact),
MultiFab(
amrex::convert(ba, IntVect::TheDimensionVector(2)), dm, NUM_SPECIES,
nGrow, MFInfo(), ebfact))};
EB_interp_CellCentroid_to_FaceCentroid(
*a_spec[lev], GetArrOfPtrs(edgstate), 0, 0, NUM_SPECIES, geom[lev],
bcRecSpec);
auto const& areafrac = ebfact.getAreaFrac();
#endif

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*a_spec[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
const Box& ebx = mfi.nodaltilebox(idim);
const Box& edomain = amrex::surroundingNodes(domain, idim);
auto const& rhoY = a_spec[lev]->const_array(mfi);
auto const& flux_dir = a_spfluxes[lev][idim]->array(mfi);

const auto bc_lo = bcRecSpec[0].lo(idim);
const auto bc_hi = bcRecSpec[0].hi(idim);

#ifdef AMREX_USE_EB
auto const& flagfab = ebfact.getMultiEBCellFlagFab()[mfi];

if (flagfab.getType(amrex::grow(ebx, 0)) != FabType::covered) {
// No cut cells in tile + nghost-cell width halo -> use non-eb routine
if (flagfab.getType(amrex::grow(ebx, nGrow)) == FabType::regular) {
amrex::ParallelFor(
ebx, [idim, rhoY, flux_dir, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux(i, j, k, idim, on_lo, on_hi, rhoY, flux_dir);
});
} else {
auto const& rhoYed_ar = edgstate[idim].const_array(mfi);
auto const& areafrac_ar = areafrac[idim]->const_array(mfi);
amrex::ParallelFor(
ebx,
[idim, rhoY, flux_dir, rhoYed_ar, areafrac_ar, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux_eb(
i, j, k, idim, on_lo, on_hi, rhoY, rhoYed_ar, areafrac_ar,
flux_dir);
});
}
}
#else
amrex::ParallelFor(
ebx, [idim, rhoY, flux_dir, edomain, bc_lo,
bc_hi] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
int idx[3] = {i, j, k};
bool on_lo =
((bc_lo == amrex::BCType::ext_dir) &&
(idx[idim] <= edomain.smallEnd(idim)));
bool on_hi =
((bc_hi == amrex::BCType::ext_dir) &&
(idx[idim] >= edomain.bigEnd(idim)));
repair_flux(i, j, k, idim, on_lo, on_hi, rhoY, flux_dir);
});
#endif
}
}
}
}

void
PeleLM::computeSpeciesEnthalpyFlux(
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& a_fluxes,
Expand Down Expand Up @@ -1020,7 +1007,7 @@ PeleLM::differentialDiffusionUpdate(
fillPatchSpecies(AmrNewTime);

// Adjust species diffusion fluxes to ensure their sum is zero
adjustSpeciesFluxes<pele::physics::PhysicsType::eos_type>(
adjustSpeciesFluxes(
GetVecOfArrOfPtrs(fluxes), GetVecOfConstPtrs(getSpeciesVect(AmrNewTime)));

// Average down fluxes^{np1,kp1}
Expand Down

0 comments on commit 472c78d

Please sign in to comment.