From 16ab7dff5185f999fd58f83e76e16b53f0dad29d Mon Sep 17 00:00:00 2001 From: Thomas Fiala Date: Mon, 2 Nov 2015 11:31:16 +0100 Subject: [PATCH] [Reactor] Implement functions getState for reactors and reactor networks These are simpler forms of getInitialConditions. --- include/cantera/kinetics/ImplicitSurfChem.h | 10 ++++ include/cantera/zeroD/ConstPressureReactor.h | 1 + include/cantera/zeroD/FlowReactor.h | 1 + .../zeroD/IdealGasConstPressureReactor.h | 1 + include/cantera/zeroD/IdealGasReactor.h | 1 + include/cantera/zeroD/Reactor.h | 8 +++ include/cantera/zeroD/ReactorNet.h | 2 + interfaces/cython/cantera/_cantera.pxd | 3 + interfaces/cython/cantera/reactor.pyx | 55 +++++++++++++++++++ src/kinetics/ImplicitSurfChem.cpp | 5 ++ src/zeroD/ConstPressureReactor.cpp | 10 +++- src/zeroD/FlowReactor.cpp | 9 ++- src/zeroD/IdealGasConstPressureReactor.cpp | 7 ++- src/zeroD/IdealGasReactor.cpp | 9 ++- src/zeroD/Reactor.cpp | 9 ++- src/zeroD/ReactorNet.cpp | 11 ++-- 16 files changed, 129 insertions(+), 13 deletions(-) diff --git a/include/cantera/kinetics/ImplicitSurfChem.h b/include/cantera/kinetics/ImplicitSurfChem.h index 1447f6479..2b28f71ca 100644 --- a/include/cantera/kinetics/ImplicitSurfChem.h +++ b/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. diff --git a/include/cantera/zeroD/ConstPressureReactor.h b/include/cantera/zeroD/ConstPressureReactor.h index c5b5c86ec..5582b4d37 100644 --- a/include/cantera/zeroD/ConstPressureReactor.h +++ b/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, diff --git a/include/cantera/zeroD/FlowReactor.h b/include/cantera/zeroD/FlowReactor.h index fa1f05537..3d059a2b7 100644 --- a/include/cantera/zeroD/FlowReactor.h +++ b/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, diff --git a/include/cantera/zeroD/IdealGasConstPressureReactor.h b/include/cantera/zeroD/IdealGasConstPressureReactor.h index 2546b8f4b..f69aceb7b 100644 --- a/include/cantera/zeroD/IdealGasConstPressureReactor.h +++ b/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, diff --git a/include/cantera/zeroD/IdealGasReactor.h b/include/cantera/zeroD/IdealGasReactor.h index ada758cda..9f2304ffd 100644 --- a/include/cantera/zeroD/IdealGasReactor.h +++ b/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); diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 7eb2c0058..124bb2e6f 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/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); /*! diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 011fad489..b7616a8b8 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/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; } diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 69e709026..404196be9 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/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() diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 323254ab5..507f1ed3c 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/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') diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp index ed1685af8..773a657e5 100644 --- a/src/kinetics/ImplicitSurfChem.cpp +++ b/src/kinetics/ImplicitSurfChem.cpp @@ -101,6 +101,11 @@ int ImplicitSurfChem::checkMatch(std::vector 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++) { diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index cd237114a..2eef4c684 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/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); diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp index 69790fe50..9c6a180da 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/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); diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 2454c4150..030d3c60d 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/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); diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 88852372b..974f42351 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/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); diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 711aa6709..1f2927afd 100644 --- a/src/zeroD/Reactor.cpp +++ b/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); diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 9cfdc139d..02aba4f03 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/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]); } }