Merge pull request #15 from DLR-Institute-of-Future-Fuels/dev

Test added for the linear interpolation and benchmark scripts extended

.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

.gitignore

@@ -14,4 +14,5 @@ venv/
 .venv/
 thermo_data/combined_data.yaml
 src/gaspype/data/therm_data.bin
-*/*/build/
+*/*/build/
+examples/*.ipynb

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])

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])

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])

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()