[1D] Reset negative species concentrations after a failed time step

include/cantera/oneD/Domain1D.h

@@ -122,6 +122,13 @@ public:
     virtual void setInitialState(doublereal* xlocal = 0) {}
     virtual void setState(size_t point, const doublereal* state, doublereal* x) {}
 
+    /*!
+     * When called, this function should reset "bad" values in the state vector
+     * such as negative species concentrations. This function may be called
+     * after a failed solution attempt.
+     */
+    virtual void resetBadValues(double* xg) {}
+
     /*!
      * Resize the domain to have nv components and np grid points. This method
      * is virtual so that subclasses can perform other actions required to

include/cantera/oneD/OneDim.h

@@ -200,6 +200,8 @@ public:
     double timeStep(int nsteps, double dt, double* x,
                     double* r, int loglevel);
 
+    void resetBadValues(double* x);
+
     //! Write statistics about the number of iterations and Jacobians at each
     //! grid level
     /*!

include/cantera/oneD/StFlow.h

@@ -55,6 +55,9 @@ public:
 
     virtual void setupGrid(size_t n, const doublereal* z);
 
+    virtual void resetBadValues(double* xg);
+
+
     thermo_t& phase() {
         return *m_thermo;
     }

src/oneD/OneDim.cpp

@@ -328,6 +328,8 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
     debuglog("===============================\n", loglevel);
 
     int n = 0;
+    int successiveFailures = 0;
+
     while (n < nsteps) {
         if (loglevel > 0) {
             doublereal ss = ssnorm(x, r);
@@ -343,6 +345,7 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
         // successful time step. Copy the new solution in r to
         // the current solution in x.
         if (m >= 0) {
+            successiveFailures = 0;
             n += 1;
             debuglog("\n", loglevel);
             copy(r, r + m_size, x);
@@ -351,13 +354,20 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
             }
             dt = std::min(dt, m_tmax);
         } else {
+            successiveFailures++;
             // No solution could be found with this time step.
             // Decrease the stepsize and try again.
             debuglog("...failure.\n", loglevel);
-            dt *= m_tfactor;
-            if (dt < m_tmin) {
-                throw CanteraError("OneDim::timeStep",
-                                   "Time integration failed.");
+            if (successiveFailures > 2) {
+                //debuglog("Resetting negative species concentrations.\n", loglevel);
+                resetBadValues(x);
+                successiveFailures = 0;
+            } else {
+                dt *= m_tfactor;
+                if (dt < m_tmin) {
+                    throw CanteraError("OneDim::timeStep",
+                                       "Time integration failed.");
+                }
             }
         }
     }
@@ -371,6 +381,13 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x,
     return dt;
 }
 
+void OneDim::resetBadValues(double* x)
+{
+    for (auto dom : m_dom) {
+        dom->resetBadValues(x);
+    }
+}
+
 void OneDim::save(const std::string& fname, std::string id,
                   const std::string& desc, doublereal* sol,
                   int loglevel)

src/oneD/StFlow.cpp

@@ -139,6 +139,16 @@ void StFlow::setupGrid(size_t n, const doublereal* z)
     }
 }
 
+void StFlow::resetBadValues(double* xg) {
+    double* x = xg + loc();
+    for (size_t j = 0; j < m_points; j++) {
+        double* Y = x + m_nv*j + c_offset_Y;
+        m_thermo->setMassFractions(Y);
+        m_thermo->getMassFractions(Y);
+    }
+}
+
+
 void StFlow::setTransport(Transport& trans, bool withSoret)
 {
     m_trans = &trans;