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,7 +9,7 @@ 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

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/_numerics.py

@@ -8,7 +8,7 @@ def null_space(A: FloatArray) -> FloatArray:
 
     Args:
         A: Input matrix of shape (m, n)
-
+    
     Return:
         Null space vectors as columns, shape (n, n - rank)
     """
@@ -17,5 +17,5 @@ def null_space(A: FloatArray) -> FloatArray:
     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()
+    Q = vh[num:,:].T.conj()
-    return Q
+    return Q

src/gaspype/_solver.py

@@ -1,18 +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:
-    from scipy.optimize import minimize
-else:
-    try:
-        from scipy.optimize import minimize
-    except ImportError:
-        def minimize(*a: Any, **b: Any) -> Any:
-            raise ImportError('scipy is required for the "gibs minimization" solver')
-
 
 def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
     """
@@ -21,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
@@ -79,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
@@ -88,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 = np.min(element_norm / (fs.array_species_elements + epsy), axis=1)
-    species_max_log = np.log(species_max + epsy)
+    logn_start = np.log(species_max + epsy)
 
-    # Prepare constant arrays
+    # global count
-    j_eq_eye = np.eye(len(species_max))
+    # count = 0
-    j_eq_ones = np.ones((len(species_max), 1))
+
+    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:
+        # Jacobian:
-        j_eq = a @ (j_eq_eye - j_eq_ones * n / np.sum(n))
+        j_eq = a - a_sum * n / n_sum
 
         # Residuals from elemental balance:
-        el_sum_norm = np.dot(el_matrix, n) + epsy
+        el_sum = np.dot(el_matrix, n)
-        resid_ab = np.log(el_sum_norm) - element_norm_log
+        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
-        #print(f'* resid_eq: {resid_eq}      resid_ab: {resid_ab}            {element_norm}')
 
-        # Jacobian for elemental balance:
+        # print(el_sum, element_norm)
-        j_ab = el_matrix * n / el_sum_norm[:, None]
+
+        # 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
+    ret = root(residuals, logn_start, jac=True, tol=1e-10)
-
+    n = np.exp(np.array(ret['x'], dtype=NDFloat))
-    for i in range(30):
+    # print(ret)
-        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)
 
     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

@@ -50,7 +50,7 @@ def test_equilibrium_cantera():
         print(f"Inp. Composition: {comp_dict}")
         print(f"Res. Composition: {gp_comp_result}")
         print(f"Ref. Composition: {ct_comp_result}")
-        print("---")
+        print(f"---")
         assert np.all(dev < 0.04), f"Deviateion: {dev}"
 
     assert np.mean(deviations) < 2e-4
@@ -74,4 +74,4 @@ def test_cantera():
 
 
 if __name__ == "__main__":
-    test_equilibrium_cantera()
+    test_cantera()