[1D] Add general support for calculating adjoint sensitivities

include/cantera/oneD/Sim1D.h

@@ -120,6 +120,11 @@ public:
         OneDim::eval(npos, m_x.data(), m_xnew.data(), rdt, count);
     }
 
+    // Evaluate the governing equations and return the vector of residuals
+    void getResidual(double rdt, double* resid) {
+        OneDim::eval(npos, m_x.data(), resid, rdt, 0);
+    }
+
     /// Refine the grid in all domains.
     int refine(int loglevel=0);
 
@@ -183,6 +188,21 @@ public:
 
     void evalSSJacobian();
 
+    //! Solve the equation \f$ J^T \lambda = b \f$.
+    /**
+     * Here, \f$ J = \partial f/\partial x \f$ is the Jacobian matrix of the
+     * system of equations \f$ f(x,p)=0 \f$. This can be used to efficiently
+     * solve for the sensitivities of a scalar objective function \f$ g(x,p) \f$
+     * to a vector of parameters \f$ p \f$ by solving:
+     * \f[ J^T \lambda = \left( \frac{\partial g}{\partial x} \right)^T \f]
+     * for \f$ \lambda \f$ and then computing:
+     * \f[
+     *     \left.\frac{dg}{dp}\right|_{f=0} = \frac{\partial g}{\partial p}
+     *         - \lambda^T \frac{\partial f}{\partial p}
+     * \f]
+     */
+    void solveAdjoint(const double* b, double* lambda);
+
     virtual void resize();
 
     //! Set a function that will be called after each successful steady-state

interfaces/cython/cantera/_cantera.pxd

@@ -723,7 +723,9 @@ cdef extern from "cantera/oneD/Sim1D.h":
         int domainIndex(string) except +translate_exception
         double value(size_t, size_t, size_t) except +translate_exception
         double workValue(size_t, size_t, size_t) except +translate_exception
-        void eval(double, int) except +translate_exception
+        size_t size()
+        void solveAdjoint(const double*, double*) except +translate_exception
+        void getResidual(double, double*) except +translate_exception
         void setJacAge(int, int)
         void setTimeStepFactor(double)
         void setMinTimeStep(double)

interfaces/cython/cantera/onedim.pyx

@@ -1067,6 +1067,68 @@ cdef class Sim1D:
         """
         self.sim.clearStats()
 
+    def solve_adjoint(self, perturb, n_params, dgdx, g=None, dp=1e-5):
+        r"""
+        Find the sensitivities of an objective function using an adjoint method.
+
+        For an objective function :math:`g(x, p)` where :math:`x` is the state
+        vector of the system and :math:`p` is a vector of parameters, this
+        computes the vector of sensitivities :math:`dg/dp`. This assumes that
+        the system of equations has already been solved to find :math:`x`.
+
+        :param perturb:
+            A function with the signature ``perturb(sim, i, dp)`` which
+            perturbs parameter ``i`` by a relative factor of ``dp``. To
+            perturb a reaction rate constant, this function could be defined
+            as::
+                def perturb(sim, i, dp):
+                    sim.gas.set_multiplier(1+dp, i)
+            Calling ``perturb(sim, i, 0)`` should restore that parameter to its
+            default value.
+        :param n_params:
+            The length of the vector of sensitivity parameters
+        :param dgdx:
+            The vector of partial derivatives of the function :math:`g(x, p)`
+            with respect to the system state :math:`x`.
+        :param g:
+            A function with the signature ``value = g(sim)`` which computes the
+            value of :math:`g(x,p)` at the current system state. This is used to
+            compute :math:`\partial g/\partial p`. If this is identically zero
+            (i.e. :math:`g` is independent of :math:`p`) then this argument may
+            be omitted.
+        :param dp:
+            A relative value by which to perturb each parameter
+        """
+        n_vars = self.sim.size()
+        cdef np.ndarray[np.double_t, ndim=1] L = np.empty(n_vars)
+        cdef np.ndarray[np.double_t, ndim=1] gg = \
+                np.ascontiguousarray(dgdx, dtype=np.double)
+
+        self.sim.solveAdjoint(&gg[0], &L[0])
+
+        cdef np.ndarray[np.double_t, ndim=1] dgdp = np.empty(n_params)
+        cdef np.ndarray[np.double_t, ndim=2] dfdp = np.empty((n_vars, n_params))
+        cdef np.ndarray[np.double_t, ndim=1] fplus = np.empty(n_vars)
+        cdef np.ndarray[np.double_t, ndim=1] fminus = np.empty(n_vars)
+        gplus = gminus = 0
+
+        for i in range(n_params):
+            perturb(self, i, dp)
+            if g:
+                gplus = g(self)
+            self.sim.getResidual(0, &fplus[0])
+
+            perturb(self, i, -dp)
+            if g:
+                gminus = g(self)
+            self.sim.getResidual(0, &fminus[0])
+
+            perturb(self, i, 0)
+            dgdp[i] = (gplus - gminus)/(2*dp)
+            dfdp[:,i] = (fplus - fminus) / (2*dp)
+
+        return dgdp - np.dot(L, dfdp)
+
     property grid_size_stats:
         """Return total grid size in each call to solve()"""
         def __get__(self):

src/oneD/Sim1D.cpp

@@ -567,6 +567,24 @@ void Sim1D::evalSSJacobian()
     OneDim::evalSSJacobian(m_x.data(), m_xnew.data());
 }
 
+void Sim1D::solveAdjoint(const double* b, double* lambda)
+{
+    evalSSJacobian();
+
+    // Form J^T
+    size_t bw = bandwidth();
+    BandMatrix Jt(size(), bw, bw);
+    for (size_t i = 0; i < size(); i++) {
+        size_t j1 = (i > bw) ? i - bw : 0;
+        size_t j2 = (i + bw >= size()) ? size() - 1: i + bw;
+        for (size_t j = j1; j <= j2; j++) {
+            Jt(j,i) = m_jac->value(i,j);
+        }
+    }
+
+    Jt.solve(b, lambda);
+}
+
 void Sim1D::resize()
 {
     OneDim::resize();