[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.

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)))

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."""

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'