From db4d0c0606d1b3b9b8e26c1b40dda667f0b5f82c Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 13 Sep 2024 17:37:23 -0600 Subject: [PATCH 01/27] basics of having eos_parm in PeleLMeX --- Source/PeleLMeX.H | 6 ++++++ Source/PeleLMeX.cpp | 5 +++++ Source/PeleLMeX_Diffusion.cpp | 17 +++++++++++++++-- Source/PeleLMeX_Setup.cpp | 8 ++++++++ Submodules/PelePhysics | 2 +- 5 files changed, 35 insertions(+), 3 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index c0c13f673..c2ec7ddb1 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -606,6 +606,7 @@ 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 void adjustSpeciesFluxes( const amrex::Vector>& a_spfluxes, @@ -1753,6 +1754,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 3dc3a2a39..d93059fb1 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_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 18c341f58..8d19f047a 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -290,7 +290,8 @@ PeleLM::computeDifferentialDiffusionFluxes( } // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes(a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); + adjustSpeciesFluxes( + a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); //---------------------------------------------------------------- //---------------------------------------------------------------- @@ -637,6 +638,7 @@ PeleLM::addSoretTerm( } } +template void PeleLM::adjustSpeciesFluxes( const Vector>& a_spfluxes, @@ -745,6 +747,17 @@ PeleLM::adjustSpeciesFluxes( } } +template <> +void +PeleLM::adjustSpeciesFluxes( + const Vector>& /*a_spfluxes*/, + Vector const& /*a_spec*/) +{ + // Manifold Model: "Species" don't sum to unity, so no need to adjust the + // fluxes + BL_PROFILE("PeleLM::adjustSpeciesFluxes()"); +} + void PeleLM::computeSpeciesEnthalpyFlux( const Vector>& a_fluxes, @@ -1004,7 +1017,7 @@ PeleLM::differentialDiffusionUpdate( fillPatchSpecies(AmrNewTime); // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes( + adjustSpeciesFluxes( GetVecOfArrOfPtrs(fluxes), GetVecOfConstPtrs(getSpeciesVect(AmrNewTime))); // Average down fluxes^{np1,kp1} diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 1a76b5956..e07a43de5 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/Submodules/PelePhysics b/Submodules/PelePhysics index 060a7aeb9..809d138e1 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 060a7aeb97e6bfd3703646c479c9ed10b2bbe152 +Subproject commit 809d138e13aa80e3853535976b2e6ed45aef56cf From d09004327f1430ec9ea2b9ec0e0e5353fc59a18d Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 16 Sep 2024 14:07:38 -0600 Subject: [PATCH 02/27] add eosparm to PeleLMeX_K functions --- Source/PeleLMeX_Advection.cpp | 20 +++++----- Source/PeleLMeX_DeriveFunc.cpp | 21 +++++----- Source/PeleLMeX_Diffusion.cpp | 26 ++++++------ Source/PeleLMeX_Eos.cpp | 27 ++++++++----- Source/PeleLMeX_K.H | 66 ++++++++++++++++++++----------- Source/PeleLMeX_Reactions.cpp | 19 +++++---- Source/PeleLMeX_TransportProp.cpp | 6 ++- 7 files changed, 114 insertions(+), 71 deletions(-) diff --git a/Source/PeleLMeX_Advection.cpp b/Source/PeleLMeX_Advection.cpp index 45d9d326d..534d0c927 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 bf3100753..c19ba0e40 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); @@ -1371,17 +1372,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 +1420,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 8d19f047a..8882f6397 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -367,6 +367,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)); @@ -380,9 +381,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); }); } } @@ -768,6 +769,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) { @@ -799,21 +801,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); }); } } @@ -1239,6 +1241,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); @@ -1276,7 +1279,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 @@ -1383,6 +1386,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( @@ -1392,7 +1396,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 0d5afd6f7..628e5d702 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_K.H b/Source/PeleLMeX_K.H index 20404845a..eece3cdaf 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_Reactions.cpp b/Source/PeleLMeX_Reactions.cpp index 6ac639c62..5ef3760c7 100644 --- a/Source/PeleLMeX_Reactions.cpp +++ b/Source/PeleLMeX_Reactions.cpp @@ -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_TransportProp.cpp b/Source/PeleLMeX_TransportProp.cpp index 420670b8d..d95b7838b 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(); @@ -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), From 0b22bbfa1d61838cf4b3988116af2df2ab2413d6 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 16 Sep 2024 14:42:28 -0600 Subject: [PATCH 03/27] more eosparm: now everywhere except Efield & BCs --- Source/PeleLMeX_DeriveFunc.cpp | 3 +- Source/PeleLMeX_Diffusion.cpp | 9 +++--- Source/PeleLMeX_Forces.cpp | 4 ++- Source/PeleLMeX_Plot.cpp | 53 ++++++++++++++++--------------- Source/PeleLMeX_Reactions.cpp | 4 +-- Source/PeleLMeX_TransportProp.cpp | 2 +- Source/PeleLMeX_Utils.cpp | 12 ++++--- 7 files changed, 49 insertions(+), 38 deletions(-) diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index c19ba0e40..2ff984df2 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -147,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}; @@ -154,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]; diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 8882f6397..bf983ba5e 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -466,10 +466,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++) { diff --git a/Source/PeleLMeX_Forces.cpp b/Source/PeleLMeX_Forces.cpp index 1c22d278c..033855a31 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_Plot.cpp b/Source/PeleLMeX_Plot.cpp index ebf9fe58f..35b008b21 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 5ef3760c7..e862e0253 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( diff --git a/Source/PeleLMeX_TransportProp.cpp b/Source/PeleLMeX_TransportProp.cpp index d95b7838b..e4d36c9c6 100644 --- a/Source/PeleLMeX_TransportProp.cpp +++ b/Source/PeleLMeX_TransportProp.cpp @@ -302,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 diff --git a/Source/PeleLMeX_Utils.cpp b/Source/PeleLMeX_Utils.cpp index c35f82e01..0be2abcce 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'"); From a50c5523739216b586cab55b7374422ad08730b9 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 16 Sep 2024 16:13:57 -0600 Subject: [PATCH 04/27] eosparm in BC stuff --- Exec/RegTests/FlameSheet/pelelmex_prob.H | 4 ++-- Exec/RegTests/FlameSheet/pelelmex_prob.cpp | 1 + Exec/RegTests/FlameSheet/pelelmex_prob_parm.H | 3 +++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Exec/RegTests/FlameSheet/pelelmex_prob.H b/Exec/RegTests/FlameSheet/pelelmex_prob.H index f5c48e40a..c540c16d6 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 d861e2363..779cd18d2 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 906b3d8e4..27d610c7c 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 From af7813dda1ce38a98e60cd8d5be4fe274ef4cdd2 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 13:15:31 -0600 Subject: [PATCH 05/27] re-order functions for template --- Source/PeleLMeX_Diffusion.cpp | 240 +++++++++++++++++----------------- 1 file changed, 120 insertions(+), 120 deletions(-) diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index bf983ba5e..ba2ef2a45 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -204,6 +204,126 @@ PeleLM::computeDifferentialDiffusionTerms( #endif } +template +void +PeleLM::adjustSpeciesFluxes( + const Vector>& a_spfluxes, + Vector 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 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( + const Vector>& /*a_spfluxes*/, + Vector 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, @@ -640,126 +760,6 @@ PeleLM::addSoretTerm( } } -template -void -PeleLM::adjustSpeciesFluxes( - const Vector>& a_spfluxes, - Vector 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 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( - const Vector>& /*a_spfluxes*/, - Vector const& /*a_spec*/) -{ - // Manifold Model: "Species" don't sum to unity, so no need to adjust the - // fluxes - BL_PROFILE("PeleLM::adjustSpeciesFluxes()"); -} - void PeleLM::computeSpeciesEnthalpyFlux( const Vector>& a_fluxes, From 472c78d670cea1f89df7842498511cbee1d26e08 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 13:24:39 -0600 Subject: [PATCH 06/27] remove manifold-specfic changes, will go in later PR --- Source/PeleLMeX.H | 1 - Source/PeleLMeX_Diffusion.cpp | 233 ++++++++++++++++------------------ 2 files changed, 110 insertions(+), 124 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index c2ec7ddb1..2afdf9a54 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -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 void adjustSpeciesFluxes( const amrex::Vector>& a_spfluxes, diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index ba2ef2a45..162d44506 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -204,126 +204,6 @@ PeleLM::computeDifferentialDiffusionTerms( #endif } -template -void -PeleLM::adjustSpeciesFluxes( - const Vector>& a_spfluxes, - Vector 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 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( - const Vector>& /*a_spfluxes*/, - Vector 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, @@ -410,8 +290,7 @@ PeleLM::computeDifferentialDiffusionFluxes( } // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes( - a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); + adjustSpeciesFluxes(a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); //---------------------------------------------------------------- //---------------------------------------------------------------- @@ -760,6 +639,114 @@ PeleLM::addSoretTerm( } } +void +PeleLM::adjustSpeciesFluxes( + const Vector>& a_spfluxes, + Vector 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 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>& a_fluxes, @@ -1020,7 +1007,7 @@ PeleLM::differentialDiffusionUpdate( fillPatchSpecies(AmrNewTime); // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes( + adjustSpeciesFluxes( GetVecOfArrOfPtrs(fluxes), GetVecOfConstPtrs(getSpeciesVect(AmrNewTime))); // Average down fluxes^{np1,kp1} From f36161c73ae17f26caddb0e3d6ad7f549710656f Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 13:32:24 -0600 Subject: [PATCH 07/27] gnu make changes for manifold --- Exec/Make.PeleLMeX | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 3580e0fa3..844d1ff5b 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -73,6 +73,20 @@ endif ifeq ($(Eos_Model),$(filter $(Eos_Model),Soave-Redlich-Kwong)) DEFINES += -DUSE_SRK_EOS endif +ifeq ($(Eos_Model),$(filter $(Eos_Model),Manifold)) + DEFINES += -DUSE_MANIFOLD_EOS + DEFINES += -DUSE_MANIFOLD_EOS + DEFINES += -DMANIFOLD_DIM=$(Manifold_Dim) + ifeq ($(Manifold_Type),Table) + DEFINES += -DMANIFOLD_EOS_TYPE=1 + else + ifeq ($(Manifold_Type),Network) + DEFINES += -DMANIFOLD_EOS_TYPE=2 + else + DEFINES += -DMANIFOLD_EOS_TYPE=0 + endif + endif +endif # Transport model switches ifeq ($(Transport_Model), Simple) @@ -88,6 +102,9 @@ endif ifeq ($(Transport_Model), Sutherland) DEFINES += -DUSE_SUTHERLAND_TRANSPORT endif +ifeq ($(Transport_Model), Manifold) + DEFINES += -DUSE_MANIFOLD_TRANSPORT +endif ifeq ($(PELE_USE_KLU), TRUE) DEFINES += -DPELE_USE_KLU From 2c569863ea2d6c78cdee75ece76ffa85fbb73d8d Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 13:34:21 -0600 Subject: [PATCH 08/27] create manifold version of adjust fluxes function --- Source/PeleLMeX.H | 1 + Source/PeleLMeX_Diffusion.cpp | 233 ++++++++++++++++++---------------- 2 files changed, 124 insertions(+), 110 deletions(-) diff --git a/Source/PeleLMeX.H b/Source/PeleLMeX.H index 2afdf9a54..c2ec7ddb1 100644 --- a/Source/PeleLMeX.H +++ b/Source/PeleLMeX.H @@ -606,6 +606,7 @@ 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 void adjustSpeciesFluxes( const amrex::Vector>& a_spfluxes, diff --git a/Source/PeleLMeX_Diffusion.cpp b/Source/PeleLMeX_Diffusion.cpp index 162d44506..ba2ef2a45 100644 --- a/Source/PeleLMeX_Diffusion.cpp +++ b/Source/PeleLMeX_Diffusion.cpp @@ -204,6 +204,126 @@ PeleLM::computeDifferentialDiffusionTerms( #endif } +template +void +PeleLM::adjustSpeciesFluxes( + const Vector>& a_spfluxes, + Vector 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 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( + const Vector>& /*a_spfluxes*/, + Vector 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, @@ -290,7 +410,8 @@ PeleLM::computeDifferentialDiffusionFluxes( } // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes(a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); + adjustSpeciesFluxes( + a_fluxes, GetVecOfConstPtrs(getSpeciesVect(a_time))); //---------------------------------------------------------------- //---------------------------------------------------------------- @@ -639,114 +760,6 @@ PeleLM::addSoretTerm( } } -void -PeleLM::adjustSpeciesFluxes( - const Vector>& a_spfluxes, - Vector 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 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>& a_fluxes, @@ -1007,7 +1020,7 @@ PeleLM::differentialDiffusionUpdate( fillPatchSpecies(AmrNewTime); // Adjust species diffusion fluxes to ensure their sum is zero - adjustSpeciesFluxes( + adjustSpeciesFluxes( GetVecOfArrOfPtrs(fluxes), GetVecOfConstPtrs(getSpeciesVect(AmrNewTime))); // Average down fluxes^{np1,kp1} From f26906a7290bbe140c48fb3a39043f50c41ebbe2 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 14:26:08 -0600 Subject: [PATCH 09/27] Add abort for Closed Chamber + Manifold --- Source/PeleLMeX_Setup.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 4bb3c7a0e..b23d5e855 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -238,6 +238,12 @@ PeleLM::readParameters() if ((verbose != 0) && (m_closed_chamber != 0)) { Print() << " Simulation performed with the closed chamber algorithm \n"; } +#ifdef USE_MANIFOLD_EOS + if (m_closed_chamber) { + amrex::Abort( + "Simulation with closed chamber not supported for Manifold EOS"); + } +#endif #ifdef PELE_USE_EFIELD ParmParse ppef("ef"); From a8dc2fce1f5554d0669ab079e0f2db5cafefa0fc Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 14:40:59 -0600 Subject: [PATCH 10/27] Better setup chackes for Manifold + wbar, closed chamber, mixfrac/prog --- Source/PeleLMeX_Setup.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index b23d5e855..534d3e421 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -170,8 +170,10 @@ PeleLM::Setup() } // Mixture fraction & Progress variable - initMixtureFraction(); - initProgressVariable(); + if (pele::physics::PhysicsType::eos_type::identifier() != "Manifold") { + initMixtureFraction(); + initProgressVariable(); + } // Initialize turbulence injection turb_inflow.init(Geom(0)); @@ -238,12 +240,12 @@ PeleLM::readParameters() if ((verbose != 0) && (m_closed_chamber != 0)) { Print() << " Simulation performed with the closed chamber algorithm \n"; } -#ifdef USE_MANIFOLD_EOS - if (m_closed_chamber) { + if ( + m_closed_chamber && + pele::physics::PhysicsType::eos_type::identifier() == "Manifold") { amrex::Abort( "Simulation with closed chamber not supported for Manifold EOS"); } -#endif #ifdef PELE_USE_EFIELD ParmParse ppef("ef"); @@ -462,6 +464,11 @@ PeleLM::readParameters() "fixed_Pr or fixed_Le is true" << std::endl; } + if ( + m_use_wbar != 0 && + pele::physics::PhysicsType::eos_type::identifier() == "Manifold") { + amrex::Abort("Use of Wbar fluxes is not compatible with Manifold EOS"); + } pp.query("deltaT_verbose", m_deltaT_verbose); pp.query("deltaT_iterMax", m_deltaTIterMax); From 1a0a7e2cb8d20174a6db3950de5a102313068a20 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 17 Sep 2024 16:10:49 -0600 Subject: [PATCH 11/27] advection of rho corrections for manifold models --- Source/PeleLMeX_Advection.cpp | 19 +++++++------------ Source/PeleLMeX_Utils.H | 13 +++++++++++++ 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/Source/PeleLMeX_Advection.cpp b/Source/PeleLMeX_Advection.cpp index 534d0c927..da5f9520e 100644 --- a/Source/PeleLMeX_Advection.cpp +++ b/Source/PeleLMeX_Advection.cpp @@ -488,9 +488,8 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) afrac] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { rho_ed(i, j, k) = 0.0; if (afrac(i, j, k) > 0.0) { // Uncovered faces - for (int n = 0; n < NUM_SPECIES; n++) { - rho_ed(i, j, k) += rhoY_ed(i, j, k, n); - } + pele::physics::PhysicsType::eos_type::RY2R( + array4_to_array(i, j, k, rhoY_ed).data(), rho_ed(i, j, k)); } }); } else // Regular boxes @@ -499,10 +498,8 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) amrex::ParallelFor( ebx, [rho_ed, rhoY_ed] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - rho_ed(i, j, k) = 0.0; - for (int n = 0; n < NUM_SPECIES; n++) { - rho_ed(i, j, k) += rhoY_ed(i, j, k, n); - } + pele::physics::PhysicsType::eos_type::RY2R( + array4_to_array(i, j, k, rhoY_ed).data(), rho_ed(i, j, k)); }); } } @@ -751,11 +748,9 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) amrex::ParallelFor( advData->AofS[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - aofsma[box_no](i, j, k, DENSITY) = 0.0; - for (int n = 0; n < NUM_SPECIES; n++) { - aofsma[box_no](i, j, k, DENSITY) += - aofsma[box_no](i, j, k, FIRSTSPEC + n); - } + pele::physics::PhysicsType::eos_type::RY2R( + array4_to_array(i, j, k, aofsma[box_no], FIRSTSPEC).data(), + aofsma[box_no](i, j, k, DENSITY)); }); } Gpu::streamSynchronize(); diff --git a/Source/PeleLMeX_Utils.H b/Source/PeleLMeX_Utils.H index f38e182c2..169782a05 100644 --- a/Source/PeleLMeX_Utils.H +++ b/Source/PeleLMeX_Utils.H @@ -1,5 +1,6 @@ #ifndef PELELM_UTILS_H #define PELELM_UTILS_H + template amrex::Gpu::DeviceVector convertToDeviceVector(amrex::Vector v) @@ -13,4 +14,16 @@ convertToDeviceVector(amrex::Vector v) #endif return v_d; } + +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::GpuArray +array4_to_array(int i, int j, int k, amrex::Array4 array, int start_comp = 0) +{ + amrex::GpuArray celldata; + for (unsigned int n = start_comp; n < N; i++) { + celldata[n] = array(i, j, k, n); + } + return celldata; +} + #endif From d55daa0b858bb46f2d804b3d0bded48bf596d1a9 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Wed, 18 Sep 2024 10:40:57 -0600 Subject: [PATCH 12/27] manifold-specific transport, derives, divu --- Source/PeleLMeX_DeriveFunc.H | 16 ++++ Source/PeleLMeX_DeriveFunc.cpp | 55 +++++++++++- Source/PeleLMeX_Eos.cpp | 4 +- Source/PeleLMeX_K.H | 137 +++++++++++++++++++++++++++--- Source/PeleLMeX_Setup.cpp | 74 +++++++++------- Source/PeleLMeX_TransportProp.cpp | 2 +- Source/PeleLMeX_Utils.H | 7 +- 7 files changed, 244 insertions(+), 51 deletions(-) diff --git a/Source/PeleLMeX_DeriveFunc.H b/Source/PeleLMeX_DeriveFunc.H index 7d74421fc..9f60602e1 100644 --- a/Source/PeleLMeX_DeriveFunc.H +++ b/Source/PeleLMeX_DeriveFunc.H @@ -290,6 +290,22 @@ void pelelmex_deruserdef( const amrex::Vector& bcrec, int level); +#ifdef USE_MANIFOLD_EOS +void pelelmex_dermaniout( + PeleLM* a_pelelm, + const amrex::Box& bx, + amrex::FArrayBox& derfab, + int dcomp, + int ncomp, + const amrex::FArrayBox& statefab, + const amrex::FArrayBox& reactfab, + const amrex::FArrayBox& pressfab, + const amrex::Geometry& geom, + amrex::Real time, + const amrex::Vector& bcrec, + int level); +#endif + #ifdef PELE_USE_EFIELD #include #endif diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index 2ff984df2..76a1a3c12 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -2,6 +2,7 @@ #include "PeleLMeX.H" #include "PeleLMeX_K.H" #include "PeleLMeX_DeriveFunc.H" +#include "PeleLMeX_Utils.H" #include #include @@ -1382,7 +1383,7 @@ pelelmex_derdiffc( 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( + getTransportCoeff( i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD, rhotheta, lambda, mu, ltransparm, leosparm); }); @@ -1428,7 +1429,7 @@ pelelmex_derlambda( 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( + getTransportCoeff( i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, LeInv, PrInv, rhoY, T, rhoD, rhotheta, lambda, mu, ltransparm, leosparm); }); @@ -1451,7 +1452,6 @@ pelelmex_derdmap( Real /*time*/, const Vector& /*bcrec*/, int /*level*/) - { AMREX_ASSERT(derfab.box().contains(bx)); auto der = derfab.array(dcomp); @@ -1460,3 +1460,52 @@ pelelmex_derdmap( der(i, j, k) = myrank; }); } + +// +// Derive manifold output quantities +// +#ifdef USE_MANIFOLD_EOS +void +pelelmex_dermaniout( + PeleLM* a_pelelm, + const Box& bx, + FArrayBox& derfab, + int dcomp, + int ncomp, + const FArrayBox& statefab, + const FArrayBox& /*reactfab*/, + const FArrayBox& /*pressfab*/, + const Geometry& /*geom*/, + Real /*time*/, + const Vector& /*bcrec*/, + int /*level*/) +{ + auto& h_manf_data = + a_pelelm->eos_parms.host_only_parm().manfunc_par->host_parm(); + auto* d_manf_data = + a_pelelm->eos_parms.host_only_parm().manfunc_par->device_parm(); + int nmanivar = h_manf_data.Nvar; + + amrex::ignore_unused(ncomp); + AMREX_ASSERT(derfab.box().contains(bx)); + AMREX_ASSERT(statefab.box().contains(bx)); + AMREX_ASSERT(derfab.nComp() >= dcomp + ncomp); + AMREX_ASSERT(statefab.nComp() >= NUM_SPECIES + 1); + AMREX_ASSERT(ncomp == nmanivar); + AMREX_ASSERT(!a_pelelm->m_incompressible); + + auto const in_dat = statefab.array(); + auto der = derfab.array(dcomp); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::Real rho, rhoinv, maniparm[NUM_SPECIES]; + // pele::physics::PhysicsType::eos_type::RY2RRinvY(array4_to_array(i, j, k, + // in_dat, FIRSTSPEC).data(), rho, rhoinv, maniparm); + pele::physics::PhysicsType::eos_type::RY2RRinvY( + in_dat.ptr(i, j, k, FIRSTSPEC), rho, rhoinv, maniparm); + pele::physics::BlackBoxFunctionFactory< + pele::physics::eos::ManifoldFunctionType> + manfunc{d_manf_data}; + manfunc.get_func()->get_all_values(maniparm, der.ptr(i, j, k)); + }); +} +#endif diff --git a/Source/PeleLMeX_Eos.cpp b/Source/PeleLMeX_Eos.cpp index 628e5d702..83c34e7a9 100644 --- a/Source/PeleLMeX_Eos.cpp +++ b/Source/PeleLMeX_Eos.cpp @@ -139,7 +139,7 @@ PeleLM::calcDivU( if (flag(i, j, k).isCovered()) { divu(i, j, k) = 0.0; } else { - compute_divu( + compute_divu( i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, use_react, leosparm); } @@ -151,7 +151,7 @@ PeleLM::calcDivU( 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( + compute_divu( i, j, k, rhoY, T, SpecD, Fourier, DiffDiff, r, extRhoY, extRhoH, divu, use_react, leosparm); }); diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index eece3cdaf..a33aa35a2 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -4,9 +4,8 @@ #include #include -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void getTransportCoeff( int i, int j, @@ -113,6 +112,70 @@ getTransportCoeff( } } +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +getTransportCoeff( + int i, + int j, + int k, + const bool /*do_fixed_Le*/, + const bool /*do_fixed_Pr*/, + const bool do_soret, + amrex::Real /*LeInv*/, + amrex::Real /*PrInv*/, + amrex::Array4 const& rhoY, + amrex::Array4 const& T, + amrex::Array4 const& rhoDi, + amrex::Array4 const& /*rhotheta*/, + amrex::Array4 const& lambda, + amrex::Array4 const& mu, + pele::physics::transport::TransParm< + pele::physics::PhysicsType::eos_type, + pele::physics::PhysicsType::transport_type> const* trans_parm, + const pele::physics::eos::EosParm* + eosparm) noexcept +{ + using namespace amrex::literals; + + // No Soret with manifold moels + AMREX_ASSERT(!do_soret); + + // Get rho & Y from rhoY + amrex::Real rho, rhoinv; + amrex::Real y[NUM_SPECIES] = {0.0}; + amrex::Real massdens[NUM_SPECIES]; + for (int n = 0; n < NUM_SPECIES; n++) { + massdens[n] = rhoY(i, j, k, n); + } + auto eos = pele::physics::PhysicsType::eos(eosparm); + eos.RY2RRinvY(massdens, rho, rhoinv, y); + + rho *= 1.0e-3_rt; // MKS -> CGS conversion + amrex::Real rhoDi_cgs[NUM_SPECIES] = {0.0}; + amrex::Real lambda_cgs = 0.0_rt; + amrex::Real mu_cgs = 0.0_rt; + amrex::Real dummy_xi = 0.0_rt; + amrex::Real* dummy_chi = nullptr; + amrex::Real Tloc = T(i, j, k); + + bool get_xi = false; + bool get_mu = true; + bool get_lam = false; + bool get_Ddiag = true; + bool get_chi = do_soret; + auto trans = pele::physics::PhysicsType::transport(); + trans.transport( + get_xi, get_mu, get_lam, get_Ddiag, get_chi, Tloc, rho, y, rhoDi_cgs, + dummy_chi, mu_cgs, dummy_xi, lambda_cgs, trans_parm); + + // Do CGS -> MKS conversions + for (int n = 0; n < NUM_SPECIES; n++) { + rhoDi(i, j, k, n) = rhoDi_cgs[n] * 1.0e-1_rt; + } + lambda(i, j, k) = 0.0; // No need to carry lambda for manifold + mu(i, j, k) = mu_cgs * 1.0e-1_rt; +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -130,15 +193,14 @@ getVelViscosity( using namespace amrex::literals; // Get rho & Y from rhoY - amrex::Real rho = 0.0_rt; + amrex::Real massdens[NUM_SPECIES]; for (int n = 0; n < NUM_SPECIES; n++) { - rho += rhoY(i, j, k, n); + massdens[n] = rhoY(i, j, k, n); } + amrex::Real rho = 0.0_rt; amrex::Real rhoinv = 1.0_rt / rho; amrex::Real y[NUM_SPECIES] = {0.0}; - for (int n = 0; n < NUM_SPECIES; n++) { - y[n] = rhoY(i, j, k, n) * rhoinv; - } + pele::physics::PhysicsType::eos_type::RY2RRinvY(massdens, rho, rhoinv, y); rho *= 1.0e-3_rt; // MKS -> CGS conversion amrex::Real temp = T(i, j, k); @@ -191,9 +253,8 @@ getPGivenRTY( P(i, j, k) = P(i, j, k) * 0.1_rt; // CGS -> MKS conversion } -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -void +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void compute_divu( int i, int j, @@ -251,6 +312,57 @@ compute_divu( } } +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void +compute_divu( + int i, + int j, + int k, + amrex::Array4 const& rhoY, + amrex::Array4 const& /*T*/, + amrex::Array4 const& specDiff, + amrex::Array4 const& /*tempDiff*/, + amrex::Array4 const& /*specEnthDiff*/, + amrex::Array4 const& rhoYdot, + amrex::Array4 const& extRhoY, + amrex::Array4 const& /*extRhoH*/, + amrex::Array4 const& divu, + int do_react, + pele::physics::eos::EosParm const* + eosparm) noexcept +{ + using namespace amrex::literals; + + // Get rho & Ys from rhoYs. + amrex::Real rho, rhoinv; + amrex::Real y[NUM_SPECIES] = {0.0_rt}; + amrex::Real massdens[NUM_SPECIES] = {0.0_rt}; + for (int n = 0; n < NUM_SPECIES; n++) { + massdens[n] = rhoY(i, j, k, n); + } + + auto eos = pele::physics::PhysicsType::eos(eosparm); + eos.RY2RRinvY(massdens, rho, rhoinv, y); + amrex::Real derivs[NUM_SPECIES - 1] = {0.0}; + rho /= 1.0e3; // Unit conversion MKS -> CGS + eos.RY2dRdY(rho, y, derivs); + + divu(i, j, k) = 0.0; + + for (int n = 0; n < NUM_SPECIES - 1; n++) { + derivs[n] *= 1.0e3; // Unit conversion CGS -> MKS + divu(i, j, k) += derivs[n] * specDiff(i, j, k, n); + if (do_react != 0) { + divu(i, j, k) += derivs[n] * (rhoYdot(i, j, k, n)); + } + divu(i, j, k) += derivs[n] * (extRhoY(i, j, k, n)); + } + // First -rhoinv: from divu = -1/rho Drho/Dt + // Second rhoinv: convert from rho*Dxi/Dt (which RHS terms are evaluated for) + // to Dxi/Dt + divu(i, j, k) *= -rhoinv * rhoinv; +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -624,7 +736,8 @@ makeVelForce( const amrex::Real& rho_incomp, int pseudo_gravity, int pseudo_gravity_dir, - const amrex::Real& /*time*/, + const amrex::Real& /*time*/ + , amrex::GpuArray const gravity, amrex::GpuArray const gp0, const amrex::Real& dV_control, diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 534d3e421..578fbef5d 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -75,34 +75,6 @@ PeleLM::Setup() // Setup the state variables variablesSetup(); - // Derived variables - derivedSetup(); - - // Evaluate variables - evaluateSetup(); - - // Tagging setup - taggingSetup(); - -#ifdef PELE_USE_SPRAY - SpraySetup(); -#endif -#ifdef PELE_USE_SOOT - if (do_soot_solve) { - soot_model->define(); - } -#endif - // Diagnostics setup - createDiagnostics(); - - // Boundary Patch Setup - if (m_do_patch_mfr != 0) { - initBPatches(Geom(0)); - } - - // Initialize Level Hierarchy data - resizeArray(); - // Initialize EOS and others if (m_incompressible == 0) { amrex::Print() << " Initialization of Eos ... \n"; @@ -169,6 +141,34 @@ PeleLM::Setup() #endif } + // Derived variables + derivedSetup(); + + // Evaluate variables + evaluateSetup(); + + // Tagging setup + taggingSetup(); + +#ifdef PELE_USE_SPRAY + SpraySetup(); +#endif +#ifdef PELE_USE_SOOT + if (do_soot_solve) { + soot_model->define(); + } +#endif + // Diagnostics setup + createDiagnostics(); + + // Boundary Patch Setup + if (m_do_patch_mfr != 0) { + initBPatches(Geom(0)); + } + + // Initialize Level Hierarchy data + resizeArray(); + // Mixture fraction & Progress variable if (pele::physics::PhysicsType::eos_type::identifier() != "Manifold") { initMixtureFraction(); @@ -241,8 +241,8 @@ PeleLM::readParameters() Print() << " Simulation performed with the closed chamber algorithm \n"; } if ( - m_closed_chamber && - pele::physics::PhysicsType::eos_type::identifier() == "Manifold") { + (m_closed_chamber != 0) && + (pele::physics::PhysicsType::eos_type::identifier() == "Manifold")) { amrex::Abort( "Simulation with closed chamber not supported for Manifold EOS"); } @@ -1114,6 +1114,20 @@ PeleLM::derivedSetup() "enstrophy", IndexType::TheCellType(), 1, pelelmex_derenstrophy, grow_box_by_two); +#ifdef USE_MANIFOLD_EOS + auto& mani_data = eos_parms.host_only_parm().manfunc_par->host_parm(); + int nmanivar = mani_data.Nvar; + Vector var_names_maniout(nmanivar); + for (int n = 0; n < nmanivar; n++) { + std::string nametmp = std::string( + &(mani_data.varnames)[n * mani_data.len_str], mani_data.len_str); + var_names_maniout[n] = "MANI_" + amrex::trim(nametmp); + } + derive_lst.add( + "maniout", IndexType::TheCellType(), nmanivar, var_names_maniout, + pelelmex_dermaniout, the_same_box); +#endif + #ifdef PELE_USE_EFIELD // Charge distribution derive_lst.add( diff --git a/Source/PeleLMeX_TransportProp.cpp b/Source/PeleLMeX_TransportProp.cpp index e4d36c9c6..23a77f96e 100644 --- a/Source/PeleLMeX_TransportProp.cpp +++ b/Source/PeleLMeX_TransportProp.cpp @@ -318,7 +318,7 @@ PeleLM::calcDiffusivity(const TimeStamp& a_time) amrex::ParallelFor( ldata_p->diff_cc, ldata_p->diff_cc.nGrowVect(), [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { - getTransportCoeff( + getTransportCoeff( i, j, k, do_fixed_Le, do_fixed_Pr, do_soret, Le_inv, Pr_inv, Array4(sma[box_no], FIRSTSPEC), Array4(sma[box_no], TEMP), Array4(dma[box_no], 0), diff --git a/Source/PeleLMeX_Utils.H b/Source/PeleLMeX_Utils.H index 169782a05..e45210962 100644 --- a/Source/PeleLMeX_Utils.H +++ b/Source/PeleLMeX_Utils.H @@ -17,11 +17,12 @@ convertToDeviceVector(amrex::Vector v) template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::GpuArray -array4_to_array(int i, int j, int k, amrex::Array4 array, int start_comp = 0) +array4_to_array( + int i, int j, int k, const amrex::Array4& array, int start_comp = 0) { amrex::GpuArray celldata; - for (unsigned int n = start_comp; n < N; i++) { - celldata[n] = array(i, j, k, n); + for (unsigned int n = 0; n < N; n++) { + celldata[n] = array(i, j, k, start_comp + n); } return celldata; } From e5a3b115f5bb02620aeb920b01b931546b063384 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 18 Oct 2024 14:27:18 -0600 Subject: [PATCH 13/27] remove divide by zero --- Source/PeleLMeX_K.H | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index a33aa35a2..5385fb32a 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -197,9 +197,7 @@ getVelViscosity( for (int n = 0; n < NUM_SPECIES; n++) { massdens[n] = rhoY(i, j, k, n); } - amrex::Real rho = 0.0_rt; - amrex::Real rhoinv = 1.0_rt / rho; - amrex::Real y[NUM_SPECIES] = {0.0}; + amrex::Real rho = 0.0_rt, rhoinv = 0.0_rt, y[NUM_SPECIES] = {0.0}; pele::physics::PhysicsType::eos_type::RY2RRinvY(massdens, rho, rhoinv, y); rho *= 1.0e-3_rt; // MKS -> CGS conversion From 7f48aec0f706c1cfe2113f89e77f92c693309ba9 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 18 Oct 2024 14:36:22 -0600 Subject: [PATCH 14/27] go to templated EOS functions to use cellData --- Source/PeleLMeX_Advection.cpp | 8 ++++---- Source/PeleLMeX_DeriveFunc.cpp | 7 ++----- Source/PeleLMeX_Utils.H | 12 ------------ Submodules/PelePhysics | 2 +- 4 files changed, 7 insertions(+), 22 deletions(-) diff --git a/Source/PeleLMeX_Advection.cpp b/Source/PeleLMeX_Advection.cpp index da5f9520e..9cf524187 100644 --- a/Source/PeleLMeX_Advection.cpp +++ b/Source/PeleLMeX_Advection.cpp @@ -489,7 +489,7 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) rho_ed(i, j, k) = 0.0; if (afrac(i, j, k) > 0.0) { // Uncovered faces pele::physics::PhysicsType::eos_type::RY2R( - array4_to_array(i, j, k, rhoY_ed).data(), rho_ed(i, j, k)); + rhoY_ed.cellData(i, j, k), rho_ed(i, j, k)); } }); } else // Regular boxes @@ -499,7 +499,7 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) ebx, [rho_ed, rhoY_ed] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { pele::physics::PhysicsType::eos_type::RY2R( - array4_to_array(i, j, k, rhoY_ed).data(), rho_ed(i, j, k)); + rhoY_ed.cellData(i, j, k), rho_ed(i, j, k)); }); } } @@ -749,8 +749,8 @@ PeleLM::computeScalarAdvTerms(std::unique_ptr& advData) advData->AofS[lev], [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept { pele::physics::PhysicsType::eos_type::RY2R( - array4_to_array(i, j, k, aofsma[box_no], FIRSTSPEC).data(), - aofsma[box_no](i, j, k, DENSITY)); + aofsma[box_no].cellData(i, j, k), aofsma[box_no](i, j, k, DENSITY), + FIRSTSPEC); }); } Gpu::streamSynchronize(); diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index 76a1a3c12..87ac8e24a 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -1494,14 +1494,11 @@ pelelmex_dermaniout( AMREX_ASSERT(ncomp == nmanivar); AMREX_ASSERT(!a_pelelm->m_incompressible); - auto const in_dat = statefab.array(); + auto const in_spec = statefab.array(FIRST_SPEC); auto der = derfab.array(dcomp); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real rho, rhoinv, maniparm[NUM_SPECIES]; - // pele::physics::PhysicsType::eos_type::RY2RRinvY(array4_to_array(i, j, k, - // in_dat, FIRSTSPEC).data(), rho, rhoinv, maniparm); - pele::physics::PhysicsType::eos_type::RY2RRinvY( - in_dat.ptr(i, j, k, FIRSTSPEC), rho, rhoinv, maniparm); + pele::physics::PhysicsType::eos_type::RY2RRinvY(in_spec.cellData(i, j, k) FIRSTSPEC), rho, rhoinv, maniparm); pele::physics::BlackBoxFunctionFactory< pele::physics::eos::ManifoldFunctionType> manfunc{d_manf_data}; diff --git a/Source/PeleLMeX_Utils.H b/Source/PeleLMeX_Utils.H index e45210962..66fb0f18a 100644 --- a/Source/PeleLMeX_Utils.H +++ b/Source/PeleLMeX_Utils.H @@ -15,16 +15,4 @@ convertToDeviceVector(amrex::Vector v) return v_d; } -template -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::GpuArray -array4_to_array( - int i, int j, int k, const amrex::Array4& array, int start_comp = 0) -{ - amrex::GpuArray celldata; - for (unsigned int n = 0; n < N; n++) { - celldata[n] = array(i, j, k, start_comp + n); - } - return celldata; -} - #endif diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index 45ae8ba14..563cb3c19 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 45ae8ba14f2a14e8f692c83a3b675883f7018bfd +Subproject commit 563cb3c198d46b40c9d5955381beee637d575d68 From 5acd662380503f535fcbcac4cd765e39acbb820c Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Fri, 18 Oct 2024 15:06:36 -0600 Subject: [PATCH 15/27] fix oops in pelephysics --- Submodules/PelePhysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index 563cb3c19..bef80356f 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 563cb3c198d46b40c9d5955381beee637d575d68 +Subproject commit bef80356f2cf7001b659e764c710a42743cb2381 From b3491d660ab3e390190ccb37e4dc78983a09b035 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Wed, 23 Oct 2024 13:06:22 -0600 Subject: [PATCH 16/27] fix compile errors for Manifold derives --- Source/PeleLMeX_DeriveFunc.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index 87ac8e24a..0837cc5a6 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -1494,11 +1494,12 @@ pelelmex_dermaniout( AMREX_ASSERT(ncomp == nmanivar); AMREX_ASSERT(!a_pelelm->m_incompressible); - auto const in_spec = statefab.array(FIRST_SPEC); + auto const in_spec = statefab.array(FIRSTSPEC); auto der = derfab.array(dcomp); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real rho, rhoinv, maniparm[NUM_SPECIES]; - pele::physics::PhysicsType::eos_type::RY2RRinvY(in_spec.cellData(i, j, k) FIRSTSPEC), rho, rhoinv, maniparm); + pele::physics::PhysicsType::eos_type::RY2RRinvY( + in_spec.cellData(i, j, k), rho, rhoinv, maniparm); pele::physics::BlackBoxFunctionFactory< pele::physics::eos::ManifoldFunctionType> manfunc{d_manf_data}; From 424b17b8181e50f1706fce88fa86c8b35e91238e Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 13:34:06 -0600 Subject: [PATCH 17/27] fix compile condition for transparm initialization --- Source/PeleLMeX_Setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 578fbef5d..35db10fe7 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -81,7 +81,7 @@ PeleLM::Setup() eos_parms.initialize(); amrex::Print() << " Initialization of Transport ... \n"; -#ifdef USE_MANIFOLD_EOS +#ifdef USE_MANIFOLD_TRANSPORT trans_parms.host_only_parm().manfunc_par = eos_parms.host_only_parm().manfunc_par; #endif From 94cd7ff3f38b8d6069f6e050f1f180397a168ee8 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 13:36:21 -0600 Subject: [PATCH 18/27] fix missing RY2RRinvY conversion --- Source/PeleLMeX_K.H | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index 5385fb32a..bd91b4c2c 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -973,15 +973,8 @@ reactionRateRhoY( auto eos = pele::physics::PhysicsType::eos(eosparm); // Get rho & Ys from rhoYs. - amrex::Real rho = 0.0_rt; - for (int n = 0; n < NUM_SPECIES; n++) { - rho += rhoY(i, j, k, n); - } - amrex::Real rhoinv = 1.0_rt / rho; - amrex::Real y[NUM_SPECIES] = {0.0_rt}; - for (int n = 0; n < NUM_SPECIES; n++) { - y[n] = rhoY(i, j, k, n) * rhoinv; - } + amrex::Real rho = 0.0, rhoinv = 0.0, y[NUM_SPECIES] = {0.0}; + eos.RY2RRinvY(rhoY.cellData(i,j,k), rho, rhoinv, y); // Get T from Y/H. amrex::Real Tloc = T(i, j, k); From d7fce1a29847db5cbd464af4c674c0c2671e3223 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 13:37:58 -0600 Subject: [PATCH 19/27] use PelePhysics that free manfunc shared pointer --- Submodules/PelePhysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index bef80356f..683ca7998 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit bef80356f2cf7001b659e764c710a42743cb2381 +Subproject commit 683ca7998d121ca8cc31c294ece942e55424f2f8 From 7ac6731317da3a4df127d3f5f04782a17987786f Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 13:44:40 -0600 Subject: [PATCH 20/27] clang-formatting --- Source/PeleLMeX_K.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index bd91b4c2c..0a515de2b 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -974,7 +974,7 @@ reactionRateRhoY( // Get rho & Ys from rhoYs. amrex::Real rho = 0.0, rhoinv = 0.0, y[NUM_SPECIES] = {0.0}; - eos.RY2RRinvY(rhoY.cellData(i,j,k), rho, rhoinv, y); + eos.RY2RRinvY(rhoY.cellData(i, j, k), rho, rhoinv, y); // Get T from Y/H. amrex::Real Tloc = T(i, j, k); From c9f511822df4c9a058bd44fa635d2bef671aa70b Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 14:19:12 -0600 Subject: [PATCH 21/27] fux derives for manifold --- Source/PeleLMeX_DeriveFunc.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Source/PeleLMeX_DeriveFunc.cpp b/Source/PeleLMeX_DeriveFunc.cpp index 0837cc5a6..edf605a19 100644 --- a/Source/PeleLMeX_DeriveFunc.cpp +++ b/Source/PeleLMeX_DeriveFunc.cpp @@ -1503,7 +1503,11 @@ pelelmex_dermaniout( pele::physics::BlackBoxFunctionFactory< pele::physics::eos::ManifoldFunctionType> manfunc{d_manf_data}; - manfunc.get_func()->get_all_values(maniparm, der.ptr(i, j, k)); + + // TODO: use get_all_values instead + for (int n = 0; n < nmanivar; ++n) { + manfunc.get_func()->get_value(n, maniparm, der(i, j, k, n)); + } }); } #endif From 962fcd15bd83f1077cabebc7645bd327a3d6b53e Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 14:23:04 -0600 Subject: [PATCH 22/27] add flamesheet test case for manifold --- Exec/RegTests/FlameSheet/GNUmakefile | 15 ++++- Exec/RegTests/FlameSheet/input.manifold | 85 ++++++++++++++++++++++++ Exec/RegTests/FlameSheet/pelelmex_prob.H | 60 +++++++++++------ 3 files changed, 135 insertions(+), 25 deletions(-) create mode 100644 Exec/RegTests/FlameSheet/input.manifold diff --git a/Exec/RegTests/FlameSheet/GNUmakefile b/Exec/RegTests/FlameSheet/GNUmakefile index 25af38941..b19495cf3 100644 --- a/Exec/RegTests/FlameSheet/GNUmakefile +++ b/Exec/RegTests/FlameSheet/GNUmakefile @@ -29,9 +29,18 @@ THREAD_SANITIZER = FALSE USE_EFIELD = FALSE # PelePhysics -Chemistry_Model = drm19 -Eos_Model = Fuego -Transport_Model = Simple +USE_MANIFOLD = FALSE +ifeq ($(USE_MANIFOLD), TRUE) + Eos_Model = Manifold + Chemistry_Model = Null + Transport_Model = Manifold + Manifold_Dim = 1 + # Manifold_Type = Table +else + Chemistry_Model = drm19 + Eos_Model = Fuego + Transport_Model = Simple +endif PELE_HOME ?= ../../.. include $(PELE_HOME)/Exec/Make.PeleLMeX diff --git a/Exec/RegTests/FlameSheet/input.manifold b/Exec/RegTests/FlameSheet/input.manifold new file mode 100644 index 000000000..d67d9d7e9 --- /dev/null +++ b/Exec/RegTests/FlameSheet/input.manifold @@ -0,0 +1,85 @@ +#----------------------DOMAIN DEFINITION------------------------ +geometry.is_periodic = 0 0 # For each dir, 0: non-perio, 1: periodic +geometry.coord_sys = 0 # 0 => cart, 1 => RZ +geometry.prob_lo = 0.0 0.0 0.0 # x_lo y_lo (z_lo) +geometry.prob_hi = 0.003 0.012 0.016 # x_hi y_hi (z_hi) + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# Interior, Inflow, Outflow, Symmetry, +# SlipWallAdiab, NoSlipWallAdiab, SlipWallIsotherm, NoSlipWallIsotherm +peleLM.lo_bc = Interior Inflow +peleLM.hi_bc = Interior Outflow +peleLM.lo_bc = SlipWallAdiab Inflow +peleLM.hi_bc = SlipWallAdiab Outflow + + +#-------------------------AMR CONTROL---------------------------- +amr.n_cell = 64 256 32 # Level 0 number of cells in each direction +amr.v = 1 # AMR verbose +amr.max_level = 0 # maximum level number allowed +amr.ref_ratio = 2 2 2 2 # refinement ratio +amr.regrid_int = 5 # how often to regrid +amr.n_error_buf = 1 1 2 2 # number of buffer cells in error est +amr.grid_eff = 0.7 # what constitutes an efficient grid +amr.blocking_factor = 16 # block factor in grid generation (min box size) +amr.max_grid_size = 64 # max box size + + +#--------------------------- Problem ------------------------------- +prob.P_mean = 101325.0 +prob.standoff = -.03 +prob.pertmag = 0.0000 +pmf.datafile = "prem_drm19_phi1_p1_t298_mani.dat" + +#-------------------------PeleLM CONTROL---------------------------- +peleLM.v = 3 +peleLM.incompressible = 0 +peleLM.rho = 1.17 +peleLM.mu = 0.0 +peleLM.sdc_iterMax = 1 +peleLM.floor_species = 0 + +peleLM.do_temporals = 1 +peleLM.temporal_int = 2 +peleLM.mass_balance = 1 + +#amr.restart = chk00005 +#amr.check_int = 20 +amr.plot_int = 1 +amr.max_step = 10 +amr.dt_shrink = 0.01 +amr.stop_time = 0.01 +#amr.stop_time = 1.00 +amr.cfl = 0.5 +amr.derive_plot_vars = avg_pressure mag_vort mass_fractions maniout + +# ------------------- INPUTS DERIVED DIAGS ------------------ +peleLM.fuel_name = CH4 + +# --------------- INPUTS TO CHEMISTRY REACTOR --------------- +peleLM.chem_integrator = "ReactorRK64" + +mac_proj.verbose = 0 +nodal_proj.verbose = 0 + +#--------------------REFINEMENT CONTROL------------------------ +amr.refinement_indicators = temp +amr.temp.max_level = 3 +amr.temp.adjacent_difference_greater = 10 +amr.temp.field_name = temp + +#amrex.fpe_trap_invalid = 1 +#amrex.fpe_trap_zero = 1 +#amrex.fpe_trap_overflow = 1 + +#----------------------- Manifold Stuff ---------------------------- +manifold.model = Table +manifold.table.v = 1 +manifold.table.filename = prem_drm19_phi1_p1_t298.ctb +manifold.nominal_pressure_cgs = 1013250.0 +manifold.compute_temperature = 1 +manifold.density_lookup_type = log +peleLM.use_wbar = 0 +peleLM.chi_correction_type = DivuFirstIter +#peleLM.chi_correction_type = NoDivu +peleLM.print_chi_convergence = 1 \ No newline at end of file diff --git a/Exec/RegTests/FlameSheet/pelelmex_prob.H b/Exec/RegTests/FlameSheet/pelelmex_prob.H index c540c16d6..40537ef34 100644 --- a/Exec/RegTests/FlameSheet/pelelmex_prob.H +++ b/Exec/RegTests/FlameSheet/pelelmex_prob.H @@ -45,7 +45,6 @@ pelelmex_initdata( amrex::Real pert = 0.0; if (prob_parm.pertmag > 0.0) { - #if (AMREX_SPACEDIM == 2) amrex::GpuArray lpert{Lx}; if (prob_parm.pertlength > 0.0) { @@ -58,10 +57,8 @@ pelelmex_initdata( 1.017 * std::sin(2 * Pi * 5 * (x - 0.0033) / lpert[0]) + 0.982 * std::sin(2 * Pi * 5 * (x - 0.014234) / lpert[0])); } - amrex::Real y1 = (y - prob_parm.standoff - 0.5 * dx[1] + pert) * 100; amrex::Real y2 = (y - prob_parm.standoff + 0.5 * dx[1] + pert) * 100; - #elif (AMREX_SPACEDIM == 3) amrex::GpuArray lpert{Lx, Ly}; if (prob_parm.pertlength > 0.0) { @@ -80,29 +77,36 @@ pelelmex_initdata( std::sin(2 * Pi * 6 * (y - 0.018) / lpert[1]) + 0.982 * std::sin(2 * Pi * 5 * (x - 0.014234) / lpert[0])); } - amrex::Real y1 = (z - prob_parm.standoff - 0.5 * dx[2] + pert) * 100; amrex::Real y2 = (z - prob_parm.standoff + 0.5 * dx[2] + pert) * 100; #endif pele::physics::PMF::pmf(pmf_data, y1, y2, pmf_vals); - state(i, j, k, TEMP) = pmf_vals[0]; - ; + for (int dim = 0; dim < AMREX_SPACEDIM - 1; dim++) { + state(i, j, k, VELX + dim) = 0.0; + } + state(i, j, k, VELX + AMREX_SPACEDIM - 1) = pmf_vals[1] * 1e-2; +#ifdef USE_MANIFOLD_EOS + for (int n = 0; n < NUM_SPECIES - 1; n++) { + massfrac[n] = pmf_vals[3 + n]; + } + amrex::Real rho_temp, tmp1, tmp2; + eos.PYT2R(tmp1, massfrac, tmp2, rho_temp); + rho_temp *= 1.0e3; + state(i, j, k, DENSITY) = rho_temp; // CGS -> MKS conversion + state(i, j, k, RHOH) = 0.0; // No RhoH transport in manifold model + for (int n = 0; n < NUM_SPECIES - 1; n++) { + state(i, j, k, FIRSTSPEC + n) = massfrac[n] * rho_temp; + } + state(i, j, k, FIRSTSPEC + NUM_SPECIES - 1) = rho_temp; +#else for (int n = 0; n < NUM_SPECIES; n++) { molefrac[n] = pmf_vals[3 + n]; } eos.X2Y(molefrac, massfrac); - state(i, j, k, VELX) = 0.0; -#if (AMREX_SPACEDIM == 2) - state(i, j, k, VELY) = pmf_vals[1] * 1e-2; -#elif (AMREX_SPACEDIM == 3) - state(i, j, k, VELY) = 0.0; - state(i, j, k, VELZ) = pmf_vals[1] * 1e-2; -#endif - amrex::Real P_cgs = prob_parm.P_mean * 10.0; // Density @@ -119,6 +123,7 @@ pelelmex_initdata( for (int n = 0; n < NUM_SPECIES; n++) { state(i, j, k, FIRSTSPEC + n) = massfrac[n] * state(i, j, k, DENSITY); } +#endif } AMREX_GPU_DEVICE @@ -144,15 +149,25 @@ bcnormal( if (sgn == 1) { pele::physics::PMF::pmf(pmf_data, prob_lo[idir], prob_lo[idir], pmf_vals); - s_ext[VELX] = 0.0; -#if (AMREX_SPACEDIM == 2) - s_ext[VELY] = pmf_vals[1] * 1e-2; -#elif (AMREX_SPACEDIM == 3) - s_ext[VELY] = 0.0; - s_ext[VELZ] = pmf_vals[1] * 1e-2; -#endif - + for (int dim = 0; dim < AMREX_SPACEDIM - 1; dim++) { + s_ext[VELX + dim] = 0.0; + } + s_ext[VELX + AMREX_SPACEDIM - 1] = pmf_vals[1] * 1e-2; s_ext[TEMP] = pmf_vals[0]; + +#ifdef USE_MANIFOLD_EOS + for (int n = 0; n < NUM_SPECIES; n++) { + massfrac[n] = pmf_vals[3 + n]; + } + amrex::Real P_dummy, rho_cgs; + eos.PYT2R(P_dummy, massfrac, s_ext[TEMP], rho_cgs); + s_ext[DENSITY] = rho_cgs * 1.0e3; + s_ext[RHOH] = 0.0; + for (int n = 0; n < NUM_SPECIES - 1; n++) { + s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY]; + } + s_ext[FIRSTSPEC + NUM_SPECIES - 1] = s_ext[DENSITY]; +#else for (int n = 0; n < NUM_SPECIES; n++) { molefrac[n] = pmf_vals[3 + n]; } @@ -170,6 +185,7 @@ bcnormal( for (int n = 0; n < NUM_SPECIES; n++) { s_ext[FIRSTSPEC + n] = massfrac[n] * s_ext[DENSITY]; } +#endif } } From 75cbe6776010c4320cd43f2838f8204b19dbf2ff Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Thu, 24 Oct 2024 14:32:52 -0600 Subject: [PATCH 23/27] remove extraneous changes in file --- Source/PeleLMeX_Utils.H | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/PeleLMeX_Utils.H b/Source/PeleLMeX_Utils.H index 66fb0f18a..f38e182c2 100644 --- a/Source/PeleLMeX_Utils.H +++ b/Source/PeleLMeX_Utils.H @@ -1,6 +1,5 @@ #ifndef PELELM_UTILS_H #define PELELM_UTILS_H - template amrex::Gpu::DeviceVector convertToDeviceVector(amrex::Vector v) @@ -14,5 +13,4 @@ convertToDeviceVector(amrex::Vector v) #endif return v_d; } - #endif From ad8e37b0029ee3d99ef18c8a55ccfddec704eaaa Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 4 Nov 2024 15:46:45 -0700 Subject: [PATCH 24/27] Update Source/PeleLMeX_Setup.cpp Co-authored-by: Marc T. Henry de Frahan --- Source/PeleLMeX_Setup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 0334d4121..31c5c2184 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -1133,7 +1133,7 @@ PeleLM::derivedSetup() #ifdef USE_MANIFOLD_EOS auto& mani_data = eos_parms.host_only_parm().manfunc_par->host_parm(); - int nmanivar = mani_data.Nvar; + const int nmanivar = mani_data.Nvar; Vector var_names_maniout(nmanivar); for (int n = 0; n < nmanivar; n++) { std::string nametmp = std::string( From f5e7637f9927c5922f66339a0a6688f05a516e32 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 4 Nov 2024 15:51:24 -0700 Subject: [PATCH 25/27] address Marc HdF's comments with minor formatting --- Source/PeleLMeX_K.H | 3 +-- Source/PeleLMeX_Setup.cpp | 4 ++-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Source/PeleLMeX_K.H b/Source/PeleLMeX_K.H index 826bac6af..c8c2e6e42 100644 --- a/Source/PeleLMeX_K.H +++ b/Source/PeleLMeX_K.H @@ -735,8 +735,7 @@ makeVelForce( const amrex::Real& rho_incomp, int pseudo_gravity, int pseudo_gravity_dir, - const amrex::Real& /*time*/ - , + const amrex::Real& /*time*/, amrex::GpuArray const gravity, amrex::GpuArray const gp0, const amrex::Real& dV_control, diff --git a/Source/PeleLMeX_Setup.cpp b/Source/PeleLMeX_Setup.cpp index 31c5c2184..dfcc24c7c 100644 --- a/Source/PeleLMeX_Setup.cpp +++ b/Source/PeleLMeX_Setup.cpp @@ -468,8 +468,8 @@ PeleLM::readParameters() << std::endl; } if ( - m_use_wbar != 0 && - pele::physics::PhysicsType::eos_type::identifier() == "Manifold") { + (m_use_wbar != 0) && + (pele::physics::PhysicsType::eos_type::identifier() == "Manifold")) { amrex::Abort("Use of Wbar fluxes is not compatible with Manifold EOS"); } From e3ec5841e5bcea77e51db1570151b9b5b01d2fe5 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Mon, 4 Nov 2024 15:58:23 -0700 Subject: [PATCH 26/27] remove copied line --- Exec/Make.PeleLMeX | 1 - 1 file changed, 1 deletion(-) diff --git a/Exec/Make.PeleLMeX b/Exec/Make.PeleLMeX index 5760d9681..e31ee739f 100644 --- a/Exec/Make.PeleLMeX +++ b/Exec/Make.PeleLMeX @@ -78,7 +78,6 @@ ifeq ($(Eos_Model),$(filter $(Eos_Model),Soave-Redlich-Kwong)) DEFINES += -DUSE_SRK_EOS endif ifeq ($(Eos_Model),$(filter $(Eos_Model),Manifold)) - DEFINES += -DUSE_MANIFOLD_EOS DEFINES += -DUSE_MANIFOLD_EOS DEFINES += -DMANIFOLD_DIM=$(Manifold_Dim) ifeq ($(Manifold_Type),Table) From 136c4b8ebd0c3c894943ca5af2b0f1cf8faa266b Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Tue, 5 Nov 2024 15:36:27 -0700 Subject: [PATCH 27/27] update PP for merged RY2R stuff --- Submodules/PelePhysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Submodules/PelePhysics b/Submodules/PelePhysics index 683ca7998..067456088 160000 --- a/Submodules/PelePhysics +++ b/Submodules/PelePhysics @@ -1 +1 @@ -Subproject commit 683ca7998d121ca8cc31c294ece942e55424f2f8 +Subproject commit 06745608874aea5230088aadeb7571a5e9afb8f1