SOC examples updated to use new slicing

CI: missing "jobs:" readded

CITATION.cff fixed and CI updated for CITATION.cff validation

.flake8

@@ -14,7 +14,7 @@ exclude =
     dist,
     .conda,
     tests/autogenerated_*,
-    docs/source/_autogenerated
+    docs/source/api
     .venv,
     venv
 

.github/workflows/ci.yml

@@ -27,6 +27,13 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         python -m pip install -e .[dev]
+        if [ "${{ matrix.python-version }}" = "3.13" ]; then
+          python -m pip install cffconvert
+        fi
+
+    - name: Validate CITATION.cff
+      if: ${{ matrix.python-version == '3.13' }}
+      run: cffconvert --validate
 
     - name: Prepare data and compile thermo database
       run: |

.github/workflows/docs.yml

@@ -9,7 +9,7 @@ permissions:
   contents: write
 
 jobs:
-  build-and-deploy:
+  build:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
@@ -41,8 +41,24 @@ jobs:
           rm ./source/*.rst
           make html
           touch ./build/html/.nojekyll
-      - name: Deploy to GitHub Pages
-        uses: JamesIves/github-pages-deploy-action@v4
+          mkdir -p ./build/html/_autogenerated
+          cp ./build/html/api/* ./build/html/_autogenerated/
+      - name: Upload artifact
+        uses: actions/upload-pages-artifact@v3
         with:
-          branch: gh-pages
-          folder: docs/build/html
+          path: docs/build/html
+
+  deploy:
+    needs: build
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read
+      pages: write
+      id-token: write
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    steps:
+      - name: Deploy to GitHub Pages
+        id: deployment
+        uses: actions/deploy-pages@v4

.github/workflows/publish.yml

@@ -10,6 +10,10 @@ jobs:
     name: Build and publish
     runs-on: ubuntu-latest
 
+    environment:
+      name: pypi
+      url: https://pypi.org/project/${{ github.event.repository.name }}/
+
     steps:
       - uses: actions/checkout@v3
 

.gitignore

@@ -9,7 +9,7 @@ __pycache__
 .pytest_cache
 tests/autogenerated_*.py
 docs/build/
-docs/source/_autogenerated/
+docs/source/api/
 venv/
 .venv/
 thermo_data/combined_data.yaml

CITATION.cff

@@ -1,4 +1,5 @@
-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
 authors:
@@ -8,12 +9,11 @@ authors:
     affiliation: "German Aerospace Center (DLR)"
     address: "Linder Höhe"
     city: Köln
-version: 1.1.1
-date-released: "2025-06-17"
+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: ""
-license: MIT License
+license: MIT
 repository-code: "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"
-documentation: "https://dlr-institute-of-future-fuels.github.io/gaspype"

README.md

@@ -1,20 +1,21 @@
 # Gaspype
-The python package provides an performant library for thermodynamic calculations
+The Python package provides a performant library for thermodynamic calculations
 like equilibrium reactions for several hundred gas species and their mixtures -
-written in Python/Numpy.
+written in Python/NumPy.
 
 Species are treated as ideal gases. Therefore the application is limited to moderate
 pressures or high temperature applications.
 
-Its designed with goal to be portable to Numpy-style GPU frameworks like JAX and PyTorch.
+It is designed with goal to be portable to NumPy-style GPU frameworks like JAX and PyTorch.
 
-## Key features
-- Tensor operations to prevent bottlenecks by the python interpreter
-- Immutable types
-- Elegant pythonic interface
-- Great readable and compact syntax when using this package
-- Good usability in Jupyter Notebook
-- High performance for multidimensional fluid arrays
+## Key Features
+
+- Pure Python implementation with NumPy vectorization for high performance
+- Immutable types and comprehensive type hints for reliability
+- Intuitive, Pythonic API for both rapid prototyping and complex multidimensional models
+- Ready for Jupyter Notebook and educational use
+- Designed for future GPU support (JAX, PyTorch)
+- Ships with a comprehensive NASA9-based species database
 
 ## Installation
 Installation with pip:
@@ -52,7 +53,7 @@ mass = fl.get_mass()
 gas_volume = fl.get_v(t=800+273.15, p=1e5)
 ```
 
-The arguments can be provided as numpy-arrays:
+The arguments can be provided as NumPy-arrays:
 
 ``` python
 import numpy as np
@@ -63,7 +64,7 @@ fl.get_density(t=t_range, p=1e5)
 array([0.10122906, 0.09574625, 0.09082685, 0.08638827, 0.08236328])
 ```
 A ```fluid``` object can have multiple compositions. A multidimensional ```fluid``` object
-can be created for example by multiplication with a numpy array:
+can be created for example by multiplication with a NumPy array:
 
 ``` python
 fl2 = gp.fluid({'H2O': 1, 'N2': 2}) + \
@@ -220,8 +221,8 @@ cd gaspype
 It's recommended to setup an venv:
 
 ```bash
-python -m venv venv
-source venv/bin/activate  # On Windows use `venv\Scripts\activate`
+python -m venv .venv
+source .venv/bin/activate  # On Windows use `.venv\Scripts\activate`
 ```
 
 Install the package and dev-dependencies while keeping the package files

docs/source/conf.py

@@ -10,7 +10,7 @@ import os
 import sys
 sys.path.insert(0, os.path.abspath("../src/"))
 
-project = 'gaspype'
+project = 'Gaspype'
 copyright = '2025, Nicolas Kruse'
 author = 'Nicolas Kruse'
 

docs/source/generate_class_list.py

@@ -27,7 +27,7 @@ def write_classes(f: TextIOWrapper, patterns: list[str], module_name: str, title
     write_dochtree(f, title, classes)
 
     for cls in classes:
-        with open(f'docs/source/_autogenerated/{cls}.md', 'w') as f2:
+        with open(f'docs/source/api/{cls}.md', 'w') as f2:
             f2.write(f'# {module_name}.{cls}\n')
             f2.write('```{eval-rst}\n')
             f2.write(f'.. autoclass:: {module_name}.{cls}\n')
@@ -43,7 +43,7 @@ def write_functions(f: TextIOWrapper, patterns: list[str], module_name: str, tit
     module = importlib.import_module(module_name)
 
     functions = [
-        name for name, obj in inspect.getmembers(module, inspect.isfunction)
+        name for name, _ in inspect.getmembers(module, inspect.isfunction)
         if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude))
     ]
 
@@ -54,7 +54,7 @@ def write_functions(f: TextIOWrapper, patterns: list[str], module_name: str, tit
 
     for func in functions:
         if not func.startswith('_'):
-            with open(f'docs/source/_autogenerated/{func}.md', 'w') as f2:
+            with open(f'docs/source/api/{func}.md', 'w') as f2:
                 f2.write(f'# {module_name}.{func}\n')
                 f2.write('```{eval-rst}\n')
                 f2.write(f'.. autofunction:: {module_name}.{func}\n')
@@ -74,9 +74,9 @@ def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
 
 if __name__ == "__main__":
     # Ensure the output directory exists
-    os.makedirs('docs/source/_autogenerated', exist_ok=True)
+    os.makedirs('docs/source/api', exist_ok=True)
 
-    with open('docs/source/_autogenerated/index.md', 'w') as f:
+    with open('docs/source/api/index.md', 'w') as f:
         f.write('# Classes and functions\n\n')
 
         write_classes(f, ['*'], 'gaspype', title='Classes')

docs/source/index.md

@@ -1,8 +1,9 @@
 ```{toctree}
 :maxdepth: 1
 :hidden:
-_autogenerated/index
-_autogenerated/examples
+api/index
+api/examples
+repo
 ```
 
 ```{include} ../../README.md

docs/source/render_examples.py

@@ -51,7 +51,7 @@ def render_examples(filter: str, example_file: str):
         f.write('## Download Jupyter Notebooks\n\n')
         for path, name in zip(files, names):
             if name.lower() != 'readme':
-                run_rendering(path, 'docs/source/_autogenerated')
+                run_rendering(path, 'docs/source/api')
                 notebook = name + '.ipynb'
                 f.write(f'- [{notebook}]({notebook})\n\n')
 
@@ -63,4 +63,4 @@ def render_examples(filter: str, example_file: str):
 
 
 if __name__ == "__main__":
-    render_examples('examples/*.md', 'docs/source/_autogenerated/examples.md')
+    render_examples('examples/*.md', 'docs/source/api/examples.md')

docs/source/repo.md

@@ -0,0 +1,3 @@
+# Code repository
+
+Code repository is on GitHub: [github.com/DLR-Institute-of-Future-Fuels/gaspype](https://github.com/DLR-Institute-of-Future-Fuels/gaspype).

examples/README.md

@@ -11,11 +11,11 @@ The conversion is done like the following automated by the
 [docs/source/render_examples.py](../docs/source/render_examples.py) script:
 ``` bash
 # Converting markdown with code sections to Jupyter Notebook and run it:
-notedown examples/soec_methane.md --to notebook --output docs/source/_autogenerated/soec_methane.ipynb --run
+notedown examples/soec_methane.md --to notebook --output docs/source/api/soec_methane.ipynb --run
 
 # Converting the Jupyter Notebook to Markdown and a folder with image
-# files placed in docs/source/_autogenerated/:
-jupyter nbconvert --to markdown docs/source/_autogenerated/soec_methane.ipynb --output soec_methane.md
+# files placed in docs/source/api/:
+jupyter nbconvert --to markdown docs/source/api/soec_methane.ipynb --output soec_methane.md
 ```
 
 A new example Markdown file can be created from a Jupyter Notebook running

examples/hydrogen_production_thermodynamics.md

@@ -0,0 +1,45 @@
+# Thermodynamics of Hydrogen Production
+
+The minimal energy required to produce hydrogen from liquid water is given by the
+Higher Heating Value (HHV). The HHV is the sum of the difference
+between the enthalpies of products and educts (LHV: Lower Heating Value) and
+the Heat of Evaporation for water.
+
+```python
+import gaspype as gp
+
+lhv = gp.fluid({'H2': 1, 'O2': 1/2, 'H2O': -1}).get_H(25 + 273.15)
+dh_v = 43990  # J/mol  (heat of evaporation for water @ 25 °C)
+
+hhv = lhv + dh_v
+print(f'LHV: {lhv/1e3:.1f} kJ/mol')
+print(f'HHV: {hhv/1e3:.1f} kJ/mol')
+```
+
+Thermodynamics also defines which part of the energy must be
+provided as work (e.g., electric power) and which part can be supplied
+as heat. This depends on temperature and pressure. For generating 1 bar
+of hydrogen the temperature dependency can be calculated as follows:
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+
+t = np.linspace(0, 2000, 128)  # 0 to 2000 °C
+p = 1e5  # Pa (=1 bar)
+
+g_products = gp.fluid({'H2': 1, 'O2': 1/2, 'H2O': 0}).get_G(t + 273.15, p)
+g_educts =   gp.fluid({'H2': 0, 'O2': 0,   'H2O': 1}).get_G(t + 273.15, p)
+
+work = g_products - g_educts  # J/mol
+heat = lhv - work  # J/mol
+
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("Temperature / °C")
+ax.set_ylabel("Energy / kWh/kg")
+k = 1e-3 / 3600 / 0.002  # Conversion factor from J/mol to kWh/kg for hydrogen
+ax.stackplot(t, k * work, k * heat, k * dh_v * np.ones_like(t))
+```
+
+Green is the heat of evaporation, orange the additional heat provided at
+the given temperature and blue the work.

examples/methane_mixtures.md

@@ -25,8 +25,8 @@ Equilibrium calculation for methane steam mixtures:
 ```python
 ratio = np.linspace(0.01, 1.5, num=64)
 
-el = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
-equilibrium_h2o = gp.equilibrium(el, t, p)
+fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
+equilibrium_h2o = gp.equilibrium(fl, t, p)
 ```
 
 
@@ -44,8 +44,8 @@ Equilibrium calculation for methane CO2 mixtures:
 
 
 ```python
-el = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
-equilibrium_co2 = gp.equilibrium(el, t, p)
+fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
+equilibrium_co2 = gp.equilibrium(fl, t, p)
 ```
 
 

examples/soec_syngas.md

@@ -1,8 +1,10 @@
-# SOEC with Methane
+# SOEC Co-Electrolysis
 
-This example shows a 1D isothermal SOEC (Solid oxide electrolyzer cell) model.
+This example shows a 1D isothermal SOEC (Solid oxide electrolyzer cell) model for
+converting carbon dioxide and steam into syngas.
 
-The operating parameters chosen here are not necessary realistic
+The operating parameters chosen here are not necessarily realistic. For example,
+a utilization of 0.95 causes issues with the formation of solid carbon.
 
 ```python
 import gaspype as gp
@@ -11,26 +13,26 @@ import numpy as np
 import matplotlib.pyplot as plt
 ```
 
-Calculation of the local equilibrium compositions on the fuel and air
-side in counter flow along the fuel flow direction:
+Calculate equilibrium compositions for fuel and air sides in counter flow
+along the fuel flow direction:
 ```python
-fuel_utilization = 0.90
-air_utilization = 0.5
-t = 800 + 273.15 #K
+utilization = 0.95
+air_dilution = 0.2
+t = 700 + 273.15  # K
 p = 1e5  # Pa
 
 fs = gp.fluid_system('H2, H2O, O2, CH4, CO, CO2')
-feed_fuel = gp.fluid({'CH4': 1, 'H2O': 0.1}, fs)
+feed_fuel = gp.fluid({'H2O': 2, 'CO2': 1}, fs)
 
-o2_full_conv = np.sum(gp.elements(feed_fuel)[['H', 'C' ,'O']] * [1/4, 1, -1/2])
+o2_full_conv = np.sum(gp.elements(feed_fuel).get_n(['C' ,'O']) * [-1/2, 1/2])
 
-feed_air = gp.fluid({'O2': 1, 'N2': 4}) * o2_full_conv / air_utilization
+feed_air = gp.fluid({'O2': 1, 'N2': 4}) * o2_full_conv * utilization * air_dilution
 
-conversion = np.linspace(0, fuel_utilization, 32)
+conversion = np.linspace(0, utilization, 128)
 perm_oxygen = o2_full_conv * conversion * gp.fluid({'O2': 1})
 
-fuel_side = gp.equilibrium(feed_fuel + perm_oxygen, t, p)
-air_side  = gp.equilibrium(feed_air  - perm_oxygen, t, p)
+fuel_side = gp.equilibrium(feed_fuel - perm_oxygen, t, p)
+air_side  = gp.equilibrium(feed_air  + perm_oxygen, t, p)
 ```
 
 Plot compositions of the fuel and air side:
@@ -64,13 +66,20 @@ ax.plot(conversion, np.stack([o2_fuel_side, o2_air_side], axis=1), '-')
 ax.legend(['o2_fuel_side', 'o2_air_side'])
 ```
 
+The high oxygen partial pressure at the inlet is in reality lower.
+The assumption that gas inter-diffusion in the flow direction is slower
+than the gas velocity does not hold at this very high gradient. However
+often the oxygen partial pressure is still to high to prevent oxidation of the
+cell/electrode. This can be effectively prevented by recycling small amounts of
+the hydrogen riche output gas.
+
 Calculation of the local nernst potential between fuel and air side:
 ```python
 z_O2 = 4
 nernst_voltage = R*t / (z_O2*F) * np.log(o2_air_side/o2_fuel_side)
 ```
 
-#Plot nernst potential:
+Plot nernst potential:
 ```python
 fig, ax = plt.subplots()
 ax.set_xlabel("Conversion")
@@ -89,16 +98,16 @@ To calculate the local current density (**node_current**) as well
 as the total cell current (**terminal_current**) the (relative)
 physical distance between the nodes (**dz**) must be calculated:
 ```python
-cell_voltage = 0.77 #V
+cell_voltage = 1.3  # V
 ASR = 0.2  # Ohm*cm²
 
-node_current = (nernst_voltage - cell_voltage) / ASR  # mA/cm² (Current density at each node)
+node_current = (nernst_voltage - cell_voltage) / ASR  # A/cm² (Current density at each node)
 
-current = (node_current[1:] + node_current[:-1]) / 2  # mA/cm² (Average current density between the nodes)
+current = (node_current[1:] + node_current[:-1]) / 2  # A/cm² (Average current density between the nodes)
 
 dz = 1/current / np.sum(1/current)  # Relative distance between each node
 
-terminal_current = np.sum(current * dz) # mA/cm² (Total cell current per cell area)
+terminal_current = np.sum(current * dz) # A/cm² (Total cell current per cell area)
 
 print(f'Terminal current: {terminal_current:.2f} A/cm²')
 ```

examples/sofc_methane.md

@@ -0,0 +1,141 @@
+# SOFC with Methane
+
+This example shows a 1D isothermal SOFC (Solid oxide fuel cell) model.
+
+The operating parameters chosen here are not necessarily realistic due to
+constraints not included in this model. Under the example condition the
+formation of solid carbon is for example very likely. A pre-reforming step
+is typically applied when operating on methane.
+
+```python
+import gaspype as gp
+from gaspype.constants import R, F
+import numpy as np
+import matplotlib.pyplot as plt
+```
+
+Calculate equilibrium compositions for fuel and air sides
+in counter flow along the fuel flow direction:
+```python
+fuel_utilization = 0.90
+air_utilization = 0.5
+t = 800 + 273.15  # K
+p = 1e5  # Pa
+
+fs = gp.fluid_system('H2, H2O, O2, CH4, CO, CO2')
+feed_fuel = gp.fluid({'CH4': 1, 'H2O': 0.1}, fs)
+
+o2_full_conv = np.sum(gp.elements(feed_fuel).get_n(['H', 'C' ,'O']) * [1/4, 1, -1/2])
+
+feed_air = gp.fluid({'O2': 1, 'N2': 4}) * o2_full_conv * fuel_utilization / air_utilization
+
+conversion = np.linspace(0, fuel_utilization, 32)
+perm_oxygen = o2_full_conv * conversion * gp.fluid({'O2': 1})
+
+fuel_side = gp.equilibrium(feed_fuel + perm_oxygen, t, p)
+air_side  = gp.equilibrium(feed_air  - perm_oxygen, t, p)
+```
+
+Plot compositions of the fuel and air side:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, fuel_side.get_x(), '-')
+ax.legend(fuel_side.species)
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Molar fraction")
+ax.plot(conversion, air_side.get_x(), '-')
+ax.legend(air_side.species)
+```
+
+Calculation of the oxygen partial pressures:
+```python
+o2_fuel_side = gp.oxygen_partial_pressure(fuel_side, t, p)
+o2_air_side = air_side.get_x('O2') * p
+```
+
+Plot oxygen partial pressures:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Oxygen partial pressure / Pa")
+ax.set_yscale('log')
+ax.plot(conversion, np.stack([o2_fuel_side, o2_air_side], axis=1), '-')
+ax.legend(['o2_fuel_side', 'o2_air_side'])
+```
+
+Calculation of the local nernst potential between fuel and air side:
+```python
+z_O2 = 4  # Number of electrons transferred per O2 molecule
+nernst_voltage = R*t / (z_O2*F) * np.log(o2_air_side/o2_fuel_side)
+```
+
+Plot nernst potential:
+```python
+fig, ax = plt.subplots()
+ax.set_xlabel("Conversion")
+ax.set_ylabel("Voltage / V")
+ax.plot(conversion, nernst_voltage, '-')
+print(np.min(nernst_voltage))
+```
+
+The model uses between each node a constant conversion. Because
+current density depends strongly on the position along the cell
+the constant conversion does not relate to a constant distance.
+
+![Alt text](../../media/soc_inverted.svg)
+
+To calculate the local current density (**node_current**) as well
+as the total cell current (**terminal_current**) the (relative)
+physical distance between the nodes (**dz**) must be calculated:
+```python
+cell_voltage = 0.77  # V
+ASR = 0.2  # Ohm*cm²
+
+node_current = (nernst_voltage - cell_voltage) / ASR  # A/cm² (Current density at each node)
+
+current = (node_current[1:] + node_current[:-1]) / 2  # A/cm² (Average current density between the nodes)
+
+dz = 1/current / np.sum(1/current)  # Relative distance between each node
+
+terminal_current = np.sum(current * dz) # A/cm² (Total cell current per cell area)
+
+print(f'Terminal current: {terminal_current:.2f} A/cm²')
+```
+
+Plot the local current density:
+```python
+z_position = np.concatenate([[0], np.cumsum(dz)])  # Relative position of each node
+
+fig, ax = plt.subplots()
+ax.set_xlabel("Relative cell position")
+ax.set_ylabel("Current density / A/cm²")
+ax.plot(z_position, node_current, '-')
+```
+
+Based on the cell current and voltage the energy balance can be calculated.
+In the following the electric cell output power (often referred to as "DC power")
+and lower heating value (LHV) are calculated. The numbers here are per cell area for
+being cell and stack size independent. The quotient of both is often referred to as
+LHV based DC efficiency.
+
+```python
+dc_power = cell_voltage * terminal_current  # W/cm²
+print(f"DC power: {dc_power:.2f} W/cm²")
+
+lhv = gp.fluid({'CH4': 1, 'H2O': -2, 'CO2': -1}).get_H(25 + 273.15)  # J/mol (LHV of methane)
+
+# LHV based chemical input power:
+lhv_power = lhv * terminal_current / (2 * z_O2 * F)  # W/cm² (two O2 per CH4 for full oxidation)
+efficiency = dc_power / lhv_power
+print(f"LHV based DC efficiency: {efficiency*100:.1f} %")
+
+# Or by shortening the therms:
+lhv_voltage = lhv / (2 * z_O2 * F)  # V
+print(f"LHV voltage: {lhv_voltage:.2f} V")
+efficiency = cell_voltage / lhv_voltage  # LHV based DC efficiency
+print(f"LHV based DC efficiency: {efficiency*100:.1f} %")
+```

examples/sulfur_oxygen_equalibrium.md

@@ -32,8 +32,11 @@ el = gp.elements({'S': 1}, fs) + oxygen_ratio * gp.elements({'O': 1}, fs)
 
 composition = gp.equilibrium(el, 800+273.15, 1e4)
 
-plt.plot(oxygen_ratio, composition.get_x())
-plt.legend(composition.species)
+fig, ax = plt.subplots()
+ax.set_xlabel("Oxygen to sulfur ratio")
+ax.set_ylabel("Molar fraction")
+ax.plot(oxygen_ratio, composition.get_x(), '-')
+ax.legend(composition.species)
 ```
 
 Calculation of the molar equilibrium fractions for sulfur and oxygen depending on temperature in °C:
@@ -47,6 +50,9 @@ el = gp.elements({'S': 1, 'O':2.5}, fs)
 t_range = np.linspace(500, 1300, num=32)
 composition = gp.equilibrium(el, t_range+273.15, 1e4)
 
-plt.plot(t_range, composition.get_x())
-plt.legend(composition.species)
+fig, ax = plt.subplots()
+ax.set_xlabel("Temperature / °C")
+ax.set_ylabel("Molar fraction")
+ax.plot(t_range, composition.get_x(), '-')
+ax.legend(composition.species)
 ```

examples/thermodynamic_carbon_activity.md

@@ -1,7 +1,7 @@
 # Carbon Activity
 
-This example shows the equilibrium calculation for solid carbon.
-
+This example shows the calculation of the carbon activity for methane mixtures
+in thermodynamic equilibrium.
 
 ```python
 import gaspype as gp
@@ -66,3 +66,25 @@ ax.plot(ratio, carbon_activity.T)
 ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed')
 ax.legend([f'{tc} °C' for tc in t_range])
 ```
+
+Let's do the equilibrium calculation for methane CO2 mixtures as well:
+
+
+```python
+fl_co2 = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
+carbon_activity_co2 = np.vstack([partial_c_activity(fl_co2, tc + 273.15, p) for tc in t_range])
+```
+
+And plot carbon activities over the CO2 to CH4 ratio:
+
+
+```python
+fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
+ax.set_xlabel("CO2/CH4")
+ax.set_ylabel("carbon activity")
+ax.set_ylim(1e-1, 1e3)
+ax.set_yscale('log')
+ax.plot(ratio, carbon_activity_co2.T)
+ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed')
+ax.legend([f'{tc} °C' for tc in t_range])
+```

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "gaspype"
-version = "1.1.1"
+version = "1.1.3"
 authors = [
   { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
 ]
@@ -18,9 +18,10 @@ dependencies = [
 ]
 
 [project.urls]
-Homepage = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"
+Homepage = "https://dlr-institute-of-future-fuels.github.io/gaspype/"
+documentation = "https://dlr-institute-of-future-fuels.github.io/gaspype/api/"
+Repository = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype.git"
 Issues = "https://github.com/DLR-Institute-of-Future-Fuels/gaspype/issues"
-documentation = "https://dlr-institute-of-future-fuels.github.io/gaspype/"
 
 [build-system]
 requires = ["setuptools>=61.0", "wheel"]
@@ -30,7 +31,7 @@ build-backend = "setuptools.build_meta"
 where = ["src"]
 
 [tool.setuptools.package-data]
-gaspype = ["data/therm_data.bin"]
+gaspype = ["data/therm_data.bin", "py.typed"]
 
 [project.optional-dependencies]
 dev = [

src/gaspype/__init__.py

@@ -13,8 +13,8 @@ Example usage:
 """
 
 from ._main import species, fluid_system, fluid, elements
-from ._operations import set_solver, get_solver, equilibrium
 from ._operations import stack, concat, carbon_activity, oxygen_partial_pressure
+from ._solver import set_solver, get_solver, equilibrium
 
 __all__ = [
     'species', 'fluid_system', 'fluid', 'elements',

src/gaspype/_main.py

@@ -7,7 +7,7 @@ from gaspype._phys_data import atomic_weights, db_reader
 import re
 import pkgutil
 from .constants import R, epsy, p0
-from .typing import FloatArray, NDFloat, Shape
+from .typing import FloatArray, NDFloat, Shape, ArrayIndices
 
 T = TypeVar('T', 'fluid', 'elements')
 
@@ -120,9 +120,9 @@ class fluid_system:
             # print(f'Warning: temperature ({T}) out of range for {s}')
 
         self.elements: list[str] = sorted(list(set(k for ac in element_compositions for k in ac.keys())))
-        self.array_species_elements = np.array([[ec[el] if el in ec else 0.0 for el in self.elements] for ec in element_compositions])
+        self.array_species_elements: FloatArray = np.array([[ec[el] if el in ec else 0.0 for el in self.elements] for ec in element_compositions])
 
-        self.array_atomic_mass = np.array([atomic_weights[el] for el in self.elements]) * 1e-3  # kg/mol
+        self.array_atomic_mass: FloatArray = np.array([atomic_weights[el] for el in self.elements]) * 1e-3  # kg/mol
         self.array_molar_mass: FloatArray = np.sum(self.array_atomic_mass * self.array_species_elements, axis=-1)  # kg/mol
 
         self.array_stoichiometric_coefficients: FloatArray = np.array(null_space(self.array_species_elements.T), dtype=NDFloat).T
@@ -483,6 +483,28 @@ class fluid:
             assert set(species) <= set(self.fs.species), f'Species {", ".join([s for s in species if s not in self.fs.species])} is/are not part of the fluid system'
             return self.array_fractions[..., [self.fs.species.index(k) for k in species]]
 
+    def get_n(self, species: str | list[str] | None = None) -> FloatArray:
+        """Get molar amount of fluid species
+
+        Args:
+            species: A single species name, a list of species names or None for
+                returning the amount of all species
+
+        Returns:
+            Returns an array of floats with the molar amount of the species.
+            If the a single species name is provided the return float array has
+            the same dimensions as the fluid type. If a list or None is provided
+            the return array has an additional dimension for the species.
+        """
+        if not species:
+            return self.array_composition
+        elif isinstance(species, str):
+            assert species in self.fs.species, f'Species {species} is not part of the fluid system'
+            return self.array_composition[..., self.fs.species.index(species)]
+        else:
+            assert set(species) <= set(self.fs.species), f'Species {", ".join([s for s in species if s not in self.fs.species])} is/are not part of the fluid system'
+            return self.array_composition[..., [self.fs.species.index(k) for k in species]]
+
     def __add__(self, other: T) -> T:
         return array_operation(self, other, np.add)
 
@@ -510,16 +532,21 @@ class fluid:
     # def __array__(self) -> FloatArray:
     #     return self.array_composition
 
-    def __getitem__(self, key: str | int | list[str] | list[int] | slice) -> FloatArray:
+    @overload
+    def __getitem__(self, key: str) -> FloatArray:
+        pass
+
+    @overload
+    def __getitem__(self, key: ArrayIndices) -> 'fluid':
+        pass
+
+    def __getitem__(self, key: str | ArrayIndices) -> Any:
         if isinstance(key, str):
             assert key in self.fs.species, f'Species {key} is not part of the fluid system'
             return self.array_composition[..., self.fs.species.index(key)]
-        elif isinstance(key, (slice, int)):
-            return self.array_composition[..., key]
         else:
-            mset = set(self.fs.species) | set(range(len(self.fs.species)))
-            assert set(key) <= mset, f'Species {", ".join([str(s) for s in key if s not in mset])} is/are not part of the fluid system'
-            return self.array_composition[..., [self.fs.species.index(k) if isinstance(k, str) else k for k in key]]
+            key_tuple = key if isinstance(key, tuple) else (key,)
+            return fluid(self.array_composition[(*key_tuple, slice(None))], self.fs)
 
     def __iter__(self) -> Iterator[dict[str, float]]:
         assert len(self.shape) < 2, 'Cannot iterate over species with more than one dimension'
@@ -614,6 +641,28 @@ class elements:
         """
         return np.sum(self.array_elemental_composition * self.fs.array_atomic_mass, axis=-1, dtype=NDFloat)
 
+    def get_n(self, elemental_species: str | list[str] | None = None) -> FloatArray:
+        """Get molar amount of elements
+
+        Args:
+            elemental_species: A single element name, a list of element names or None for
+                returning the amount of all element
+
+        Returns:
+            Returns an array of floats with the molar amount of the elements.
+            If the a single element name is provided the return float array has
+            the same dimensions as the fluid type. If a list or None is provided
+            the return array has an additional dimension for the elements.
+        """
+        if not elemental_species:
+            return self.array_elemental_composition
+        elif isinstance(elemental_species, str):
+            assert elemental_species in self.fs.elements, f'Element {elemental_species} is not part of the fluid system'
+            return self.array_elemental_composition[..., self.fs.elements.index(elemental_species)]
+        else:
+            assert set(elemental_species) <= set(self.fs.elements), f'Elements {", ".join([s for s in elemental_species if s not in self.fs.elements])} is/are not part of the fluid system'
+            return self.array_elemental_composition[..., [self.fs.elements.index(k) for k in elemental_species]]
+
     def __add__(self, other: 'fluid | elements') -> 'elements':
         return array_operation(self, other, np.add)
 
@@ -639,16 +688,21 @@ class elements:
     def __array__(self) -> FloatArray:
         return self.array_elemental_composition
 
-    def __getitem__(self, key: str | int | list[str] | list[int] | slice) -> FloatArray:
+    @overload
+    def __getitem__(self, key: str) -> FloatArray:
+        pass
+
+    @overload
+    def __getitem__(self, key: ArrayIndices) -> 'elements':
+        pass
+
+    def __getitem__(self, key: str | ArrayIndices) -> Any:
         if isinstance(key, str):
             assert key in self.fs.elements, f'Element {key} is not part of the fluid system'
             return self.array_elemental_composition[..., self.fs.elements.index(key)]
-        elif isinstance(key, (slice, int)):
-            return self.array_elemental_composition[..., key]
         else:
-            mset = set(self.fs.elements) | set(range(len(self.fs.elements)))
-            assert set(key) <= mset, f'Elements {", ".join([str(s) for s in key if s not in mset])} is/are not part of the fluid system'
-            return self.array_elemental_composition[..., [self.fs.elements.index(k) if isinstance(k, str) else k for k in key]]
+            key_tuple = key if isinstance(key, tuple) else (key,)
+            return elements(self.array_elemental_composition[(*key_tuple, slice(None))], self.fs)
 
     def __iter__(self) -> Iterator[dict[str, float]]:
         assert len(self.shape) < 2, 'Cannot iterate over elements with more than one dimension'

src/gaspype/_operations.py

@@ -1,48 +1,9 @@
-from typing import Literal, Any
 from math import exp
-from scipy.optimize import minimize, root
 import numpy as np
-from ._main import T, elements, fluid, fluid_system
-from .typing import NDFloat, FloatArray
-from .constants import p0, epsy, R
-
-
-def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
-    """
-    Select a solver for chemical equilibrium.
-
-    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.
-
-        - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
-          the elemental composition constant using the SLSQP implementation of scipy
-
-    Args:
-        solver: Name of the solver
-    """
-    global _equilibrium_solver
-    if solver == 'gibs minimization':
-        _equilibrium_solver = equilibrium_gmin
-    elif solver == 'system of equations':
-        _equilibrium_solver = equilibrium_eq
-    else:
-        raise ValueError('Unknown solver')
-
-
-def get_solver() -> Literal['gibs minimization', 'system of equations']:
-    """Returns the selected solver name.
-
-    Returns:
-        Solver name
-    """
-    if _equilibrium_solver == equilibrium_gmin:
-        return 'gibs minimization'
-    else:
-        assert _equilibrium_solver == equilibrium_eq
-        return 'system of equations'
+from ._main import T, elements, fluid
+from .typing import FloatArray
+from .constants import p0, R
+from ._solver import equilibrium
 
 
 def stack(arrays: list[T], axis: int = 0) -> T:
@@ -81,111 +42,6 @@ def concat(arrays: list[T], axis: int = 0) -> T:
         axis=axis), a0.fs)
 
 
-def equilibrium_gmin(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
-    """Calculate the equilibrium composition of a fluid based on minimizing the Gibbs free energy"""
-    def element_balance(n: FloatArray, fs: fluid_system, ref: FloatArray) -> FloatArray:
-        return np.dot(n, fs.array_species_elements) - ref  # type: ignore
-
-    def gibbs_rt(n: FloatArray, grt: FloatArray, p_rel: float):  # type: ignore
-        # Calculate G/(R*T)
-        return np.sum(n * (grt + np.log(p_rel * n / np.sum(n) + epsy)))
-
-    cons: dict[str, Any] = {'type': 'eq', 'fun': element_balance, 'args': [fs, element_composition]}
-    bnds = [(0, None) for _ in fs.species]
-    grt = fs.get_species_g_rt(t)
-    p_rel = p / p0
-
-    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)  # type: ignore
-
-    return sol
-
-
-def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
-    """Calculate the equilibrium composition of a fluid based on equilibrium equations"""
-    el_max = np.max(element_composition)
-    element_norm = element_composition / el_max
-
-    a = fs.array_stoichiometric_coefficients
-    a_sum = np.sum(a)
-    el_matrix = fs.array_species_elements.T
-
-    # Log equilibrium constants for each reaction equation
-    b = -np.sum(fs.get_species_g_rt(t) * a, axis=1)
-
-    # Pressure corrected log equilibrium constants
-    bp = b - np.sum(a * np.log(p / p0), axis=1)
-
-    logn_start = np.ones(el_matrix.shape[1]) * 0.1
-
-    def residuals(logn: FloatArray):  # type: ignore
-        n = np.exp(logn)
-        n_sum = np.sum(n)
-
-        # Residuals from equilibrium equations:
-        eq_resid = np.dot(a, logn - np.log(n_sum)) - bp
-
-        # Derivative:
-        j_eq = a - a_sum * n / n_sum
-
-        # Residuals from elemental balance:
-        el_error = np.dot(el_matrix, n) - element_norm
-        ab_resid = np.log1p(el_error)
-
-        # Derivative:
-        j_ab = el_matrix * n / np.expand_dims(el_error + 1, axis=1)
-
-        return (np.hstack([eq_resid, ab_resid]), np.concatenate([j_eq, j_ab], axis=0))
-
-    ret = root(residuals, logn_start, jac=True, tol=1e-30)
-    n = np.exp(np.array(ret['x'], dtype=NDFloat))
-
-    return n * el_max
-
-
-def equilibrium(f: fluid | elements, t: float | FloatArray, p: float = 1e5) -> fluid:
-    """Calculate the equilibrium composition of a fluid at a given temperature and pressure"
-
-    Args:
-        f: Fluid or elements object
-        t: Temperature in Kelvin
-        p: Pressure in Pascal
-
-    Returns:
-        A new fluid object with the equilibrium composition
-    """
-    assert isinstance(f, (fluid, elements)), 'Argument f must be a fluid or elements'
-    m_shape: int = f.fs.array_stoichiometric_coefficients.shape[0]
-    if isinstance(f, fluid):
-        if not m_shape:
-            return f
-    else:
-        if not m_shape:
-            def linalg_lstsq(array_elemental_composition: FloatArray, matrix: FloatArray) -> Any:
-                # TODO: np.dot(np.linalg.pinv(a), b) is eqivalent to lstsq(a, b).
-                # the constant np.linalg.pinv(a) can be precomputed for each fs.
-                return np.dot(np.linalg.pinv(matrix), array_elemental_composition)
-
-            # print('-->', f.array_elemental_composition.shape, f.fs.array_species_elements.transpose().shape)
-            composition = np.apply_along_axis(linalg_lstsq, -1, f.array_elemental_composition, f.fs.array_species_elements.transpose())
-            return fluid(composition, f.fs)
-
-    assert np.min(f.array_elemental_composition) >= 0, 'Input element fractions must be 0 or positive'
-    if isinstance(t, np.ndarray):
-        assert f.shape == tuple(), 'Multidimensional temperature can currently only used for 0D fluids'
-        t_composition = np.zeros(t.shape + (f.fs.array_species_elements.shape[0],))
-        for t_index in np.ndindex(t.shape):
-            t_composition[t_index] = _equilibrium_solver(f.fs, f.array_elemental_composition, float(t[t_index]), p)
-        return fluid(t_composition, f.fs)
-    else:
-        composition = np.ones(f.shape + (len(f.fs.species),), dtype=float)
-        for index in np.ndindex(f.shape):
-            # print(composition.shape, index, _equilibrium(f.fs, f._element_composition[index], t, p))
-            composition[index] = _equilibrium_solver(f.fs, f.array_elemental_composition[index], t, p)
-        return fluid(composition, f.fs)
-
-
 def carbon_activity(f: fluid | elements, t: float, p: float) -> float:
     """Calculate the activity of carbon in a fluid at a given temperature and pressure.
     At a value of 1 the fluid is in equilibrium with solid graphite. At a value > 1
@@ -292,6 +148,3 @@ def oxygen_partial_pressure(f: fluid | elements, t: float, p: float) -> FloatArr
         return np.apply_along_axis(get_oxygen, -1, x)
     else:
         return get_oxygen(x)
-
-
-_equilibrium_solver = equilibrium_eq

src/gaspype/_solver.py

@@ -0,0 +1,167 @@
+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
+
+
+def set_solver(solver: Literal['gibs minimization', 'system of equations']) -> None:
+    """
+    Select a solver for chemical equilibrium.
+
+    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.
+
+        - **gibs minimization**: Minimizes the total Gibbs Enthalpy while keeping
+          the elemental composition constant using the SLSQP implementation of scipy
+
+    Args:
+        solver: Name of the solver
+    """
+    global _equilibrium_solver
+    if solver == 'gibs minimization':
+        _equilibrium_solver = equilibrium_gmin
+    elif solver == 'system of equations':
+        _equilibrium_solver = equilibrium_eq
+    else:
+        raise ValueError('Unknown solver')
+
+
+def get_solver() -> Literal['gibs minimization', 'system of equations']:
+    """Returns the selected solver name.
+
+    Returns:
+        Solver name
+    """
+    if _equilibrium_solver == equilibrium_gmin:
+        return 'gibs minimization'
+    else:
+        assert _equilibrium_solver == equilibrium_eq
+        return 'system of equations'
+
+
+def equilibrium_gmin(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
+    """Calculate the equilibrium composition of a fluid based on minimizing the Gibbs free energy"""
+    def element_balance(n: FloatArray, fs: fluid_system, ref: FloatArray) -> FloatArray:
+        return np.dot(n, fs.array_species_elements) - ref  # type: ignore
+
+    def gibbs_rt(n: FloatArray, grt: FloatArray, p_rel: float):  # type: ignore
+        # Calculate G/(R*T)
+        return np.sum(n * (grt + np.log(p_rel * n / np.sum(n) + epsy)))
+
+    cons: dict[str, Any] = {'type': 'eq', 'fun': element_balance, 'args': [fs, element_composition]}
+    bnds = [(0, None) for _ in fs.species]
+    grt = fs.get_species_g_rt(t)
+    p_rel = p / p0
+
+    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)  # type: ignore
+
+    return sol
+
+
+def equilibrium_eq(fs: fluid_system, element_composition: FloatArray, t: float, p: float) -> FloatArray:
+    """Calculate the equilibrium composition of a fluid based on equilibrium equations"""
+    el_max = np.max(element_composition)
+    element_norm = element_composition / el_max
+    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
+    b = -np.sum(fs.get_species_g_rt(t) * a, axis=1)
+
+    # Pressure corrected log equilibrium constants
+    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)
+
+    # global count
+    # count = 0
+
+    weighting = 100
+
+    def residuals(logn: FloatArray) -> tuple[FloatArray, FloatArray]:
+        # global count
+        # count += 1
+        # print('------', count)
+        # assert count < 100
+        n = np.exp(logn)
+        n_sum = np.sum(n)
+
+        # Residuals from equilibrium equations:
+        resid_eq = np.dot(a, logn - np.log(n_sum)) - bp
+
+        # Jacobian:
+        j_eq = a - a_sum * n / n_sum
+
+        # Residuals from elemental balance:
+        el_sum = np.dot(el_matrix, n)
+        resid_ab = weighting * (np.log(el_sum) - element_norm_log)
+
+        # 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))
+
+    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
+
+
+def equilibrium(f: fluid | elements, t: float | FloatArray, p: float = 1e5) -> fluid:
+    """Calculate the isobaric equilibrium composition of a fluid at a given temperature and pressure"
+
+    Args:
+        f: Fluid or elements object
+        t: Temperature in Kelvin
+        p: Pressure in Pascal
+
+    Returns:
+        A new fluid object with the equilibrium composition
+    """
+    assert isinstance(f, (fluid, elements)), 'Argument f must be a fluid or elements'
+    m_shape: int = f.fs.array_stoichiometric_coefficients.shape[0]
+    if isinstance(f, fluid):
+        if not m_shape:
+            return f
+    else:
+        if not m_shape:
+            def linalg_lstsq(array_elemental_composition: FloatArray, matrix: FloatArray) -> Any:
+                # TODO: np.dot(np.linalg.pinv(a), b) is eqivalent to lstsq(a, b).
+                # the constant np.linalg.pinv(a) can be precomputed for each fs.
+                return np.dot(np.linalg.pinv(matrix), array_elemental_composition)
+
+            # print('-->', f.array_elemental_composition.shape, f.fs.array_species_elements.transpose().shape)
+            composition = np.apply_along_axis(linalg_lstsq, -1, f.array_elemental_composition, f.fs.array_species_elements.transpose())
+            return fluid(composition, f.fs)
+
+    assert np.min(f.array_elemental_composition) >= 0, 'Input element fractions must be 0 or positive'
+    if isinstance(t, np.ndarray):
+        assert f.shape == tuple(), 'Multidimensional temperature can currently only used for 0D fluids'
+        t_composition = np.zeros(t.shape + (f.fs.array_species_elements.shape[0],))
+        for t_index in np.ndindex(t.shape):
+            t_composition[t_index] = _equilibrium_solver(f.fs, f.array_elemental_composition, float(t[t_index]), p)
+        return fluid(t_composition, f.fs)
+    else:
+        composition = np.ones(f.shape + (len(f.fs.species),), dtype=float)
+        for index in np.ndindex(f.shape):
+            # print(composition.shape, index, _equilibrium(f.fs, f._element_composition[index], t, p))
+            composition[index] = _equilibrium_solver(f.fs, f.array_elemental_composition[index], t, p)
+        return fluid(composition, f.fs)
+
+
+_equilibrium_solver = equilibrium_eq

src/gaspype/constants.py

@@ -21,4 +21,4 @@ p0 = 1e5  # Pa
 t0 = 273.15 + 25  # K
 p_atm = 101325  # Pa
 
-epsy = 1e-18
+epsy = 1e-30

src/gaspype/typing.py

@@ -1,6 +1,10 @@
 from numpy import float64
 from numpy.typing import NDArray
+from typing import Sequence
+from types import EllipsisType
 
 Shape = tuple[int, ...]
 NDFloat = float64
 FloatArray = NDArray[NDFloat]
+ArrayIndex = int | slice | None | EllipsisType | Sequence[int]
+ArrayIndices = ArrayIndex | tuple[ArrayIndex, ...]

tests/benchmark_cp.py

@@ -0,0 +1,42 @@
+import cantera as ct
+import numpy as np
+import time
+import gaspype as gp
+
+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
+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])

tests/benchmark_cp_comp.py

@@ -0,0 +1,55 @@
+import cantera as ct
+import numpy as np
+import time
+import gaspype as gp
+
+gas = ct.Solution("gri30.yaml")
+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)
+
+# Generate random compositions for H2, H2O, N2
+rng = np.random.default_rng(seed=42)
+fractions = rng.random((n_states, 3))
+fractions /= fractions.sum(axis=1)[:, None]  # normalize
+
+# Convert to full 53-species mole fraction array
+X = np.zeros((n_states, n_species))
+X[:, gas.species_index('H2')] = fractions[:, 0]
+X[:, gas.species_index('H2O')] = fractions[:, 1]
+X[:, gas.species_index('N2')] = fractions[:, 2]
+
+# Build SolutionArray
+states = ct.SolutionArray(gas, n_states)
+
+time.sleep(0.5)
+
+# Vectorized assignment
+t0 = time.perf_counter()
+states.TPX = temperatures, pressures, X
+cp_values = states.cp_mole
+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({'H2': 1}) * fractions[:, 0]
+    + gp.fluid({'H2O': 1}) * fractions[:, 1]
+    + gp.fluid({'N2': 1}) * fractions[:, 2]
+)
+
+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])

tests/benchmark_equalibrium.py

@@ -0,0 +1,54 @@
+import cantera as ct
+import gaspype as gp
+import numpy as np
+import time
+from gaspype import fluid_system
+
+# -----------------------
+# Settings
+# -----------------------
+n_temps = 1000
+temps_C = np.linspace(300, 1000, n_temps)     # °C
+temperatures = temps_C + 273.15               # K
+pressure = 101325                             # Pa (1 atm)
+
+composition = {"CH4": 0.8, "H2O": 0.2}
+species_to_track = ["CH4", "H2O", "CO", "CO2", "H2", "O2"]
+
+# -----------------------
+# Cantera benchmark
+# -----------------------
+gas = ct.Solution("gri30.yaml")
+gas.TPX = temperatures[0], pressure, composition
+
+eq_cantera = np.zeros((n_temps, len(species_to_track)))
+
+time.sleep(0.5)
+t0 = time.perf_counter()
+for i, T in enumerate(temperatures):
+    gas.TP = T, pressure
+    gas.equilibrate('TP')
+    eq_cantera[i, :] = [gas.X[gas.species_index(s)] for s in species_to_track]
+elapsed_cantera = time.perf_counter() - t0
+print(f"Cantera: {elapsed_cantera:.4f} s")
+
+# -----------------------
+# Gaspype benchmark
+# -----------------------
+# Construct the fluid with composition and tracked species
+fluid = gp.fluid(composition, fs=fluid_system(species_to_track))
+
+time.sleep(0.5)
+t0 = time.perf_counter()
+eq_gaspype = gp.equilibrium(fluid, t=temperatures, p=pressure)
+elapsed_gaspype = time.perf_counter() - t0
+print(f"Gaspype: {elapsed_gaspype:.4f} s")
+
+# -----------------------
+# Compare first 5 results
+# -----------------------
+print("First 5 equilibrium compositions (mole fractions):")
+for i in range(5):
+    print(f"T = {temperatures[i]:.1f} K")
+    print("  Cantera:", eq_cantera[i])
+    print("  Gaspype :", eq_gaspype.array_composition[i])

tests/md_to_code.py

@@ -96,6 +96,7 @@ def segments_to_test(segments: Iterable[markdown_segment], script_language: str
                 ret_block_flag = lines[-1] if (not re.match(r'^[^(]*=', lines[-1]) and
                                                not lines[-1].startswith('import ') and
                                                not lines[-1].startswith('from ') and
+                                               not lines[-1].startswith('print(') and
                                                not lines[-1].startswith(' ')) else None
                 # print('Last line: ', ret_block_flag, '-----------', lines[-1])
 

tests/test_edge_case_results.py

@@ -0,0 +1,20 @@
+import pytest
+import gaspype as gp
+
+
+def test_calculation_cold():
+    # Testing a non-equilibrium case. Result is only
+    # determined by stoichiometry.
+    fs = gp.fluid_system(['CH4', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+
+    el = gp.elements(gp.fluid({'H2': 1, 'CH4': 0.08, 'O2': 0.05}), fs=fs)
+
+    composition = gp.equilibrium(el, 300, 1e5)
+
+    print(el)
+    print(composition)
+
+    ref_result = [8.00000000e-02, 1.00000000e-01, 9.00000000e-01, 3.01246781e-23,
+                  2.90783583e-27, 3.60487456e-82]
+
+    assert composition.array_composition == pytest.approx(ref_result, abs=0.001)

tests/test_results_cantera.py

@@ -1,7 +1,8 @@
 import gaspype as gp
 import numpy as np
 # import pytest
-import cantera as ct  # type: ignore
+import cantera as ct
+from gaspype.typing import NDFloat
 
 
 def test_equilibrium_cantera():
@@ -16,7 +17,7 @@ def test_equilibrium_cantera():
 
     composition = gp.fluid({'H2': 1}, fs) +\
         gp.fluid({'CH4': 1}, fs) * np.linspace(0, 0.05, 30) +\
-        gp.fluid({'O2': 1}, fs) * np.expand_dims(np.linspace(0, 0.5, 30), axis=1)
+        gp.fluid({'O2': 1}, fs) * np.linspace(0, 0.5, 30)[:, None]
 
     t = 1495 + 273.15  # K
     p = 1e5  # Pa
@@ -40,10 +41,12 @@ def test_equilibrium_cantera():
         #if flow1.X[indeces][0] > 0.01:
         #    print(flow1.X[indeces])
 
-    ct_result_array = np.stack(ct_results)  # type: ignore
+    ct_result_array = np.stack(ct_results, dtype=NDFloat)  # type: ignore
 
-    max_deviation = np.max(np.abs((gp_result_array - ct_result_array)))
-    mean_deviation = np.mean(np.abs((gp_result_array - ct_result_array)))
+    deviations = np.abs(gp_result_array - ct_result_array)
 
-    assert max_deviation < 0.04
-    assert mean_deviation < 2e-4
+    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

tests/test_slicing.py

@@ -12,28 +12,28 @@ def test_str_index():
     assert el['C'].shape == (2, 3, 4)
 
 
-def test_str_list_index():
-    assert fl[['CO2', 'H2', 'CO']].shape == (2, 3, 4, 3)
-    assert el[['C', 'H', 'O']].shape == (2, 3, 4, 3)
+def test_single_axis_int_index():
+    assert fl[0].shape == (3, 4)
+    assert fl[1].shape == (3, 4)
+    assert el[1].shape == (3, 4)
+    assert el[0].shape == (3, 4)
 
 
-def test_int_list_index():
-    assert fl[[1, 2, 0, 5]].shape == (2, 3, 4, 4)
-    assert el[[1, 2, 0, 3]].shape == (2, 3, 4, 4)
+def test_single_axis_int_list():
+    assert fl[:, [0, 1]].shape == (2, 2, 4)
+    assert el[:, [0, 1]].shape == (2, 2, 4)
 
 
-def test_mixed_list_index():
-    assert el[[1, 'H', 0, 'O']].shape == (2, 3, 4, 4)
-
-
-def test_int_index():
-    assert fl[5].shape == (2, 3, 4)
-    assert el[-1].shape == (2, 3, 4)
-
-
-def test_slice_index():
-    assert fl[0:3].shape == (2, 3, 4, 3)
-    assert fl[:].shape == (2, 3, 4, 6)
-
-    assert el[0:3].shape == (2, 3, 4, 3)
-    assert el[:].shape == (2, 3, 4, 4)
+def test_multi_axis_int_index():
+    assert fl[0, 1].shape == (4,)
+    assert fl[0, 1, 2].shape == tuple()
+    assert fl[0, 2].shape == (4,)
+    assert fl[:, 2, :].shape == (2, 4)
+    assert fl[0, [1, 2]].shape == (2, 4)
+    assert fl[..., 0].shape == (2, 3)
+    assert el[0, 1].shape == (4,)
+    assert el[0, 1, 2].shape == tuple()
+    assert el[0, 2].shape == (4,)
+    assert el[:, 2, :].shape == (2, 4)
+    assert el[0, [1, 2]].shape == (2, 4)
+    assert el[..., 0].shape == (2, 3)