Compare commits

..

No commits in common. "main" and "v1.1.2" have entirely different histories.
main ... v1.1.2

28 changed files with 106 additions and 572 deletions

View File

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

View File

@ -27,13 +27,6 @@ jobs:
run: | run: |
python -m pip install --upgrade pip python -m pip install --upgrade pip
python -m pip install -e .[dev] 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 - name: Prepare data and compile thermo database
run: | run: |

View File

@ -9,7 +9,7 @@ permissions:
contents: write contents: write
jobs: jobs:
build: build-and-deploy:
runs-on: ubuntu-latest runs-on: ubuntu-latest
steps: steps:
- uses: actions/checkout@v4 - uses: actions/checkout@v4
@ -41,24 +41,8 @@ jobs:
rm ./source/*.rst rm ./source/*.rst
make html make html
touch ./build/html/.nojekyll touch ./build/html/.nojekyll
mkdir -p ./build/html/_autogenerated
cp ./build/html/api/* ./build/html/_autogenerated/
- name: Upload artifact
uses: actions/upload-pages-artifact@v3
with:
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 - name: Deploy to GitHub Pages
id: deployment uses: JamesIves/github-pages-deploy-action@v4
uses: actions/deploy-pages@v4 with:
branch: gh-pages
folder: docs/build/html

View File

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

2
.gitignore vendored
View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -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: [docs/source/render_examples.py](../docs/source/render_examples.py) script:
``` bash ``` bash
# Converting markdown with code sections to Jupyter Notebook and run it: # Converting markdown with code sections to Jupyter Notebook and run it:
notedown examples/soec_methane.md --to notebook --output docs/source/api/soec_methane.ipynb --run notedown examples/soec_methane.md --to notebook --output docs/source/_autogenerated/soec_methane.ipynb --run
# Converting the Jupyter Notebook to Markdown and a folder with image # Converting the Jupyter Notebook to Markdown and a folder with image
# files placed in docs/source/api/: # files placed in docs/source/_autogenerated/:
jupyter nbconvert --to markdown docs/source/api/soec_methane.ipynb --output soec_methane.md jupyter nbconvert --to markdown docs/source/_autogenerated/soec_methane.ipynb --output soec_methane.md
``` ```
A new example Markdown file can be created from a Jupyter Notebook running A new example Markdown file can be created from a Jupyter Notebook running

View File

@ -1,45 +0,0 @@
# 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.

View File

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

View File

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

View File

@ -1,141 +0,0 @@
# 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} %")
```

View File

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

View File

@ -1,7 +1,7 @@
# Carbon Activity # Carbon Activity
This example shows the calculation of the carbon activity for methane mixtures This example shows the equilibrium calculation for solid carbon.
in thermodynamic equilibrium.
```python ```python
import gaspype as gp import gaspype as gp
@ -66,25 +66,3 @@ ax.plot(ratio, carbon_activity.T)
ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed') ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed')
ax.legend([f'{tc} °C' for tc in t_range]) 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])
```

View File

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

View File

@ -7,7 +7,7 @@ from gaspype._phys_data import atomic_weights, db_reader
import re import re
import pkgutil import pkgutil
from .constants import R, epsy, p0 from .constants import R, epsy, p0
from .typing import FloatArray, NDFloat, Shape, ArrayIndices from .typing import FloatArray, NDFloat, Shape
T = TypeVar('T', 'fluid', 'elements') T = TypeVar('T', 'fluid', 'elements')
@ -483,28 +483,6 @@ 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' 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]] 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: def __add__(self, other: T) -> T:
return array_operation(self, other, np.add) return array_operation(self, other, np.add)
@ -532,21 +510,16 @@ class fluid:
# def __array__(self) -> FloatArray: # def __array__(self) -> FloatArray:
# return self.array_composition # return self.array_composition
@overload def __getitem__(self, key: str | int | list[str] | list[int] | slice) -> FloatArray:
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): if isinstance(key, str):
assert key in self.fs.species, f'Species {key} is not part of the fluid system' 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)] return self.array_composition[..., self.fs.species.index(key)]
elif isinstance(key, (slice, int)):
return self.array_composition[..., key]
else: else:
key_tuple = key if isinstance(key, tuple) else (key,) mset = set(self.fs.species) | set(range(len(self.fs.species)))
return fluid(self.array_composition[(*key_tuple, slice(None))], self.fs) 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]]
def __iter__(self) -> Iterator[dict[str, float]]: def __iter__(self) -> Iterator[dict[str, float]]:
assert len(self.shape) < 2, 'Cannot iterate over species with more than one dimension' assert len(self.shape) < 2, 'Cannot iterate over species with more than one dimension'
@ -641,28 +614,6 @@ class elements:
""" """
return np.sum(self.array_elemental_composition * self.fs.array_atomic_mass, axis=-1, dtype=NDFloat) 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': def __add__(self, other: 'fluid | elements') -> 'elements':
return array_operation(self, other, np.add) return array_operation(self, other, np.add)
@ -688,21 +639,16 @@ class elements:
def __array__(self) -> FloatArray: def __array__(self) -> FloatArray:
return self.array_elemental_composition return self.array_elemental_composition
@overload def __getitem__(self, key: str | int | list[str] | list[int] | slice) -> FloatArray:
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): if isinstance(key, str):
assert key in self.fs.elements, f'Element {key} is not part of the fluid system' 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)] return self.array_elemental_composition[..., self.fs.elements.index(key)]
elif isinstance(key, (slice, int)):
return self.array_elemental_composition[..., key]
else: else:
key_tuple = key if isinstance(key, tuple) else (key,) mset = set(self.fs.elements) | set(range(len(self.fs.elements)))
return elements(self.array_elemental_composition[(*key_tuple, slice(None))], self.fs) 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]]
def __iter__(self) -> Iterator[dict[str, float]]: def __iter__(self) -> Iterator[dict[str, float]]:
assert len(self.shape) < 2, 'Cannot iterate over elements with more than one dimension' assert len(self.shape) < 2, 'Cannot iterate over elements with more than one dimension'

View File

View File

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

View File

@ -1,42 +0,0 @@
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])

View File

@ -1,55 +0,0 @@
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])

View File

@ -1,54 +0,0 @@
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])

View File

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

View File

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