[Reactor] Implement enthalpy of formation sensitivity analysis

include/cantera/numerics/FuncEval.h

@@ -65,6 +65,9 @@ public:
     //! This is the array which is perturbed and passed back as the fourth
     //! argument to eval().
     vector_fp m_sens_params;
+
+    //! Scaling factors for each sensitivity parameter
+    vector_fp m_paramScales;
 };
 
 }

include/cantera/zeroD/Reactor.h

@@ -126,6 +126,10 @@ public:
     //! (in the homogeneous phase).
     virtual void addSensitivityReaction(size_t rxn);
 
+    //! Add a sensitivity parameter associated with the enthalpy formation of
+    //! species *k* (in the homogeneous phase)
+    virtual void addSensitivitySpeciesEnthalpy(size_t k);
+
     //! 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

include/cantera/zeroD/ReactorBase.h

@@ -21,11 +21,17 @@ const int ConstPressureReactorType = 4;
 const int IdealGasReactorType = 5;
 const int IdealGasConstPressureReactorType = 6;
 
+enum class SensParameterType {
+    reaction,
+    enthalpy
+};
+
 struct SensitivityParameter
 {
     size_t local; //!< local parameter index
     size_t global; //!< global parameter index
     double value; //!< nominal value of the parameter
+    SensParameterType type; //!< type of sensitivity parameter
 };
 
 /**

include/cantera/zeroD/ReactorNet.h

@@ -118,15 +118,19 @@ public:
     //! Return the sensitivity of the *k*-th solution component with respect to
     //! the *p*-th sensitivity parameter.
     /*!
-     *  The normalized sensitivity coefficient \f$ S_{ki} \f$ of solution
+     *  The sensitivity coefficient \f$ S_{ki} \f$ of solution variable \f$ y_k
-     *  variable \f$ y_k \f$ with respect to sensitivity parameter \f$ p_i \f$
+     *  \f$ with respect to sensitivity parameter \f$ p_i \f$ is defined as:
-     *  is defined as:
      *
-     *  \f[ S_{ki} = \frac{p_i}{y_k} \frac{\partial y_k}{\partial p_i} \f]
+     *  \f[ S_{ki} = \frac{1}{y_k} \frac{\partial y_k}{\partial p_i} \f]
      *
      *  For reaction sensitivities, the parameter is a multiplier on the forward
      *  rate constant (and implicitly on the reverse rate constant for
-     *  reversible reactions).
+     *  reversible reactions) which has a nominal value of 1.0, and the
+     *  sensitivity is nondimensional.
+     *
+     *  For species enthalpy sensitivities, the parameter is a perturbation to
+     *  the molar enthalpy of formation, such that the dimensions of the
+     *  sensitivity are kmol/J.
      */
     double sensitivity(size_t k, size_t p);
 
