Add 2 new run-time diagnostics: PDF and conditional stats (#173)

* Add PDF and conditionals. WIP. * Fix mask in DiagPDF. * More WIP on Conditional. * Add Barrier in FramePlane. * Update diags to have better control over the outputed variables. PDF and conditional write plain ASCII. * Add a light descriptor of the reaction state. * Add coordinate as a derived. * Manage diagnostic variables. * Add reactVariables checks. One can now tag based on reaction rate of a species. * Add diagnostic example in FlameSheet 3D. * Trailing whitespace. * More trailing whitespace. * Add documentation on the new diagnostics. * Trailing whitespace in manual. * Fix undeclared in derived. * Fix capture of field indices in lambda. * Update Plane definition in TurbInflow regtest. * Fix geom vs. geomdata inconsistency in dervied. * Typo. * Add user defined derive. * Add header and formatting to PDF and conditionals diagnostic outputs. * Remove debugging statement. * Enable prescribed range with Conditionals. * use index for derived components * dont divide by zero * varnames for user defined derive * Missing update on intel CI. --------- Co-authored-by: Bruce Perry <[email protected]> Co-authored-by: Bruce Perry <[email protected]>
AMReX-Combustion · Feb 28, 2023 · 3fd1187 · 3fd1187
1 parent 54f7076
commit 3fd1187
Show file tree

Hide file tree

Showing 24 changed files with 1,279 additions and 133 deletions.
diff --git a/.github/workflows/intel.yml b/.github/workflows/intel.yml
@@ -32,8 +32,8 @@ jobs:
         set +e
         source /opt/intel/oneapi/setvars.sh
         set -e
-        export CXX=$(which dpcpp)
-        export CC=$(which icc)
+        export CXX=$(which icpx)
+        export CC=$(which icx)
         cd Exec/RegTests/FlameSheet/
         make TPL COMP=intel DIM=3 USE_SYCL=TRUE USE_MPI=FALSE
         make -j $(nproc) COMP=intel DIM=3 USE_SYCL=TRUE USE_MPI=FALSE
@@ -58,9 +58,8 @@ jobs:
         set +e
         source /opt/intel/oneapi/setvars.sh
         set -e
-        dpcpp -v
-        export CXX=$(which dpcpp)
-        export CC=$(which icc)
+        export CXX=$(which icpx)
+        export CC=$(which icx)
         cd Exec/RegTests/EB_EnclosedFlame/
         make TPL COMP=intel DIM=3 USE_SYCL=TRUE USE_MPI=FALSE
         make -j $(nproc) COMP=intel DIM=3 USE_SYCL=TRUE USE_MPI=FALSE
diff --git a/Docs/source/manual/LMeXControls.rst b/Docs/source/manual/LMeXControls.rst
@@ -148,6 +148,9 @@ The following list of derived variables are available in PeleLMeX:
     * - `rhominsumrhoY`
       - 1
       - Rho minus sum of rhoYs, for debug purposes
+    * - `coordinates`
+      - AMREX_SPACEDIM
+      - Cell-center coordinates
 
 Note that `mixture_fraction` and `progress_variable` requires additional inputs from the users as described below.
 
@@ -331,28 +334,58 @@ state:
     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:
+Analysing the data a-posteriori can become extremely cumbersome when dealing with extreme datasets.
+PeleLMeX offers a set of diagnostics available at runtime and more are under development.
+Currently, the list of diagnostic contains:
+
+* `DiagFramePlane` : extract a plane aligned in the 'x','y' or 'z' direction across the AMR hierarchy, writing
+  a 2D plotfile compatible with Amrvis, Paraview or yt. Only available for 3D simulations.
+* `DiagPDF` : extract the PDF of a given variable and write it to an ASCII file.
+* `DiagConditional` : extract statistics (average and standard deviation, integral or sum) of a
+  set of variables conditioned on the value of given variable and write it to an ASCII file.
+
+When using `DiagPDF` or `DiagConditional`, it is possible to narrow down the diagnostic to a region of interest
+by specifying a set of filters, defining a range of interest for a variable. Note also the for these two diagnostics,
+fine-covered regions are masked. The following provide examples for each diagnostic:
 
 ::
 
     #--------------------------DIAGNOSTICS------------------------
 
-    peleLM.diagnostics = xnormal ynormal
-    peleLM.xnormal.type = DiagFramePlane
-    peleLM.xnormal.file = xNorm5mm
-    peleLM.xnormal.normal = 0
-    peleLM.xnormal.center = 0.005
-    peleLM.xnormal.int    = 5
-    peleLM.xnormal.interpolation = Linear
-
-    peleLM.ynormal.type = DiagFramePlane
-    peleLM.ynormal.file = yNormCent
-    peleLM.ynormal.normal = 1
-    peleLM.ynormal.center = 0.0
-    peleLM.ynormal.int    = 10
-    peleLM.ynormal.interpolation = Quadratic
-
+    peleLM.diagnostics = xnormP condT pdfTest
+
+    peleLM.xnormP.type = DiagFramePlane                             # Diagnostic type
+    peleLM.xnormP.file = xNorm5mm                                   # Output file prefix
+    peleLM.xnormP.normal = 0                                        # Plane normal (0, 1 or 2 for x, y or z)
+    peleLM.xnormP.center = 0.005                                    # Coordinate in the normal direction
+    peleLM.xnormP.int    = 5                                        # Frequency (as step #) for performing the diagnostic
+    peleLM.xnormP.interpolation = Linear                            # [OPT, DEF=Linear] Interpolation type : Linear or Quadratic
+    peleLM.xnormP.field_names = x_velocity mag_vort density         # List of variables outputed to the 2D pltfile
+
+    peleLM.condT.type = DiagConditional                             # Diagnostic type
+    peleLM.condT.file = condTest                                    # Output file prefix
+    peleLM.condT.int  = 5                                           # Frequency (as step #) for performing the diagnostic
+    peleLM.condT.filters = xHigh stoich                             # [OPT, DEF=None] List of filters
+    peleLM.condT.xHigh.field_name = x                               # Filter field
+    peleLM.condT.xHigh.value_greater = 0.006                        # Filter definition : value_greater, value_less, value_inrange
+    peleLM.condT.stoich.field_name = mixture_fraction               # Filter field
+    peleLM.condT.stoich.value_inrange = 0.053 0.055                 # Filter definition : value_greater, value_less, value_inrange
+    peleLM.condT.conditional_type = Average                         # Conditional type : Average, Integral or Sum
+    peleLM.condT.nBins = 50                                         # Number of bins for the conditioning variable
+    peleLM.condT.condition_field_name = temp                        # Conditioning variable name
+    peleLM.condT.field_names = HeatRelease I_R(CH4) I_R(H2)         # List of variables to be treated
+
+    peleLM.pdfTest.type = DiagPDF                                   # Diagnostic type
+    peleLM.pdfTest.file = PDFTest                                   # Output file prefix
+    peleLM.pdfTest.int  = 5                                         # Frequency (as step #) for performing the diagnostic
+    peleLM.pdfTest.filters = innerFlame                             # [OPT, DEF=None] List of filters
+    peleLM.pdfTest.innerFlame.field_name = temp                     # Filter field
+    peleLM.pdfTest.innerFlame.value_inrange = 450.0 1500.0          # Filter definition : value_greater, value_less, value_inrange
+    peleLM.pdfTest.nBins = 50                                       # Number of bins for the PDF
+    peleLM.pdfTest.normalized = 1                                   # [OPT, DEF=1] PDF is normalized (i.e. integral is unity) ?
+    peleLM.pdfTest.volume_weighted = 1                              # [OPT, DEF=1] Computation of the PDF is volume weighted ?
+    peleLM.pdfTest.range = 0.0 2.0                                  # [OPT, DEF=data min/max] Specify the range of the PDF
+    peleLM.pdfTest.field_name = x_velocity                          # Variable of interest
 
 Run-time control
 --------------------

diff --git a/Exec/RegTests/FlameSheet/input.3d-regt b/Exec/RegTests/FlameSheet/input.3d-regt
@@ -10,9 +10,8 @@ geometry.prob_hi     = 0.016 0.016 0.032        # x_hi y_hi (z_hi)
 peleLM.lo_bc = Interior Interior Inflow
 peleLM.hi_bc = Interior Interior Outflow
 
-
 #-------------------------AMR CONTROL----------------------------
-amr.n_cell          = 32 32 64      # Level 0 number of cells in each direction   
+amr.n_cell          = 32 32 64      # Level 0 number of cells in each direction
 amr.v               = 1                # AMR verbose
 amr.max_level       = 1                # maximum level number allowed
 amr.ref_ratio       = 2 2 2 2          # refinement ratio
@@ -57,8 +56,46 @@ peleLM.chem_integrator = "ReactorCvode"
 peleLM.use_typ_vals_chem = 1          # Use species/temp typical values in CVODE
 ode.rtol = 1.0e-6                     # Relative tolerance of the chemical solve
 ode.atol = 1.0e-5                     # Absolute tolerance factor applied on typical values
-cvode.solve_type = GMRES               # CVODE Linear solve type (for Newton direction) 
-cvode.max_order  = 4                   # CVODE max BDF order. 
+cvode.solve_type = GMRES               # CVODE Linear solve type (for Newton direction)
+cvode.max_order  = 4                   # CVODE max BDF order.
+
+#-------------------DIAGNOSTIC CONTROL-----------------------
+peleLM.diagnostics = xnormal ynormal condT pdfVel
+peleLM.xnormal.type = DiagFramePlane
+peleLM.xnormal.file = xNorm5mm
+peleLM.xnormal.normal = 0
+peleLM.xnormal.center = 0.008
+peleLM.xnormal.int    = 5
+peleLM.xnormal.interpolation = Linear
+peleLM.xnormal.field_names = Y(O2) x_velocity temp
+
+peleLM.ynormal.type = DiagFramePlane
+peleLM.ynormal.file = yNormCent
+peleLM.ynormal.normal = 1
+peleLM.ynormal.center = 0.008
+peleLM.ynormal.int    = 5
+peleLM.ynormal.interpolation = Quadratic
+peleLM.ynormal.field_names = density Y(H2) Y(CH4) mag_vort
+
+peleLM.condT.type = DiagConditional
+peleLM.condT.file = condTest
+peleLM.condT.int  = 2
+peleLM.condT.filters = innerXCore
+peleLM.condT.innerXCore.field_name = x
+peleLM.condT.innerXCore.value_inrange = 0.006 0.010
+peleLM.condT.conditional_type = Average
+peleLM.condT.nBins = 50
+peleLM.condT.condition_field_name = temp
+peleLM.condT.field_names = HeatRelease I_R(CH4) I_R(H2)
+
+peleLM.pdfVel.type = DiagPDF
+peleLM.pdfVel.file = PDFTest
+peleLM.pdfVel.int  = 2
+peleLM.pdfVel.filters = innerFlame
+peleLM.pdfVel.innerFlame.field_name = temp
+peleLM.pdfVel.innerFlame.value_inrange = 350 1500
+peleLM.pdfVel.nBins = 50
+peleLM.pdfVel.field_name = x_velocity
 
 #--------------------REFINEMENT CONTROL------------------------
 #amr.refinement_indicators = temp

diff --git a/Exec/RegTests/TurbInflow/input.3d b/Exec/RegTests/TurbInflow/input.3d
@@ -71,6 +71,7 @@ peleLM.testPlane.type = DiagFramePlane
 peleLM.testPlane.normal = 0
 peleLM.testPlane.center = 0.005
 peleLM.testPlane.int    = 5
+peleLM.testPlane.field_names = x_velocity y_velocity z_velocity mag_vort
 
 #--------------------REFINEMENT CONTROL------------------------
 amr.refinement_indicators = O2

diff --git a/Source/DeriveUserDefined.cpp b/Source/DeriveUserDefined.cpp
@@ -0,0 +1,27 @@
+#include <PeleLMDeriveFunc.H>
+#include <PeleLM_Index.H>
+#include <PelePhysics.H>
+#include <mechanism.H>
+#include <PeleLM.H>
+#include <PeleLM_K.H>
+
+using namespace amrex;
+
+//
+// User-defined derived variables list
+//
+Vector<std::string> pelelm_setuserderives()
+{
+  //Vector<std::string> var_names({"derUserDefine_null"});
+  return {"derUserDefine_null"}; //var_names;
+}
+
+//
+// User-defined derived definition
+//
+void pelelm_deruserdef (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>& /*bcrec*/, int /*level*/)
+{
+    Abort("Using derUserDefine derived requires providing a definition in local DeriveUserDefined.cpp");
+}
diff --git a/Source/Diagnostics/DiagBase.H b/Source/Diagnostics/DiagBase.H
@@ -3,14 +3,16 @@
 
 #include <AMReX_ParmParse.H>
 #include <AMReX_MultiFab.H>
+#include "DiagFilter.H"
 #include "Factory.H"
 
 class DiagBase : public pele::physics::Factory<DiagBase>
 {
 public:
   static std::string base_identifier() { return "DiagBase"; }
 
-  virtual void init(const std::string &a_prefix) = 0;
+  virtual void init(const std::string &a_prefix,
+                    std::string_view a_diagName);
 
   virtual void close() = 0;
 
@@ -21,20 +23,28 @@ public:
   virtual void prepare(int a_nlevels,
                        const amrex::Vector<amrex::Geometry> &a_geoms,
                        const amrex::Vector<amrex::BoxArray> &a_grids,
-                       const amrex::Vector<amrex::DistributionMapping> &a_dmap) = 0;
+                       const amrex::Vector<amrex::DistributionMapping> &a_dmap,
+                       const amrex::Vector<std::string> &a_varNames);
 
   virtual void processDiag(int a_nstep,
                            const amrex::Real &a_time,
                            const amrex::Vector<const amrex::MultiFab*> &a_state,
-                           const amrex::Vector<std::string> &a_varNames) =0;
+                           const amrex::Vector<std::string> &a_varNames) = 0;
+
+  virtual void addVars(amrex::Vector<std::string> &a_varList);
+
+  int getFieldIndex(const std::string &a_field,
+                    const amrex::Vector<std::string> &a_varList);
 
 protected:
