Skip to content

Commit

Permalink
Eosparm support in PeleLMeX (AMReX-Combustion#414)
Browse files Browse the repository at this point in the history 
* basics of having eos_parm in PeleLMeX

* add eosparm to PeleLMeX_K functions

* more eosparm: now everywhere except Efield & BCs

* eosparm in BC stuff

* re-order functions for template

* remove manifold-specfic changes, will go in later PR
  • Loading branch information
@baperry2
baperry2 authored Sep 23, 2024
1 parent 93fecb4 commit af2853a
Show file tree
Hide file tree
Showing 17 changed files with 188 additions and 112 deletions.
4 changes: 2 additions & 2 deletions Exec/RegTests/FlameSheet/pelelmex_prob.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ pelelmex_initdata(

constexpr amrex::Real Pi = 3.14159265358979323846264338327950288;

auto eos = pele::physics::PhysicsType::eos();
auto eos = pele::physics::PhysicsType::eos(prob_parm.eosparm);
amrex::GpuArray<amrex::Real, NUM_SPECIES + 4> pmf_vals = {0.0};
amrex::Real molefrac[NUM_SPECIES] = {0.0};
amrex::Real massfrac[NUM_SPECIES] = {0.0};
Expand Down Expand Up @@ -140,7 +140,7 @@ bcnormal(
amrex::Real molefrac[NUM_SPECIES] = {0.0};
amrex::Real massfrac[NUM_SPECIES] = {0.0};

auto eos = pele::physics::PhysicsType::eos();
auto eos = pele::physics::PhysicsType::eos(prob_parm.eosparm);
if (sgn == 1) {
pele::physics::PMF::pmf(pmf_data, prob_lo[idir], prob_lo[idir], pmf_vals);

Expand Down
1 change: 1 addition & 0 deletions Exec/RegTests/FlameSheet/pelelmex_prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,6 @@ PeleLM::readProbParm()
pp.query("pertmag", PeleLM::prob_parm->pertmag);
pp.query("pertlength", PeleLM::prob_parm->pertlength);

PeleLM::prob_parm->eosparm = PeleLM::eos_parms.device_parm();
PeleLM::pmf_data.initialize();
}
3 changes: 3 additions & 0 deletions Exec/RegTests/FlameSheet/pelelmex_prob_parm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include <AMReX_REAL.H>
#include <AMReX_GpuMemory.H>
#include <PelePhysics.H>

using namespace amrex::literals;

Expand All @@ -13,5 +14,7 @@ struct ProbParm
amrex::Real pertmag = 0.0004_rt;
amrex::Real pertlength = 0.008_rt;
int meanFlowDir = 1;
pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type> const*
eosparm;
};
#endif
5 changes: 5 additions & 0 deletions Source/PeleLMeX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1753,6 +1753,11 @@ public:
pele::physics::PhysicsType::transport_type>>
trans_parms;

// EOS pointer
static pele::physics::PeleParams<
pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>>
eos_parms;

// Reactor pointer
std::string m_chem_integrator;
std::unique_ptr<pele::physics::reactions::ReactorBase> m_reactor;
Expand Down
5 changes: 5 additions & 0 deletions Source/PeleLMeX.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,17 @@ pele::physics::PeleParams<pele::physics::transport::TransParm<
pele::physics::PhysicsType::transport_type>>
PeleLM::trans_parms;

pele::physics::PeleParams<
pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>>
PeleLM::eos_parms;

PeleLM::PeleLM() = default;

PeleLM::~PeleLM()
{
if (m_incompressible == 0) {
trans_parms.deallocate();
eos_parms.deallocate();
m_reactor->close();
}

Expand Down
20 changes: 11 additions & 9 deletions Source/PeleLMeX_Advection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,7 @@ PeleLM::getScalarAdvForce(
// Get t^{n} data pointer
auto* ldata_p = getLevelDataPtr(lev, AmrOldTime);
auto* ldataR_p = getLevelDataReactPtr(lev);
auto const* leosparm = eos_parms.device_parm();

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
Expand All @@ -276,11 +277,11 @@ PeleLM::getScalarAdvForce(
amrex::ParallelFor(
bx,
[rho, rhoY, T, dn, ddn, r, fY, fT, extRhoY, extRhoH, dp0dt = m_dp0dt,
is_closed_ch = m_closed_chamber,
do_react = m_do_react] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
is_closed_ch = m_closed_chamber, do_react = m_do_react,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
buildAdvectionForcing(
i, j, k, rho, rhoY, T, dn, ddn, r, extRhoY, extRhoH, dp0dt,
is_closed_ch, do_react, fY, fT);
is_closed_ch, do_react, fY, fT, leosparm);
});
}
}
Expand Down Expand Up @@ -552,6 +553,7 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
for (MFIter mfi(ldata_p->state, TilingIfNotGPU()); mfi.isValid(); ++mfi) {

Box const& bx = mfi.tilebox();
auto const* leosparm = eos_parms.device_parm();

#ifdef AMREX_USE_EB
auto const& flagfab = ebfact.getMultiEBCellFlagFab()[mfi];
Expand All @@ -573,21 +575,21 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr<AdvanceAdvData>& advData)
// boxes
const auto& afrac = areafrac[idim]->array(mfi);
amrex::ParallelFor(
ebx, [rho, rhoY, T, rhoHm,
afrac] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
ebx, [rho, rhoY, T, rhoHm, afrac,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
if (afrac(i, j, k) <= 0.0) { // Covered faces
rhoHm(i, j, k) = 0.0;
} else {
getRHmixGivenTY(i, j, k, rho, rhoY, T, rhoHm);
getRHmixGivenTY(i, j, k, rho, rhoY, T, rhoHm, leosparm);
}
});
} else // Regular boxes
#endif
{
amrex::ParallelFor(
ebx, [rho, rhoY, T,
rhoHm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getRHmixGivenTY(i, j, k, rho, rhoY, T, rhoHm);
ebx, [rho, rhoY, T, rhoHm,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getRHmixGivenTY(i, j, k, rho, rhoY, T, rhoHm, leosparm);
});
}
}
Expand Down
24 changes: 14 additions & 10 deletions Source/PeleLMeX_DeriveFunc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,9 @@ pelelmex_derheatrelease(
auto const react = reactfab.const_array(0);
auto const& Hi = EnthFab.array();
auto HRR = derfab.array(dcomp);
auto const* leosparm = a_pelelm->eos_parms.device_parm();
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getHGivenT(i, j, k, temp, Hi);
getHGivenT(i, j, k, temp, Hi, leosparm);
HRR(i, j, k) = 0.0;
for (int n = 0; n < NUM_SPECIES; n++) {
HRR(i, j, k) -= Hi(i, j, k, n) * react(i, j, k, n);
Expand Down Expand Up @@ -146,14 +147,15 @@ pelelmex_dermolefrac(
AMREX_ASSERT(!a_pelelm->m_incompressible);
auto const in_dat = statefab.array();
auto der = derfab.array(dcomp);
auto const* leosparm = a_pelelm->eos_parms.device_parm();
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
amrex::Real Yt[NUM_SPECIES] = {0.0};
amrex::Real Xt[NUM_SPECIES] = {0.0};
amrex::Real rhoinv = 1.0 / in_dat(i, j, k, DENSITY);
for (int n = 0; n < NUM_SPECIES; n++) {
Yt[n] = in_dat(i, j, k, FIRSTSPEC + n) * rhoinv;
}
auto eos = pele::physics::PhysicsType::eos();
auto eos = pele::physics::PhysicsType::eos(leosparm);
eos.Y2X(Yt, Xt);
for (int n = 0; n < NUM_SPECIES; n++) {
der(i, j, k, n) = Xt[n];
Expand Down Expand Up @@ -1371,17 +1373,18 @@ pelelmex_derdiffc(
auto lambda = dummies.array(0);
auto mu = dummies.array(1);
auto const* ltransparm = a_pelelm->trans_parms.device_parm();
auto const* leosparm = a_pelelm->eos_parms.device_parm();
auto rhotheta = do_soret ? derfab.array(dcomp + NUM_SPECIES)
: dummies.array(2); // dummy for no soret
amrex::Real LeInv = a_pelelm->m_Lewis_inv;
amrex::Real PrInv = a_pelelm->m_Prandtl_inv;
amrex::ParallelFor(
bx,
[do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD, rhotheta,
lambda, mu, ltransparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
bx, [do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD,
rhotheta, lambda, mu, ltransparm,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getTransportCoeff(
i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T,
rhoD, rhotheta, lambda, mu, ltransparm);
rhoD, rhotheta, lambda, mu, ltransparm, leosparm);
});
}

Expand Down Expand Up @@ -1418,15 +1421,16 @@ pelelmex_derlambda(
auto mu = dummies.array(0);
auto rhotheta = dummies.array(NUM_SPECIES + 1);
auto const* ltransparm = a_pelelm->trans_parms.device_parm();
auto const* leosparm = a_pelelm->eos_parms.device_parm();
amrex::Real LeInv = a_pelelm->m_Lewis_inv;
amrex::Real PrInv = a_pelelm->m_Prandtl_inv;
amrex::ParallelFor(
bx,
[do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD, rhotheta,
lambda, mu, ltransparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
bx, [do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD,
rhotheta, lambda, mu, ltransparm,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getTransportCoeff(
i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T,
rhoD, rhotheta, lambda, mu, ltransparm);
rhoD, rhotheta, lambda, mu, ltransparm, leosparm);
});
}

Expand Down
35 changes: 20 additions & 15 deletions Source/PeleLMeX_Diffusion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,6 +366,7 @@ PeleLM::addWbarTerm(
// Compute Wbar on all the levels
int nGrow = 1; // Need one ghost cell to compute gradWbar
Vector<MultiFab> Wbar(finest_level + 1);
auto const* leosparm = eos_parms.device_parm();
for (int lev = 0; lev <= finest_level; ++lev) {

Wbar[lev].define(grids[lev], dmap[lev], 1, nGrow, MFInfo(), Factory(lev));
Expand All @@ -379,9 +380,9 @@ PeleLM::addWbarTerm(
auto const& rhoY_arr = a_spec[lev]->const_array(mfi);
auto const& Wbar_arr = Wbar[lev].array(mfi);
amrex::ParallelFor(
gbx, [rho_arr, rhoY_arr,
Wbar_arr] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr);
gbx, [rho_arr, rhoY_arr, Wbar_arr,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getMwmixGivenRY(i, j, k, rho_arr, rhoY_arr, Wbar_arr, leosparm);
});
}
}
Expand Down Expand Up @@ -464,10 +465,11 @@ PeleLM::addWbarTerm(
// Wbar flux is : - \rho Y_m / \overline{W} * D_m * \nabla
// \overline{W} with beta_m = \rho * D_m below
amrex::ParallelFor(
ebx,
[need_wbar_fluxes, gradWbar_ar, beta_ar, rhoY, spFlux_ar,
spwbarFlux_ar] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
auto eos = pele::physics::PhysicsType::eos();
ebx, [need_wbar_fluxes, gradWbar_ar, beta_ar, rhoY, spFlux_ar,
spwbarFlux_ar,
eosparm =
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
auto eos = pele::physics::PhysicsType::eos(eosparm);
// Get Wbar from rhoYs
amrex::Real rho = 0.0;
for (int n = 0; n < NUM_SPECIES; n++) {
Expand Down Expand Up @@ -755,6 +757,7 @@ PeleLM::computeSpeciesEnthalpyFlux(

// Get the species BCRec
auto bcRecSpec = fetchBCRecArray(FIRSTSPEC, NUM_SPECIES);
auto const* leosparm = eos_parms.device_parm();

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

Expand Down Expand Up @@ -786,21 +789,21 @@ PeleLM::computeSpeciesEnthalpyFlux(
} else if (flagfab.getType(gbx) != FabType::regular) { // EB containing
// boxes
amrex::ParallelFor(
gbx, [Temp_arr, Hi_arr,
flag] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
gbx, [Temp_arr, Hi_arr, flag,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
if (flag(i, j, k).isCovered()) {
Hi_arr(i, j, k) = 0.0;
} else {
getHGivenT(i, j, k, Temp_arr, Hi_arr);
getHGivenT(i, j, k, Temp_arr, Hi_arr, leosparm);
}
});
} else
#endif
{
amrex::ParallelFor(
gbx,
[Temp_arr, Hi_arr] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getHGivenT(i, j, k, Temp_arr, Hi_arr);
gbx, [Temp_arr, Hi_arr,
leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
getHGivenT(i, j, k, Temp_arr, Hi_arr, leosparm);
});
}
}
Expand Down Expand Up @@ -1226,6 +1229,7 @@ PeleLM::deltaTIter_prepare(
std::unique_ptr<AdvanceAdvData>& advData,
std::unique_ptr<AdvanceDiffData>& diffData)
{
auto const* leosparm = eos_parms.device_parm();
for (int lev = 0; lev <= finest_level; ++lev) {

auto* ldataOld_p = getLevelDataPtr(lev, AmrOldTime);
Expand Down Expand Up @@ -1263,7 +1267,7 @@ PeleLM::deltaTIter_prepare(
fourier(i, j, k) + diffDiff(i, j, k));

// Get \rho * Cp_{mix}
getCpmixGivenRYT(i, j, k, rho, rhoY, T, rhocp);
getCpmixGivenRYT(i, j, k, rho, rhoY, T, rhocp, leosparm);
rhocp(i, j, k) *= rho(i, j, k);

// Save T
Expand Down Expand Up @@ -1370,6 +1374,7 @@ PeleLM::deltaTIter_update(
//------------------------------------------------------------------------
// Recompute RhoH
for (int lev = 0; lev <= finest_level; ++lev) {
auto const* leosparm = eos_parms.device_parm();
auto* ldata_p = getLevelDataPtr(lev, AmrNewTime);
auto const& sma = ldata_p->state.arrays();
amrex::ParallelFor(
Expand All @@ -1379,7 +1384,7 @@ PeleLM::deltaTIter_update(
i, j, k, Array4<Real const>(sma[box_no], DENSITY),
Array4<Real const>(sma[box_no], FIRSTSPEC),
Array4<Real const>(sma[box_no], TEMP),
Array4<Real>(sma[box_no], RHOH));
Array4<Real>(sma[box_no], RHOH), leosparm);
});
}
Gpu::streamSynchronize();
Expand Down
Loading

0 comments on commit af2853a

Please sign in to comment.