Support Manifold EoS in RK64 Reactor

.github/workflows/ci.yml

@@ -271,9 +271,15 @@ jobs:
               pip install numpy
               ccache -z
               make -j ${{env.NPROCS}} Eos_Model=Fuego Chemistry_Model=dodecane_lu TINY_PROFILE=TRUE USE_CCACHE=TRUE ${{matrix.amrex_build_args}}
-              bash exec_ignDelay.sh
+              bash exec_ignDelay.sh firstpass
               python check_ignDelay.py
               if [ $? -ne 0 ]; then exit 1; fi; \
+              rm log PPreaction.txt inputs/inputs.0d_refine
+              make realclean
+              make -j ${{env.NPROCS}} Eos_Model=Manifold Transport_Model=Manifold Manifold_Dim=1 Chemistry_Model=Null TINY_PROFILE=TRUE USE_CCACHE=TRUE ${{matrix.amrex_build_args}}
+              bash exec_ignDelay.sh manifold
+              python check_ignDelay.py 0.0766 10
+              if [ $? -ne 0 ]; then exit 1; fi; \
           fi;
           make realclean
       - name: Ignition delay ccache report

Mechanisms/Null/mechanism.cpp

@@ -14,14 +14,14 @@ void
 CKSYME_STR(amrex::Vector<std::string>& ename)
 {
   ename.resize(1);
-  ename[0] = "XO";
+  ename[0] = "X0";
 }
 void
 CKSYMS_STR(amrex::Vector<std::string>& kname)
 {
   kname.resize(NUM_SPECIES);
   for (int i = 0; i < NUM_SPECIES; ++i) {
-    kname[0] = "X" + std::to_string(i);
+    kname[i] = "X" + std::to_string(i);
   }
 }
 

README.md

@@ -15,7 +15,7 @@ PelePhysics contains C++ source code for the following physics and other modules
 * *Reactions*: Routines for substepped integration of stiff chemical systems, including with CVODE
 * *Spray*: Lagrangian spray droplet library, formerly part of [PeleMP](https://github.com/AMReX-Combustion/PeleMP)
 * *Soot*: An implementation of the Hybrid Method of Moments soot model, formerly part of [PeleMP](https://github.com/AMReX-Combustion/PeleMP)
-* *Radiation*: Raidative heat transfer model based on the spherical harmonics method, formerly [PeleRad](https://github.com/AMReX-Combustion/PeleRad)
+* *Radiation*: Radiative heat transfer model based on the spherical harmonics method, formerly [PeleRad](https://github.com/AMReX-Combustion/PeleRad)
 * *Utility*: Several handy modules for data management across other Pele codes
 
 Additionally, PelePhysics contains a variety of stand-alone tools to aid in Pele simulation workflows, found in the `Support` directory.

Source/Eos/Manifold.H

@@ -92,7 +92,7 @@ struct Manifold
     const amrex::Real /*R*/,
     const amrex::Real /*T*/,
     const amrex::Real* /*Y[]*/,
-    amrex::Real* /*E*/)
+    amrex::Real /*E*/)
   {
     amrex::Error("RTY2E does not have significance for Manifold EOS");
   }

Source/Reactions/ReactorBase.H

@@ -81,11 +81,28 @@ public:
 
   void set_typ_vals_ode(const std::vector<amrex::Real>& ExtTypVals);
 
+  // Manifold EOS needs an eosparm - right now that is only propagated through
+  // the RK64 reactor
+  virtual void set_eos_parm(
+    const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+      eosparm)
+  {
+    auto eos = pele::physics::PhysicsType::eos(eosparm);
+    if (eos.identifier() == "Manifold") {
+      amrex::Abort(
+        "Manifold EOS only supported with ReactorRK64 and ReactorNull for now");
+    } else {
+      m_eosparm = eosparm;
+    }
+  }
+
   ~ReactorBase() override = default;
 
 protected:
   int verbose{0};
   amrex::GpuArray<amrex::Real, NUM_SPECIES + 1> m_typ_vals = {0.0};
+  const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+    m_eosparm;
 };
 } // namespace pele::physics::reactions
 #endif

Source/Reactions/ReactorNull.H

@@ -86,6 +86,13 @@ public:
       rhoE, frcEExt, FC_in, y_vect, vect_energy, FCunt, dt);
   }
 
