Convection discretization rewrite

python/ember/_ember.pxd

@@ -57,6 +57,7 @@ cdef extern from "readConfig.h":
         double globalTimestep, diffusionTimestepMultiplier
 
         string splittingMethod, chemistryIntegrator, rateMultiplierFunctionType
+        string convectionScheme
         double integratorRelTol, integratorMomentumAbsTol, integratorEnergyAbsTol
         double integratorSpeciesAbsTol, integratorMinTimestep
 
@@ -87,7 +88,7 @@ cdef extern from "readConfig.h":
         double ignition_center, ignition_stddev
         cbool alwaysUpdateHeatFlux
 
-        double vtol, dvtol, rmTol, dampConst, gridMin, gridMax
+        double errTol, rmTol, dampConst, gridMin, gridMax
         double uniformityTol, absvtol
         double boundaryTol, boundaryTolRm, unstrainedDownstreamWidth
         int addPointCount

python/ember/_ember.pyx

@@ -201,6 +201,7 @@ cdef class ConfigOptions:
         opts.twinFlame = G.twinFlame
         opts.chemistryIntegrator = stringify(G.chemistryIntegrator)
         opts.splittingMethod = stringify(G.splittingMethod)
+        opts.convectionScheme = stringify(G.convectionScheme)
         opts.setContinuityBC(stringify(G.continuityBC))
         opts.errorStopCount = G.errorStopCount
         opts.stopIfError = G.errorStopCount > 0
@@ -291,8 +292,7 @@ cdef class ConfigOptions:
 
         # Grid
         opts.centerGridMin = self.grid.centerGridMin
-        opts.vtol = self.grid.vtol
-        opts.dvtol = self.grid.dvtol
+        opts.errTol = self.grid.errTol
         opts.rmTol = self.grid.rmTol
         opts.dampConst = self.grid.dampConst
         opts.gridMax = self.grid.gridMax

python/ember/examples/example_laminarFlameSpeed.py

@@ -23,7 +23,10 @@
     General(fixedLeftLocation=True,
             fixedBurnedVal=False,
             nThreads=4),
-    Grid(vtol=0.1, dvtol=0.15, gridMin=5e-6, gridMax=0.001),
+    # errTol=5e-4 is the docstring's "minimal accuracy" tier, chosen to
+    # preserve this example's historical resolution/runtime (the library
+    # default errTol=1e-4 is finer).
+    Grid(errTol=5e-4, gridMin=5e-6, gridMax=0.001),
     PositionControl(proportionalGain=2000, xInitial=0.005, xFinal=0.005),
     TerminationCondition(tolerance=1e-5),
     Times(profileStepInterval=50))

python/ember/input.py

@@ -335,6 +335,13 @@ class General(Options):
     #: terms. Options are ``strang`` and ``balanced``.
     splittingMethod = StringOption("balanced", ("strang",), level=2)
 
+    #: Discretization scheme to use for the convection term. The default,
+    #: ``secondOrderLimited``, is 2nd-order accurate in smooth regions and
+    #: falls back locally to a 1st-order limiter at extrema. ``firstOrderUpwind``
+    #: is the legacy scheme used prior to version 1.7.
+    convectionScheme = StringOption("secondOrderLimited",
+                                    ("firstOrderUpwind",), level=2)
+
     #: Number of integration failures to tolerate in the chemistry
     #: integrator before aborting.
     errorStopCount = IntegerOption(100, level=2)
@@ -387,17 +394,29 @@ class Chemistry(Options):
 
 class Grid(Options):
     """ Parameters controlling the adaptive grid """
-    #: Maximum relative scalar variation of each state vector
-    #: component between consecutive grid points. For high accuracy,
-    #: ``vtol = 0.08``; For minimal accuracy, ``vtol = 0.20``.
-    vtol = FloatOption(0.12)
-
-    #: Maximum relative variation of the gradient of each state vector
-    #: component between consecutive grid points. For high accuracy,
-    #: ``dvtol = 0.12``; For minimal accuracy, ``dvtol = 0.4``.
-    dvtol = FloatOption(0.2)
-
-    #: Relative tolerance (compared to vtol and dvtol) for grid point removal.
+    #: Target relative accuracy of the solution on the adapted grid. The
+    #: estimated local truncation error of each state-vector component is
+    #: kept below ``errTol`` times that component's range, using an error
+    #: estimate that accounts for the order of the selected
+    #: ``general.convectionScheme``. The same value yields similar solution
+    #: accuracy under either scheme for the validation cases (highly strained
+    #: flames may need a few-times tighter ``errTol`` under
+    #: ``firstOrderUpwind`` for matched accuracy); the higher-order scheme
+    #: needs fewer grid points. On the calibration study cases the default
+    #: gives consumption-speed errors of a few parts in 10^4 with grids of
+    #: roughly 100-170 points. For high accuracy, ``errTol = 2e-5``; for
+    #: minimal accuracy, ``errTol = 5e-4``.
+    errTol = FloatOption(1e-4, min=1e-12)
+
+    #: Deprecated and ignored. Grid resolution is controlled by
+    #: :attr:`errTol`.
+    vtol = FloatOption(None, level=3)
+
+    #: Deprecated and ignored. Grid resolution is controlled by
+    #: :attr:`errTol`.
+    dvtol = FloatOption(None, level=3)
+
+    #: Relative tolerance (compared to errTol) for grid point removal.
     rmTol = FloatOption(0.6, level=1)
 
     #: Parameter to limit numerical diffusion in regions with high
