Updated the version of the CITATION.cff file

CITATION.cff

@@ -1,4 +1,4 @@
-cff-version: 1.1.0
+cff-version: 1.2.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.3
+version: v1.1.4
 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.3"
+version = "1.1.4"
 authors = [
   { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
 ]
@@ -14,7 +14,6 @@ classifiers = [
 ]
 dependencies = [
   "numpy>2.0.0",
-  "scipy>1.12.0",
 ]
 
 [project.urls]
@@ -42,7 +41,6 @@ dev = [
     "cantera",
     "pyyaml>=6.0.1",
     "types-PyYAML",
-    "scipy-stubs",
     "matplotlib"
 ]
 doc_build = [

src/gaspype/_solver.py

@@ -1,10 +1,18 @@
-from typing import Literal, Any
-from scipy.optimize import minimize, root
+from typing import Literal, Any, TYPE_CHECKING
 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:
     """
@@ -13,8 +21,7 @@ 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 using the null_space
-          implementation of scipy.
+          The minimal set of equilibrium equations is derived by SVD calculating null space.
 
         - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
           the elemental composition constant using the SLSQP implementation of scipy
@@ -72,7 +79,6 @@ 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
@@ -82,42 +88,49 @@ 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 / (fs.array_species_elements + epsy), axis=1)
-    logn_start = np.log(species_max + epsy)
+    species_max = np.min((element_norm + epsy) / (fs.array_species_elements + epsy), axis=1)
+    species_max_log = np.log(species_max + epsy)
 
-    # global count
-    # count = 0
-
-    weighting = 100
+    # Prepare constant arrays
+    j_eq_eye = np.eye(len(species_max))
+    j_eq_ones = np.ones((len(species_max), 1))
 
     def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
-        # global count
-        # count += 1
-        # print('------', count)
-        # assert count < 100
-        n = np.exp(logn)
+        n: FloatArray = np.exp(logn)  # n is the molar amount normalized by el_max
         n_sum = np.sum(n)
 
         # Residuals from equilibrium equations:
-        resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
+        resid_eq = a @ (logn - np.log(n_sum)) - bp
 
-        # Jacobian:
-        j_eq = a - a_sum * n / n_sum
+        # Jacobian for equilibrium equations:
+        j_eq = a @ (j_eq_eye - j_eq_ones * n / np.sum(n))
 
         # Residuals from elemental balance:
-        el_sum = np.dot(el_matrix, n)
-        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
+        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}')
 
-        # print(el_sum, element_norm)
-
-        # Jacobian
-        j_ab = weighting * el_matrix * n / el_sum[:, np.newaxis]
+        # Jacobian for elemental balance:
+        j_ab = el_matrix * n / el_sum_norm[:, None]
 
         return (np.hstack([resid_eq, resid_ab]), np.concatenate([j_eq, j_ab], axis=0))
 
-    ret = root(residuals, logn_start, jac=True, tol=1e-10)
-    n = np.exp(np.array(ret['x'], dtype=NDFloat))
-    # print(ret)
+    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)
 
     return n * el_max
 

tests/benchmark_equalibrium.py

@@ -44,6 +44,10 @@ 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

@@ -0,0 +1,25 @@
+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,6 +55,7 @@ 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():
@@ -84,7 +85,7 @@ def test_equilibrium():
             result_values = gp.equilibrium(fl, t, p).array_fractions
 
             print(index, gp.get_solver(), '----')
-            print(molar_comp)
+            print('Species:        ' + ''.join(f"{s:14}" for s in fs.species))
             outp(result_values, 'Under test: ')
             outp(reference_values, 'Reference:  ')
 
@@ -123,3 +124,7 @@ 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(f"---")
+        print("---")
         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_cantera()
+    test_equilibrium_cantera()