[Reactor] Eliminate re-ordering of sensitivity parameters

Passing the full parameter vector to evalEqs for each reactor and wall
eliminates the need to re-order the parameter vector. Instead, each reactor and
wall just needs to know the indices of its sensitivity parameters, which are now
returned by ReactorNet::registerSensitivityReaction.

include/cantera/zeroD/Reactor.h

@@ -126,19 +126,6 @@ public:
     //! (in the homogeneous phase).
     virtual void addSensitivityReaction(size_t rxn);
 
-    //! Return a vector specifying the ordering of objects to use when
-    //! determining sensitivity parameter indices.
-    /*!
-     *  Used to construct ReactorNet::m_sensOrder.
-     *
-     *  @return A vector of pairs where the first element of each pair is a
-     *      pointer to either a Reactor object or a Wall object and the second
-     *      element is either 0 (in the case of a Reactor) or in the case of a
-     *      Wall indicates that the sensitivity parameters are associated with
-     *      surface chemistry on the left (0) or right (1) side of the wall.
-     */
-    std::vector<std::pair<void*, int> > getSensitivityOrder() const;
-
     //! Return the index in the solution vector for this reactor of the
     //! component named *nm*. Possible values for *nm* are "mass", "volume",
     //! "int_energy", the name of a homogeneous phase species, or the name of a
@@ -193,10 +180,8 @@ protected:
     bool m_energy;
     size_t m_nv;
 
-    size_t m_nsens;
+    // Data associated each sensitivity parameter
-    std::vector<size_t> m_pnum;
+    std::vector<SensitivityParameter> m_sensParams;
-    std::vector<size_t> m_nsens_wall;
-    vector_fp m_mult_save;
 };
 }
 

include/cantera/zeroD/ReactorBase.h

@@ -21,6 +21,13 @@ const int ConstPressureReactorType = 4;
 const int IdealGasReactorType = 5;
 const int IdealGasConstPressureReactorType = 6;
 
