[1D] Add function for restoring last successful steady-state solution

include/cantera/oneD/Sim1D.h

@@ -141,6 +141,11 @@ public:
     //! integration.
     void restoreTimeSteppingSolution();
 
+    //! Set the current solution vector and grid to the last successful steady-
+    //! state solution. This can be used to examine the solver progress after a
+    //! failure during grid refinement.
+    void restoreSteadySolution();
+
     void getInitialSoln();
 
     void setSolution(const doublereal* soln) {
@@ -160,7 +165,15 @@ protected:
     vector_fp m_x;
 
     //! the solution vector after the last successful timestepping
-    vector_fp m_xlast;
+    vector_fp m_xlast_ts;
+
+    //! the solution vector after the last successful steady-state solve (stored
+    //! before grid refinement)
+    vector_fp m_xlast_ss;
+
+    //! the grids for each domain after the last successful steady-state solve
+    //! (stored before grid refinement)
+    std::vector<vector_fp> m_grid_last_ss;
 
     //! a work array used to hold the residual or the new solution
     vector_fp m_xnew;

interfaces/cython/cantera/_cantera.pxd

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

interfaces/cython/cantera/onedim.pyx

@@ -811,6 +811,14 @@ cdef class Sim1D:
         """
         self.sim.restoreTimeSteppingSolution()
 
+    def restore_steady_solution(self):
+        """
+        Set the current solution vector to the last successful steady-state
+        solution. This can be used to examine the solver progress after a
+        failure during grid refinement.
+        """
+        self.sim.restoreSteadySolution()
+
     def show_stats(self, print_time=True):
         """
         Show the statistics for the last solution.

src/oneD/Sim1D.cpp

@@ -21,7 +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_xlast_ts.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)]);
@@ -131,7 +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_xlast_ts.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);
@@ -171,7 +171,16 @@ void Sim1D::showSolution()
 
 void Sim1D::restoreTimeSteppingSolution()
 {
-    m_x = m_xlast;
+    m_x = m_xlast_ts;
+}
+
+void Sim1D::restoreSteadySolution()
+{
+    m_x = m_xlast_ss;
+    for (size_t n = 0; n < nDomains(); n++) {
+        vector_fp& z = m_grid_last_ss[n];
+        domain(n).setupGrid(z.size(), z.data());
+    }
 }
 
 void Sim1D::getInitialSoln()
@@ -268,7 +277,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;
+                m_xlast_ts = m_x;
                 if (loglevel > 6) {
                     save("debug_sim1d.xml", "debug", "After timestepping");
                 }
@@ -329,10 +338,16 @@ int Sim1D::refine(int loglevel)
     doublereal xmid, zmid;
     std::vector<size_t> dsize;
 
+    m_xlast_ss = m_x;
+    m_grid_last_ss.clear();
+
     for (size_t n = 0; n < nDomains(); n++) {
         Domain1D& d = domain(n);
         Refiner& r = d.refiner();
 
+        // Save the old grid corresponding to the converged solution
+        m_grid_last_ss.push_back(d.grid());
+
         // determine where new points are needed
         ianalyze = r.analyze(d.grid().size(), d.grid().data(), &m_x[start(n)]);
         if (ianalyze < 0) {