Add mixture fraction and progress variable derived variables (#76)

* Pull mixture fraction into LM -> LMeX
Add integral of fuel consumption and HR to temporals.

Conflicts:
	Source/PeleLM.H
	Source/PeleLMSetup.cpp
	Source/PeleLMUtils.cpp

* Add progress variable pieces to .H

* Derived now takes in the pelelm object to access to members.
Add progress variable kernel.

* Declare and init mixfrac and prog var during Setup

* Progress variable init functions.

* FlameSheet input.2d-regt updated to show usage of mixture fraction and progress
variable.

* Add doc for mixture fraction and progress variable.

Docs/source/LMeXControls.rst

@@ -147,6 +147,30 @@ Linear solvers are a key component of PeleLMeX algorithm, separate controls are
 
 ### Run-time diagnostics
 
+Combustion diagnostics often involve the use of a mixture fraction and/or a progress variable, both of which can be defined
+at run time and added to the derived variables included in the plotfile. If `mixture_fraction` or `progress_variable` is
+added to the `amr.derive_plot_vars` list, one need to provide input for defining those. The mixture fraction is based on
+Bilger's element definition and one needs to provide the composition of the 'fuel' and 'oxidizer' tanks using a Cantera-like
+format (<species>:<value>) which assumes unspecified species at zero, or a list of floats, in which case all the species must
+be specified in the order they appear in the mechanism file.
+The progress variable definition in based on a linear combination of the species mass fractions and temperature, and can be
+specified in a manner similar to the mixture fraction, providing a list of weights and the prescription of a 'cold' and 'hot'
+state:
+
+::
+
+    # ------------------- INPUTS DERIVED DIAGS ------------------
+    peleLM.fuel_name = CH4 
+    peleLM.mixtureFraction.format = Cantera
+    peleLM.mixtureFraction.type   = mass
+    peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767
+    peleLM.mixtureFraction.fuelTank = H2:0.5 CH4:0.5
+    peleLM.progressVariable.format = Cantera
+    peleLM.progressVariable.weights = CO:1.0 CO2:1.0
+    peleLM.progressVariable.coldState = CO:0.0 CO2:0.0
+    peleLM.progressVariable.hotState = CO:0.000002 CO2:0.0666
+
+
 A set of diagnostics available at runtime are currently under development. The following provide an example for extracting
 the state variables on a 'x','y' or 'z' aligned plane and writting a 2D plotfile compatible with Amrvis, Paraview or yt:
 

Exec/RegTests/FlameSheet/input.2d-regt

@@ -50,7 +50,18 @@ amr.dt_shrink = 0.0001
 amr.stop_time = 0.001
 #amr.stop_time = 1.00
 amr.cfl = 0.5
-amr.derive_plot_vars = avg_pressure mag_vort mass_fractions
+amr.derive_plot_vars = avg_pressure mag_vort mass_fractions mixture_fraction progress_variable
+
+# ------------------- INPUTS DERIVED DIAGS ------------------
+peleLM.fuel_name = CH4
+peleLM.mixtureFraction.format = Cantera
+peleLM.mixtureFraction.type   = mass
+peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767
+peleLM.mixtureFraction.fuelTank = H2:0.5 CH4:0.5
+peleLM.progressVariable.format = Cantera
+peleLM.progressVariable.weights = CO:1.0 CO2:1.0
+peleLM.progressVariable.coldState = CO:0.0 CO2:0.0
+peleLM.progressVariable.hotState = CO:0.000002 CO2:0.0666
 
 # --------------- INPUTS TO CHEMISTRY REACTOR ---------------
 peleLM.chem_integrator = "ReactorCvode"

Source/PeleLM.H

@@ -723,6 +723,16 @@ class PeleLM : public amrex::AmrCore {
    void updateTypicalValuesChem();
 
    void checkMemory(const std::string &a_message);
+
+   // Mixture fraction & Progress variable
+   void initMixtureFraction();
+   void parseComposition(amrex::Vector<std::string> compositionIn,
+                         std::string                compositionType,
+                         amrex::Real               *massFrac);
+   void parseVars(const amrex::Vector<std::string> &a_varNames,
+                  const amrex::Vector<std::string> &a_stringIn,
+                  amrex::Vector<amrex::Real> &a_rVars);
+   void initProgressVariable();
    //-----------------------------------------------------------------------------
 
 #ifdef AMREX_USE_EB
@@ -945,6 +955,7 @@ class PeleLM : public amrex::AmrCore {
    amrex::Vector<amrex::MultiFab* > getDivUVect(const PeleLM::TimeStamp &a_time);
    amrex::Vector<amrex::MultiFab* > getDiffusivityVect(const PeleLM::TimeStamp &a_time);
    amrex::Vector<amrex::MultiFab* > getViscosityVect(const PeleLM::TimeStamp &a_time);
+   amrex::Vector<amrex::MultiFab* > getIRVect();
 #ifdef PELE_USE_EFIELD
    amrex::Vector<amrex::MultiFab* > getPhiVVect(const PeleLM::TimeStamp &a_time);
    amrex::Vector<amrex::MultiFab* > getnEVect(const PeleLM::TimeStamp &a_time);
@@ -1216,6 +1227,24 @@ class PeleLM : public amrex::AmrCore {
 #endif
    //-----------------------------------------------------------------------------
 
+   //-----------------------------------------------------------------------------
+   // Combustion physics
+   int fuelID = -1;
+
+   // Mixture fraction
+   amrex::Array<amrex::Real, 4> Beta_mix;
+   amrex::Array<amrex::Real, NUM_SPECIES> spec_Bilger_fact;
+   amrex::Real Zfu = -1.0;
+   amrex::Real Zox = -1.0;
+
+   // Progress variable
+   // C is a linear combination of Ys and T
+   amrex::GpuArray<amrex::Real,NUM_SPECIES+1> m_Cweights;
+   amrex::Real m_C0 = -1.0;
+   amrex::Real m_C1 = -1.0;
+   int m_Crevert = 0;           // Flip the definition with 1 - ∑_k wgt*v_k
+   //-----------------------------------------------------------------------------
+
    //-----------------------------------------------------------------------------
    // Linear Solvers
 

Source/PeleLM.cpp

@@ -239,6 +239,16 @@ PeleLM::getViscosityVect(const TimeStamp &a_time) {
    return r;
 }
 
+Vector<MultiFab *>
+PeleLM::getIRVect() {
+   Vector<MultiFab*> r;
+   r.reserve(finest_level+1);
+   for (int lev = 0; lev <= finest_level; ++lev) {
+      r.push_back(&(m_leveldatareact[lev]->I_R));
+   }
+   return r;
+}
+
 void
 PeleLM::averageDownState(const PeleLM::TimeStamp &a_time)
 {

Source/PeleLMDerive.H

@@ -8,8 +8,10 @@
 #include <AMReX_MultiFabUtil.H>
 
 // Mimic AMReX Derive behavior but without StateDescriptor
+// Forward declaration of PeleLM
+class PeleLM;
 
-typedef void (*PeleLMDeriveFunc) (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+typedef void (*PeleLMDeriveFunc) (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                                   const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                                   const amrex::Geometry& geomdata,
                                   amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);

Source/PeleLMDeriveFunc.H

@@ -8,28 +8,41 @@
 #include <PeleLMEFDeriveFunc.H>
 #endif
 
-void pelelm_dermassfrac (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+// Forward declaration of PeleLM
+class PeleLM;
+
+void pelelm_dermassfrac (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                          const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                          const amrex::Geometry& geomdata,
                          amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
 
-void pelelm_dertemp (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_dertemp (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                      const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                      const amrex::Geometry& geomdata,
                      amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
 
-void pelelm_derkineticenergy (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_derkineticenergy (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                               const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                               const amrex::Geometry& geomdata,
                               amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
 
-void pelelm_deravgpress (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_deravgpress (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                          const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                          const amrex::Geometry& geomdata,
                          amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
 
-void pelelm_dermgvort (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_dermgvort (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
                        const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
                        const amrex::Geometry& geomdata,
                        amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
+
+void pelelm_dermixfrac (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+                        const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
+                        const amrex::Geometry& geomdata,
+                        amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
+
+void pelelm_derprogvar (PeleLM* a_pelelm, const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int ncomp,
+                        const amrex::FArrayBox& statefab, const amrex::FArrayBox& pressfab,
+                        const amrex::Geometry& geomdata,
+                        amrex::Real time, const amrex::Vector<amrex::BCRec> &bcrec, int level);
 #endif

Source/PeleLMDeriveFunc.cpp

@@ -2,13 +2,14 @@
 #include <PeleLM_Index.H>
 #include <PelePhysics.H>
 #include <mechanism.H>
+#include <PeleLM.H>
 
 using namespace amrex;
 
 //
 // Extract temp
 //
-void pelelm_dertemp (const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_dertemp (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
                      const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
                      const Geometry& /*geomdata*/,
                      Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
@@ -29,7 +30,7 @@ void pelelm_dertemp (const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
 //
 // Extract species mass fractions Y_n
 //
-void pelelm_dermassfrac (const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
+void pelelm_dermassfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
                          const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
                          const Geometry& /*geomdata*/,
                          Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
@@ -53,7 +54,7 @@ void pelelm_dermassfrac (const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
 //
 // Compute cell averaged pressure from nodes
 //
-void pelelm_deravgpress (const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
+void pelelm_deravgpress (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
                          const FArrayBox& /*statefab*/, const FArrayBox& pressfab,
                          const Geometry& /*geomdata*/,
                          Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
@@ -81,7 +82,7 @@ void pelelm_deravgpress (const Box& bx, FArrayBox& derfab, int dcomp, int /*ncom
 //
 // Compute vorticity magnitude
 //
-void pelelm_dermgvort (const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
+void pelelm_dermgvort (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
                        const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
                        const Geometry& geomdata,
                        Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
@@ -126,7 +127,7 @@ void pelelm_dermgvort (const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*
 //
 // Compute the kinetic energy
 //
-void pelelm_derkineticenergy (const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
+void pelelm_derkineticenergy (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
                               const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
                               const Geometry& /*geomdata*/,
                               Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
@@ -146,3 +147,88 @@ void pelelm_derkineticenergy (const Box& bx, FArrayBox& derfab, int dcomp, int /
                                           + vel(i,j,k,2)*vel(i,j,k,2)) );
     });
 }
+
+//
+// Compute mixture fraction
+//
+void pelelm_dermixfrac (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
+                        const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
+                        const Geometry& /*geomdata*/,
+                        Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
+
+{
+    AMREX_ASSERT(derfab.box().contains(bx));
+    AMREX_ASSERT(statefab.box().contains(bx));
+    AMREX_ASSERT(ncomp == 1);
+
+    if (a_pelelm->Zfu < 0.0) amrex::Abort("Mixture fraction not initialized");
+
+    auto const density   = statefab.array(DENSITY);
+    auto const rhoY      = statefab.array(FIRSTSPEC);
+    auto       mixt_frac = derfab.array(dcomp);
+
+    amrex::Real Zox_lcl   = a_pelelm->Zox;
+    amrex::Real Zfu_lcl   = a_pelelm->Zfu;
+    amrex::Real denom_inv = 1.0 / (Zfu_lcl - Zox_lcl);
+    amrex::GpuArray<amrex::Real,NUM_SPECIES> fact_Bilger;
+    for (int n=0; n<NUM_SPECIES; ++n) {
+        fact_Bilger[n] = a_pelelm->spec_Bilger_fact[n];
+    }
+
+    amrex::ParallelFor(bx,
+    [density, rhoY, mixt_frac, fact_Bilger, Zox_lcl, denom_inv] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+        amrex::Real rho_inv = 1.0_rt / density(i,j,k);
+        mixt_frac(i,j,k) = 0.0_rt;
+        for (int n = 0; n<NUM_SPECIES; ++n) {
+            mixt_frac(i,j,k) += ( rhoY(i,j,k,n) * fact_Bilger[n] ) * rho_inv;
+        }
+        mixt_frac(i,j,k) = ( mixt_frac(i,j,k) - Zox_lcl ) * denom_inv;
+    });
+}
+
+//
+// Compute progress variable
+//
+void pelelm_derprogvar (PeleLM* a_pelelm, const Box& bx, FArrayBox& derfab, int dcomp, int ncomp,
+                        const FArrayBox& statefab, const FArrayBox& /*pressfab*/,
+                        const Geometry& /*geomdata*/,
+                        Real /*time*/, const Vector<BCRec>& /*bcrec*/, int /*level*/)
+
+{
+    AMREX_ASSERT(derfab.box().contains(bx));
+    AMREX_ASSERT(statefab.box().contains(bx));
+    AMREX_ASSERT(ncomp == 1); 
+
+    if (a_pelelm->m_C0 < 0.0) amrex::Abort("Progress variable not initialized");
+
+    auto const density  = statefab.array(DENSITY);
+    auto const rhoY     = statefab.array(FIRSTSPEC);
+    auto const temp     = statefab.array(TEMP);
+    auto       prog_var = derfab.array(dcomp);
+
+    amrex::Real C0_lcl   = a_pelelm->m_C0;
+    amrex::Real C1_lcl   = a_pelelm->m_C1;
+    amrex::Real denom_inv = 1.0 / (C1_lcl - C0_lcl);
+    amrex::GpuArray<amrex::Real,NUM_SPECIES+1> Cweights;
+    for (int n=0; n<NUM_SPECIES+1; ++n) {
+        Cweights[n] = a_pelelm->m_Cweights[n];
+    }   
+
+    amrex::ParallelFor(bx, [=,revert=a_pelelm->m_Crevert]
+    AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {   
+        amrex::Real rho_inv = 1.0_rt / density(i,j,k);
+        prog_var(i,j,k) = 0.0_rt;
+        for (int n = 0; n<NUM_SPECIES; ++n) {
+            prog_var(i,j,k) += ( rhoY(i,j,k,n) * Cweights[n] ) * rho_inv;
+        }   
+        prog_var(i,j,k) += temp(i,j,k) * Cweights[NUM_SPECIES];
+        if (revert) {
+           prog_var(i,j,k) = 1.0 - ( prog_var(i,j,k) - C0_lcl ) * denom_inv;
+        } else {
+           prog_var(i,j,k) = ( prog_var(i,j,k) - C0_lcl ) * denom_inv;
+        }   
+    });
+}
+

Source/PeleLMSetup.cpp

@@ -105,6 +105,10 @@ void PeleLM::Setup() {
 #endif
    }
 
+   // Mixture fraction & Progress variable
+   initMixtureFraction();
+   initProgressVariable();
+
    // Initiliaze turbulence injection
    turb_inflow.init(Geom(0));
 
@@ -606,6 +610,20 @@ void PeleLM::variablesSetup() {
    } else {
       typical_values.resize(NVAR,-1.0);
    }
+
+   if ( !m_incompressible ) {
+      // -----------------------------------------
+      // Combustion
+      // -----------------------------------------
+      ParmParse pp("peleLM");
+      std::string fuel_name = "";
+      pp.query("fuel_name",fuel_name);
+      fuel_name = "rho.Y("+fuel_name+")";
+      if (isStateVariable(fuel_name)) {
+         fuelID = stateVariableIndex(fuel_name) - FIRSTSPEC;
+      }
+   }
+
    if (max_level > 0 && !m_initial_grid_file.empty()) {
      readGridFile(m_initial_grid_file, m_initial_ba);
      if (verbose > 0) {
@@ -685,6 +703,12 @@ void PeleLM::derivedSetup()
    // Kinetic energy
    derive_lst.add("kinetic_energy",IndexType::TheCellType(),1,pelelm_derkineticenergy,the_same_box);
 
+   // Mixture fraction
+   derive_lst.add("mixture_fraction",IndexType::TheCellType(),1,pelelm_dermixfrac,the_same_box);
+
+   // Progress variable
+   derive_lst.add("progress_variable",IndexType::TheCellType(),1,pelelm_derprogvar,the_same_box);
+
 #ifdef PELE_USE_EFIELD
    // Charge distribution
    derive_lst.add("chargedistrib",IndexType::TheCellType(),1,pelelm_derchargedist,the_same_box);

Source/PeleLMTemporals.cpp

@@ -355,10 +355,26 @@ void PeleLM::writeTemporals()
    }
    Real kinetic_energy = MFSum(GetVecOfConstPtrs(kinEnergy),0);
 
+   // Combustion
+   Real fuelConsumptionInt = 0.0;
+   Real heatReleaseRateInt = 0.0;
+   if (fuelID > 0 && !(m_chem_integrator == "ReactorNull")) {
+       fuelConsumptionInt =  MFSum(GetVecOfConstPtrs(getIRVect()),fuelID);
+       for (int lev = 0; lev <= finest_level; ++lev) { 
+          getHeatRelease(lev, kinEnergy[lev].get());  // Re-use kinEnergy container
+       }
+       heatReleaseRateInt = MFSum(GetVecOfConstPtrs(kinEnergy),0);
+   }
+
    // Get min/max/mean for non-species state components
 
-   tmpStateFile << m_nstep << " " << m_cur_time << " " << kinetic_energy << " " << m_pNew << " \n";
-   tmpStateFile.flush();
+   tmpStateFile << m_nstep << " " << m_cur_time                 // Time
+                << " " << kinetic_energy                        // Kinetic energy
+                << " " << m_pNew                                // Thermo. pressure
+                << " " << fuelConsumptionInt                    // Integ fuel burning rate
+                << " " << heatReleaseRateInt                    // Integ heat release rate
+                << " \n";
+   tmpStateFile.flush(); 
 }
 
 void PeleLM::openTempFile()