CITATION.cff

@@ -1,4 +1,4 @@
-cff-version: 1.2.0
+cff-version: 1.1.0
 message: "If you use this software, please cite it as below."
 title: Gaspype
 abstract: Gaspype is a performant library for thermodynamic calculations with ideal gases
@@ -9,11 +9,11 @@ authors:
     affiliation: "German Aerospace Center (DLR)"
     address: "Linder Höhe"
     city: Köln
-version: v1.1.4
+version: v1.1.3
 date-released: "2025-06-24"
-identifiers:
-  - description: This is the collection of archived snapshots of all versions of Gaspype
-    type: doi
-    value: "10.5281/zenodo.17047601"
+#identifiers:
+#  - description: This is the collection of archived snapshots of all versions of Gaspype
+#    type: doi
+#    value: ""
 license: MIT
 repository-code: "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"

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.species)
+ax.legend(fs.active_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.species)
+ax.legend(fs.active_species)
 ```

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "gaspype"
-version = "1.1.4"
+version = "1.1.3"
 authors = [
   { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
 ]
@@ -14,6 +14,7 @@ classifiers = [
 ]
 dependencies = [
   "numpy>2.0.0",
+  "scipy>1.12.0",
 ]
 
 [project.urls]
@@ -41,6 +42,7 @@ dev = [
     "cantera",
     "pyyaml>=6.0.1",
     "types-PyYAML",
+    "scipy-stubs",
     "matplotlib"
 ]
 doc_build = [

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 ._numerics import null_space
+from scipy.linalg 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  # for backward compatibility
+        self.active_species = species
         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 species list of the fluid_system.
+                the active_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 species list of the fluid_system.
+                the active_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

@@ -1,21 +0,0 @@
-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

src/gaspype/_solver.py

@@ -1,19 +1,10 @@
-from typing import Literal, Any, TYPE_CHECKING
+from typing import Literal, Any
+from scipy.optimize import minimize, root
 import numpy as np
 from ._main import elements, fluid, fluid_system
 from .typing import NDFloat, FloatArray
 from .constants import p0, epsy
 
-if TYPE_CHECKING:
-    def minimize(*a: Any, **b: Any) -> dict[str, FloatArray]:
-        ...
-else:
-    try:
-        from scipy.optimize import minimize
-    except ImportError:
-        def minimize(*a, **b):
-            raise ImportError('scipy is required for the "gibs minimization" solver')
-
 
 def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
     """
@@ -22,7 +13,8 @@ def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> N
     Solvers:
         - **system of equations** (default): Finds the root for a system of
           equations covering a minimal set of equilibrium equations and elemental balance.
-          The minimal set of equilibrium equations is derived by SVD calculating null space.
+          The minimal set of equilibrium equations is derived by SVD using the null_space
+          implementation of scipy.
 
         - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
           the elemental composition constant using the SLSQP implementation of scipy
@@ -68,7 +60,7 @@ def equilibrium_gmin(fs: fluid_system, element_composition: FloatArray, t: float
 
     start_composition_array = np.ones_like(fs.species, dtype=float)
     sol = np.array(minimize(gibbs_rt, start_composition_array, args=(grt, p_rel), method='SLSQP',
-                   bounds=bnds, constraints=cons, options={'maxiter': 2000, 'ftol': 1e-12})['x'], dtype=NDFloat)
+                   bounds=bnds, constraints=cons, options={'maxiter': 2000, 'ftol': 1e-12})['x'], dtype=NDFloat)  # type: ignore
 
     return sol
 
@@ -80,6 +72,7 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
     element_norm_log = np.log(element_norm + epsy)
 
     a = fs.array_stoichiometric_coefficients
+    a_sum = np.sum(a)
     el_matrix = fs.array_species_elements.T
 
     # Log equilibrium constants for each reaction equation
@@ -89,49 +82,42 @@ def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float,
     bp = b - np.sum(a * np.log(p / p0), axis=1)
 
     # Calculating the maximum possible amount for each species based on the elements
-    species_max = np.min((element_norm + epsy) / (fs.array_species_elements + epsy), axis=1)
-    species_max_log = np.log(species_max + epsy)
+    species_max = np.min(element_norm / (fs.array_species_elements + epsy), axis=1)
+    logn_start = np.log(species_max + epsy)
 
-    # Prepare constant arrays
-    j_eq_eye = np.eye(len(species_max))
-    j_eq_ones = np.ones((len(species_max), 1))
+    # global count
+    # count = 0
+
+    weighting = 100
 
     def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
-        n: FloatArray = np.exp(logn)  # n is the molar amount normalized by el_max
+        # global count
+        # count += 1
+        # print('------', count)
+        # assert count < 100
+        n = np.exp(logn)
         n_sum = np.sum(n)
 
         # Residuals from equilibrium equations:
-        resid_eq = a @ (logn - np.log(n_sum)) - bp
+        resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
 
-        # Jacobian for equilibrium equations:
-        j_eq = a @ (j_eq_eye - j_eq_ones * n / np.sum(n))
+        # Jacobian:
+        j_eq = a - a_sum * n / n_sum
 
         # Residuals from elemental balance:
-        el_sum_norm = np.dot(el_matrix, n) + epsy
-        resid_ab = np.log(el_sum_norm) - element_norm_log
-        #print(f'* resid_eq: {resid_eq}      resid_ab: {resid_ab}            {element_norm}')
+        el_sum = np.dot(el_matrix, n)
+        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
 
-        # Jacobian for elemental balance:
-        j_ab = el_matrix * n / el_sum_norm[:, None]
+        # print(el_sum, element_norm)
+
+        # Jacobian
+        j_ab = weighting * el_matrix * n / el_sum[:, np.newaxis]
 
         return (np.hstack([resid_eq, resid_ab]), np.concatenate([j_eq, j_ab], axis=0))
 
-    logn: FloatArray = species_max_log  # Set start values
-
-    for i in range(30):
-        rF, J = residuals(logn)
-
-        delta = np.linalg.solve(J, -rF)
-
-        logn = logn + delta
-        logn = np.minimum(logn, species_max_log + 1)
-
-        #print(f'{i} F: {np.linalg.norm(rF):.5f} lognmin={np.min(logn):.3f}, lognmax={np.max(logn):.3f}, delta={np.linalg.norm(delta):.3f} logn=')
-        if np.linalg.norm(rF) < 1e-10:
-            #print(f'Converged in {i} iterations')
-            break
-
-    n = np.exp(logn)
+    ret = root(residuals, logn_start, jac=True, tol=1e-10)
+    n = np.exp(np.array(ret['x'], dtype=NDFloat))
+    # print(ret)
 
     return n * el_max
 

tests/benchmark_equalibrium.py

@@ -44,10 +44,6 @@ eq_gaspype = gp.equilibrium(fluid, t=temperatures, p=pressure)
 elapsed_gaspype = time.perf_counter() - t0
 print(f"Gaspype: {elapsed_gaspype:.4f} s")
 
-# Check if elemental balance of result is correct
-el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2)
-assert np.all(el_err < 1e-20)
-
 # -----------------------
 # Compare first 5 results
 # -----------------------

tests/test_equalibrium.py

@@ -1,25 +0,0 @@
-import gaspype as gp
-
-
-def test_single_equilibrium():
-    # Compare equilibrium calculations to Cantera results
-
-    # gp.set_solver('system of equations')
-    # gp.set_solver('gibs minimization')
-
-    # fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
-    fs = gp.fluid_system(['CH4', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
-    # fs = gp.fluid_system([s for s in flow1.species_names if s in gps])
-
-    composition = gp.elements({'H': 2, 'O': 0, 'C': 0}, fs)
-
-    t = 1495 + 273.15  # K
-    p = 1e5  # Pa
-
-    fl = gp.equilibrium(composition, t, p)
-
-    print(fl)
-
-
-if __name__ == "__main__":
-    test_single_equilibrium()

tests/test_results.py

@@ -55,7 +55,6 @@ def test_nh3_data():
 def test_equilibrium():
     # Compare equilibrium calculations to Cycle-Tempo results
     df = pd.read_csv('tests/test_data/cycle_temp_matlab_ref.csv', sep=';', decimal=',').fillna(0)
-    #gp.set_solver('gibs minimization')
     fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'C4H10,n-butane', 'H2O', 'H2', 'CO2', 'CO'])
 
     for index, row in df.iterrows():
@@ -85,7 +84,7 @@ def test_equilibrium():
             result_values = gp.equilibrium(fl, t, p).array_fractions
 
             print(index, gp.get_solver(), '----')
-            print('Species:        ' + ''.join(f"{s:14}" for s in fs.species))
+            print(molar_comp)
             outp(result_values, 'Under test: ')
             outp(reference_values, 'Reference:  ')
 
@@ -124,7 +123,3 @@ def test_carbon():
             assert result_values > 0.9
         else:
             assert result_values < 1.1
-
-
-if __name__ == '__main__':
-    test_equilibrium()

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.species for i, sn in enumerate(flow1.species_names) if flsn == sn]  # type: ignore
+        indeces = [i for flsn in fs.active_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,33 +45,8 @@ def test_equilibrium_cantera():
 
     deviations = np.abs(gp_result_array - ct_result_array)
 
-    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("---")
+    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}")
         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_equilibrium_cantera()