@@ -171,11 +175,15 @@ public:
     size_t globalComponentIndex(const std::string& component, 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
+    //! sensitivity parameters so that the ReactorNet can keep track of the
     //! order in which sensitivity parameters are added.
+    //! @param value The nominal value of the parameter
+    //! @param scale A scaling factor to be applied to the sensitivity
+    //!     coefficient
     //! @returns the index of this parameter in the vector of sensitivity
     //!     parameters (global across all reactors)
-    size_t registerSensitivityReaction(const std::string& name);
+    size_t registerSensitivityParameter(const std::string& name, double value,
+                                        double scale);
 
     //! The name of the p-th sensitivity parameter added to this ReactorNet.
     const std::string& sensitivityParameterName(size_t p) {

interfaces/cython/cantera/_cantera.pxd

@@ -491,6 +491,7 @@ cdef extern from "cantera/zeroD/Reactor.h":
         void getState(double*)
 
         void addSensitivityReaction(size_t) except +
+        void addSensitivitySpeciesEnthalpy(size_t) except +
         size_t nSensParams()
 
 

interfaces/cython/cantera/reactor.pyx

@@ -221,6 +221,14 @@ cdef class Reactor(ReactorBase):
         """
         self.reactor.addSensitivityReaction(m)
 
+    def add_sensitivity_species_enthalpy(self, k):
+        """
+        Specifies that the sensitivity of the state variables with respect to
+        species *k* should be computed. The reactor must be part of a network
+        first.
+        """
+        self.reactor.addSensitivitySpeciesEnthalpy(self.thermo.species_index(k))
+
     def component_index(self, name):
         """
         Returns the index of the component named *name* in the system. This

interfaces/cython/cantera/test/test_reactor.py

@@ -1011,12 +1011,18 @@ class TestReactorSensitivities(utilities.CanteraTest):
         # Single reactor, changing the order in which parameters are added
         gas = ct.Solution('h2o2.xml')
 
-        def setup():
+        def setup(params):
             net = ct.ReactorNet()
             gas.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5'
 
             r = reactorClass(gas)
             net.add_reactor(r)
+
+            for kind, p in params:
+                if kind == 'r':
+                    r.add_sensitivity_reaction(p)
+                elif kind == 's':
+                    r.add_sensitivity_species_enthalpy(p)
             return r, net
 
         def integrate(r, net):
@@ -1024,26 +1030,28 @@ class TestReactorSensitivities(utilities.CanteraTest):
                 net.step()
             return net.sensitivities()
 
-        r1,net1 = setup()
+        def check_names(reactor, net, params):
-        params1 = [2,10,18,19]
+            for i,(kind,p) in enumerate(params):
-        for p in params1:
+                rname, comp = net.sensitivity_parameter_name(i).split(': ')
-            r1.add_sensitivity_reaction(p)
+                self.assertEqual(reactor.name, rname)
+                if kind == 'r':
+                    self.assertEqual(gas.reaction_equation(p), comp)
+                elif kind == 's':
+                    self.assertEqual(p + ' enthalpy', comp)
+
+        params1 = [('r', 2), ('r', 10), ('r', 18), ('r', 19), ('s', 'O2'),
+                   ('s', 'OH'), ('s', 'H2O2')]
+        r1,net1 = setup(params1)
         S1 = integrate(r1, net1)
+        check_names(r1, net1, params1)
 
-        pname = lambda r,i: '%s: %s' % (r.name, gas.reaction_equation(i))
+        params2 = [('r', 19), ('s', 'H2O2'), ('s', 'OH'), ('r', 10),
-        for i,p in enumerate(params1):
+                   ('s', 'O2'), ('r', 2), ('r', 18)]
-            self.assertEqual(pname(r1,p), net1.sensitivity_parameter_name(i))
+        r2,net2 = setup(params2)
-
-        r2,net2 = setup()
-        params2 = [19,10,2,18]
-        for p in params2:
-            r2.add_sensitivity_reaction(p)
         S2 = integrate(r2, net2)
+        check_names(r2, net2, params2)
 
-        for i,p in enumerate(params2):
+        for i,j in enumerate((5,3,6,0,4,2,1)):
-            self.assertEqual(pname(r2,p), net2.sensitivity_parameter_name(i))
-
-        for i,j in enumerate((2,1,3,0)):
             self.assertArrayNear(S1[:,i], S2[:,j])
 
     def test_parameter_order1a(self):

src/numerics/CVodesIntegrator.cpp

@@ -235,7 +235,12 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
     if (flag != CV_SUCCESS) {
         throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensInit");
     }
-    vector_fp atol(m_np, m_abstolsens);
+    vector_fp atol(m_np);
+    for (size_t n = 0; n < m_np; n++) {
+        // This scaling factor is tuned so that reaction and species enthalpy
+        // sensitivities can be computed simultaneously with the same abstol.
+        atol[n] = m_abstolsens / func.m_paramScales[n];
+    }
     flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data());
 }
 
@@ -312,7 +317,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
     if (func.nparams() > 0) {
         sensInit(t0, func);
         flag = CVodeSetSensParams(m_cvode_mem, func.m_sens_params.data(),
-                                  NULL, NULL);
+                                  func.m_paramScales.data(), NULL);
+        if (flag != CV_SUCCESS) {
+            throw CanteraError("CVodesIntegrator::initialize",
+                               "CVodeSetSensParams failed.");
+        }
     }
     applyOptions();
 }

src/zeroD/Reactor.cpp

@@ -320,8 +320,25 @@ void Reactor::addSensitivityReaction(size_t rxn)
                            "Reaction number out of range ({})", rxn);
     }
 
-    size_t p = network().registerSensitivityReaction(name()+": "+m_kin->reactionString(rxn));
+    size_t p = network().registerSensitivityParameter(
-    m_sensParams.emplace_back(SensitivityParameter{rxn, p, 1.0});
+        name()+": "+m_kin->reactionString(rxn), 1.0, 1.0);
+    m_sensParams.emplace_back(
+        SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
+}
+
+void Reactor::addSensitivitySpeciesEnthalpy(size_t k)
+{
+    if (k >= m_thermo->nSpecies()) {
+        throw CanteraError("Reactor::addSensitivitySpeciesEnthalpy",
+                           "Species index out of range ({})", k);
+    }
+
+    size_t p = network().registerSensitivityParameter(
+        name() + ": " + m_thermo->speciesName(k) + " enthalpy",
+        0.0, GasConstant * 298.15);
+    m_sensParams.emplace_back(
+        SensitivityParameter{k, p, m_thermo->Hf298SS(k),
+                             SensParameterType::enthalpy});
 }
 
 size_t Reactor::speciesIndex(const string& nm) const
