From fe6b5d3c0be7ffc8a003ee9a891677f45a9a3c49 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Fri, 2 Nov 2012 20:07:25 +0000 Subject: [PATCH] Changed the ordering of Reactor sensitivity coefficients The order now matches the order in which the corresponding sensitivity reactions are added to the ReactorNet, regardless of the order in which Reactors and Walls are added to the network. Sensitivity parameter names can be accessed using the "sensitivityParameterName" method of ReactorNet, and the "sensParamID" methods of Reactor and Wall have been removed as they no longer meaningful. --- include/cantera/zeroD/Reactor.h | 6 +----- include/cantera/zeroD/ReactorBase.h | 10 ++++++++++ include/cantera/zeroD/ReactorNet.h | 26 ++++++++++++++++++++++++- include/cantera/zeroD/Wall.h | 11 +++-------- interfaces/cython/cantera/_cantera.pxd | 3 +-- interfaces/cython/cantera/reactor.pyx | 9 +++------ src/zeroD/Reactor.cpp | 19 +++++++++++++++++- src/zeroD/ReactorBase.cpp | 18 ++++++++++++++++- src/zeroD/ReactorNet.cpp | 27 +++++++++++++++++++++++++- src/zeroD/Wall.cpp | 13 +++++++++++-- 10 files changed, 115 insertions(+), 27 deletions(-) diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 3f7924c67..5cd39b68c 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -122,10 +122,7 @@ public: virtual size_t nSensParams(); virtual void addSensitivityReaction(size_t rxn); - - virtual std::string sensParamID(int p) { - return m_pname[p]; - } + std::vector > getSensitivityOrder() const; virtual size_t componentIndex(const std::string& nm) const; @@ -143,7 +140,6 @@ protected: size_t m_nsens; std::vector m_pnum; - std::vector m_pname; std::vector m_nsens_wall; vector_fp m_mult_save; diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 3dd9b1a88..e934aaffd 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -13,6 +13,7 @@ namespace Cantera { class FlowDevice; class Wall; +class ReactorNet; const int ReservoirType = 1; const int ReactorType = 2; @@ -162,6 +163,12 @@ public: return 1; } + //! The ReactorNet that this reactor belongs to. + ReactorNet& network(); + + //! Set the ReactorNet that this reactor belongs to. + void setNetwork(ReactorNet* net); + protected: //! Number of homogeneous species in the mixture @@ -183,6 +190,9 @@ protected: std::string m_name; double m_rho0; + //! The ReactorNet that this reactor is part of + ReactorNet* m_net; + private: void tilt(const std::string& method="") const { diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 94b66cc77..bc6a8989f 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -122,7 +122,7 @@ public: void updateState(doublereal* y); double sensitivity(size_t k, size_t p) { - return m_integ->sensitivity(k, p)/m_integ->solution(k); + return m_integ->sensitivity(k, m_sensIndex[p])/m_integ->solution(k); } double sensitivity(const std::string& species, size_t p, int reactor=0) { @@ -149,6 +149,17 @@ public: size_t globalComponentIndex(const std::string& species, size_t reactor=0); + //! Used by Reactor and Wall objects to register the addition of + //! sensitivity reactions so that the ReactorNet can keep track of the + //! order in which sensitivity parameters are added. + void registerSensitivityReaction(void* reactor, size_t reactionIndex, + const std::string& name, int leftright=0); + + //! The name of the p-th sensitivity parameter added to this ReactorNet. + const std::string& sensitivityParameterName(size_t p) { + return m_paramNames.at(p); + } + void connect(size_t i, size_t j) { m_connect[j*m_nr + i] = 1; m_connect[i*m_nr + j] = 1; @@ -181,6 +192,19 @@ protected: bool m_verbose; size_t m_ntotpar; std::vector m_nparams; + + //! Names corresponding to each sensitivity parameter + std::vector m_paramNames; + + //! Structure used to determine the order of sensitivity parameters + //! m_sensOrder[Reactor or Wall, leftright][reaction number] = parameter index + std::map, std::map > m_sensOrder; + + //! Mapping from the order in which sensitivity parameters were added to + //! the ReactorNet to the order in which they occur in the integrator + //! output. + std::vector m_sensIndex; + vector_int m_connect; vector_fp m_ydot; diff --git a/include/cantera/zeroD/Wall.h b/include/cantera/zeroD/Wall.h index 1abcefc1d..e3af061af 100644 --- a/include/cantera/zeroD/Wall.h +++ b/include/cantera/zeroD/Wall.h @@ -109,6 +109,9 @@ public: /// Install the wall between two reactors or reservoirs bool install(ReactorBase& leftReactor, ReactorBase& rightReactor); + /// Called just before the start of integration + virtual void initialize(); + /// True if the wall is correctly configured and ready to use. virtual bool ready() { return (m_left != 0 && m_right != 0); @@ -171,13 +174,6 @@ public: } } void addSensitivityReaction(int leftright, size_t rxn); - std::string sensitivityParamID(int leftright, size_t p) { - if (leftright == 0) { - return m_pname_left[p]; - } else { - return m_pname_right[p]; - } - } void setSensitivityParameters(int lr, double* params); void resetSensitivityParameters(int lr); @@ -200,7 +196,6 @@ protected: std::vector m_pleft, m_pright; Cantera::vector_fp m_leftmult_save, m_rightmult_save; - std::vector m_pname_left, m_pname_right; private: diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 5bd918bd1..e821808cc 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -226,7 +226,6 @@ cdef extern from "cantera/zeroD/Reactor.h": cbool energyEnabled() void addSensitivityReaction(size_t) except + - string sensParamID(int) size_t nSensParams() @@ -260,7 +259,6 @@ cdef extern from "cantera/zeroD/Wall.h": double Q(double) void addSensitivityReaction(int, size_t) except + - string sensitivityParamID(int, size_t) size_t nSensParams(int) @@ -305,6 +303,7 @@ cdef extern from "cantera/zeroD/ReactorNet.h": double sensitivity(size_t, size_t) except + double sensitivity(string&, size_t, int) except + size_t nparams() + string sensitivityParameterName(size_t) except + cdef extern from "cantera/thermo/ThermoFactory.h" namespace "Cantera": diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 19eade726..bf2db02c8 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -190,9 +190,6 @@ cdef class Reactor(ReactorBase): def addSensitivityReaction(self, m): self.reactor.addSensitivityReaction(m) - def sensitivityParameterName(self, m): - return pystr(self.reactor.sensParamID(m)) - cdef class Reservoir(ReactorBase): """ @@ -277,9 +274,6 @@ cdef class WallSurface: def addSensitivityReaction(self, int m): self.cxxwall.addSensitivityReaction(self.side, m) - def sensitivityParameterName(self, int m): - return pystr(self.cxxwall.sensitivityParamID(m, self.side)) - cdef class Wall: r""" @@ -792,6 +786,9 @@ cdef class ReactorNet: data[k,p] = self.net.sensitivity(k,p) return data + def sensitivityParameterName(self, int p): + return pystr(self.net.sensitivityParameterName(p)) + property nSensitivityParams: def __get__(self): return self.net.nparams() diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 8233532a8..54e6973ce 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -11,6 +11,7 @@ #include "cantera/zeroD/Wall.h" #include "cantera/kinetics/InterfaceKinetics.h" #include "cantera/thermo/SurfPhase.h" +#include "cantera/zeroD/ReactorNet.h" #include @@ -86,6 +87,7 @@ void Reactor::initialize(doublereal t0) size_t nt = 0, maxnt = 0; for (size_t m = 0; m < m_nwalls; m++) { + m_wall[m]->initialize(); if (m_wall[m]->kinetics(m_lr[m])) { nt = m_wall[m]->kinetics(m_lr[m])->nTotalSpecies(); if (nt > maxnt) { @@ -102,6 +104,7 @@ void Reactor::initialize(doublereal t0) } } m_work.resize(maxnt); + std::sort(m_pnum.begin(), m_pnum.end()); m_init = true; } @@ -315,11 +318,25 @@ void Reactor::addSensitivityReaction(size_t rxn) if (rxn >= m_kin->nReactions()) throw CanteraError("Reactor::addSensitivityReaction", "Reaction number out of range ("+int2str(rxn)+")"); + + network().registerSensitivityReaction(this, rxn, + name()+": "+m_kin->reactionString(rxn)); m_pnum.push_back(rxn); - m_pname.push_back(name()+": "+m_kin->reactionString(rxn)); m_mult_save.push_back(1.0); } +std::vector > Reactor::getSensitivityOrder() const +{ + std::vector > order; + order.push_back(std::make_pair(const_cast(this), 0)); + for (size_t n = 0; n < m_nwalls; n++) { + if (m_nsens_wall[n]) { + order.push_back(std::make_pair(m_wall[n], m_lr[n])); + } + } + return order; +} + size_t Reactor::componentIndex(const string& nm) const { diff --git a/src/zeroD/ReactorBase.cpp b/src/zeroD/ReactorBase.cpp index 881a36754..0187fb20a 100644 --- a/src/zeroD/ReactorBase.cpp +++ b/src/zeroD/ReactorBase.cpp @@ -23,7 +23,8 @@ ReactorBase::ReactorBase(const string& name) : m_nsp(0), m_enthalpy(0.0), m_intEnergy(0.0), m_pressure(0.0), - m_nwalls(0) + m_nwalls(0), + m_net(0) { m_name = name; } @@ -68,6 +69,21 @@ Wall& ReactorBase::wall(size_t n) return *m_wall[n]; } +ReactorNet& ReactorBase::network() +{ + if (m_net) { + return *m_net; + } else { + throw CanteraError("ReactorBase::network", + "Reactor is not part of a ReactorNet"); + } +} + +void ReactorBase::setNetwork(ReactorNet* net) +{ + m_net = net; +} + doublereal ReactorBase::residenceTime() { int nout = static_cast(m_outlet.size()); diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index e0f1183df..96d91ce2f 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -56,6 +56,7 @@ void ReactorNet::initialize() if (m_nr == 0) throw CanteraError("ReactorNet::initialize", "no reactors in network!"); + size_t sensParamNumber = 0; for (n = 0; n < m_nr; n++) { if (m_r[n]->type() >= ReactorType) { m_r[n]->initialize(m_time); @@ -64,7 +65,16 @@ void ReactorNet::initialize() nv = r->neq(); m_size.push_back(nv); m_nparams.push_back(r->nSensParams()); - m_ntotpar += 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_nreactors++; @@ -164,6 +174,7 @@ double ReactorNet::step(doublereal time) void ReactorNet::addReactor(ReactorBase* r, bool iown) { + r->setNetwork(this); if (r->type() >= ReactorType) { m_r.push_back(r); m_iown.push_back(iown); @@ -277,4 +288,18 @@ size_t ReactorNet::globalComponentIndex(const string& species, size_t reactor) return start + m_reactors[n]->componentIndex(species); } +void ReactorNet::registerSensitivityReaction(void* reactor, + size_t reactionIndex, const std::string& name, int leftright) +{ + std::pair R = std::make_pair(reactor, leftright); + if (m_sensOrder.count(R) && + m_sensOrder[R].count(reactionIndex)) { + throw CanteraError("ReactorNet::registerSensitivityReaction", + "Attempted to register duplicate sensitivity reaction"); + } + m_paramNames.push_back(name); + m_sensOrder[R][reactionIndex] = m_ntotpar; + m_ntotpar++; +} + } diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index c394abff7..c294314ca 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -1,5 +1,6 @@ #include "cantera/zeroD/Wall.h" #include "cantera/zeroD/ReactorBase.h" +#include "cantera/zeroD/ReactorNet.h" #include "cantera/numerics/Func1.h" #include "cantera/kinetics/InterfaceKinetics.h" #include "cantera/thermo/SurfPhase.h" @@ -31,6 +32,12 @@ bool Wall::install(ReactorBase& rleft, ReactorBase& rright) return true; } +void Wall::initialize() +{ + std::sort(m_pleft.begin(), m_pleft.end()); + std::sort(m_pright.begin(), m_pright.end()); +} + /** Specify the kinetics managers for the surface mechanisms on * the left side and right side of the wall. Enter 0 if there is * no reaction mechanism. @@ -137,13 +144,15 @@ void Wall::addSensitivityReaction(int leftright, size_t rxn) throw CanteraError("Wall::addSensitivityReaction", "Reaction number out of range ("+int2str(rxn)+")"); if (leftright == 0) { + m_left->network().registerSensitivityReaction(this, rxn, + m_chem[0]->reactionString(rxn), leftright); m_pleft.push_back(rxn); m_leftmult_save.push_back(1.0); - m_pname_left.push_back(m_chem[0]->reactionString(rxn)); } else { + m_right->network().registerSensitivityReaction(this, rxn, + m_chem[1]->reactionString(rxn), leftright); m_pright.push_back(rxn); m_rightmult_save.push_back(1.0); - m_pname_right.push_back(m_chem[1]->reactionString(rxn)); } }