-  std::string m_diagfile;
-  int m_verbose = 0;
-  amrex::Real m_per = -1.0;
-  int m_interval = -1;
-  bool need_update = true;
-  bool first_time = true;
+  std::string m_diagfile{""};
+  int m_verbose{0};
+  amrex::Real m_per{-1.0};
+  int m_interval{-1};
+  bool need_update{true};
+  bool first_time{true};
+  amrex::Vector<DiagFilter> m_filters{};
+  amrex::Gpu::DeviceVector<DiagFilterData> m_filterData;
 };
 
 #endif
diff --git a/Source/Diagnostics/DiagBase.cpp b/Source/Diagnostics/DiagBase.cpp
@@ -1,4 +1,57 @@
 #include "DiagBase.H"
+#include "AMReX_ParmParse.H"
+
+void
+DiagBase::init(const std::string &a_prefix,
+               std::string_view a_diagName)
+{
+    amrex::ParmParse pp(a_prefix);
+
+    // IO
+    pp.query("int", m_interval);
+    pp.query("per", m_per);
+    m_diagfile = a_diagName;
+    pp.query("file",m_diagfile);
+    AMREX_ASSERT(m_interval>0 || m_per>0.0);
+
+    // Filters
+    int nFilters = 0;
+    nFilters = pp.countval("filters");
+    amrex::Vector<std::string> filtersName;
+    if (nFilters > 0) {
+        m_filters.resize(nFilters);
+        filtersName.resize(nFilters);
+    }
+    for (int n = 0; n < nFilters; ++n) {
+        pp.get("filters",filtersName[n],n);
+        const std::string filter_prefix = a_prefix + "." + filtersName[n];
+        m_filters[n].init(filter_prefix);
+    }
+}
+
+void
+DiagBase::prepare(int a_nlevels,
+                  const amrex::Vector<amrex::Geometry> &a_geoms,
+                  const amrex::Vector<amrex::BoxArray> &a_grids,
+                  const amrex::Vector<amrex::DistributionMapping> &a_dmap,
+                  const amrex::Vector<std::string> &a_varNames)
+{
+    if (first_time) {
+        int nFilters = m_filters.size();
+        // Move the filter data to the device
+        for (int n = 0; n < nFilters; ++n) {
+            m_filters[n].setup(a_varNames);
+        }
+        amrex::Vector<DiagFilterData> hostFilterData;
+        for (int n = 0; n < nFilters; ++n) {
+            hostFilterData.push_back(m_filters[n].m_fdata);
+        }
+        m_filterData.resize(nFilters);
+        amrex::Gpu::copy(amrex::Gpu::hostToDevice,
+                         hostFilterData.begin(),hostFilterData.end(),
+                         m_filterData.begin());
+    }
+}
 
 bool
 DiagBase::doDiag(const amrex::Real &a_time,
@@ -14,3 +67,28 @@ DiagBase::doDiag(const amrex::Real &a_time,
 
     return willDo;
 }
+
+void
+DiagBase::addVars(amrex::Vector<std::string> &a_varList) {
+    int nFilters = m_filters.size();
+    for (int n = 0; n < nFilters; ++n) {
+        a_varList.push_back(m_filters[n].m_filterVar);
+    }
+}
+
+int
+DiagBase::getFieldIndex(const std::string &a_field,
+                        const amrex::Vector<std::string> &a_varList)
+{
+    int index = -1;
+    for (int n{0}; n < a_varList.size(); ++n) {
+        if (a_varList[n] == a_field) {
+            index = n;
+            break;
+        }
+    }
+    if (index < 0) {
+        amrex::Abort("Field : "+a_field+" wasn't found in available fields");
+    }
+    return index;
+}