+struct SensitivityParameter
+{
+    size_t local; //!< local parameter index
+    size_t global; //!< global parameter index
+    double value; //!< nominal value of the parameter
+};
+
 /**
  * Base class for stirred reactors. Allows using any substance model, with
  * arbitrary inflow, outflow, heat loss/gain, surface chemistry, and volume

include/cantera/zeroD/ReactorNet.h

@@ -162,7 +162,7 @@ public:
     virtual void getState(doublereal* y);
 
     virtual size_t nparams() {
-        return m_ntotpar;
+        return m_sens_params.size();
     }
 
     //! Return the index corresponding to the component named *component* in the
@@ -173,8 +173,9 @@ public:
     //! 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,
+    //! @returns the index of this parameter in the vector of sensitivity
-                                     const std::string& name, int leftright=0);
+    //!     parameters (global across all reactors)
+    size_t registerSensitivityReaction(const std::string& name);
 
     //! The name of the p-th sensitivity parameter added to this ReactorNet.
     const std::string& sensitivityParameterName(size_t p) {
@@ -217,21 +218,10 @@ protected:
 
     int m_maxErrTestFails;
     bool m_verbose;
-    size_t m_ntotpar;
-    std::vector<size_t> m_nparams;
 
     //! Names corresponding to each sensitivity parameter
     std::vector<std::string> m_paramNames;
 
-    //! Structure used to determine the order of sensitivity parameters
-    //! m_sensOrder[Reactor or Wall, leftright][reaction number] = parameter
-    //! index
-    std::map<std::pair<void*, int>, std::map<size_t, size_t> > 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<size_t> m_sensIndex;
-
     vector_fp m_ydot;
 };
 }

include/cantera/zeroD/Wall.h

@@ -7,12 +7,11 @@
 
 #include "cantera/base/ctexceptions.h"
 #include "cantera/numerics/Func1.h"
+#include "cantera/zeroD/ReactorBase.h"
 
 namespace Cantera
 {
 
-// forward references
-class ReactorBase;
 class Kinetics;
 class SurfPhase;
 
@@ -122,7 +121,7 @@ public:
     bool install(ReactorBase& leftReactor, ReactorBase& rightReactor);
 
     //! Called just before the start of integration
-    virtual void initialize();
+    virtual void initialize() {}
 
     //! True if the wall is correctly configured and ready to use.
     virtual bool ready() {
@@ -186,8 +185,8 @@ public:
         }
     }
     void addSensitivityReaction(int leftright, size_t rxn);
-    void setSensitivityParameters(int lr, double* params);
+    void setSensitivityParameters(double* params);
-    void resetSensitivityParameters(int lr);
+    void resetSensitivityParameters();
 
 protected:
     ReactorBase* m_left;
@@ -201,7 +200,7 @@ protected:
     Func1* m_qf;
     vector_fp m_leftcov, m_rightcov;
 
-    std::vector<size_t> m_pleft, m_pright;
+    std::vector<SensitivityParameter> m_pleft, m_pright;
     vector_fp m_leftmult_save, m_rightmult_save;
 };
 

src/zeroD/FlowReactor.cpp

@@ -76,18 +76,7 @@ void FlowReactor::evalEqs(doublereal time, doublereal* y,
                           doublereal* ydot, doublereal* params)
 {
     m_thermo->restoreState(m_state);
-
+    applySensitivity(params);
-    double mult;
-    size_t n, npar;
-
-    // process sensitivity parameters
-    if (params) {
-        npar = nSensParams();
-        for (n = 0; n < npar; n++) {
-            mult = m_kin->multiplier(m_pnum[n]);
-            m_kin->setMultiplier(m_pnum[n], mult*params[n]);
-        }
-    }
 
     // distance equation
     ydot[0] = m_speed;
@@ -104,18 +93,10 @@ void FlowReactor::evalEqs(doublereal time, doublereal* y,
         fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
     }
     doublereal rrho = 1.0/m_thermo->density();
-    for (n = 0; n < m_nsp; n++) {
+    for (size_t n = 0; n < m_nsp; n++) {
         ydot[n+2] *= mw[n]*rrho;
     }
-
+    resetSensitivity(params);
-    // reset sensitivity parameters
-    if (params) {
-        npar = nSensParams();
-        for (n = 0; n < npar; n++) {
-            mult = m_kin->multiplier(m_pnum[n]);
-            m_kin->setMultiplier(m_pnum[n], mult/params[n]);
-        }
-    }
 }
 
 size_t FlowReactor::componentIndex(const string& nm) const

src/zeroD/Reactor.cpp

@@ -21,8 +21,7 @@ Reactor::Reactor() :
     m_mass(0.0),
     m_chem(false),
     m_energy(true),
-    m_nv(0),
+    m_nv(0)
-    m_nsens(npos)
 {}
 
 void Reactor::setKineticsMgr(Kinetics& kin)
@@ -117,22 +116,15 @@ void Reactor::initialize(doublereal t0)
         }
     }
     m_work.resize(maxnt);
-    std::sort(m_pnum.begin(), m_pnum.end());
 }
 
 size_t Reactor::nSensParams()
 {
-    if (m_nsens == npos) {
+    size_t ns = m_sensParams.size();
-        // determine the number of sensitivity parameters
+    for (size_t m = 0; m < m_wall.size(); m++) {
-        size_t m, ns;
+        ns += m_wall[m]->nSensParams(m_lr[m]);
-        m_nsens = m_pnum.size();
-        for (m = 0; m < m_wall.size(); m++) {
-            ns = m_wall[m]->nSensParams(m_lr[m]);
-            m_nsens_wall.push_back(ns);
-            m_nsens += ns;
-        }
     }
-    return m_nsens;
+    return ns;
 }
 
 void Reactor::syncState()
@@ -328,22 +320,8 @@ void Reactor::addSensitivityReaction(size_t rxn)
                            "Reaction number out of range ({})", rxn);
     }
 
-    network().registerSensitivityReaction(this, rxn,
+    size_t p = network().registerSensitivityReaction(name()+": "+m_kin->reactionString(rxn));
-                                          name()+": "+m_kin->reactionString(rxn));
+    m_sensParams.emplace_back(SensitivityParameter{rxn, p, 1.0});
-    m_pnum.push_back(rxn);
-    m_mult_save.push_back(1.0);
-}
-
-std::vector<std::pair<void*, int> > Reactor::getSensitivityOrder() const
-{
-    std::vector<std::pair<void*, int> > order;
-    order.emplace_back(const_cast<Reactor*>(this), 0);
-    for (size_t n = 0; n < m_wall.size(); n++) {
-        if (m_nsens_wall[n]) {
-            order.emplace_back(m_wall[n], m_lr[n]);
-        }
-    }
-    return order;
 }
 
 size_t Reactor::speciesIndex(const string& nm) const
@@ -408,19 +386,14 @@ void Reactor::applySensitivity(double* params)
     if (!params) {
         return;
     }
-    size_t npar = m_pnum.size();
+    for (auto& p : m_sensParams) {
-    for (size_t n = 0; n < npar; n++) {
+        p.value = m_kin->multiplier(p.local);
-        double mult = m_kin->multiplier(m_pnum[n]);
+        m_kin->setMultiplier(p.local, p.value*params[p.global]);
-        m_kin->setMultiplier(m_pnum[n], mult*params[n]);
     }
-    size_t ploc = npar;
     for (size_t m = 0; m < m_wall.size(); m++) {
-        if (m_nsens_wall[m] > 0) {
+        m_wall[m]->setSensitivityParameters(params);
-            m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc);
-            ploc += m_nsens_wall[m];
-        }
     }
-
+    m_kin->invalidateCache();
 }
 
 void Reactor::resetSensitivity(double* params)
@@ -428,17 +401,11 @@ void Reactor::resetSensitivity(double* params)
     if (!params) {
         return;
     }
-    size_t npar = m_pnum.size();
+    for (auto& p : m_sensParams) {
-    for (size_t n = 0; n < npar; n++) {
+        m_kin->setMultiplier(p.local, p.value);
-        double mult = m_kin->multiplier(m_pnum[n]);
-        m_kin->setMultiplier(m_pnum[n], mult/params[n]);
     }
-    size_t ploc = npar;
     for (size_t m = 0; m < m_wall.size(); m++) {
-        if (m_nsens_wall[m] > 0) {
+        m_wall[m]->resetSensitivityParameters();
-            m_wall[m]->resetSensitivityParameters(m_lr[m]);
-            ploc += m_nsens_wall[m];
-        }
     }
 }
 

src/zeroD/ReactorNet.cpp

@@ -15,7 +15,7 @@ ReactorNet::ReactorNet() :
     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(0.0), m_maxErrTestFails(0),
-    m_verbose(false), m_ntotpar(0)
+    m_verbose(false)
 {
     m_integ = newIntegrator("CVODE");
 
@@ -80,19 +80,11 @@ void ReactorNet::initialize()
         throw CanteraError("ReactorNet::initialize",
                            "no reactors in network!");
     }
-    size_t sensParamNumber = 0;
     m_start.assign(1, 0);
     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());
-        for (const auto& sens_obj : r.getSensitivityOrder()) {
-            for (const auto& order : m_sensOrder[sens_obj]) {
-                m_sensIndex.resize(std::max(order.second + 1, m_sensIndex.size()));
-                m_sensIndex[order.second] = sensParamNumber++;
-            }
-        }
         m_nv += nv;
         m_start.push_back(m_nv);
 
@@ -171,12 +163,9 @@ void ReactorNet::eval(doublereal t, doublereal* y,
                       doublereal* ydot, doublereal* p)
 {
     size_t n;
-    size_t pstart = 0;
     updateState(y);
     for (n = 0; n < m_reactors.size(); n++) {
-        m_reactors[n]->evalEqs(t, y + m_start[n],
+        m_reactors[n]->evalEqs(t, y + m_start[n], ydot + m_start[n], p);
-                               ydot + m_start[n], p + pstart);
-        pstart += m_nparams[n];
     }
     checkFinite("ydot", ydot, m_nv);
 }
@@ -186,11 +175,11 @@ double ReactorNet::sensitivity(size_t k, size_t p)
     if (!m_init) {
         initialize();
     }
-    if (p >= m_sensIndex.size()) {
+    if (p >= m_sens_params.size()) {
         throw IndexError("ReactorNet::sensitivity",
-                         "m_sensIndex", p, m_sensIndex.size()-1);
+                         "m_sens_params", p, m_sens_params.size()-1);
     }
-    return m_integ->sensitivity(k, m_sensIndex[p])/m_integ->solution(k);
+    return m_integ->sensitivity(k, p) / m_integ->solution(k);
 }
 
 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
@@ -249,24 +238,16 @@ size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
     return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
 }
 
-void ReactorNet::registerSensitivityReaction(void* reactor,
+size_t ReactorNet::registerSensitivityReaction(const std::string& name)
-        size_t reactionIndex, const std::string& name, int leftright)
 {
     if (m_integrator_init) {
         throw CanteraError("ReactorNet::registerSensitivityReaction",
                            "Sensitivity reactions cannot be added after the"
                            "integrator has been initialized.");
     }
-    std::pair<void*, int> R = {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++;
     m_sens_params.push_back(1.0);
+    return m_sens_params.size() - 1;
 }
 
 }

src/zeroD/Wall.cpp

@@ -31,12 +31,6 @@ 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());
-}
-
 void Wall::setKinetics(Kinetics* left, Kinetics* right)
 {
     m_chem[0] = left;
@@ -146,53 +140,36 @@ void Wall::addSensitivityReaction(int leftright, size_t rxn)
                            "Reaction number out of range ({})", rxn);
     }
     if (leftright == 0) {
-        m_left->network().registerSensitivityReaction(this, rxn,
+        size_t p = m_left->network().registerSensitivityReaction(
-                m_chem[0]->reactionString(rxn), leftright);
+            m_chem[0]->reactionString(rxn));
-        m_pleft.push_back(rxn);
+        m_pleft.emplace_back(SensitivityParameter{rxn, p, 1.0});
-        m_leftmult_save.push_back(1.0);
     } else {
-        m_right->network().registerSensitivityReaction(this, rxn,
+        size_t p = m_right->network().registerSensitivityReaction(
-                m_chem[1]->reactionString(rxn), leftright);
+            m_chem[1]->reactionString(rxn));
-        m_pright.push_back(rxn);
+        m_pright.emplace_back(SensitivityParameter{rxn, p, 1.0});
-        m_rightmult_save.push_back(1.0);
     }
 }
 
-void Wall::setSensitivityParameters(int lr, double* params)
+void Wall::setSensitivityParameters(double* params)
 {
     // process sensitivity parameters
-    size_t n, npar;
+    for (auto& p : m_pleft) {
-    if (lr == 0) {
+        p.value = m_chem[0]->multiplier(p.local);
-        npar = m_pleft.size();
+        m_chem[0]->setMultiplier(p.local, p.value*params[p.global]);
-        for (n = 0; n < npar; n++) {
+    }
-            m_leftmult_save[n] = m_chem[0]->multiplier(m_pleft[n]);
+    for (auto& p : m_pright) {
-            m_chem[0]->setMultiplier(m_pleft[n],
+        p.value = m_chem[1]->multiplier(p.local);
-                                     m_leftmult_save[n]*params[n]);
+        m_chem[1]->setMultiplier(p.local, p.value*params[p.global]);
-        }
-    } else {
-        npar = m_pright.size();
-        for (n = 0; n < npar; n++) {
-            m_rightmult_save[n] = m_chem[1]->multiplier(m_pright[n]);
-            m_chem[1]->setMultiplier(m_pright[n],
-                                     m_rightmult_save[n]*params[n]);
-        }
     }
 }
 
-void Wall::resetSensitivityParameters(int lr)
+void Wall::resetSensitivityParameters()
 {
-    size_t n, npar;
+    for (auto& p : m_pleft) {
-    if (lr == 0) {
+        m_chem[0]->setMultiplier(p.local, p.value);
-        npar = m_pleft.size();
+    }
-        for (n = 0; n < npar; n++) {
+    for (auto& p : m_pright) {
-            m_chem[0]->setMultiplier(m_pleft[n], m_leftmult_save[n]);
+        m_chem[1]->setMultiplier(p.local, p.value);
-        }
-    } else {
-        npar = m_pright.size();
-        for (n = 0; n < npar; n++) {
-            m_chem[1]->setMultiplier(m_pright[n],
-                                     m_rightmult_save[n]);
-        }
     }
 }
 }