From d08223f0426c3e7134a883ce81885c8623cb08ed Mon Sep 17 00:00:00 2001 From: Marina Date: Sun, 21 Sep 2025 13:18:08 +0100 Subject: [PATCH 1/2] Full update for example_data submodule pr nakamura mech --- AUTHORS.md | 1 + interfaces/cython/cantera/onedim.py | 34 +++++++- samples/python/onedim/species_sensitivity.py | 90 ++++++++++++++++++++ test/python/test_onedim.py | 47 ++++++++++ 4 files changed, 169 insertions(+), 3 deletions(-) create mode 100644 samples/python/onedim/species_sensitivity.py 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..191faeee672 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) + + i_X = self.flame.component_index(species) + self.domains[1].n_components * ix + + dgdx = np.zeros(n_vars) + dgdx[i_X] = 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) @@ -1539,4 +1567,4 @@ def set_initial_guess(self, data=None, group=None): self.set_profile('velocity', [0.0, 1.0], [uu, 0]) self.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) - self.set_profile("lambda", [0.0, 1.0], [L, L]) + self.set_profile("lambda", [0.0, 1.0], [L, L]) \ No newline at end of file diff --git a/samples/python/onedim/species_sensitivity.py b/samples/python/onedim/species_sensitivity.py new file mode 100644 index 00000000000..26cc454081e --- /dev/null +++ b/samples/python/onedim/species_sensitivity.py @@ -0,0 +1,90 @@ +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" +sens.sensitivity = flame.get_species_reaction_sensitivities(species, distance) + +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 Nakamura mechanism") +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 + +# get index in the grid at distance +ix = np.argmin(np.abs(flame.grid - distance)) + +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 Nakamura Mech") +ax.set_xlabel(r"Sensitivity: $\frac{\partial\:\ln X_{i}}{\partial\:\ln k}$") +ax.grid(axis='x') +plt.tight_layout() +plt.show() \ No newline at end of file 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) From 36e8004a8f2b4c28fea7b5f5bb4a26ce598e4b93 Mon Sep 17 00:00:00 2001 From: Marina Date: Sat, 4 Oct 2025 13:44:39 +0100 Subject: [PATCH 2/2] added newline and updated distance on example file --- interfaces/cython/cantera/onedim.py | 6 +++--- samples/python/onedim/species_sensitivity.py | 14 +++++--------- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/interfaces/cython/cantera/onedim.py b/interfaces/cython/cantera/onedim.py index 191faeee672..8701ed19824 100644 --- a/interfaces/cython/cantera/onedim.py +++ b/interfaces/cython/cantera/onedim.py @@ -189,10 +189,10 @@ def g(sim): n_vars = sum(dd.n_components * dd.n_points for dd in self.domains) - i_X = self.flame.component_index(species) + self.domains[1].n_components * ix + ix_domain = self.flame.component_index(species) + self.domains[1].n_components * ix dgdx = np.zeros(n_vars) - dgdx[i_X] = 1 + dgdx[ix_domain] = 1 X0 = g(self) def perturb(sim, i, dp): @@ -1567,4 +1567,4 @@ def set_initial_guess(self, data=None, group=None): self.set_profile('velocity', [0.0, 1.0], [uu, 0]) self.set_profile('spread_rate', [0.0, 1.0], [0.0, a]) - self.set_profile("lambda", [0.0, 1.0], [L, L]) \ No newline at end of file + self.set_profile("lambda", [0.0, 1.0], [L, L]) diff --git a/samples/python/onedim/species_sensitivity.py b/samples/python/onedim/species_sensitivity.py index 26cc454081e..4132c0232cf 100644 --- a/samples/python/onedim/species_sensitivity.py +++ b/samples/python/onedim/species_sensitivity.py @@ -22,7 +22,6 @@ 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) @@ -46,26 +45,23 @@ # Use the adjoint method to calculate species sensitivities at a set distance in the flame domain distance = 0.02 species = "N2O" -sens.sensitivity = flame.get_species_reaction_sensitivities(species, distance) +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 Nakamura mechanism") +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 -# get index in the grid at distance -ix = np.argmin(np.abs(flame.grid - distance)) - Su0 = flame.X[gas.species_index(species), ix] fwd = [] for m in range(flame.gas.n_reactions): @@ -83,8 +79,8 @@ 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 Nakamura Mech") +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() \ No newline at end of file +plt.show()