@@ -387,12 +404,17 @@ void Reactor::applySensitivity(double* params)
         return;
     }
     for (auto& p : m_sensParams) {
-        p.value = m_kin->multiplier(p.local);
+        if (p.type == SensParameterType::reaction) {
-        m_kin->setMultiplier(p.local, p.value*params[p.global]);
+            p.value = m_kin->multiplier(p.local);
+            m_kin->setMultiplier(p.local, p.value*params[p.global]);
+        } else if (p.type == SensParameterType::enthalpy) {
+            m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]);
+        }
     }
     for (size_t m = 0; m < m_wall.size(); m++) {
         m_wall[m]->setSensitivityParameters(params);
     }
+    m_thermo->invalidateCache();
     m_kin->invalidateCache();
 }
 
@@ -402,11 +424,17 @@ void Reactor::resetSensitivity(double* params)
         return;
     }
     for (auto& p : m_sensParams) {
-        m_kin->setMultiplier(p.local, p.value);
+        if (p.type == SensParameterType::reaction) {
+            m_kin->setMultiplier(p.local, p.value);
+        } else if (p.type == SensParameterType::enthalpy) {
+            m_thermo->resetHf298(p.local);
+        }
     }
     for (size_t m = 0; m < m_wall.size(); m++) {
         m_wall[m]->resetSensitivityParameters();
     }
+    m_thermo->invalidateCache();
+    m_kin->invalidateCache();
 }
 
 }

src/zeroD/ReactorNet.cpp

@@ -13,7 +13,7 @@ namespace Cantera
 ReactorNet::ReactorNet() :
     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_atols(1.0e-15), m_atolsens(1.0e-6),
     m_maxstep(0.0), m_maxErrTestFails(0),
     m_verbose(false)
 {
@@ -179,7 +179,11 @@ double ReactorNet::sensitivity(size_t k, size_t p)
         throw IndexError("ReactorNet::sensitivity",
                          "m_sens_params", p, m_sens_params.size()-1);
     }
-    return m_integ->sensitivity(k, p) / m_integ->solution(k);
+    double denom = m_integ->solution(k);
+    if (denom == 0.0) {
+        denom = SmallNumber;
+    }
+    return m_integ->sensitivity(k, p) / denom;
 }
 
 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
@@ -238,15 +242,17 @@ size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor)
     return m_start[reactor] + m_reactors[reactor]->componentIndex(component);
 }
 
-size_t ReactorNet::registerSensitivityReaction(const std::string& name)
+size_t ReactorNet::registerSensitivityParameter(
+    const std::string& name, double value, double scale)
 {
     if (m_integrator_init) {
-        throw CanteraError("ReactorNet::registerSensitivityReaction",
+        throw CanteraError("ReactorNet::registerSensitivityParameter",
-                           "Sensitivity reactions cannot be added after the"
+                           "Sensitivity parameters cannot be added after the"
                            "integrator has been initialized.");
     }
     m_paramNames.push_back(name);
-    m_sens_params.push_back(1.0);
+    m_sens_params.push_back(value);
+    m_paramScales.push_back(scale);
     return m_sens_params.size() - 1;
 }
 

src/zeroD/Wall.cpp

@@ -140,13 +140,15 @@ void Wall::addSensitivityReaction(int leftright, size_t rxn)
                            "Reaction number out of range ({})", rxn);
     }
     if (leftright == 0) {
-        size_t p = m_left->network().registerSensitivityReaction(
+        size_t p = m_left->network().registerSensitivityParameter(
-            m_chem[0]->reactionString(rxn));
+            m_chem[0]->reactionString(rxn), 1.0, 1.0);
-        m_pleft.emplace_back(SensitivityParameter{rxn, p, 1.0});
+        m_pleft.emplace_back(
+            SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
     } else {
-        size_t p = m_right->network().registerSensitivityReaction(
+        size_t p = m_right->network().registerSensitivityParameter(
-            m_chem[1]->reactionString(rxn));
+            m_chem[1]->reactionString(rxn), 1.0, 1.0);
-        m_pright.emplace_back(SensitivityParameter{rxn, p, 1.0});
+        m_pright.emplace_back(
+            SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction});
     }
 }