first commit for feature

Author: wandadars

include/cantera/numerics/SteadyStateSystem.h

@@ -72,6 +72,11 @@ class SteadyStateSystem
     //! x. On return, array r contains the steady-state residual values.
     double ssnorm(double* x, double* r);
 
+    //! Transient max norm (infinity norm) of the residual evaluated using solution
+    //! x and the current timestep (rdt). On return, array r contains the
+    //! transient residual values.
+    double tsnorm(double* x, double* r);
+
     //! Total solution vector length;
     size_t size() const {
         return m_size;
@@ -266,6 +271,9 @@ class SteadyStateSystem
     //! @param[out] rsd  Storage for the residual, length size()
     void evalSSJacobian(double* x, double* rsd);
 
+    //! Determine the timestep growth factor after a successful step.
+    double timeStepIncreaseFactor(double* x_before, double* x_after);
+
     //! Array of number of steps to take after each unsuccessful steady-state solve
     //! before re-attempting the steady-state solution. For subsequent time stepping
     //! calls, the final value is reused. See setTimeStep().
@@ -325,6 +333,7 @@ class SteadyStateSystem
     double m_jacobianRelPerturb = 1e-5;
     //! Absolute perturbation of each component in finite difference Jacobian
     double m_jacobianAbsPerturb = 1e-10;
+
 };
 
 }

include/cantera/oneD/MultiNewton.h

@@ -35,6 +35,11 @@ class MultiNewton
         return m_n;
     }
 
+    //! Number of Newton iterations taken in the most recent solve() call.
+    int lastIterations() const {
+        return m_lastIterations;
+    }
+
     //! Compute the undamped Newton step.  The residual function is evaluated
     //! at `x`, but the Jacobian is not recomputed.
     //! @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object.
@@ -164,13 +169,16 @@ class MultiNewton
     double m_dampFactor = sqrt(2.0);
 
     //! Maximum number of damping iterations
-    size_t m_maxDampIter = 7;
+    size_t m_maxDampIter = 14;
 
     //! number of variables
     size_t m_n;
 
     //! Elapsed CPU time spent computing the Jacobian.
     double m_elapsed = 0.0;
+
+    //! Number of Newton iterations taken in the last solve() call.
+    int m_lastIterations = 0;
 };
 }
 

include/cantera/oneD/Sim1D.h

@@ -337,6 +337,17 @@ class Sim1D : public OneDim
         m_steady_callback = callback;
     }
 
+    //! Set the maximum number of regrid attempts after a timestep failure.
+    //! Set to 0 to disable regrid-on-timestep-failure retries.
+    void setTimeStepRegridMax(int nmax) {
+        m_ts_regrid_max = nmax;
+    }
+
+    //! Get the maximum number of regrid attempts after a timestep failure.
+    int timeStepRegridMax() const {
+        return m_ts_regrid_max;
+    }
+
 protected:
     //! the solution vector after the last successful steady-state solve (stored
     //! before grid refinement)
@@ -349,6 +360,9 @@ class Sim1D : public OneDim
     //! User-supplied function called after a successful steady-state solve.
     Func1* m_steady_callback;
 
+    //! 0 disables regrid-on-timestep-failure retries
+    int m_ts_regrid_max = 10;
+
 private:
     //! Calls method _finalize in each domain.
     void finalize();

interfaces/cython/cantera/_onedim.pxd

@@ -151,6 +151,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void setTimeStepFactor(double)
         void setMinTimeStep(double)
         void setMaxTimeStep(double)
+        void setTimeStepRegridMax(int)
+        int timeStepRegridMax()
         void setMaxGridPoints(int, size_t) except +translate_exception
         size_t maxGridPoints(size_t) except +translate_exception
         void setGridMin(int, double) except +translate_exception

interfaces/cython/cantera/_onedim.pyi

@@ -297,6 +297,7 @@ class Sim1D:
     def max_time_step_count(self) -> int: ...
     @max_time_step_count.setter
     def max_time_step_count(self, nmax: int) -> None: ...
+    def set_time_step_regrid(self, max_tries: int = 10) -> None: ...
     def set_jacobian_perturbation(
         self, relative: float, absolute: float, threshold: float
     ) -> None: ...