@@ -910,9 +929,6 @@ def __init__(self, *args):
         self.terminationCondition = get(TerminationCondition)
         self.extinction = get(Extinction)
 
-    def evaluate(self):
-        return ConcreteConfig(self)
-
     def __iter__(self):
         for item in self.__dict__.values():
             if isinstance(item, Options):
@@ -927,7 +943,29 @@ def stringify(self):
 
         return 'conf = Config(\n' + ',\n'.join(ans) + ')\n'
 
+    def _warnDeprecated(self):
+        if getattr(self, '_deprecationWarned', False):
+            return
+        self._deprecationWarned = True
+        if self.grid.vtol.isSet or self.grid.dvtol.isSet:
+            print("WARNING: 'grid.vtol' and 'grid.dvtol' are deprecated and have"
+                  " no effect.\n"
+                  "         Grid resolution is now controlled by the local-error"
+                  " tolerance 'grid.errTol';\n"
+                  "         the same errTol gives similar accuracy for either"
+                  " convection scheme (for strongly strained flames the"
+                  " first-order scheme may need a tighter errTol for matched"
+                  " accuracy).")
+        if self.grid.errTol.value is not None and self.grid.errTol.value > 0.5:
+            print("WARNING: 'grid.errTol' > 0.5 effectively disables grid"
+                  " refinement.")
+
+    def evaluate(self):
+        self._warnDeprecated()
+        return ConcreteConfig(self)
+
     def validate(self):
+        self._warnDeprecated()
         error = False
         cylindricalFlame = True if self.general.flameGeometry == 'cylindrical' else False
         discFlame = True if self.general.flameGeometry == 'disc' else False

src/convectionDifferencer.cpp

@@ -0,0 +1,81 @@
+#include "convectionDifferencer.h"
+
+#include <cassert>
+
+void ConvectionDifferencer::setScheme(Scheme s)
+{
+    scheme = s;
+}
+
+void ConvectionDifferencer::resize(size_t nPoints)
+{
+    sigma.resize(nPoints);
+}
+
+void ConvectionDifferencer::computeDerivatives(const dvec& y, const dvec& v,
+                                               const OneDimGrid& grid, dvec& dydx)
+{
+    const size_t jj = grid.jj;
+    const dvec& hh = grid.hh;
+
+    if (scheme == Scheme::FirstOrderUpwind) {
+        // Legacy 1st-order upwind, reproduced bit-identically from the original
+        // convection loops: forward difference where the advecting velocity is
+        // negative or at the left boundary node, backward difference otherwise.
+        for (size_t j = 0; j < jj; j++) {
+            if (v[j] < 0 || j == 0) {
+                dydx[j] = (y[j+1] - y[j]) / hh[j];
+            } else {
+                dydx[j] = (y[j] - y[j-1]) / hh[j-1];
+            }
+        }
+        return;
+    }
+
+    // *** SecondOrderLimited ***
+    const size_t nPoints = grid.nPoints;
+    assert(nPoints >= 3);
+    if (static_cast<size_t>(sigma.size()) != nPoints) {
+        sigma.resize(nPoints);
+    }
+    const dvec& dlj = grid.dlj;
+
+    // Pass 1: van Albada limited node slopes. The outermost nodes fall back to
+    // a zero slope, making the outermost face reconstructions locally 1st-order
+    // (these are small-gradient far-field regions).
+    sigma[0] = 0.0;
+    sigma[jj] = 0.0;
+    for (size_t j = 1; j < jj; j++) {
+        double sm = (y[j] - y[j-1]) / hh[j-1]; // one-sided slope from the left
+        double sp = (y[j+1] - y[j]) / hh[j];   // one-sided slope from the right
+        if (sm * sp > 0) {
+            // van Albada limiter. Returns exactly the common slope when
+            // sm == sp (so linear data yields the exact derivative). The
+            // sm*sp > 0 guard both clips at extrema and avoids the 0/0 form on
+            // flat data (sm == sp == 0), without an epsilon that would break
+            // the linear-exactness property.
+            sigma[j] = (sm * sp) * (sm + sp) / (sm * sm + sp * sp);
+        } else {
+            sigma[j] = 0.0;
+        }
+    }
+
+    // Pass 2: upwind-biased face reconstruction. The advective derivative at
+    // node j is the difference of the reconstructed face values divided by the
+    // node spacing dlj[j] = (hh[j-1] + hh[j])/2. Each face value is
+    // reconstructed from the upwind side selected by sign(v[j]). dydx[0] is
+    // unused by the callers (node 0 uses a boundary closure).
+    dydx[0] = 0.0;
+    for (size_t j = 1; j < jj; j++) {
+        double yFaceR; // reconstruction at the right face x_{j+1/2}
+        double yFaceL; // reconstruction at the left face  x_{j-1/2}
+        if (v[j] < 0) {
+            yFaceR = y[j+1] - sigma[j+1] * hh[j] / 2;
+            yFaceL = y[j]   - sigma[j]   * hh[j-1] / 2;
+        } else {
+            yFaceR = y[j]   + sigma[j]   * hh[j] / 2;
+            yFaceL = y[j-1] + sigma[j-1] * hh[j-1] / 2;
+        }
+        dydx[j] = (yFaceR - yFaceL) / dlj[j];
+    }
+}

