diff --git a/AUTHORS.md b/AUTHORS.md index 0c6ea987fad..942fe698c10 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -44,6 +44,7 @@ update, please report on Cantera's - **Cory Kinney** [@corykinney](https://github.com/corykinney) - **Gandhali Kogekar** [@gkogekar](https://github.com/gkogekar) - **Daniel Korff** [@korffdm](https://github.com/korffdm) - Colorado School of Mines +- **Marina Kovaleva** [@marina8888](https://github.com/marina8888) - Tohoku University - **Ashwin Kumar** [@mgashwinkumar](https://github.com/mgashwinkumar) - Virginia Tech - **Jon Kristofer** - **Samesh Lakothia** [@sameshl](https://github.com/sameshl) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index b5467776c3f..8701ed19824 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -172,6 +172,34 @@ def set_profile(self, component, positions, values): """ super().set_profile(self.flame, component, positions, values) + def get_species_reaction_sensitivities(self, species: str, ix: int) -> np.ndarray: + r""" + Compute the normalized sensitivities of a species, + :math:`X_i` with respect to the reaction rate constants :math:`k_i`: + .. math:: + s_i = \frac{k_i}{X_i} \frac{dX_i}{dk_i} + :param species: + Species of interest for sensitivity calculation (i.e "NO") + :param distance: + Index of flame grid point for which the sensitivity analysis is undertaken + """ + + def g(sim): + return sim.X[sim.gas.species_index(species), ix] + + n_vars = sum(dd.n_components * dd.n_points for dd in self.domains) + + ix_domain = self.flame.component_index(species) + self.domains[1].n_components * ix + + dgdx = np.zeros(n_vars) + dgdx[ix_domain] = 1 + X0 = g(self) + + def perturb(sim, i, dp): + sim.gas.set_multiplier(1 + dp, i) + + return self.solve_adjoint(perturb, self.gas.n_reactions, dgdx) / X0 + @property def max_grid_points(self): """ @@ -773,12 +801,12 @@ def get_flame_speed_reaction_sensitivities(self): def g(sim): return sim.velocity[0] - Nvars = sum(D.n_components * D.n_points for D in self.domains) + n_vars = sum(dd.n_components * dd.n_points for dd in self.domains) # Index of u[0] in the global solution vector i_Su = self.inlet.n_components + self.flame.component_index('velocity') - dgdx = np.zeros(Nvars) + dgdx = np.zeros(n_vars) dgdx[i_Su] = 1 Su0 = g(self) diff --git a/samples/python/onedim/species_sensitivity.py b/samples/python/onedim/species_sensitivity.py new file mode 100644 index 00000000000..4132c0232cf --- /dev/null +++ b/samples/python/onedim/species_sensitivity.py @@ -0,0 +1,86 @@ +r""" +Sensitivity analysis +======================================== +In this example we simulate an impinging jet flame, premixed ammonia/hydrogen-air flame, +calculating the sensitivity of N2O to the A-factor reaction rate constants. +Requires: cantera >= 3.0.0, pandas +.. tags:: Python, combustion, 1D flow, species reaction, premixed flame, + sensitivity analysis, plotting +""" +import cantera as ct +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +# Define the gas mixture +gas = ct.Solution("example_data/nakamura.yaml") # Use the Nakamura mechanism for ammonia combustion blends +gas.TP = 295, ct.one_atm +air = "O2:0.21,N2:0.79" +gas.set_equivalence_ratio(phi=0.6, fuel="NH3:0.7, H2:0.3", oxidizer=air) +flame = ct.ImpingingJet(gas=gas, width=0.02) +flame.inlet.mdot = 0.255 * gas.density +flame.surface.T = 493.5 +flame.set_initial_guess("equil") + +# Refine grid to improve accuracy +flame.set_refine_criteria(ratio=3, slope=0.025, curve=0.05) + +# Solve the flame +flame.solve(loglevel=1, auto=True) # Increase loglevel for more output + +# Plot temperature profile +plt.figure(figsize=(8, 6)) +plt.plot(flame.grid * 1e3, flame.T, label="Flame Temperature", color="red") +plt.xlabel("Distance (mm)") +plt.ylabel("Temperature (K)") +plt.title("Temperature Profile of a Flame") +plt.grid(True, linestyle='--', alpha=0.7) +plt.legend() +plt.tight_layout() +plt.show() + +# Create a DataFrame to store sensitivity-analysis data +sens = pd.DataFrame(index=gas.reaction_equations(), columns=["sensitivity"]) + +# Use the adjoint method to calculate species sensitivities at a set distance in the flame domain +distance = 0.02 +species = "N2O" +ix = np.argmin(np.abs(flame.grid - distance)) +sens.sensitivity = flame.get_species_reaction_sensitivities(species, ix) + +sens = sens.iloc[(-sens['sensitivity'].abs()).argsort()] +fig, ax = plt.subplots() + +sens.head(15).plot.barh(ax=ax, legend=None) +ax.invert_yaxis() # put the largest sensitivity on top +ax.set_title(f"Sensitivities for {species} Using the Adjoint-Based Method") +ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln X_{N2O}}{\partial\:\ln k}$") +ax.grid(axis='x') +plt.tight_layout() +plt.show() + +# Forward sensitivities +dk = 3e-4 + +Su0 = flame.X[gas.species_index(species), ix] +fwd = [] +for m in range(flame.gas.n_reactions): + flame.gas.set_multiplier(1.0) # reset all multipliers + flame.gas.set_multiplier(1 + dk, m) # perturb reaction m + flame.solve(loglevel=0, refine_grid=False) + Suplus = flame.X[gas.species_index(species), ix] + flame.gas.set_multiplier(1 - dk, m) # perturb reaction m + flame.solve(loglevel=0, refine_grid=False) + Suminus = flame.X[gas.species_index(species), ix] + fwd.append((Suplus - Suminus) / (2 * Su0 * dk)) +sens = pd.DataFrame(index=gas.reaction_equations(), columns=["sensitivity with forward"], data=fwd) +sens = sens.iloc[(-sens['sensitivity with forward'].abs()).argsort()] +fig, ax = plt.subplots() + +sens.head(15).plot.barh(ax=ax, legend=None) +ax.invert_yaxis() # put the largest sensitivity on top +ax.set_title(f"Sensitivities for {species} Using Brute Force Method") +ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln X_{i}}{\partial\:\ln k}$") +ax.grid(axis='x') +plt.tight_layout() +plt.show() diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 3bb891423dd..3043e779335 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -530,6 +530,53 @@ def test_adjoint_sensitivities(self): fwd = (Suplus-Suminus)/(2*Su0*dk) assert fwd == approx(dSdk_adj[m], rel=5e-3, abs=1e-7) + def test_adjoint_sensitivities2(self): + self.gas = ct.Solution('gri30.yaml') + self.gas.TP = 300, 101325 + self.gas.set_equivalence_ratio(1.0, 'CH4', {"O2": 1.0, "N2": 3.76}) + + self.flame = ct.FreeFlame(self.gas, width=0.014) + self.flame.set_refine_criteria(ratio=3, slope=0.07, curve=0.14) + self.flame.solve(loglevel=0, refine_grid=False) + + # Adjoint sensitivities + ix = -1 + species = "NO" + dk = 1e-4 + + adjoint_sensitivities = self.flame.get_species_reaction_sensitivities(species, ix) + reaction_equations = self.gas.reaction_equations() + adjoint_sens = [(reaction, sensitivity) for reaction, sensitivity in + zip(reaction_equations, adjoint_sensitivities)] + adjoint_sens_sorted = sorted(adjoint_sens, key=lambda x: abs(x[1]), reverse=True) + + Su0 = self.flame.X[self.gas.species_index(species), ix] + fwd_sensitivities = [] + for m in range(self.flame.gas.n_reactions): + self.flame.gas.set_multiplier(1.0) # reset all multipliers + self.flame.gas.set_multiplier(1 + dk, m) # perturb reaction m + self.flame.solve(loglevel=0, refine_grid=False) + Suplus = self.flame.X[self.gas.species_index(species), ix] + self.flame.gas.set_multiplier(1 - dk, m) # perturb reaction m + self.flame.solve(loglevel=0, refine_grid=False) + Suminus = self.flame.X[self.gas.species_index(species), ix] + fwd_sensitivities.append((Suplus - Suminus) / (2 * Su0 * dk)) + fwd_sens = [(reaction, sensitivity) for reaction, sensitivity in zip(reaction_equations, fwd_sensitivities)] + fwd_sens_sorted = sorted(fwd_sens, key=lambda x: abs(x[1]), reverse=True) + + # Extract top 10 reactions from both lists + top_reactions_adjoint = [reaction for reaction, _ in adjoint_sens_sorted[:10]] + top_reactions_fwd = [reaction for reaction, _ in fwd_sens_sorted[:10]] + + # Count common reactions + common_reactions = list(set(top_reactions_fwd) & set(top_reactions_adjoint)) + + # Assert that at least 8 reactions match + assert len( + common_reactions) >= 8, f"Only {len(common_reactions)} Fwd and adjoint based species sensitivities do not match" + + + def test_jacobian_options(self): reactants = {'H2': 0.65, 'O2': 0.5, 'AR': 2} self.create_sim(p=ct.one_atm, Tin=300, reactants=reactants, width=0.03)