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)); } }