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'