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 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_cp.py b/tests/benchmark_cp.py index 330c71a..5147a68 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..966f7f4 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]) 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]) 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()