extended the test_results_cantera

examples/methane_mixtures.md

@@ -37,7 +37,7 @@ ax.set_ylabel("molar fraction")
 ax.set_ylim(0, 1.1)
 #ax.set_xlim(0, 100)
 ax.plot(ratio, equilibrium_h2o.get_x())
-ax.legend(fs.active_species)
+ax.legend(fs.species)
 ```
 
 Equilibrium calculation for methane CO2 mixtures:
@@ -56,5 +56,5 @@ ax.set_ylabel("molar fraction")
 ax.set_ylim(0, 1.1)
 #ax.set_xlim(0, 100)
 ax.plot(ratio, equilibrium_co2.get_x())
-ax.legend(fs.active_species)
+ax.legend(fs.species)
 ```

src/gaspype/_main.py

@@ -2,7 +2,7 @@ import numpy as np
 from numpy.typing import NDArray
 from typing import Sequence, Any, TypeVar, Iterator, overload, Callable
 from math import log as ln, ceil
-from scipy.linalg import null_space
+from ._numerics import null_space
 from gaspype._phys_data import atomic_weights, db_reader
 import re
 import pkgutil
@@ -90,7 +90,7 @@ class fluid_system:
 
         self._t_offset = int(t_min)
         self.species = species
-        self.active_species = species
+        self.active_species = species  # for backward compatibility
         element_compositions: list[dict[str, int]] = list()
 
         for i, s in enumerate(species):
@@ -220,7 +220,7 @@ class fluid:
                 The array can be multidimensional, the size of the last dimension
                 must match the number of species defined for the fluid_system.
                 The indices of the last dimension correspond to the indices in
-                the active_species list of the fluid_system.
+                the species list of the fluid_system.
             fs: Reference to a fluid_system. Is optional if composition is
                 defined by a dict. If not specified a new fluid_system with
                 the components from the dict is created.
@@ -585,7 +585,7 @@ class elements:
                 The array can be multidimensional, the size of the last dimension
                 must match the number of elements used in the fluid_system.
                 The indices of the last dimension correspond to the indices in
-                the active_species list of the fluid_system.
+                the species list of the fluid_system.
             fs: Reference to a fluid_system.
             shape: Tuple or list for the dimensions the fluid array. Can
                 only be used if composition argument is a dict. Otherwise

src/gaspype/_numerics.py

@@ -0,0 +1,21 @@
+import numpy as np
+from .typing import FloatArray
+
+
+def null_space(A: FloatArray) -> FloatArray:
+    """
+    Compute an orthonormal basis for the null space of A using NumPy SVD.
+
+    Args:
+        A: Input matrix of shape (m, n)
+    
+    Return:
+        Null space vectors as columns, shape (n, n - rank)
+    """
+    u, s, vh = np.linalg.svd(A, full_matrices=True)
+    M, N = u.shape[0], vh.shape[1]
+    rcond = np.finfo(s.dtype).eps * max(M, N)
+    tol = np.amax(s, initial=0.) * rcond
+    num = np.sum(s > tol, dtype=int)
+    Q = vh[num:,:].T.conj()
+    return Q

tests/test_results_cantera.py

@@ -28,14 +28,14 @@ def test_equilibrium_cantera():
 
     flow1 = ct.Solution('gri30.yaml')  # type: ignore
     ct_results = []
-    comp = (composition.array_fractions[i, j] for i in range(composition.shape[0]) for j in range(composition.shape[1]))
+    comp = [composition.array_fractions[i, j] for i in range(composition.shape[0]) for j in range(composition.shape[1])]
     for c in comp:
         comp_dict = {s: v for s, v in zip(fs.species, c)}
 
         flow1.TP = t, p
         flow1.X = comp_dict
         flow1.equilibrate('TP')  # type: ignore
-        indeces = [i for flsn in fs.active_species for i, sn in enumerate(flow1.species_names) if flsn == sn]  # type: ignore
+        indeces = [i for flsn in fs.species for i, sn in enumerate(flow1.species_names) if flsn == sn]  # type: ignore
         ct_results.append(flow1.X[indeces])  # type: ignore
 
         #if flow1.X[indeces][0] > 0.01:
@@ -45,8 +45,33 @@ def test_equilibrium_cantera():
 
     deviations = np.abs(gp_result_array - ct_result_array)
 
-    for dev, gp_comp_result, ct_comp_result in zip(deviations, gp_result_array, ct_result_array):
-        print(f"Composition: {gp_comp_result} / {ct_comp_result}")
+    for dev, gp_comp_result, ct_comp_result, c in zip(deviations, gp_result_array, ct_result_array, comp):
+        comp_dict = {s: v for s, v in zip(fs.species, c)}
+        print(f"Inp. Composition: {comp_dict}")
+        print(f"Res. Composition: {gp_comp_result}")
+        print(f"Ref. Composition: {ct_comp_result}")
+        print(f"---")
         assert np.all(dev < 0.04), f"Deviateion: {dev}"
 
     assert np.mean(deviations) < 2e-4
+
+
+def test_cantera():
+
+    t = 1495 + 273.15  # K
+    p = 1e5  # Pa
+
+    flow1 = ct.Solution('gri30.yaml')  # type: ignore
+    flow1.TP = t, p
+    inp_comp = {'CH4': 0.0, 'H2O': 0.0, 'H2': 0.9508196721311476, 'CO2': 0.0, 'CO': 0.0, 'O2': 0.04918032786885246}
+    flow1.X = inp_comp
+    flow1.equilibrate('TP')  # type: ignore
+
+    results: dict[str, float] = {sn: float(flow1.X[i]) for flsn in inp_comp for i, sn in enumerate(flow1.species_names) if flsn == sn}  # type: ignore
+
+    print(inp_comp)
+    print(results)
+
+
+if __name__ == "__main__":
+    test_cantera()