From e34ac55f4ff610e55eafa8c75ebb460d2ac24435 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Fri, 2 Nov 2012 20:07:15 +0000
Subject: [PATCH] [Cython] Added sensitivity analysis to Reactor / ReactorNet

---
 include/cantera/zeroD/ReactorNet.h     | 12 ++++++
 interfaces/cython/cantera/_cantera.pxd | 16 ++++++++
 interfaces/cython/cantera/reactor.pyx  | 53 ++++++++++++++++++++++++++
 3 files changed, 81 insertions(+)

diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h
index dc70c4f15..94b66cc77 100644
--- a/include/cantera/zeroD/ReactorNet.h
+++ b/include/cantera/zeroD/ReactorNet.h
@@ -75,10 +75,22 @@ public:
     doublereal rtol() {
         return m_rtol;
     }
+
+    /// Absolute integration tolerance
     doublereal atol() {
         return m_atols;
     }
 
+    /// Relative sensitivity tolerance
+    doublereal rtolSensitivity() const {
+        return m_rtolsens;
+    }
+
+    /// Absolute sensitivity tolerance
+    doublereal atolSensitivity() const {
+        return m_atolsens;
+    }
+
     /**
      * Advance the state of all reactors in time.
      * @param time Time to advance to (s).
diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd
index 7c58fdeac..5bd918bd1 100644
--- a/interfaces/cython/cantera/_cantera.pxd
+++ b/interfaces/cython/cantera/_cantera.pxd
@@ -225,6 +225,10 @@ cdef extern from "cantera/zeroD/Reactor.h":
         void setEnergy(int)
         cbool energyEnabled()
 
+        void addSensitivityReaction(size_t) except +
+        string sensParamID(int)
+        size_t nSensParams()
+
 
 cdef extern from "cantera/zeroD/FlowReactor.h":
     cdef cppclass CxxFlowReactor "Cantera::FlowReactor" (CxxReactor):
@@ -255,6 +259,10 @@ cdef extern from "cantera/zeroD/Wall.h":
         double vdot(double)
         double Q(double)
 
+        void addSensitivityReaction(int, size_t) except +
+        string sensitivityParamID(int, size_t)
+        size_t nSensParams(int)
+
 
 cdef extern from "cantera/zeroD/flowControllers.h":
     cdef cppclass CxxFlowDevice "Cantera::FlowDevice":
@@ -289,6 +297,14 @@ cdef extern from "cantera/zeroD/ReactorNet.h":
         void setMaxTimeStep(double)
         cbool verbose()
         void setVerbose(cbool)
+        size_t neq()
+
+        void setSensitivityTolerances(double, double)
+        double rtolSensitivity()
+        double atolSensitivity()
+        double sensitivity(size_t, size_t) except +
+        double sensitivity(string&, size_t, int) except +
+        size_t nparams()
 
 
 cdef extern from "cantera/thermo/ThermoFactory.h" namespace "Cantera":
diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx
index 58135f7fb..19eade726 100644
--- a/interfaces/cython/cantera/reactor.pyx
+++ b/interfaces/cython/cantera/reactor.pyx
@@ -187,6 +187,12 @@ cdef class Reactor(ReactorBase):
         def __set__(self, pybool value):
             self.reactor.setEnergy(int(value))
 
+    def addSensitivityReaction(self, m):
+        self.reactor.addSensitivityReaction(m)
+
+    def sensitivityParameterName(self, m):
+        return pystr(self.reactor.sensParamID(m))
+
 
 cdef class Reservoir(ReactorBase):
     """
@@ -268,6 +274,12 @@ cdef class WallSurface:
                     np.ascontiguousarray(coverages, dtype=np.double)
             self.cxxwall.setCoverages(self.side, &data[0])
 
+    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"""
@@ -737,6 +749,24 @@ cdef class ReactorNet:
         def __set__(self, tol):
             self.net.setTolerances(-1, tol)
 
+    property rtolSensitivity:
+        """
+        The relative error tolerance for sensitivity analysis.
+        """
+        def __get__(self):
+            return self.net.rtolSensitivity()
+        def __set__(self, tol):
+            self.net.setSensitivityTolerances(tol, -1)
+
+    property atolSensitivity:
+        """
+        The absolute error tolerance for sensitivity analysis.
+        """
+        def __get__(self):
+            return self.net.atolSensitivity()
+        def __set__(self, tol):
+            self.net.setSensitivityTolerances(-1, tol)
+
     property verbose:
         """
         If *True*, verbose debug information will be printed during
@@ -746,3 +776,26 @@ cdef class ReactorNet:
             return pybool(self.verbose())
         def __set__(self, pybool v):
             self.net.setVerbose(v)
+
+    def sensitivity(self, species, int p, int r=0):
+        if isinstance(species, int):
+            return self.net.sensitivity(species,p)
+        elif isinstance(species, str):
+            return self.net.sensitivity(stringify(species), p, r)
+
+    def sensitivities(self):
+        cdef np.ndarray[np.double_t, ndim=2] data = \
+                np.empty((self.nVars, self.nSensitivityParams))
+        cdef int p,k
+        for p in range(self.nSensitivityParams):
+            for k in range(self.nVars):
+                data[k,p] = self.net.sensitivity(k,p)
+        return data
+
+    property nSensitivityParams:
+        def __get__(self):
+            return self.net.nparams()
+
+    property nVars:
+        def __get__(self):
+            return self.net.neq()