From 318fa04cf78863a150793b91928072ce9652357f Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sun, 21 Feb 2016 22:40:14 -0500
Subject: [PATCH] [1D] Reset negative species concentrations after a failed
 time step

---
 include/cantera/oneD/Domain1D.h |  7 +++++++
 include/cantera/oneD/OneDim.h   |  2 ++
 include/cantera/oneD/StFlow.h   |  3 +++
 src/oneD/OneDim.cpp             | 25 +++++++++++++++++++++----
 src/oneD/StFlow.cpp             | 10 ++++++++++
 5 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h
index 5dd93e06f..1846611d1 100644
--- a/include/cantera/oneD/Domain1D.h
+++ b/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
diff --git a/include/cantera/oneD/OneDim.h b/include/cantera/oneD/OneDim.h
index 95068972d..016fb4038 100644
--- a/include/cantera/oneD/OneDim.h
+++ b/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
     /*!
diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h
index 53ac8e320..9509d842b 100644
--- a/include/cantera/oneD/StFlow.h
+++ b/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;
     }
diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp
index ec19fd4a4..58dd015a9 100644
--- a/src/oneD/OneDim.cpp
+++ b/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)
diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
index f2e8151dd..77209488d 100644
--- a/src/oneD/StFlow.cpp
+++ b/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;