+  void set_eos_parm(
+    const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+      eosparm) override
+  {
+    m_eosparm = eosparm;
+  }
+
 private:
   utils::FlattenOps<Ordering> flatten_ops;
   int m_reactor_type{0};

Source/Reactions/ReactorNull.cpp

@@ -40,6 +40,7 @@ ReactorNull::react(
 #endif
 
   int captured_reactor_type = m_reactor_type;
+  const auto* leosparm = m_eosparm;
 
   ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
     amrex::Real renergy_loc =
@@ -60,7 +61,7 @@ ReactorNull::react(
     }
     amrex::Real energy_loc = renergy_loc / rho_loc;
     amrex::Real T_loc = T_in(i, j, k, 0);
-    auto eos = pele::physics::PhysicsType::eos();
+    auto eos = pele::physics::PhysicsType::eos(leosparm);
     if (captured_reactor_type == ReactorTypes::e_reactor_type) {
       eos.REY2T(rho_loc, energy_loc, Y_loc, T_loc);
     } else if (captured_reactor_type == ReactorTypes::h_reactor_type) {

Source/Reactions/ReactorRK64.H

@@ -120,6 +120,13 @@ public:
       rhoE, frcEExt, FC_in, y_vect, vect_energy, FCunt, dt);
   }
 
+  void set_eos_parm(
+    const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+      eosparm) override
+  {
+    m_eosparm = eosparm;
+  }
+
 private:
   amrex::Real absTol{1e-10};
   int rk64_nsubsteps_guess{10};

Source/Reactions/ReactorRK64.cpp

