From d55daa0b858bb46f2d804b3d0bded48bf596d1a9 Mon Sep 17 00:00:00 2001 From: Bruce Perry Date: Wed, 18 Sep 2024 10:40:57 -0600 Subject: [PATCH] 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 7d74421f..9f60602e 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 2ff984df..76a1a3c1 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 628e5d70..83c34e7a 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 eece3cda..a33aa35a 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 534d3e42..578fbef5 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 e4d36c9c..23a77f96 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 169782a0..e4521096 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; }