[1D] Set limit on number of timesteps without reaching steady state

include/cantera/oneD/OneDim.h

@@ -224,6 +224,19 @@ public:
     void setTimeStepFactor(doublereal tfactor) {
         m_tfactor = tfactor;
     }
+
+    //! Set the maximum number of timeteps allowed before successful
+    //! steady-state solve
+    void setMaxTimeStepCount(int nmax) {
+        m_nsteps_max = nmax;
+    }
+
+    //! Return the maximum number of timeteps allowed before successful
+    //! steady-state solve
+    int maxTimeStepCount() const {
+        return m_nsteps_max;
+    }
+
     void setJacAge(int ss_age, int ts_age=-1) {
         m_ss_jac_age = ss_age;
         if (ts_age > 0) {
@@ -340,6 +353,9 @@ private:
     //! Number of time steps taken in the current call to solve()
     int m_nsteps;
 
+    //! Maximum number of timesteps allowed per call to solve()
+    int m_nsteps_max;
+
     //! Number of time steps taken in each call to solve() (e.g. for each
     //! successive grid refinement)
     vector_int m_timeSteps;

interfaces/cython/cantera/_cantera.pxd

@@ -679,6 +679,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void showSolution() except +
         void setTimeStep(double, size_t, int*) except +
         void restoreTimeSteppingSolution() except +
+        void setMaxTimeStepCount(int)
+        int maxTimeStepCount()
         void getInitialSoln() except +
         void solve(int, cbool) except +translate_exception
         void refine(int) except +

interfaces/cython/cantera/onedim.pyx

@@ -647,6 +647,16 @@ cdef class Sim1D:
             data.push_back(n)
         self.sim.setTimeStep(stepsize, data.size(), &data[0])
 
+    property max_time_step_count:
+        """
+        Get/Set the maximum number of time steps allowed before reaching the
+        steady-state solution
+        """
+        def __get__(self):
+            return self.sim.maxTimeStepCount()
+        def __set__(self, nmax):
+            self.sim.setMaxTimeStepCount(nmax)
+
     def set_initial_guess(self):
         """
         Set the initial guess for the solution. Derived classes extend this

src/oneD/OneDim.cpp

@@ -18,7 +18,8 @@ OneDim::OneDim()
       m_bw(0), m_size(0),
       m_init(false), m_pts(0), m_solve_time(0.0),
       m_ss_jac_age(10), m_ts_jac_age(20),
-      m_interrupt(0), m_nevals(0), m_evaltime(0.0), m_nsteps(0)
+      m_interrupt(0), m_nevals(0), m_evaltime(0.0), m_nsteps(0),
+      m_nsteps_max(500)
 {
     m_newt.reset(new MultiNewton(1));
 }
@@ -29,7 +30,8 @@ OneDim::OneDim(vector<Domain1D*> domains) :
     m_bw(0), m_size(0),
     m_init(false), m_solve_time(0.0),
     m_ss_jac_age(10), m_ts_jac_age(20),
-    m_interrupt(0), m_nevals(0), m_evaltime(0.0), m_nsteps(0)
+    m_interrupt(0), m_nevals(0), m_evaltime(0.0), m_nsteps(0),
+    m_nsteps_max(500)
 {
     // create a Newton iterator, and add each domain.
     m_newt.reset(new MultiNewton(1));
@@ -358,6 +360,11 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
                 dt *= 1.5;
             }
             dt = std::min(dt, m_tmax);
+            if (m_nsteps == m_nsteps_max) {
+                throw CanteraError("OneDim::timeStep",
+                    "Took maximum number of timesteps allowed ({}) without "
+                    "reaching steady-state solution.", m_nsteps_max);
+            }
         } else {
             successiveFailures++;
             // No solution could be found with this time step.