[1D] Add ability to access last successful timestepped solution

This can be helpful when debugging 1D solver failures.

include/cantera/oneD/Sim1D.h

@@ -136,6 +136,11 @@ public:
     //! Initialize the solution with a previously-saved solution.
     void restore(const std::string& fname, const std::string& id, int loglevel=2);
 
+    //! Set the current solution vector to the last successful time-stepping
+    //! solution. This can be used to examine the solver progress after a failed
+    //! integration.
+    void restoreTimeSteppingSolution();
+
     void getInitialSoln();
 
     void setSolution(const doublereal* soln) {
@@ -154,6 +159,9 @@ protected:
     //! the solution vector
     vector_fp m_x;
 
+    //! the solution vector after the last successful timestepping
+    vector_fp m_xlast;
+
     //! a work array used to hold the residual or the new solution
     vector_fp m_xnew;
 

interfaces/cython/cantera/_cantera.pxd

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

interfaces/cython/cantera/onedim.pyx

@@ -793,7 +793,16 @@ cdef class Sim1D:
         self.sim.restore(stringify(filename), stringify(name), loglevel)
         self._initialized = True
 
+    def restore_time_stepping_solution(self):
+        """
+        Set the current solution vector to the last successful time-stepping
+        solution. This can be used to examine the solver progress after a failed
+        integration.
+        """
+        self.sim.restoreTimeSteppingSolution()
+
     def show_stats(self, print_time=True):
+
         """
         Show the statistics for the last solution.
 

src/oneD/Sim1D.cpp

@@ -21,6 +21,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
     // resize the internal solution vector and the work array, and perform
     // domain-specific initialization of the solution vector.
     m_x.resize(size(), 0.0);
+    m_xlast.resize(size(), 0.0);
     m_xnew.resize(size(), 0.0);
     for (size_t n = 0; n < nDomains(); n++) {
         domain(n)._getInitialSoln(&m_x[start(n)]);
@@ -130,6 +131,7 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
         sz += domain(m).nComponents() * intValue((*xd[m])["points"]);
     }
     m_x.resize(sz);
+    m_xlast.resize(sz, 0.0);
     m_xnew.resize(sz);
     for (size_t m = 0; m < nDomains(); m++) {
         domain(m).restore(*xd[m], &m_x[domain(m).loc()], loglevel);
@@ -167,6 +169,11 @@ void Sim1D::showSolution()
     }
 }
 
+void Sim1D::restoreTimeSteppingSolution()
+{
+    m_x = m_xlast;
+}
+
 void Sim1D::getInitialSoln()
 {
     for (size_t n = 0; n < nDomains(); n++) {
@@ -261,6 +268,7 @@ void Sim1D::solve(int loglevel, bool refine_grid)
                 }
                 dt = timeStep(nsteps, dt, m_x.data(), m_xnew.data(),
                               loglevel-1);
+                m_xlast = m_x;
                 if (loglevel > 6) {
                     save("debug_sim1d.xml", "debug", "After timestepping");
                 }