From 0e9ca0fc4b9004da6863458047b850ae3a750ac0 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 3 Jun 2014 16:52:00 +0000 Subject: [PATCH] [Reactor] Eliminate tracking of Reservoirs by ReactorNet Since Reservoirs don't have any equations to solve, there is no reason to add them to a ReactorNet. --- include/cantera/zeroD/ReactorNet.h | 17 ++--- interfaces/cython/cantera/_cantera.pxd | 2 +- interfaces/cython/cantera/reactor.pyx | 4 +- src/clib/ctreactor.cpp | 3 +- src/zeroD/ReactorNet.cpp | 89 ++++++++++++-------------- 5 files changed, 52 insertions(+), 63 deletions(-) diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 965e63fe2..92ce111d2 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -114,13 +114,13 @@ public: //@} //! Add the reactor *r* to this reactor network. - void addReactor(ReactorBase* r, bool iown = false); + void addReactor(Reactor* r, bool iown = false); //! Return a reference to the *n*-th reactor in this network. The reactor //! indices are determined by the order in which the reactors were added //! to the reactor network. - ReactorBase& reactor(int n) { - return *m_r[n]; + Reactor& reactor(int n) { + return *m_reactors[n]; } //! Returns `true` if verbose logging output is enabled. @@ -198,26 +198,23 @@ public: return m_paramNames.at(p); } +protected: void connect(size_t i, size_t j) { - m_connect[j*m_nr + i] = 1; - m_connect[i*m_nr + j] = 1; + m_connect[j*m_reactors.size() + i] = 1; + m_connect[i*m_reactors.size() + j] = 1; } bool connected(size_t i, size_t j) { - return (m_connect[m_nr*i + j] == 1); + return (m_connect[m_reactors.size()*i + j] == 1); } -protected: /** * Initialize the reactor network. Called automatically the first time * advance or step is called. */ void initialize(); - std::vector m_r; std::vector m_reactors; - size_t m_nr; - size_t m_nreactors; Integrator* m_integ; doublereal m_time; bool m_init; diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 932a67eca..d58b6616c 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -324,7 +324,7 @@ cdef extern from "cantera/zeroD/flowControllers.h": cdef extern from "cantera/zeroD/ReactorNet.h": cdef cppclass CxxReactorNet "Cantera::ReactorNet": CxxReactorNet() - void addReactor(CxxReactorBase*) + void addReactor(CxxReactor*) void advance(double) except + double step(double) except + double time() diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index e95d0af21..e5a21b0aa 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -729,10 +729,10 @@ cdef class ReactorNet: for R in reactors: self.add_reactor(R) - def add_reactor(self, ReactorBase r): + def add_reactor(self, Reactor r): """Add a reactor to the network.""" self._reactors.append(r) - self.net.addReactor(r.rbase) + self.net.addReactor(r.reactor) def advance(self, double t): """ diff --git a/src/clib/ctreactor.cpp b/src/clib/ctreactor.cpp index d0d49803f..05ecb318a 100644 --- a/src/clib/ctreactor.cpp +++ b/src/clib/ctreactor.cpp @@ -303,7 +303,8 @@ extern "C" { int reactornet_addreactor(int i, int n) { try { - NetworkCabinet::item(i).addReactor(&ReactorCabinet::item(n)); + NetworkCabinet::item(i).addReactor( + &dynamic_cast(ReactorCabinet::item(n))); return 0; } catch (...) { return handleAllExceptions(-1, ERR); diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index d68e6c7db..5202cd6a5 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -11,7 +11,7 @@ using namespace std; namespace Cantera { -ReactorNet::ReactorNet() : Cantera::FuncEval(), m_nr(0), m_nreactors(0), +ReactorNet::ReactorNet() : Cantera::FuncEval(), m_integ(0), m_time(0.0), m_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), @@ -33,14 +33,12 @@ ReactorNet::ReactorNet() : Cantera::FuncEval(), m_nr(0), m_nreactors(0), ReactorNet::~ReactorNet() { - for (size_t n = 0; n < m_nr; n++) { + for (size_t n = 0; n < m_reactors.size(); n++) { if (m_iown[n]) { - delete m_r[n]; + delete m_reactors[n]; } - m_r[n] = 0; + m_reactors[n] = 0; } - m_r.clear(); - m_reactors.clear(); delete m_integ; } @@ -49,57 +47,51 @@ void ReactorNet::initialize() size_t n, nv; char buf[100]; m_nv = 0; - m_reactors.clear(); - m_nreactors = 0; writelog("Initializing reactor network.\n", m_verbose); - if (m_nr == 0) + if (m_reactors.empty()) throw CanteraError("ReactorNet::initialize", "no reactors in network!"); size_t sensParamNumber = 0; m_start.assign(1, 0); - for (n = 0; n < m_nr; n++) { - if (m_r[n]->type() >= ReactorType) { - m_r[n]->initialize(m_time); - Reactor* r = (Reactor*)m_r[n]; - m_reactors.push_back(r); - nv = r->neq(); - m_nparams.push_back(r->nSensParams()); - std::vector > sens_objs = r->getSensitivityOrder(); - for (size_t i = 0; i < sens_objs.size(); i++) { - std::map& s = m_sensOrder[sens_objs[i]]; - for (std::map::iterator iter = s.begin(); - iter != s.end(); - ++iter) { - m_sensIndex.resize(std::max(iter->second + 1, m_sensIndex.size())); - m_sensIndex[iter->second] = sensParamNumber++; - } + for (n = 0; n < m_reactors.size(); n++) { + Reactor& r = *m_reactors[n]; + r.initialize(m_time); + nv = r.neq(); + m_nparams.push_back(r.nSensParams()); + std::vector > sens_objs = r.getSensitivityOrder(); + for (size_t i = 0; i < sens_objs.size(); i++) { + std::map& s = m_sensOrder[sens_objs[i]]; + for (std::map::iterator iter = s.begin(); + iter != s.end(); + ++iter) { + m_sensIndex.resize(std::max(iter->second + 1, m_sensIndex.size())); + m_sensIndex[iter->second] = sensParamNumber++; } - m_nv += nv; - m_start.push_back(m_nv); - m_nreactors++; + } + m_nv += nv; + m_start.push_back(m_nv); - if (m_verbose) { - sprintf(buf,"Reactor %s: %s variables.\n", - int2str(n).c_str(), int2str(nv).c_str()); - writelog(buf); - sprintf(buf," %s sensitivity params.\n", - int2str(r->nSensParams()).c_str()); - writelog(buf); - } - if (m_r[n]->type() == FlowReactorType && m_nr > 1) { - throw CanteraError("ReactorNet::initialize", - "FlowReactors must be used alone."); - } + if (m_verbose) { + sprintf(buf,"Reactor %s: %s variables.\n", + int2str(n).c_str(), int2str(nv).c_str()); + writelog(buf); + sprintf(buf," %s sensitivity params.\n", + int2str(r.nSensParams()).c_str()); + writelog(buf); + } + if (r.type() == FlowReactorType && m_reactors.size() > 1) { + throw CanteraError("ReactorNet::initialize", + "FlowReactors must be used alone."); } } - m_connect.resize(m_nr*m_nr,0); + m_connect.resize(m_reactors.size()*m_reactors.size(), 0); m_ydot.resize(m_nv,0.0); size_t i, j, nin, nout, nw; ReactorBase* r, *rj; - for (i = 0; i < m_nr; i++) { + for (i = 0; i < m_reactors.size(); i++) { r = m_reactors[i]; - for (j = 0; j < m_nr; j++) { + for (j = 0; j < m_reactors.size(); j++) { if (i == j) { connect(i,j); } else { @@ -172,13 +164,12 @@ double ReactorNet::step(doublereal time) return m_time; } -void ReactorNet::addReactor(ReactorBase* r, bool iown) +void ReactorNet::addReactor(Reactor* r, bool iown) { r->setNetwork(this); if (r->type() >= ReactorType) { - m_r.push_back(r); + m_reactors.push_back(r); m_iown.push_back(iown); - m_nr++; writelog("Adding reactor "+r->name()+"\n", m_verbose); } else { writelog("Not adding reactor "+r->name()+ @@ -193,7 +184,7 @@ void ReactorNet::eval(doublereal t, doublereal* y, size_t pstart = 0; updateState(y); - for (n = 0; n < m_nreactors; n++) { + for (n = 0; n < m_reactors.size(); n++) { m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p + pstart); pstart += m_nparams[n]; @@ -229,7 +220,7 @@ void ReactorNet::evalJacobian(doublereal t, doublereal* y, void ReactorNet::updateState(doublereal* y) { - for (size_t n = 0; n < m_nreactors; n++) { + for (size_t n = 0; n < m_reactors.size(); n++) { m_reactors[n]->updateState(y + m_start[n]); } } @@ -237,7 +228,7 @@ void ReactorNet::updateState(doublereal* y) void ReactorNet::getInitialConditions(doublereal t0, size_t leny, doublereal* y) { - for (size_t n = 0; n < m_nreactors; n++) { + 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]); }