From 317bdf73851fdb63e2e5d0f75dac834b00a4ad46 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Fri, 10 Jul 2026 09:59:12 +0200 Subject: [PATCH 1/4] Equilibrium benchmark updated: NASA CEA added --- .gitignore | 3 ++- tests/benchmark_equalibrium.py | 38 +++++++++++++++++++++++++++++++++- 2 files changed, 39 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 98fb789..6a9551c 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,5 @@ venv/ .venv/ thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin -*/*/build/ \ No newline at end of file +*/*/build/ +examples/*.ipynb diff --git a/tests/benchmark_equalibrium.py b/tests/benchmark_equalibrium.py index 4c0fe3a..8378b1f 100644 --- a/tests/benchmark_equalibrium.py +++ b/tests/benchmark_equalibrium.py @@ -4,9 +4,16 @@ import numpy as np import time from gaspype import fluid_system +try: + import cea + CEA_AVAILABLE = True +except ImportError: + CEA_AVAILABLE = False + # ----------------------- # Settings # ----------------------- + n_temps = 1000 temps_C = np.linspace(300, 1000, n_temps) # °C temperatures = temps_C + 273.15 # K @@ -48,11 +55,40 @@ print(f"Gaspype: {elapsed_gaspype:.4f} s") el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2) assert np.all(el_err < 1e-20) +if CEA_AVAILABLE: + reac_names = ["CH4", "H2O"] + prod_names = ["CH4", "H2O", "CO", "CO2", "H2", "O2", "H", "O", "OH"] + + reac = cea.Mixture(reac_names) + prod = cea.Mixture(prod_names) + + solver = cea.EqSolver(prod, reactants=reac) + solution = cea.EqSolution(solver) + + input_moles = np.array([8.0, 2.0]) + input_weights = reac.moles_to_weights(input_moles) + + p_bar = pressure * 1e-5 # Convert to bar + + eq_cea = np.zeros((n_temps, len(species_to_track))) + + time.sleep(0.5) + t0 = time.perf_counter() + for i, T in enumerate(temperatures): + solver.solve(solution, cea.TP, T, p_bar, input_weights) + if solution.converged: + for j, s in enumerate(species_to_track): + eq_cea[i, j] = solution.mole_fractions.get(s, 0.0) + elapsed_cea = time.perf_counter() - t0 + print(f"CEA: {elapsed_cea:.4f} s") + # ----------------------- # Compare first 5 results # ----------------------- -print("First 5 equilibrium compositions (mole fractions):") +print("\nFirst 5 equilibrium compositions:") for i in range(5): print(f"T = {temperatures[i]:.1f} K") print(" Cantera:", eq_cantera[i]) print(" Gaspype :", eq_gaspype.array_composition[i]) + if CEA_AVAILABLE: + print(" CEA :", eq_cea[i]) From 662239e49716c2e09a09fc7e8b163d5e6b62c618 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Fri, 10 Jul 2026 10:31:02 +0200 Subject: [PATCH 2/4] Lookup properties benchmark updated: NASA CEA added --- tests/benchmark_cp.py | 35 +++++++++++++++++++++++++++++------ tests/benchmark_cp_comp.py | 29 +++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 6 deletions(-) diff --git a/tests/benchmark_cp.py b/tests/benchmark_cp.py index 330c71a..49fe1ae 100644 --- a/tests/benchmark_cp.py +++ b/tests/benchmark_cp.py @@ -3,22 +3,25 @@ import numpy as np import time import gaspype as gp +try: + import cea + CEA_AVAILABLE = True +except ImportError: + CEA_AVAILABLE = False + gas = ct.Solution("gri30.yaml") composition = {"H2": 0.3, "H2O": 0.3, "N2": 0.4} n_species = gas.n_species n_states = 1_000_000 -# Random temperatures and pressures temperatures = np.linspace(300.0, 2500.0, n_states) pressures = np.full(n_states, ct.one_atm) -# Create a SolutionArray with many states at once states = ct.SolutionArray(gas, len(temperatures)) time.sleep(0.5) -# Vectorized assignment t0 = time.perf_counter() states.TPX = temperatures, pressures, composition cp_values = states.cp_mole @@ -27,16 +30,36 @@ elapsed = time.perf_counter() - t0 print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized cantera)") print("First 5 Cp values (J/mol-K):", cp_values[:5] / 1000) - -# Vectorized fluid creation fluid = gp.fluid(composition) time.sleep(0.5) -# Benchmark: calculate Cp for all states at once t0 = time.perf_counter() cp_values = fluid.get_cp(t=temperatures) elapsed = time.perf_counter() - t0 print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)") print("First 5 Cp values (J/mol·K):", cp_values[:5]) + +if CEA_AVAILABLE: + cea_mix = cea.Mixture(['H2', 'H2O', 'N2']) + mole_fracs = np.array([0.3, 0.3, 0.4]) + MW = np.array([2.016, 18.015, 28.014]) + mass_weights = mole_fracs * MW / (mole_fracs * MW).sum() + avg_MW = np.sum(mole_fracs * MW) + p_bar = cea.units.atm_to_bar(1.0) + + time.sleep(0.5) + + # the current NASA CEA Python API does not provide a NumPy-style + # vectorized interface for thermodynamic property lookups + t0 = time.perf_counter() + cea_cp = np.zeros(n_states) + for i in range(n_states): + cea_cp[i] = cea_mix.calc_property(cea.FROZEN_CP, mass_weights, temperatures[i], p_bar) + elapsed = time.perf_counter() - t0 + + cea_cp_molar = cea_cp * avg_MW / 1000 + + print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (CEA)") + print("First 5 Cp values (J/mol-K):", cea_cp_molar[:5]) diff --git a/tests/benchmark_cp_comp.py b/tests/benchmark_cp_comp.py index 7a01615..6c6c46f 100644 --- a/tests/benchmark_cp_comp.py +++ b/tests/benchmark_cp_comp.py @@ -3,6 +3,12 @@ import numpy as np import time import gaspype as gp +try: + import cea + CEA_AVAILABLE = True +except ImportError: + CEA_AVAILABLE = False + gas = ct.Solution("gri30.yaml") n_species = gas.n_species n_states = 1_000_000 @@ -53,3 +59,26 @@ elapsed = time.perf_counter() - t0 print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)") print("First 5 Cp values (J/mol·K):", cp_values[:5]) + + +if CEA_AVAILABLE: + MW = np.array([2.016, 18.015, 28.014]) + mass_weights = fractions * MW / (fractions * MW).sum(axis=1)[:, None] + avg_MW = np.sum(fractions * MW, axis=1) + p_bar = cea.units.atm_to_bar(1.0) + cea_mix = cea.Mixture(['H2', 'H2O', 'N2']) + + time.sleep(0.5) + + # the current NASA CEA Python API does not provide a NumPy-style + # vectorized interface for thermodynamic property lookups + t0 = time.perf_counter() + cea_cp = np.zeros(n_states) + for i in range(n_states): + cea_cp[i] = cea_mix.calc_property(cea.FROZEN_CP, mass_weights[i], temperatures[i], p_bar) + elapsed = time.perf_counter() - t0 + + cea_cp_molar = cea_cp * avg_MW / 1000 + + print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (CEA)") + print("First 5 Cp values (J/mol-K):", cea_cp_molar[:5]) From 0ae78ee6317b2769b7c25e41cabed26e53444924 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Tue, 14 Jul 2026 08:54:59 +0200 Subject: [PATCH 3/4] Test added for interpolation error of the lookup table --- tests/benchmark_cp.py | 2 +- tests/benchmark_cp_comp.py | 2 +- tests/test_lookup_table_error.py | 54 ++++++++++++++++++++++++++++++++ 3 files changed, 56 insertions(+), 2 deletions(-) create mode 100644 tests/test_lookup_table_error.py diff --git a/tests/benchmark_cp.py b/tests/benchmark_cp.py index 49fe1ae..5147a68 100644 --- a/tests/benchmark_cp.py +++ b/tests/benchmark_cp.py @@ -51,7 +51,7 @@ if CEA_AVAILABLE: time.sleep(0.5) - # the current NASA CEA Python API does not provide a NumPy-style + # the current NASA CEA Python API does not provide a NumPy-style # vectorized interface for thermodynamic property lookups t0 = time.perf_counter() cea_cp = np.zeros(n_states) diff --git a/tests/benchmark_cp_comp.py b/tests/benchmark_cp_comp.py index 6c6c46f..966f7f4 100644 --- a/tests/benchmark_cp_comp.py +++ b/tests/benchmark_cp_comp.py @@ -70,7 +70,7 @@ if CEA_AVAILABLE: time.sleep(0.5) - # the current NASA CEA Python API does not provide a NumPy-style + # the current NASA CEA Python API does not provide a NumPy-style # vectorized interface for thermodynamic property lookups t0 = time.perf_counter() cea_cp = np.zeros(n_states) diff --git a/tests/test_lookup_table_error.py b/tests/test_lookup_table_error.py new file mode 100644 index 0000000..a5cc3a7 --- /dev/null +++ b/tests/test_lookup_table_error.py @@ -0,0 +1,54 @@ +""" +Testing the interpolation error of the lookup table for Cp, G, and H for multiple species. +""" + +import gaspype as gp +import numpy as np +from gaspype.constants import R + +fl = gp.fluid({'CH4': 1, 'CO2': 1, 'H2O': 1, 'SO2': 1, 'NH3': 1, 'NO2': 1, 'H2S': 1, 'C2H6': 1, 'C3H8': 1, 'H2': 1, 'O2': 1}) +t_values = np.arange(200, 2000, step=10, dtype=int) + + +def test_cp_continuity_multi_species(): + cp_actual = fl.fs.get_species_cp(t_values) + cp_neighbor_mean = (fl.fs.get_species_cp(t_values - 1) + fl.fs.get_species_cp(t_values + 1)) / 2 + diff = (cp_actual - cp_neighbor_mean) / cp_actual / 2 + avg_diff = np.mean(np.abs(diff)) + max_diff = np.max(np.abs(diff)) + print(f'Average difference of Cp: {avg_diff} J/mol/K') + print(f'Max difference of Cp: {max_diff} J/mol/K') + assert avg_diff < 1e-5 + assert max_diff < 1e-5 + + +def test_g_continuity_multi_species(): + g_actual = fl.fs.get_species_g_rt(t_values) * R * np.reshape(t_values, (-1, 1)) + g_neighbor_mean = (fl.fs.get_species_g_rt(t_values - 1) + fl.fs.get_species_g_rt(t_values + 1)) / 2 * R * np.reshape(t_values, (-1, 1)) + + diff = (g_actual - g_neighbor_mean) / g_actual / 2 + avg_diff = np.mean(np.abs(diff)) + max_diff = np.max(np.abs(diff)) + print(f'Average difference of G: {avg_diff} J/mol') + print(f'Max difference of G: {max_diff} J/mol') + assert avg_diff < 1e-4 + assert max_diff < 1e-4 + + +def test_h_continuity_multi_species(): + h_actual = fl.fs.get_species_h(t_values) + h_neighbor_mean = (fl.fs.get_species_h(t_values - 1) + fl.fs.get_species_h(t_values + 1)) / 2 + + rel_diff = (h_actual - h_neighbor_mean) / h_actual / 2 + avg_diff = np.mean(np.abs(rel_diff)) + max_diff = np.max(np.abs(rel_diff)) + print(f'Average difference of H: {avg_diff} J/mol') + print(f'Max difference of H: {max_diff} J/mol') + assert avg_diff < 1e-4 + assert max_diff < 1e-4 + + +if __name__ == "__main__": + test_cp_continuity_multi_species() + test_g_continuity_multi_species() + test_h_continuity_multi_species() From e0ddfb85bbca9103730c65d79adca8153d378552 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Tue, 14 Jul 2026 08:59:29 +0200 Subject: [PATCH 4/4] CI extended for Python 3.14 --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6b3d005..36f6112 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: - python-version: ["3.10", 3.11, 3.12, 3.13] + python-version: ["3.10", 3.11, 3.12, 3.13, 3.14] steps: - name: Check out code