@@ -64,6 +64,7 @@ ReactorRK64::react(
   const int captured_nsubsteps_min = rk64_nsubsteps_min;
   const int captured_nsubsteps_max = rk64_nsubsteps_max;
   const amrex::Real captured_abstol = absTol;
+  const auto* leosparm = m_eosparm;
   RK64Params rkp;
 
   amrex::Gpu::DeviceVector<int> v_nsteps(ncells, 0);
@@ -103,7 +104,7 @@ ReactorRK64::react(
       for (int stage = 0; stage < rkp.nstages_rk64; stage++) {
         utils::fKernelSpec<Ordering>(
           0, 1, current_time - time_init, captured_reactor_type, soln_reg, ydot,
-          rhoe_init, rhoesrc_ext, rYsrc_ext);
+          rhoe_init, rhoesrc_ext, rYsrc_ext, leosparm);
 
         for (int sp = 0; sp < neq; sp++) {
           error_reg[sp] += rkp.err_rk64[stage] * dt_rk * ydot[sp];
@@ -131,9 +132,9 @@ ReactorRK64::react(
           rkp.betaerr_rk64 * pow((captured_abstol / max_err), rkp.exp2_rk64);
         dt_rk = amrex::max<amrex::Real>(dt_rk_min, dt_rk * change_factor);
       }
+      // Don't overstep the integration time
+      dt_rk = amrex::min<amrex::Real>(dt_rk, time_out - current_time);
     }
-    // Don't overstep the integration time
-    dt_rk = amrex::min<amrex::Real>(dt_rk, time_out - current_time);
     d_nsteps[icell] = nsteps;
     // copy data back
     for (int sp = 0; sp < neq; sp++) {
@@ -192,6 +193,7 @@ ReactorRK64::react(
   const int captured_nsubsteps_min = rk64_nsubsteps_min;
   const int captured_nsubsteps_max = rk64_nsubsteps_max;
   const amrex::Real captured_abstol = absTol;
+  const auto* leosparm = m_eosparm;
   RK64Params rkp;
 
   int ncells = static_cast<int>(box.numPts());
@@ -210,21 +212,18 @@ ReactorRK64::react(
     amrex::Real current_time = time_init;
     const int neq = (NUM_SPECIES + 1);
 
-    amrex::Real rho = 0.0;
+    auto eos = pele::physics::PhysicsType::eos(leosparm);
     for (int sp = 0; sp < NUM_SPECIES; sp++) {
       soln_reg[sp] = rY_in(i, j, k, sp);
       carryover_reg[sp] = soln_reg[sp];
-      rho += rY_in(i, j, k, sp);
     }
-    amrex::Real rho_inv = 1.0 / rho;
+    amrex::Real rho = 0.0, rho_inv = 0.0;
     amrex::Real mass_frac[NUM_SPECIES] = {0.0};
-    for (int sp = 0; sp < NUM_SPECIES; sp++) {
-      mass_frac[sp] = rY_in(i, j, k, sp) * rho_inv;
-    }
+    eos.RY2RRinvY(soln_reg, rho, rho_inv, mass_frac);
+
     amrex::Real temp = T_in(i, j, k, 0);
 
     amrex::Real Enrg_loc = rEner_in(i, j, k, 0) * rho_inv;
-    auto eos = pele::physics::PhysicsType::eos();
     if (captured_reactor_type == ReactorTypes::e_reactor_type) {
       eos.REY2T(rho, Enrg_loc, mass_frac, temp);
     } else if (captured_reactor_type == ReactorTypes::h_reactor_type) {
@@ -255,7 +254,7 @@ ReactorRK64::react(
       for (int stage = 0; stage < rkp.nstages_rk64; stage++) {
         utils::fKernelSpec<Ordering>(
           0, 1, current_time - time_init, captured_reactor_type, soln_reg, ydot,
-          rhoe_init, rhoesrc_ext, rYsrc_ext);
+          rhoe_init, rhoesrc_ext, rYsrc_ext, leosparm);
 
         for (int sp = 0; sp < neq; sp++) {
           error_reg[sp] += rkp.err_rk64[stage] * dt_rk * ydot[sp];
@@ -290,15 +289,11 @@ ReactorRK64::react(
     // copy data back
     int icell = (k - lo.z) * len.x * len.y + (j - lo.y) * len.x + (i - lo.x);
     d_nsteps[icell] = nsteps;
-    rho = 0.0;
     for (int sp = 0; sp < NUM_SPECIES; sp++) {
       rY_in(i, j, k, sp) = soln_reg[sp];
-      rho += rY_in(i, j, k, sp);
-    }
-    rho_inv = 1.0 / rho;
-    for (int sp = 0; sp < NUM_SPECIES; sp++) {
-      mass_frac[sp] = rY_in(i, j, k, sp) * rho_inv;
     }
+    eos.RY2RRinvY(soln_reg, rho, rho_inv, mass_frac);
+
     temp = soln_reg[NUM_SPECIES];
     rEner_in(i, j, k, 0) = rhoe_init[0] + dt_react * rhoesrc_ext[0];
     Enrg_loc = rEner_in(i, j, k, 0) * rho_inv;

Source/Reactions/ReactorUtils.H

@@ -369,29 +369,27 @@ fKernelSpec(
   amrex::Real* ydot_d, // NOLINT(readability-non-const-parameter)
   const amrex::Real* rhoe_init,
   const amrex::Real* rhoesrc_ext,
-  const amrex::Real* rYs)
+  const amrex::Real* rYs,
+  const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+    eosparm = nullptr)
 {
-  amrex::Real rho_pt = 0.0;
-  amrex::GpuArray<amrex::Real, NUM_SPECIES> massfrac = {0.0};
+  auto eos = pele::physics::PhysicsType::eos(eosparm);
+  amrex::Real rho_pt = 0.0, rho_pt_inv = 0.0;
+  amrex::GpuArray<amrex::Real, NUM_SPECIES> massfrac = {0.0}, massdens = {0.0};
   for (int n = 0; n < NUM_SPECIES; n++) {
-    massfrac[n] = yvec_d[vec_index<OrderType>(n, icell, ncells)];
-    rho_pt += massfrac[n];
-  }
-  const amrex::Real rho_pt_inv = 1.0 / rho_pt;
-
-  for (int n = 0; n < NUM_SPECIES; n++) {
-    massfrac[n] *= rho_pt_inv;
+    massdens[n] = yvec_d[vec_index<OrderType>(n, icell, ncells)];
   }
+  eos.RY2RRinvY(massdens.data(), rho_pt, rho_pt_inv, massfrac.data());
 
   const amrex::Real nrg_pt =
     (rhoe_init[icell] + rhoesrc_ext[icell] * dt_save) * rho_pt_inv;
 
   amrex::Real temp_pt =
     yvec_d[vec_index<OrderType>(NUM_SPECIES, icell, ncells)];
 
-  amrex::Real Cv_pt = 0.0;
+  amrex::Real Cv_pt = 1.0;
   amrex::GpuArray<amrex::Real, NUM_SPECIES> ei_pt = {0.0};
-  auto eos = pele::physics::PhysicsType::eos();
+#ifndef USE_MANIFOLD_EOS
   if (reactor_type == ReactorTypes::e_reactor_type) {
     eos.REY2T(rho_pt, nrg_pt, massfrac.arr, temp_pt);
     eos.RTY2Ei(rho_pt, temp_pt, massfrac.arr, ei_pt.arr);
@@ -403,6 +401,7 @@ fKernelSpec(
   } else {
     amrex::Abort("Wrong reactor type. Choose between 1 (e) or 2 (h).");
   }
+#endif
 
   amrex::GpuArray<amrex::Real, NUM_SPECIES> cdots_pt = {0.0};
   eos.RTY2WDOT(rho_pt, temp_pt, massfrac.arr, cdots_pt.arr);

Testing/Exec/EosEval/main.cpp

@@ -190,6 +190,8 @@ main(int argc, char* argv[])
       amrex::WriteSingleLevelPlotfile(
         outfile, VarPlt, plt_VarsName, geom, 0.0, 0);
     }
+
+    eos_parms.deallocate();
   }
 
   amrex::Finalize();

Testing/Exec/IgnitionDelay/GPU_misc.H

@@ -23,7 +23,9 @@ initialize_data(
   amrex::Real t0,
   amrex::Real equiv_ratio,
   amrex::Real p,
-  amrex::GeometryData const& geomdata) noexcept
+  amrex::GeometryData const& geomdata,
+  const pele::physics::eos::EosParm<pele::physics::PhysicsType::eos_type>*
+    eosparm) noexcept
 {
   const amrex::Real* plo = geomdata.ProbLo();
   const amrex::Real* phi = geomdata.ProbHi();
@@ -50,26 +52,30 @@ initialize_data(
     P[n] = L[n] / 4.0;
   }
   // Y
+  auto eos = pele::physics::PhysicsType::eos(eosparm);
+
+#ifndef USE_MANIFOLD_EOS
   for (int n = 0; n < NUM_SPECIES; n++) {
     X[n] = 0.0;
   }
-  // X[O2_ID]   = 0.7;
-  // X[fuel_id] = 0.3;
+  // Assumes n-Dodecane
   X[O2_ID] = 1.0 / (equiv_ratio / 18.5 + 1 + 0.79 / 0.21);
   X[N2_ID] = 0.79 * X[O2_ID] / 0.21;
   X[fuel_id] = 1.0 - X[O2_ID] - X[N2_ID];
-
-  auto eos = pele::physics::PhysicsType::eos();
-
   eos.X2Y(&X[0], &Y[0]);
+#else
+  Y[0] = 0.0; // Progress variable
+  Y[1] = 1.0; // Density
+#endif
 
   // T
   temp = t0;
   // temp =  Temp_lo + (Temp_hi-Temp_lo)*y/L[1] + dTemp *
   // std::sin(2.0*pi*y/P[1]);
   if (energy_type == 1) {
     // get rho and E
-    eos.PYT2RE(pressure, &Y[0], temp, density, energy);
+    eos.PYT2R(pressure, &Y[0], temp, density);
+    eos.RTY2E(density, temp, &Y[0], energy);
   } else {
     // get rho and H
     eos.PYT2R(pressure, &Y[0], temp, density);

Testing/Exec/IgnitionDelay/README.md

@@ -19,4 +19,4 @@ Ignition delay will call the Pele executable find in this directory. If you have
 
 3. Adjust your parameters in `inputs/inputs.0d_firstpass`
 
-4. Execute `bash exec_ignDelay.sh`. The result is outputted in `log`. 
+4. Execute `bash exec_ignDelay.sh firstpass`. The result is outputted in `log`.

Testing/Exec/IgnitionDelay/check_ignDelay.py

@@ -14,11 +14,19 @@ def readResult(filename):
 
 
 if __name__=="__main__":
+    # Two arguments: intended correct value and uncertainty factor
+    if len(sys.argv) > 1:
+        correct_value = float(sys.argv[1])
+        uncertainty_factor = float(sys.argv[2])
+    else:
+        correct_value = 0.0763
+        uncertainty_factor = 10.0
+
     # Read ignition delay output
     tIgn,uncertainty = readResult('log')
 
     # Compare to theoretical value
-    if abs(tIgn-0.0763)>10*uncertainty:
+    if abs(tIgn-correct_value)>uncertainty_factor*uncertainty:
          sys.exit(1)
     else:
          sys.exit(0)