src/convectionDifferencer.h

@@ -0,0 +1,59 @@
+#pragma once
+
+#include "mathUtils.h"
+#include "grid.h"
+
+//! Computes upwind-biased advective derivatives \f$ \partial y / \partial x \f$
+//! on a nonuniform grid for the convection sub-problem.
+//!
+//! Two schemes are provided:
+//!   - #Scheme::FirstOrderUpwind reproduces the legacy 1st-order upwind
+//!     stencil bit-identically (forward difference where the advecting velocity
+//!     is negative or at the left boundary node, backward difference
+//!     otherwise).
+//!   - #Scheme::SecondOrderLimited applies a van Albada limited, upwind-biased
+//!     face reconstruction that is 2nd-order in smooth regions regardless of
+//!     grid nonuniformity and falls back to local 1st-order upwind at extrema.
+//!
+//! The kernel computes derivatives at nodes `0 .. jj-1`. The right-boundary
+//! node `jj` uses a boundary-condition-specific closure and is handled by the
+//! callers.
+class ConvectionDifferencer
+{
+public:
+    //! Selects the discretization used by computeDerivatives().
+    enum class Scheme {
+        FirstOrderUpwind, //!< Legacy 1st-order upwind (bit-identical)
+        SecondOrderLimited //!< van Albada limited ~2nd-order upwind-biased
+    };
+
+    //! Select the active discretization scheme.
+    void setScheme(Scheme s);
+
+    //! Size the internal limiter scratch vector. This is the intended
+    //! pre-allocation entry point: call it with the current number of grid
+    //! points before computeDerivatives() to avoid a reallocation on first
+    //! use. computeDerivatives() also self-sizes the scratch vector
+    //! defensively if it is ever found to be the wrong size (e.g. after a
+    //! regrid where resize() was not called again), so calling it is not
+    //! strictly required for correctness, only for avoiding that fallback.
+    void resize(size_t nPoints);
+
+    //! Compute the advective derivative \f$ dy/dx \f$ at nodes `0 .. jj-1`.
+    //!
+    //! The upwind branch is chosen per node from the sign of the advecting
+    //! velocity `v[j]` (forward difference where `v[j] < 0`). The right-boundary
+    //! node `jj` is not written; callers apply their own boundary closure there.
+    //!
+    //! @param y    Field values at each grid point.
+    //! @param v    Advecting velocity at each grid point (only its sign is used
+    //!             for branch selection).
+    //! @param grid Grid providing the spacing arrays `hh` and `dlj`.
+    //! @param dydx Output derivative, written for indices `0 .. jj-1`.
+    void computeDerivatives(const dvec& y, const dvec& v,
+                            const OneDimGrid& grid, dvec& dydx);
+
+private:
+    Scheme scheme = Scheme::SecondOrderLimited;
+    dvec sigma; //!< limited node slopes (van Albada), one per grid point
+};