From 8b5ca80303dfd3e0d28f7f34a2ffb27303a5e922 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 10 Jun 2014 16:23:03 +0000 Subject: [PATCH] [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 | 2 ++ include/cantera/zeroD/ReactorBase.h | 5 +++++ include/cantera/zeroD/ReactorNet.h | 15 ++++++++++++- interfaces/cython/cantera/_cantera.pxd | 2 ++ interfaces/cython/cantera/reactor.pyx | 16 ++++++++++++++ .../cython/cantera/test/test_reactor.py | 22 +++++++++++++++++++ src/zeroD/Reactor.cpp | 6 +++++ src/zeroD/ReactorBase.cpp | 12 ++++++++++ src/zeroD/ReactorNet.cpp | 18 ++++++++++++++- 9 files changed, 96 insertions(+), 2 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 94dde02c5..6c31d5de2 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/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); diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 9fbe77c9b..650127015 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/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; diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 40186d87f..5a3d8b86e 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/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 diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 00a68cbf9..b7cdf9ce8 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/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) diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index cf4324a47..437674f15 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/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): diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index 3a91e4b6c..86b7b2c7c 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/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 diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index ef1a4357f..3af7ffc4f 100644 --- a/src/zeroD/Reactor.cpp +++ b/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++) { diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index 47dc0f06b..9c60e43c7 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/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); diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 49c31da1f..769dd2b86 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/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());