[Reactor] Implement functions getState for reactors and reactor networks

These are simpler forms of getInitialConditions.

include/cantera/kinetics/ImplicitSurfChem.h

@@ -134,6 +134,8 @@ public:
 
     //! Set the initial conditions for the solution vector
     /*!
+     *  Essentially calls getState()
+     *
      *  @param t0  Initial time
      *  @param leny  Length of the solution vector
      *  @param y   Value of the solution vector to be used.
@@ -143,6 +145,14 @@ public:
     virtual void getInitialConditions(doublereal t0,
                                       size_t leny, doublereal* y);
 
+    //! Get the current state of the solution vector
+    /*!
+     *  @param y   Value of the solution vector to be used.
+     *            On output, this contains the initial value
+     *           of the solution.
+     */
+    virtual void getState(doublereal* y);
+
     /*!
      * Get the specifications for the problem from the values
      * in the ThermoPhase objects for all phases.

include/cantera/zeroD/ConstPressureReactor.h

@@ -31,6 +31,7 @@ public:
 
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
+    virtual void getState(doublereal* y);
 
     virtual void initialize(doublereal t0 = 0.0);
     virtual void evalEqs(doublereal t, doublereal* y,

include/cantera/zeroD/FlowReactor.h

@@ -26,6 +26,7 @@ public:
 
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
+    virtual void getState(doublereal* y);
 
     virtual void initialize(doublereal t0 = 0.0);
     virtual void evalEqs(doublereal t, doublereal* y,

include/cantera/zeroD/IdealGasConstPressureReactor.h

@@ -33,6 +33,7 @@ public:
 
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
+    virtual void getState(doublereal* y);
 
     virtual void initialize(doublereal t0 = 0.0);
     virtual void evalEqs(doublereal t, doublereal* y,

include/cantera/zeroD/IdealGasReactor.h

@@ -30,6 +30,7 @@ public:
 
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
+    virtual void getState(doublereal* y);
 
     virtual void initialize(doublereal t0 = 0.0);
 

include/cantera/zeroD/Reactor.h

@@ -95,6 +95,8 @@ public:
 
     //! Called by ReactorNet to get the initial conditions.
     /*!
+     *  Essentially calls function getState()
+     *
      *  @param[in] t0 Time at which initial conditions are determined
      *  @param[in] leny Length of *y* (unused)
      *  @param[out] y state vector representing the initial state of the reactor
@@ -102,6 +104,12 @@ public:
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
 
+    //! Get the the current state of the reactor.
+    /*!
+     *  @param[out] y state vector representing the initial state of the reactor
+     */
+    virtual void getState(doublereal* y);
+
     virtual void initialize(doublereal t0 = 0.0);
 
     /*!

include/cantera/zeroD/ReactorNet.h

@@ -195,6 +195,8 @@ public:
                       doublereal* ydot, doublereal* p);
     virtual void getInitialConditions(doublereal t0, size_t leny,
                                       doublereal* y);
+    virtual void getState(doublereal* y);
+
     virtual size_t nparams() {
         return m_ntotpar;
     }

interfaces/cython/cantera/_cantera.pxd

@@ -475,6 +475,8 @@ cdef extern from "cantera/zeroD/Reactor.h":
         void setEnergy(int)
         cbool energyEnabled()
         size_t componentIndex(string&)
+        size_t neq()
+        void getState(double*)
 
         void addSensitivityReaction(size_t) except +
         size_t nSensParams()
@@ -550,6 +552,7 @@ cdef extern from "cantera/zeroD/ReactorNet.h":
         cbool verbose()
         void setVerbose(cbool)
         size_t neq()
+        void getState(double*)
 
         void setSensitivityTolerances(double, double)
         double rtolSensitivity()

interfaces/cython/cantera/reactor.pyx

@@ -235,6 +235,48 @@ cdef class Reactor(ReactorBase):
             raise IndexError('No such component: {!r}'.format(name))
         return k
 
+    property n_vars:
+        """
+        The number of state variables in the reactor.
+        Equal to:
+
+        `Reactor` and `IdealGasReactor`: `n_species` + 3 (mass, volume,
+        internal energy or temperature).
+
+        `ConstPressureReactor` and `IdealGasConstPressureReactor`:
+        `n_species` + 2 (mass, enthalpy or temperature).
+        """
+        def __get__(self):
+            return self.reactor.neq()
+
+    def get_state(self):
+        """
+        Get the state vector of the reactor.
+
+        The order of the variables (i.e. rows) is:
+
+        `Reactor` or `IdealGasReactor`:
+
+          - 0  - mass
+          - 1  - volume
+          - 2  - internal energy or temperature
+          - 3+ - mass fractions of the species
+
+        `ConstPressureReactor` or `IdealGasConstPressureReactor`:
+
+          - 0  - mass
+          - 1  - enthalpy or temperature
+          - 2+ - mass fractions of the species
+
+        You can use the function `component_index` to determine the location
+        of a specific component
+        """
+        if not self.n_vars:
+            raise Exception('Reactor empty or network not initialized.')
+        cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars)
+        self.reactor.getState(&y[0])
+        return y
+
 
 cdef class Reservoir(ReactorBase):
     """
@@ -941,6 +983,19 @@ cdef class ReactorNet:
         def __get__(self):
             return self.net.neq()
 
+    def get_state(self):
+        """
+        Get the combined state vector of the reactor network.
+
+        The combined state vector consists of the concatenated state vectors of
+        all entities contained.
+        """
+        if not self.n_vars:
+            raise Exception('ReactorNet empty or not initialized.')
+        cdef np.ndarray[np.double_t, ndim=1] y = np.zeros(self.n_vars)
+        self.net.getState(&y[0])
+        return y
+
     def __reduce__(self):
         raise NotImplementedError('ReactorNet object is not picklable')
 

src/kinetics/ImplicitSurfChem.cpp

@@ -101,6 +101,11 @@ int ImplicitSurfChem::checkMatch(std::vector<ThermoPhase*> m_vec, ThermoPhase* t
 
 void ImplicitSurfChem::getInitialConditions(doublereal t0, size_t lenc,
         doublereal* c)
+{
+    getState(c);
+}
+
+void ImplicitSurfChem::getState(doublereal* c)
 {
     size_t loc = 0;
     for (size_t n = 0; n < m_nsurf; n++) {

src/zeroD/ConstPressureReactor.cpp

@@ -13,10 +13,16 @@ using namespace std;
 namespace Cantera
 {
 
-void ConstPressureReactor::getInitialConditions(double t0, size_t leny, double* y)
+void ConstPressureReactor::getInitialConditions(double t0, size_t leny,
+                                                double* y)
+{
+    getState(y);
+}
+
+void ConstPressureReactor::getState(double* y)
 {
     if (m_thermo == 0) {
-        throw CanteraError("getInitialConditions",
+        throw CanteraError("getState",
                            "Error: reactor is empty.");
     }
     m_thermo->restoreState(m_state);

src/zeroD/FlowReactor.cpp

@@ -25,10 +25,15 @@ FlowReactor::FlowReactor() :
 }
 
 void FlowReactor::getInitialConditions(double t0, size_t leny, double* y)
+{
+    getState(y);
+}
+
+void FlowReactor::getState(double* y)
 {
     if (m_thermo == 0) {
-        writelog("Error: reactor is empty.\n");
-        return;
+        throw CanteraError("getState",
+                           "Error: reactor is empty.");
     }
     m_thermo->restoreState(m_state);
     m_thermo->getMassFractions(y+2);

src/zeroD/IdealGasConstPressureReactor.cpp

@@ -27,9 +27,14 @@ void IdealGasConstPressureReactor::setThermoMgr(ThermoPhase& thermo)
 
 void IdealGasConstPressureReactor::getInitialConditions(double t0, size_t leny,
                                                         double* y)
+{
+    getState(y);
+}
+
+void IdealGasConstPressureReactor::getState(double* y)
 {
     if (m_thermo == 0) {
-        throw CanteraError("getInitialConditions",
+        throw CanteraError("getState",
                            "Error: reactor is empty.");
     }
     m_thermo->restoreState(m_state);

src/zeroD/IdealGasReactor.cpp

@@ -23,10 +23,15 @@ void IdealGasReactor::setThermoMgr(ThermoPhase& thermo)
 }
 
 void IdealGasReactor::getInitialConditions(double t0, size_t leny, double* y)
+{
+    getState(y);
+}
+
+void IdealGasReactor::getState(double* y)
 {
     if (m_thermo == 0) {
-        cout << "Error: reactor is empty." << endl;
-        return;
+        throw CanteraError("getState",
+                           "Error: reactor is empty.");
     }
     m_thermo->restoreState(m_state);
 

src/zeroD/Reactor.cpp

@@ -28,10 +28,15 @@ Reactor::Reactor() :
 {}
 
 void Reactor::getInitialConditions(double t0, size_t leny, double* y)
+{
+    getState(y);
+}
+
+void Reactor::getState(double* y)
 {
     if (m_thermo == 0) {
-        cout << "Error: reactor is empty." << endl;
-        return;
+        throw CanteraError("getState",
+                           "Error: reactor is empty.");
     }
     m_thermo->restoreState(m_state);
 

src/zeroD/ReactorNet.cpp

@@ -175,12 +175,15 @@ void ReactorNet::updateState(doublereal* y)
     }
 }
 
-void ReactorNet::getInitialConditions(doublereal t0,
-                                      size_t leny, doublereal* y)
+void ReactorNet::getInitialConditions(double t0, size_t leny, double* y)
+{
+    getState(y);
+}
+
+void ReactorNet::getState(double* y)
 {
     for (size_t n = 0; n < m_reactors.size(); n++) {
-        m_reactors[n]->getInitialConditions(t0, m_start[n+1]-m_start[n],
-                                            y + m_start[n]);
+        m_reactors[n]->getState(y + m_start[n]);
     }
 }