[1D] Allow replacing flow domain grid after creating Sim1D

This makes it possible to try solving on a different grid if the initial
solution attempt fails without needing to recreate the Sim1D object.

include/cantera/oneD/OneDim.h

@@ -180,8 +180,8 @@ public:
      */
     Domain1D* pointDomain(size_t i);
 
-    //! Call after one or more grids has been refined.
-    void resize();
+    //! Call after one or more grids has changed size, e.g. after being refined.
+    virtual void resize();
 
     vector_int& transientMask() {
         return m_mask;

include/cantera/oneD/Sim1D.h

@@ -160,6 +160,8 @@ public:
 
     void evalSSJacobian();
 
+    virtual void resize();
+
 protected:
     //! the solution vector
     vector_fp m_x;

interfaces/cython/cantera/_cantera.pxd

@@ -695,6 +695,7 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void restore(string, string, int) except +
         void writeStats(int) except +
         void clearStats()
+        void resize() except +
         vector[size_t]& gridSizeStats()
         vector[double]& jacobianTimeStats()
         vector[double]& evalTimeStats()

interfaces/cython/cantera/onedim.pyx

@@ -736,6 +736,7 @@ cdef class Sim1D:
         Load the initial solution from each domain into the global solution
         vector.
         """
+        self.sim.resize()
         self.sim.getInitialSoln()
 
     def solve(self, loglevel=1, refine_grid=True, auto=False):

interfaces/cython/cantera/test/test_onedim.py

@@ -424,6 +424,21 @@ class TestFreeFlame(utilities.CanteraTest):
                 vals[i] = bad[i]
                 self.sim.set_refine_criteria(*vals)
 
+    def test_replace_grid(self):
+        self.create_sim(ct.one_atm, 300.0, 'H2:1.1, O2:1, AR:5')
+        self.solve_fixed_T()
+        ub = self.sim.u[-1]
+
+        # add some points to the grid
+        grid = list(self.sim.grid)
+        for i in (7,5,4,2):
+            grid.insert(i, 0.5 * (grid[i-1] + grid[i]))
+        self.sim.flame.grid = grid
+        self.sim.set_initial_guess()
+
+        self.solve_fixed_T()
+        self.assertNear(self.sim.u[-1], ub, 1e-2)
+
 
 class TestDiffusionFlame(utilities.CanteraTest):
     # Note: to re-create the reference file:

src/oneD/Sim1D.cpp

@@ -21,8 +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_xnew.resize(size(), 0.0);
+    resize();
     for (size_t n = 0; n < nDomains(); n++) {
         domain(n)._getInitialSoln(&m_x[start(n)]);
     }
@@ -122,21 +121,19 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
             " correct number of domains. Found {} expected {}.\n",
             xd.size(), nDomains());
     }
-    size_t sz = 0;
     for (size_t m = 0; m < nDomains(); m++) {
-        if (loglevel > 0 && xd[m]->attrib("id") != domain(m).id()) {
+        Domain1D& dom = domain(m);
+        if (loglevel > 0 && xd[m]->attrib("id") != dom.id()) {
             writelog("Warning: domain names do not match: '" +
-                     (*xd[m])["id"] + + "' and '" + domain(m).id() + "'\n");
+                     (*xd[m])["id"] + + "' and '" + dom.id() + "'\n");
         }
-        sz += domain(m).nComponents() * intValue((*xd[m])["points"]);
+        dom.resize(domain(m).nComponents(), intValue((*xd[m])["points"]));
     }
-    m_x.resize(sz);
+    resize();
     m_xlast_ts.clear();
-    m_xnew.resize(sz);
     for (size_t m = 0; m < nDomains(); m++) {
         domain(m).restore(*xd[m], &m_x[domain(m).loc()], loglevel);
     }
-    resize();
     finalize();
 }
 
@@ -424,10 +421,6 @@ int Sim1D::refine(int loglevel)
 
     // Replace the current solution vector with the new one
     m_x = xnew;
-
-    // resize the work array
-    m_xnew.resize(xnew.size());
-
     resize();
     finalize();
     return np;
@@ -517,8 +510,6 @@ int Sim1D::setFixedTemperature(doublereal t)
     // Replace the current solution vector with the new one
     m_x = xnew;
 
-    // resize the work array
-    m_xnew = xnew;
     resize();
     finalize();
     return np;
@@ -573,4 +564,12 @@ void Sim1D::evalSSJacobian()
 {
     OneDim::evalSSJacobian(m_x.data(), m_xnew.data());
 }
+
+void Sim1D::resize()
+{
+    OneDim::resize();
+    m_x.resize(size(), 0.0);
+    m_xnew.resize(size(), 0.0);
+}
+
 }