[Reactor] Implement reactor network reinitialization

Adds ReactorNet::reinitialize, which skips all one-time initialization and
re-uses the same CVODES integrator. The Reactor::syncState() method is
introduced for applying new initial conditions for individual Reactor objects.

This approach increases efficiency when solving many similar problems with short
integration times, for example when being used as the chemistry term integrator
in an operator-split CFD code.

include/cantera/zeroD/Reactor.h

@@ -115,6 +115,8 @@ public:
     virtual void evalEqs(doublereal t, doublereal* y,
                          doublereal* ydot, doublereal* params);
 
+    virtual void syncState();
+
     //! Set the state of the reactor to correspond to the state vector *y*.
     virtual void updateState(doublereal* y);
 

include/cantera/zeroD/ReactorBase.h

@@ -130,6 +130,11 @@ public:
         m_thermo->restoreState(m_state);
     }
 
+    //! Set the state of the reactor to correspond to the state of the
+    //! associated ThermoPhase object. This is the inverse of restoreState().
+    //! Calling this will trigger integrator reinitialization.
+    virtual void syncState();
+
     //! return a reference to the contents.
     thermo_t& contents() {
         return *m_thermo;

include/cantera/zeroD/ReactorNet.h

@@ -36,7 +36,7 @@ public:
      */
     void setInitialTime(doublereal time) {
         m_time = time;
-        m_init = false;
+        m_integrator_init = false;
     }
 
     //! Set the maximum time step.
@@ -205,6 +205,18 @@ public:
         return m_paramNames.at(p);
     }
 
