From 4753ae4a867e4accdb7f1053fecb32a355a73a60 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Fri, 11 Nov 2016 22:47:05 -0500
Subject: [PATCH] [1D/Python] Add adjoint laminar flame speed sensitivity
 calculation

This is approximately an order of magnitude faster than the 'forward' method for
calculating these sensitivities. It also eliminates the need to adjust the
solver tolerances.
---
 .../examples/onedim/flamespeed_sensitivity.py | 26 ++---------------
 interfaces/cython/cantera/onedim.py           | 28 +++++++++++++++++++
 interfaces/cython/cantera/test/test_onedim.py | 22 +++++++++++++++
 3 files changed, 53 insertions(+), 23 deletions(-)

diff --git a/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py b/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
index 5436e70af..f5ebae78b 100644
--- a/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
+++ b/interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
@@ -15,8 +15,6 @@ Tin = 300.0  # unburned gas temperature [K]
 reactants = 'CH4:0.45, O2:1.0, N2:3.76'
 
 width = 0.03  # m
-tol_ss = [1.0e-9, 1.0e-14]  # [rtol atol] for steady-state problem
-tol_ts = [1.0e-5, 1.0e-14]  # [rtol atol] for time stepping
 
 # IdealGasMix object used to compute mixture properties
 gas = ct.Solution('gri30.xml', 'gri30_mix')
@@ -24,35 +22,17 @@ gas.TPX = Tin, p, reactants
 
 # Flame object
 f = ct.FreeFlame(gas, width=width)
-f.flame.set_steady_tolerances(default=tol_ss)
-f.flame.set_transient_tolerances(default=tol_ts)
 f.set_refine_criteria(ratio=3, slope=0.07, curve=0.14)
 
 f.solve(loglevel=1, auto=True)
-
-Su0 = f.u[0]
 print('\nmixture-averaged flamespeed = {:7f} m/s\n'.format(f.u[0]))
 
-print('Initial Solution:')
-f.show_stats()
-
-# Perturbation size. This must be large compared to the steady-state relative
-# tolerance (tol_ss[0]. Sensitivities less than approximately tol_ss[0] / dk
-# are not reliable.
-dk = 1e-2
+# Use the adjoint method to calculate sensitivities
+sens = f.get_flame_speed_reaction_sensitivities()
 
 print()
 print('Rxn #   k/S*dS/dk    Reaction Equation')
 print('-----   ----------   ----------------------------------')
 for m in range(gas.n_reactions):
-    gas.set_multiplier(1.0) # reset all multipliers
-    gas.set_multiplier(1+dk, m) # perturb reaction m
-    f.solve(loglevel=0, refine_grid=False)
-    Su = f.u[0]
     print('{: 5d}   {: 10.3e}   {}'.format(
-          m, (Su-Su0)/(Su0*dk), gas.reaction_equation(m)))
-
-# Sensitivity analysis requires additional function evaluations on the final
-# grid, but no additional Jacobian evaluations.
-print('\nInitial Solution + Sensitivity calculations:')
-f.show_stats()
+          m, sens[m], gas.reaction_equation(m)))
diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py
index e41f6ba19..9b26a8f29 100644
--- a/interfaces/cython/cantera/onedim.py
+++ b/interfaces/cython/cantera/onedim.py
@@ -453,6 +453,34 @@ class FreeFlame(FlameBase):
             self.set_profile(self.gas.species_name(n),
                              locs, [Y0[n], Y0[n], Yeq[n], Yeq[n]])
 
+    def get_flame_speed_reaction_sensitivities(self):
+        r"""
+        Compute the normalized sensitivities of the laminar flame speed
+        :math:`S_u` with respect to the reaction rate constants :math:`k_i`:
+
+        .. math::
+
+            s_i = \frac{k_i}{S_u} \frac{dS_u}{dk_i}
+        """
+
+        def g(sim):
+            return sim.u[0]
+
+        Nvars = sum(D.n_components * D.n_points for D in self.domains)
+
+        # Index of u[0] in the global solution vector
+        i_Su = self.inlet.n_components + self.flame.component_index('u')
+
+        dgdx = np.zeros(Nvars)
+        dgdx[i_Su] = 1
+
+        Su0 = g(self)
+
+        def perturb(sim, i, dp):
+            sim.gas.set_multiplier(1+dp, i)
+
+        return self.solve_adjoint(perturb, self.gas.n_reactions, dgdx) / Su0
+
 
 class BurnerFlame(FlameBase):
     """A burner-stabilized flat flame."""
diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py
index d667a2910..37c7e9dc8 100644
--- a/interfaces/cython/cantera/test/test_onedim.py
+++ b/interfaces/cython/cantera/test/test_onedim.py
@@ -227,6 +227,28 @@ class TestFreeFlame(utilities.CanteraTest):
     def test_mixture_averaged_case8(self):
         self.run_mix(phi=2.0, T=400, width=2.0, p=5.0, refine=False)
 
+    def test_adjoint_sensitivities(self):
+        self.run_mix(phi=0.5, T=300, width=0.1, p=1.0, refine=True)
+        self.sim.flame.set_steady_tolerances(default=(1e-10, 1e-15))
+        self.sim.solve(loglevel=0, refine_grid=False)
+
+        # Adjoint sensitivities
+        dSdk_adj = self.sim.get_flame_speed_reaction_sensitivities()
+
+        # Forward sensitivities
+        dk = 1e-4
+        Su0 = self.sim.u[0]
+        for m in range(self.gas.n_reactions):
+            self.gas.set_multiplier(1.0) # reset all multipliers
+            self.gas.set_multiplier(1+dk, m) # perturb reaction m
+            self.sim.solve(loglevel=0, refine_grid=False)
+            Suplus = self.sim.u[0]
+            self.gas.set_multiplier(1-dk, m) # perturb reaction m
+            self.sim.solve(loglevel=0, refine_grid=False)
+            Suminus = self.sim.u[0]
+            fwd = (Suplus-Suminus)/(2*Su0*dk)
+            self.assertNear(fwd, dSdk_adj[m], 5e-3)
+
     # @utilities.unittest.skip('sometimes slow')
     def test_multicomponent(self):
         reactants = 'H2:1.1, O2:1, AR:5.3'