From af2853a55804b9d81fa016e5781c24fe7c5f07f2 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 23 Sep 2024 15:45:52 -0600 Subject: [PATCH] Eosparm support in PeleLMeX (#414) * 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 --- Exec/RegTests/FlameSheet/pelelmex_prob.H | 4 +- Exec/RegTests/FlameSheet/pelelmex_prob.cpp | 1 + Exec/RegTests/FlameSheet/pelelmex_prob_parm.H | 3 + Source/PeleLMeX.H | 5 ++ Source/PeleLMeX.cpp | 5 ++ Source/PeleLMeX_Advection.cpp | 20 +++--- Source/PeleLMeX_DeriveFunc.cpp | 24 ++++--- Source/PeleLMeX_Diffusion.cpp | 35 +++++----- Source/PeleLMeX_Eos.cpp | 27 +++++--- Source/PeleLMeX_Forces.cpp | 4 +- Source/PeleLMeX_K.H | 66 ++++++++++++------- Source/PeleLMeX_Plot.cpp | 53 ++++++++------- Source/PeleLMeX_Reactions.cpp | 23 ++++--- Source/PeleLMeX_Setup.cpp | 8 +++ Source/PeleLMeX_TransportProp.cpp | 8 ++- Source/PeleLMeX_Utils.cpp | 12 ++-- Submodules/PelePhysics | 2 +- 17 files changed, 188 insertions(+), 112 deletions(-) diff --git a/Exec/RegTests/FlameSheet/pelelmex_prob.H b/Exec/RegTests/FlameSheet/pelelmex_prob.H index f5c48e40..c540c16d 100644 --- a/Exec/RegTests/FlameSheet/pelelmex_prob.H +++ b/Exec/RegTests/FlameSheet/pelelmex_prob.H @@ -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 pmf_vals = {0.0}; amrex::Real molefrac[NUM_SPECIES] = {0.0}; amrex::Real massfrac[NUM_SPECIES] = {0.0}; @@ -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); diff --git a/Exec/RegTests/FlameSheet/pelelmex_prob.cpp b/Exec/RegTests/FlameSheet/pelelmex_prob.cpp index d861e236..779cd18d 100644 --- a/Exec/RegTests/FlameSheet/pelelmex_prob.cpp +++ b/Exec/RegTests/FlameSheet/pelelmex_prob.cpp @@ -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(); } diff --git a/Exec/RegTests/FlameSheet/pelelmex_prob_parm.H b/Exec/RegTests/FlameSheet/pelelmex_prob_parm.H index 906b3d8e..27d610c7 100644 --- a/Exec/RegTests/FlameSheet/pelelmex_prob_parm.H +++ b/Exec/RegTests/FlameSheet/pelelmex_prob_parm.H @@ -3,6 +3,7 @@ #include #include +#include using namespace amrex::literals; @@ -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 const* + eosparm; }; #endif diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index c0c13f67..2afdf9a5 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -1753,6 +1753,11 @@ public: pele::physics::PhysicsType::transport_type>> trans_parms; + // EOS pointer + static pele::physics::PeleParams< + pele::physics::eos::EosParm> + eos_parms; + // Reactor pointer std::string m_chem_integrator; std::unique_ptr m_reactor; diff --git a/Source/PeleLMeX.cpp b/Source/PeleLMeX.cpp index 3dc3a2a3..d93059fb 100644 --- a/Source/PeleLMeX.cpp +++ b/Source/PeleLMeX.cpp @@ -12,12 +12,17 @@ pele::physics::PeleParams> PeleLM::trans_parms; +pele::physics::PeleParams< + pele::physics::eos::EosParm> + PeleLM::eos_parms; + PeleLM::PeleLM() = default; PeleLM::~PeleLM() { if (m_incompressible == 0) { trans_parms.deallocate(); + eos_parms.deallocate(); m_reactor->close(); } diff --git a/Source/PeleLMeX_Advection.cpp b/Source/PeleLMeX_Advection.cpp index 45d9d326..534d0c92 100644 --- a/Source/PeleLMeX_Advection.cpp +++ b/Source/PeleLMeX_Advection.cpp @@ -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()) @@ -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); }); } } @@ -552,6 +553,7 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& 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]; @@ -573,21 +575,21 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& 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); }); } } diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index bf310075..2ff984df 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -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); @@ -146,6 +147,7 @@ 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}; @@ -153,7 +155,7 @@ pelelmex_dermolefrac( 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]; @@ -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); }); } @@ -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); }); } diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 18c341f5..162d4450 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -366,6 +366,7 @@ PeleLM::addWbarTerm( // Compute Wbar on all the levels int nGrow = 1; // Need one ghost cell to compute gradWbar Vector 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)); @@ -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); }); } } @@ -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++) { @@ -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) { @@ -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); }); } } @@ -1226,6 +1229,7 @@ PeleLM::deltaTIter_prepare( std::unique_ptr& advData, std::unique_ptr& diffData) { + auto const* leosparm = eos_parms.device_parm(); for (int lev = 0; lev <= finest_level; ++lev) { auto* ldataOld_p = getLevelDataPtr(lev, AmrOldTime); @@ -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 @@ -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( @@ -1379,7 +1384,7 @@ PeleLM::deltaTIter_update( i, j, k, Array4(sma[box_no], DENSITY), Array4(sma[box_no], FIRSTSPEC), Array4(sma[box_no], TEMP), - Array4(sma[box_no], RHOH)); + Array4(sma[box_no], RHOH), leosparm); }); } Gpu::streamSynchronize(); diff --git a/Source/PeleLMeX_Eos.cpp b/Source/PeleLMeX_Eos.cpp index 0d5afd6f..628e5d70 100644 --- a/Source/PeleLMeX_Eos.cpp +++ b/Source/PeleLMeX_Eos.cpp @@ -24,6 +24,7 @@ PeleLM::setThermoPress(int lev, const TimeStamp& a_time) auto* ldata_p = getLevelDataPtr(lev, a_time); auto const& sma = ldata_p->state.arrays(); + auto const* leosparm = eos_parms.device_parm(); amrex::ParallelFor( ldata_p->state, @@ -31,8 +32,8 @@ PeleLM::setThermoPress(int lev, const TimeStamp& a_time) getPGivenRTY( i, j, k, Array4(sma[box_no], DENSITY), Array4(sma[box_no], FIRSTSPEC), - Array4(sma[box_no], TEMP), - Array4(sma[box_no], RHORT)); + Array4(sma[box_no], TEMP), Array4(sma[box_no], RHORT), + leosparm); }); Gpu::streamSynchronize(); } @@ -121,6 +122,7 @@ PeleLM::calcDivU( auto const& extRhoH = m_extSource[lev]->const_array(mfi, RHOH); auto const& divu = ldata_p->divu.array(mfi); int use_react = ((m_do_react != 0) && (m_skipInstantRR == 0)) ? 1 : 0; + auto const* leosparm = eos_parms.device_parm(); #ifdef AMREX_USE_EB if (flagfab.getType(bx) == FabType::covered) { // Covered boxes @@ -132,24 +134,26 @@ PeleLM::calcDivU( // boxes amrex::ParallelFor( bx, [rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, - use_react, flag] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + use_react, flag, + leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { if (flag(i, j, k).isCovered()) { divu(i, j, k) = 0.0; } else { compute_divu( i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, - divu, use_react); + divu, use_react, leosparm); } }); } else #endif { amrex::ParallelFor( - bx, [rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, - use_react] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bx, + [rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, + use_react, leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { compute_divu( i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, - divu, use_react); + divu, use_react, leosparm); }); } } @@ -219,6 +223,7 @@ PeleLM::setTemperature(int lev, const TimeStamp& a_time) auto* ldata_p = getLevelDataPtr(lev, a_time); auto const& sma = ldata_p->state.arrays(); + auto const* leosparm = eos_parms.device_parm(); amrex::ParallelFor( ldata_p->state, @@ -226,7 +231,8 @@ PeleLM::setTemperature(int lev, const TimeStamp& a_time) getTfromHY( i, j, k, Array4(sma[box_no], DENSITY), Array4(sma[box_no], FIRSTSPEC), - Array4(sma[box_no], RHOH), Array4(sma[box_no], TEMP)); + Array4(sma[box_no], RHOH), Array4(sma[box_no], TEMP), + leosparm); }); Gpu::streamSynchronize(); } @@ -286,6 +292,7 @@ PeleLM::adjustPandDivU(std::unique_ptr& advData) auto const& tma = ThetaHalft[lev]->arrays(); auto const& sma_o = getLevelDataPtr(lev, AmrOldTime)->state.const_arrays(); auto const& sma_n = getLevelDataPtr(lev, AmrNewTime)->state.const_arrays(); + auto const* leosparm = eos_parms.device_parm(); amrex::ParallelFor( *ThetaHalft[lev], [=, pOld = m_pOld, pNew = m_pNew] AMREX_GPU_DEVICE( @@ -293,10 +300,10 @@ PeleLM::adjustPandDivU(std::unique_ptr& advData) auto theta = tma[box_no]; Real gammaInv_o = getGammaInv( i, j, k, Array4(sma_o[box_no], FIRSTSPEC), - Array4(sma_o[box_no], TEMP)); + Array4(sma_o[box_no], TEMP), leosparm); Real gammaInv_n = getGammaInv( i, j, k, Array4(sma_n[box_no], FIRSTSPEC), - Array4(sma_n[box_no], TEMP)); + Array4(sma_n[box_no], TEMP), leosparm); theta(i, j, k) = 0.5 * (gammaInv_o / pOld + gammaInv_n / pNew); }); } diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 1c22d278..033855a3 100644 --- a/Source/PeleLMeX_Forces.cpp +++ b/Source/PeleLMeX_Forces.cpp @@ -204,15 +204,17 @@ PeleLM::addSpark(const TimeStamp& a_timestamp) Print() << m_spark[n] << " active" << std::endl; } - auto eos = pele::physics::PhysicsType::eos(); auto statema = getLevelDataPtr(lev, a_timestamp)->state.const_arrays(); auto extma = m_extSource[lev]->arrays(); + auto const* leosparm = eos_parms.device_parm(); amrex::ParallelFor( *m_extSource[lev], [=, spark_duration = m_spark_duration[n], spark_temp = m_spark_temp[n], + eosparm = leosparm, spark_radius = m_spark_radius [n]] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { + auto eos = pele::physics::PhysicsType::eos(eosparm); Real dist_to_center = std::sqrt(AMREX_D_TERM( (i - spark_idx[0]) * (i - spark_idx[0]) * dx[0] * dx[0], +(j - spark_idx[1]) * (j - spark_idx[1]) * dx[1] * dx[1], diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index 20404845..eece3cda 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -24,11 +24,13 @@ getTransportCoeff( amrex::Array4 const& mu, pele::physics::transport::TransParm< pele::physics::PhysicsType::eos_type, - pele::physics::PhysicsType::transport_type> const* trans_parm) noexcept + pele::physics::PhysicsType::transport_type> const* trans_parm, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real mwtinv[NUM_SPECIES] = {0.0}; eos.inv_molecular_weight(mwtinv); @@ -170,11 +172,13 @@ getPGivenRTY( amrex::Array4 const& rho, amrex::Array4 const& rhoY, amrex::Array4 const& T, - amrex::Array4 const& P) noexcept + amrex::Array4 const& P, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real rhoinv = 1.0_rt / rho(i, j, k); amrex::Real rho_cgs = rho(i, j, k) * 0.001_rt; @@ -203,11 +207,13 @@ compute_divu( amrex::Array4 const& extRhoY, amrex::Array4 const& extRhoH, amrex::Array4 const& divu, - int do_react) noexcept + int do_react, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real mwtinv[NUM_SPECIES] = {0.0}; eos.inv_molecular_weight(mwtinv); @@ -425,11 +431,13 @@ getMwmixGivenRY( int k, amrex::Array4 const& rho, amrex::Array4 const& rhoY, - amrex::Array4 const& Mwmix) noexcept + amrex::Array4 const& Mwmix, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real rhoinv = 1.0_rt / rho(i, j, k); amrex::Real y[NUM_SPECIES] = {0.0_rt}; for (int n = 0; n < NUM_SPECIES; n++) { @@ -591,11 +599,13 @@ getHGivenT( int j, int k, amrex::Array4 const& T, - amrex::Array4 const& Hi) noexcept + amrex::Array4 const& Hi, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real hi_spec[NUM_SPECIES] = {0.0_rt}; eos.T2Hi(T(i, j, k), hi_spec); for (int n = 0; n < NUM_SPECIES; n++) { @@ -665,11 +675,13 @@ getRHmixGivenTY( amrex::Array4 const& rho, amrex::Array4 const& rhoY, amrex::Array4 const& T, - amrex::Array4 const& Hmix) noexcept + amrex::Array4 const& Hmix, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real rhoinv = 1.0_rt / rho(i, j, k); amrex::Real y[NUM_SPECIES] = {0.0_rt}; for (int n = 0; n < NUM_SPECIES; n++) { @@ -691,11 +703,13 @@ getTfromHY( amrex::Array4 const& rho, amrex::Array4 const& rhoY, amrex::Array4 const& rhoH, - amrex::Array4 const& T) noexcept + amrex::Array4 const& T, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real rhoinv = 1.0_rt / rho(i, j, k); amrex::Real y[NUM_SPECIES] = {0.0}; for (int n = 0; n < NUM_SPECIES; n++) { @@ -717,11 +731,13 @@ getCpmixGivenRYT( amrex::Array4 const& rho, amrex::Array4 const& rhoY, amrex::Array4 const& T, - amrex::Array4 const& cpmix) noexcept + amrex::Array4 const& cpmix, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); amrex::Real rhoinv = 1.0_rt / rho(i, j, k); amrex::Real y[NUM_SPECIES] = {0.0}; for (int n = 0; n < NUM_SPECIES; n++) { @@ -751,11 +767,13 @@ buildAdvectionForcing( int const& closed_chamber, int do_react, amrex::Array4 const& forceY, - amrex::Array4 const& forceT) noexcept + amrex::Array4 const& forceT, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); // Get species enthalpy amrex::Real hi_spec[NUM_SPECIES] = {0.0}; eos.T2Hi(T(i, j, k), hi_spec); @@ -835,11 +853,13 @@ reactionRateRhoY( amrex::Array4 const& rhoY, amrex::Array4 const& rhoH, amrex::Array4 const& T, - amrex::Array4 const& rhoYdot) noexcept + amrex::Array4 const& rhoYdot, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); // Get rho & Ys from rhoYs. amrex::Real rho = 0.0_rt; @@ -1005,11 +1025,13 @@ getGammaInv( int j, int k, amrex::Array4 const& rhoY, - amrex::Array4 const& T) noexcept + amrex::Array4 const& T, + pele::physics::eos::EosParm const* + eosparm) noexcept { using namespace amrex::literals; - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(eosparm); // Get rho & Y from rhoY amrex::Real rho = 0.0_rt; for (int n = 0; n < NUM_SPECIES; n++) { diff --git a/Source/PeleLMeX_Plot.cpp b/Source/PeleLMeX_Plot.cpp index ebf9fe58..35b008b2 100644 --- a/Source/PeleLMeX_Plot.cpp +++ b/Source/PeleLMeX_Plot.cpp @@ -975,6 +975,7 @@ PeleLM::initLevelDataFromPlt(int a_lev, const std::string& a_dataPltFile) ldata_p->gp.setVal(0.0); ProbParm const* lprobparm = prob_parm_d; + auto const* leosparm = eos_parms.device_parm(); // If m_do_patch_flow_variables is set as true, call user-defined function to // patch flow variables @@ -994,38 +995,40 @@ PeleLM::initLevelDataFromPlt(int a_lev, const std::string& a_dataPltFile) auto const& rhoY_arr = ldata_p->state.array(mfi, FIRSTSPEC); auto const& rhoH_arr = ldata_p->state.array(mfi, RHOH); auto const& temp_arr = ldata_p->state.array(mfi, TEMP); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - auto eos = pele::physics::PhysicsType::eos(); - Real massfrac[NUM_SPECIES] = {0.0}; - Real sumYs = 0.0; - for (int n = 0; n < NUM_SPECIES; n++) { - massfrac[n] = rhoY_arr(i, j, k, n); + amrex::ParallelFor( + bx, + [=, eosparm = leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + auto eos = pele::physics::PhysicsType::eos(eosparm); + Real massfrac[NUM_SPECIES] = {0.0}; + Real sumYs = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + massfrac[n] = rhoY_arr(i, j, k, n); #ifdef N2_ID - if (n != N2_ID) { - sumYs += massfrac[n]; - } + if (n != N2_ID) { + sumYs += massfrac[n]; + } #endif - } + } #ifdef N2_ID - massfrac[N2_ID] = 1.0 - sumYs; + massfrac[N2_ID] = 1.0 - sumYs; #endif - // Get density - Real P_cgs = lprobparm->P_mean * 10.0; - Real rho_cgs = 0.0; - eos.PYT2R(P_cgs, massfrac, temp_arr(i, j, k), rho_cgs); - rho_arr(i, j, k) = rho_cgs * 1.0e3; + // Get density + Real P_cgs = lprobparm->P_mean * 10.0; + Real rho_cgs = 0.0; + eos.PYT2R(P_cgs, massfrac, temp_arr(i, j, k), rho_cgs); + rho_arr(i, j, k) = rho_cgs * 1.0e3; - // Get enthalpy - Real h_cgs = 0.0; - eos.TY2H(temp_arr(i, j, k), massfrac, h_cgs); - rhoH_arr(i, j, k) = h_cgs * 1.0e-4 * rho_arr(i, j, k); + // Get enthalpy + Real h_cgs = 0.0; + eos.TY2H(temp_arr(i, j, k), massfrac, h_cgs); + rhoH_arr(i, j, k) = h_cgs * 1.0e-4 * rho_arr(i, j, k); - // Fill rhoYs - for (int n = 0; n < NUM_SPECIES; n++) { - rhoY_arr(i, j, k, n) = massfrac[n] * rho_arr(i, j, k); - } - }); + // Fill rhoYs + for (int n = 0; n < NUM_SPECIES; n++) { + rhoY_arr(i, j, k, n) = massfrac[n] * rho_arr(i, j, k); + } + }); } // Initialize thermodynamic pressure diff --git a/Source/PeleLMeX_Reactions.cpp b/Source/PeleLMeX_Reactions.cpp index 6ac639c6..e862e025 100644 --- a/Source/PeleLMeX_Reactions.cpp +++ b/Source/PeleLMeX_Reactions.cpp @@ -84,7 +84,7 @@ PeleLM::advanceChemistry(int lev, const Real& a_dt, MultiFab& a_extForcing) auto const& FnE = a_extForcing.array(mfi, NUM_SPECIES + 1); auto const& rhoYe_n = ldataNew_p->state.array(mfi, FIRSTSPEC + E_ID); auto const& FrhoYe = a_extForcing.array(mfi, E_ID); - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(&eos_parms.host_parm()); Real mwt[NUM_SPECIES] = {0.0}; eos.molecular_weight(mwt); ParallelFor( @@ -245,7 +245,7 @@ PeleLM::advanceChemistryBAChem( auto const& FnE = chemForcing.array(mfi, NUM_SPECIES + 1); auto const& rhoYe_o = chemState.array(mfi, E_ID); auto const& FrhoYe = chemForcing.array(mfi, E_ID); - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(&eos_parms.host_parm()); Real mwt[NUM_SPECIES] = {0.0}; eos.molecular_weight(mwt); ParallelFor( @@ -385,6 +385,7 @@ PeleLM::computeInstantaneousReactionRate( { BL_PROFILE("PeleLMeX::computeInstantaneousReactionRate()"); auto* ldata_p = getLevelDataPtr(lev, a_time); + auto const* leosparm = eos_parms.device_parm(); #ifdef AMREX_USE_EB auto const& ebfact = EBFactory(lev); @@ -411,23 +412,23 @@ PeleLM::computeInstantaneousReactionRate( }); } else if (flagfab.getType(bx) != FabType::regular) { // EB containing boxes amrex::ParallelFor( - bx, [rhoY, rhoH, T, rhoYdot, - flag] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bx, [rhoY, rhoH, T, rhoYdot, flag, + leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { if (flag(i, j, k).isCovered()) { for (int n = 0; n < NUM_SPECIES; n++) { rhoYdot(i, j, k, n) = 0.0; } } else { - reactionRateRhoY(i, j, k, rhoY, rhoH, T, rhoYdot); + reactionRateRhoY(i, j, k, rhoY, rhoH, T, rhoYdot, leosparm); } }); } else #endif { amrex::ParallelFor( - bx, [rhoY, rhoH, T, - rhoYdot] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - reactionRateRhoY(i, j, k, rhoY, rhoH, T, rhoYdot); + bx, [rhoY, rhoH, T, rhoYdot, + leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + reactionRateRhoY(i, j, k, rhoY, rhoH, T, rhoYdot, leosparm); }); } } @@ -478,6 +479,7 @@ PeleLM::getHeatRelease(int a_lev, MultiFab* a_HR) { auto* ldataNew_p = getLevelDataPtr(a_lev, AmrNewTime); auto* ldataR_p = getLevelDataReactPtr(a_lev); + auto const* leosparm = eos_parms.device_parm(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -491,8 +493,9 @@ PeleLM::getHeatRelease(int a_lev, MultiFab* a_HR) auto const& Hi = EnthFab.array(); auto const& HRR = a_HR->array(mfi); amrex::ParallelFor( - bx, [T, Hi, HRR, react] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - getHGivenT(i, j, k, T, Hi); + bx, [T, Hi, HRR, react, + leosparm] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + getHGivenT(i, j, k, T, 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); diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index b719ddfb..4bb3c7a0 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -105,7 +105,14 @@ PeleLM::Setup() // Initialize EOS and others if (m_incompressible == 0) { + amrex::Print() << " Initialization of Eos ... \n"; + eos_parms.initialize(); + amrex::Print() << " Initialization of Transport ... \n"; +#ifdef USE_MANIFOLD_EOS + trans_parms.host_only_parm().manfunc_par = + eos_parms.host_only_parm().manfunc_par; +#endif trans_parms.initialize(); if ((m_les_verbose != 0) and m_do_les) { // Say what transport model we're // going to use @@ -142,6 +149,7 @@ PeleLM::Setup() m_reactor = pele::physics::reactions::ReactorBase::create(m_chem_integrator); m_reactor->init(reactor_type, ncells_chem); + m_reactor->set_eos_parm(eos_parms.device_parm()); // For ReactorNull, we need to also skip instantaneous RR used in divU if (m_chem_integrator == "ReactorNull") { m_skipInstantRR = 1; diff --git a/Source/PeleLMeX_TransportProp.cpp b/Source/PeleLMeX_TransportProp.cpp index 420670b8..e4d36c9c 100644 --- a/Source/PeleLMeX_TransportProp.cpp +++ b/Source/PeleLMeX_TransportProp.cpp @@ -64,6 +64,7 @@ PeleLM::calcTurbViscosity(const TimeStamp& a_time) } else { // get cp_cc (valid in 1 grow cell for interpolation to FCs) int ngrow = 1; + auto const* leosparm = eos_parms.device_parm(); cp_cc.define(ba, dm, 1, ngrow, MFInfo(), factory); auto const& state_arr = ldata_p->state.const_arrays(); auto const& cp_arr = cp_cc.arrays(); @@ -74,7 +75,7 @@ PeleLM::calcTurbViscosity(const TimeStamp& a_time) i, j, k, Array4(state_arr[box_no], DENSITY), Array4(state_arr[box_no], FIRSTSPEC), Array4(state_arr[box_no], TEMP), - Array4(cp_arr[box_no])); + Array4(cp_arr[box_no]), leosparm); }); Gpu::streamSynchronize(); @@ -292,6 +293,7 @@ PeleLM::calcDiffusivity(const TimeStamp& a_time) // Transport data pointer auto const* ltransparm = trans_parms.device_parm(); + auto const* leosparm = eos_parms.device_parm(); // MultiArrays auto const& sma = ldata_p->state.const_arrays(); @@ -300,7 +302,7 @@ PeleLM::calcDiffusivity(const TimeStamp& a_time) auto const& kma = ldata_p->mob_cc.arrays(); GpuArray mwt{0.0}; { - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(leosparm); eos.molecular_weight(mwt.arr); } #endif @@ -322,7 +324,7 @@ PeleLM::calcDiffusivity(const TimeStamp& a_time) Array4(sma[box_no], TEMP), Array4(dma[box_no], 0), Array4(dma[box_no], NUM_SPECIES + 1 + soret_idx), Array4(dma[box_no], NUM_SPECIES), - Array4(dma[box_no], NUM_SPECIES + 1), ltransparm); + Array4(dma[box_no], NUM_SPECIES + 1), ltransparm, leosparm); #ifdef PELE_USE_EFIELD getKappaSp( i, j, k, mwt.arr, zk, Array4(sma[box_no], FIRSTSPEC), diff --git a/Source/PeleLMeX_Utils.cpp b/Source/PeleLMeX_Utils.cpp index c35f82e0..0be2abcc 100644 --- a/Source/PeleLMeX_Utils.cpp +++ b/Source/PeleLMeX_Utils.cpp @@ -658,6 +658,7 @@ PeleLM::floorSpecies(const TimeStamp& a_time) auto* ldata_p = getLevelDataPtr(lev, a_time); auto const& sma = ldata_p->state.arrays(); + auto const* leosparm = eos_parms.device_parm(); amrex::ParallelFor( ldata_p->state, @@ -676,7 +677,7 @@ PeleLM::floorSpecies(const TimeStamp& a_time) } // ... as well as rhoh - auto eos = pele::physics::PhysicsType::eos(); + auto eos = pele::physics::PhysicsType::eos(leosparm); Real massfrac[NUM_SPECIES] = {0.0}; Real rhoinv = Real(1.0) / sma[box_no](i, j, k, DENSITY); for (int n = 0; n < NUM_SPECIES; n++) { @@ -1447,7 +1448,8 @@ PeleLM::updateTypicalValuesChem() } typical_values_chem[NUM_SPECIES] = typical_values[TEMP]; #ifdef PELE_USE_EFIELD - auto eos = pele::physics::PhysicsType::eos(); + auto const* leosparm = &eos_parms.host_parm(); + auto eos = pele::physics::PhysicsType::eos(leosparm); Real mwt[NUM_SPECIES] = {0.0}; eos.molecular_weight(mwt); typical_values_chem[E_ID] = @@ -1727,7 +1729,8 @@ PeleLM::initMixtureFraction() } } - auto eos = pele::physics::PhysicsType::eos(); + auto const* leosparm = &eos_parms.host_parm(); + auto eos = pele::physics::PhysicsType::eos(leosparm); // Overwrite with user-defined value if provided in input file ParmParse pp("peleLM"); std::string MFformat; @@ -1881,7 +1884,8 @@ PeleLM::parseComposition( massFrac[i] = compoIn[i]; } } else if (compositionType == "mole") { // mole - auto eos = pele::physics::PhysicsType::eos(); + auto const* leosparm = &eos_parms.host_parm(); + auto eos = pele::physics::PhysicsType::eos(leosparm); eos.X2Y(compoIn, massFrac); } else { Abort("Unknown mixtureFraction.type ! Should be 'mass' or 'mole'"); diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index aaede479..b7036447 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit aaede479e96ac87b8a61b5f6be425ad17c323062 +Subproject commit b70364474894ad24bee238fada521635d5eefaae