+    //! Reinitialize the integrator. Used to solve a new problem (different
+    //! initial conditions) but with the same configuration of the reactor
+    //! network. Can be called manually, or automatically after calling
+    //! setInitialTime or modifying a reactor's contents.
+    void reinitialize();
+
+    //! Called to trigger integrator reinitialization before further
+    //! integration.
+    void setNeedsReinit() {
+        m_integrator_init = false;
+    }
+
 protected:
     void connect(size_t i, size_t j) {
         m_connect[j*m_reactors.size() + i] = 1;
@@ -225,6 +237,7 @@ protected:
     Integrator* m_integ;
     doublereal m_time;
     bool m_init;
+    bool m_integrator_init; //! True if integrator initialization is current
     size_t m_nv;
 
     //! m_start[n] is the starting point in the state vector for reactor n

interfaces/cython/cantera/_cantera.pxd

@@ -251,6 +251,7 @@ cdef extern from "cantera/zeroD/ReactorBase.h" namespace "Cantera":
         CxxReactorBase()
         void setThermoMgr(CxxThermoPhase&) except +
         void restoreState() except +
+        void syncState() except +
         double volume()
         string name()
         void setName(string)
@@ -327,6 +328,7 @@ cdef extern from "cantera/zeroD/ReactorNet.h":
         void addReactor(CxxReactor&)
         void advance(double) except +
         double step(double) except +
+        void reinitialize() except +
         double time()
         void setInitialTime(double)
         void setTolerances(double, double)

interfaces/cython/cantera/reactor.pyx

@@ -44,6 +44,14 @@ cdef class ReactorBase:
         def __set__(self, name):
             self.rbase.setName(stringify(name))
 
+    def syncState(self):
+        """
+        Set the state of the Reactor to match that of the associated
+        `ThermoPhase` object. After calling syncState(), call
+        ReactorNet.reinitialize() before further integration.
+        """
+        self.rbase.syncState()
+
     property thermo:
         """The `ThermoPhase` object representing the reactor's contents."""
         def __get__(self):
@@ -748,6 +756,14 @@ cdef class ReactorNet:
         """
         return self.net.step(t)
 
+    def reinitialize(self):
+        """
+        Reinitialize the integrator after making changing to the state of the
+        system. Changes to Reactor contents will automatically trigger
+        reinitialization.
+        """
+        self.net.reinitialize()
+
     property time:
         """The current time [s]."""
         def __get__(self):

interfaces/cython/cantera/test/test_reactor.py

@@ -474,6 +474,28 @@ class TestReactor(utilities.CanteraTest):
         self.assertNear(p1a, p1b)
         self.assertNear(p2a, p2b)
 
+    def test_reinitialize(self):
+        self.make_reactors(T1=300, T2=1000, independent=False)
+        self.add_wall(U=200, A=1.0)
+        self.net.advance(1.0)
+        T1a = self.r1.T
+        T2a = self.r2.T
+
+        self.r1.thermo.TD = 300, None
+        self.r1.syncState()
+
+        self.r2.thermo.TD = 1000, None
+        self.r2.syncState()
+
+        self.assertNear(self.r1.T, 300)
+        self.assertNear(self.r2.T, 1000)
+        self.net.advance(2.0)
+        T1b = self.r1.T
+        T2b = self.r2.T
+
+        self.assertNear(T1a, T1b)
+        self.assertNear(T2a, T2b)
+
     def test_unpicklable(self):
         self.make_reactors()
         import pickle

src/zeroD/Reactor.cpp

@@ -116,6 +116,12 @@ size_t Reactor::nSensParams()
     return m_nsens;
 }
 
+void Reactor::syncState()
+{
+    ReactorBase::syncState();
+    m_mass = m_thermo->density() * m_vol;
+}
+
 void Reactor::updateState(doublereal* y)
 {
     for (size_t i = 0; i < m_nv; i++) {

src/zeroD/ReactorBase.cpp

@@ -7,6 +7,7 @@
 #include "cantera/zeroD/ReactorBase.h"
 #include "cantera/zeroD/FlowDevice.h"
 #include "cantera/zeroD/Wall.h"
+#include "cantera/zeroD/ReactorNet.h"
 
 using namespace std;
 namespace Cantera
@@ -39,6 +40,17 @@ void ReactorBase::setThermoMgr(thermo_t& thermo)
     m_pressure = m_thermo->pressure();
 }
 
+void ReactorBase::syncState()
+{
+    m_thermo->saveState(m_state);
+    m_enthalpy = m_thermo->enthalpy_mass();
+    m_intEnergy = m_thermo->intEnergy_mass();
+    m_pressure = m_thermo->pressure();
+    if (m_net) {
+        m_net->setNeedsReinit();
+    }
+}
+
 void ReactorBase::addInlet(FlowDevice& inlet)
 {
     m_inlet.push_back(&inlet);

src/zeroD/ReactorNet.cpp

@@ -12,7 +12,7 @@ namespace Cantera
 {
 
 ReactorNet::ReactorNet() : Cantera::FuncEval(),
-    m_integ(0), m_time(0.0), m_init(false),
+    m_integ(0), m_time(0.0), m_init(false), m_integrator_init(false),
     m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4),
     m_atols(1.0e-15), m_atolsens(1.0e-4),
     m_maxstep(-1.0), m_maxErrTestFails(0),
@@ -132,9 +132,21 @@ void ReactorNet::initialize()
         writelog(buf);
     }
     m_integ->initialize(m_time, *this);
+    m_integrator_init = true;
     m_init = true;
 }
 
+void ReactorNet::reinitialize()
+{
+    if (m_init) {
+        writelog("Re-initializing reactor network.\n", m_verbose);
+        m_integ->reinitialize(m_time, *this);
+        m_integrator_init = true;
+    } else {
+        initialize();
+    }
+}
+
 void ReactorNet::advance(doublereal time)
 {
     if (!m_init) {
@@ -142,6 +154,8 @@ void ReactorNet::advance(doublereal time)
             m_maxstep = time - m_time;
         }
         initialize();
+    } else if (!m_integrator_init) {
+        reinitialize();
     }
     m_integ->integrate(time);
     m_time = time;
@@ -155,6 +169,8 @@ double ReactorNet::step(doublereal time)
             m_maxstep = time - m_time;
         }
         initialize();
+    } else if (!m_integrator_init) {
+        reinitialize();
     }
     m_time = m_integ->step(time);
     updateState(m_integ->solution());