interfaces/cython/cantera/_onedim.pyx

@@ -1425,6 +1425,14 @@ cdef class Sim1D:
         """
         self.sim.setTimeStepFactor(tfactor)
 
+    def set_time_step_regrid(self, max_tries=10):
+        """
+        Set the maximum number of regrid attempts after a time step failure.
+
+        Setting ``max_tries`` to zero disables regrid-on-timestep-failure retries.
+        """
+        self.sim.setTimeStepRegridMax(max_tries)
+
     def set_min_time_step(self, tsmin):
         """ Set the minimum time step. """
         self.sim.setMinTimeStep(tsmin)

src/numerics/SteadyStateSystem.cpp

@@ -108,6 +108,46 @@ void SteadyStateSystem::solve(int loglevel)
     }
 }
 
+double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_after)
+{
+    int heuristic = 2; // Can be 1-4
+    m_work1.resize(m_size);
+    const double grow_factor = 1.25;
+
+    if (heuristic == 1) {
+        // Steady-state residual gate. If the steady-state residual decreased, allow
+        // increase in timestep.
+        double ss_before = ssnorm(x_before, m_work1.data());
+        double ss_after = ssnorm(x_after, m_work1.data());
+        return (ss_after < ss_before) ? grow_factor : 1.0;
+    } else if (heuristic == 2) {
+        // Transient residual gate. If the transient residual decreased, allow
+        // increase in timestep.
+        double ts_before = tsnorm(x_before, m_work1.data());
+        double ts_after = tsnorm(x_after, m_work1.data());
+        return (ts_after < ts_before) ? grow_factor : 1.0;
+    } else if (heuristic == 3) {
+        // Residual-ratio scaling gate (transient residual). If the transient residual
+        // decreased significantly, allow larger increase in timestep.
+        double ts_before = tsnorm(x_before, m_work1.data());
+        double ts_after = tsnorm(x_after, m_work1.data());
+        if (!(ts_after > 0.0) || !(ts_before > ts_after)) {
+            return 1.0;
+        }
+        const double max_growth = 1.5;
+        const double exponent = 0.2;
+        double ratio = ts_before / ts_after;
+        double factor = std::pow(ratio, exponent);
+        return std::min(max_growth, std::max(1.0, factor));
+    } else if (heuristic == 4) {
+        // Newton-iteration gate. If the last Newton solve required only a few
+        // iterations, allow increase in timestep.
+        const int max_iters_for_growth = 3;
+        return (newton().lastIterations() <= max_iters_for_growth) ? grow_factor : 1.0;
+    }
+    return 1.0;
+}
+
 double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r, int loglevel)
 {
     // set the Jacobian age parameter to the transient value
@@ -123,12 +163,14 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
         writelog("============================");
     }
     while (n < nsteps) {
-        if (loglevel == 1) { // At level 1, output concise information
-            double ss = ssnorm(x, r);
-            writelog("\n{:<5d}  {:<6.4e}   {:>7.4f}", n, dt, log10(ss));
-        } else if (loglevel > 1) {
-            double ss = ssnorm(x, r);
-            writelog("\nTimestep ({}) dt= {:<11.4e}  log(ss)= {:<7.4f}", n, dt, log10(ss));
+        if (loglevel >= 1) {
+            double ss_before = ssnorm(x, r);
+            if (loglevel == 1) { // At level 1, output concise information
+                writelog("\n{:<5d}  {:<6.4e}   {:>7.4f}", n, dt, log10(ss_before));
+            } else {
+                writelog("\nTimestep ({}) dt= {:<11.4e}  log(ss)= {:<7.4f}", n, dt,
+                         log10(ss_before));
+            }
         }
 
         // set up for time stepping with stepsize dt
@@ -148,17 +190,18 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
             successiveFailures = 0;
             m_nsteps++;
             n += 1;
-            copy(r, r + m_size, x);
-            // No Jacobian evaluations were performed, so a larger timestep can be used
+            double grow_factor = 1.0;
             if (m_jac->nEvals() == j0) {
-                dt *= 1.5;
+                grow_factor = timeStepIncreaseFactor(x, r);
             }
+            copy(r, r + m_size, x);
+            dt *= grow_factor;
             if (m_time_step_callback) {
                 m_time_step_callback->eval(dt);
             }
             dt = std::min(dt, m_tmax);
             if (m_nsteps >= m_nsteps_max) {
-                throw CanteraError("OneDim::timeStep",
+                throw CanteraError("SteadyStateSystem::timeStep",
                     "Took maximum number of timesteps allowed ({}) without "
                     "reaching steady-state solution.", m_nsteps_max);
             }
@@ -181,7 +224,7 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
                 if (dt < m_tmin) {
                     string err_msg = fmt::format(
                         "Time integration failed. Minimum timestep ({}) reached.", m_tmin);
-                    throw CanteraError("OneDim::timeStep", err_msg);
+                    throw CanteraError("SteadyStateSystem::timeStep", err_msg);
                 }
             }
         }
@@ -212,6 +255,16 @@ double SteadyStateSystem::ssnorm(double* x, double* r)
     return ss;
 }
 
+double SteadyStateSystem::tsnorm(double* x, double* r)
+{
+    eval(x, r, m_rdt, 0);
+    double ts = 0.0;
+    for (size_t i = 0; i < m_size; i++) {
+        ts = std::max(fabs(r[i]), ts);
+    }
+    return ts;
+}
+
 void SteadyStateSystem::setTimeStep(double stepsize, size_t n, const int* tsteps)
 {
     m_tstep = stepsize;
@@ -261,7 +314,7 @@ void SteadyStateSystem::initTimeInteg(double dt, double* x)
     m_rdt = 1.0/dt;
 
     // if the stepsize has changed, then update the transient part of the Jacobian
-    if (fabs(rdt_old - m_rdt) > Tiny) {
+    if (fabs(rdt_old - m_rdt) > Tiny && m_jac_ok) {
         m_jac->updateTransient(m_rdt, m_mask.data());
     }
 }
@@ -273,7 +326,9 @@ void SteadyStateSystem::setSteadyMode()
     }
 
     m_rdt = 0.0;
-    m_jac->updateTransient(m_rdt, m_mask.data());
+    if (m_jac_ok) {
+        m_jac->updateTransient(m_rdt, m_mask.data());
+    }
 }
 
 void SteadyStateSystem::setJacAge(int ss_age, int ts_age)

src/oneD/MultiNewton.cpp

@@ -206,11 +206,13 @@ int MultiNewton::solve(double* x0, double* x1, SteadyStateSystem& r, int logleve
     double s1=1.e30;
 
     copy(x0, x0 + m_n, &m_x[0]);
+    m_lastIterations = 0;
 
     double rdt = r.rdt();
     int nJacReeval = 0;
     auto jac = r.linearSolver();
     while (true) {
+        m_lastIterations++;
         // Check whether the Jacobian should be re-evaluated.
         if (jac->age() > m_maxAge) {
             if (loglevel > 1) {

src/oneD/Sim1D.cpp

@@ -344,8 +344,32 @@ void Sim1D::solve(int loglevel, bool refine_grid)
         clearDebugFile();
     }
 
+    int retries = 0;
     while (new_points > 0) {
-        SteadyStateSystem::solve(loglevel);
+        if (refine_grid == true) {
+            try {
+                SteadyStateSystem::solve(loglevel);
+            } catch (CanteraError& err) {
+                if (err.getMethod() == "SteadyStateSystem::timeStep" && retries < m_ts_regrid_max) {
+                    writelog("\nTime stepping failed; attempting to refine the grid and retry "
+                            "({}/{})...\n", retries+1, m_ts_regrid_max);
+                    new_points = refine(loglevel);
+                    if (new_points > 0) {
+                        retries++;
+                        continue;
+                    }
+                    if (loglevel > 0) {
+                        if (new_points == 0) {
+                            writelog("Regrid retry aborted: no new points added.\n");
+                        }
+                    }
+                }
+                throw;
+            }
+        } else {
+            SteadyStateSystem::solve(loglevel);
+        }
+
         if (loglevel > 0) {
             writelog("\nNewton steady-state solve succeeded.\n\n");
             writelog("Problem solved on [");