From 662239e49716c2e09a09fc7e8b163d5e6b62c618 Mon Sep 17 00:00:00 2001
From: Nicolas <Nicolas@nonan.net>
Date: Fri, 10 Jul 2026 10:31:02 +0200
Subject: [PATCH] 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])