nicolas
/
gaspype
mirror of https://github.com/DLR-Institute-of-Future-Fuels/gaspype.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Compare commits

base: nicolas:main
Branches Tags
nicolas:main
nicolas:dev
nicolas:v1.1.3
nicolas:v1.1.2
nicolas:v1.1.1
nicolas:v1.1
nicolas:v1.0.1
nicolas:v1.0
..
compare: nicolas:v1.0.1
Branches Tags
nicolas:dev
nicolas:main
nicolas:v1.1.3
nicolas:v1.1.2
nicolas:v1.1.1
nicolas:v1.1
nicolas:v1.0.1
nicolas:v1.0

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

55 changed files with 6324 additions and 9699 deletions
Show Stats Expand all files Collapse all files

7
.flake8
View File

@ -13,10 +13,7 @@ exclude =
build, build,
dist, dist,
.conda, .conda,
tests/autogenerated_*, tests/autogenerated_*
docs/source/api
.venv,
venv
# Enable specific plugins or options # Enable specific plugins or options
# Example: Enabling flake8-docstrings # Example: Enabling flake8-docstrings
@ -24,4 +21,4 @@ select = C,E,F,W,D
# Specify custom error codes to ignore or enable # Specify custom error codes to ignore or enable
per-file-ignores = per-file-ignores =
tests/*: D tests/*: D

50
.github/workflows/ci.yml vendored
View File

@ -1,50 +0,0 @@
name: CI Pipeline
on:
push:
branches: [main]
pull_request:
branches: [main]
jobs:
build:
runs-on: ubuntu-latest
strategy:
matrix:
python-version: ["3.10", 3.11, 3.12, 3.13]
steps:
- name: Check out code
uses: actions/checkout@v4
- name: Set up Python ${{ matrix.python-version }}
uses: actions/setup-python@v5
with:
python-version: ${{ matrix.python-version }}
- name: Install dependencies
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: |
python thermo_data/combine_data.py thermo_data/combined_data.yaml thermo_data/nasa9*.yaml thermo_data/nasa9*.xml
python thermo_data/compile_to_bin.py thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin
- name: Lint code with flake8
run: flake8
- name: Type checking with mypy
run: mypy
- name: Run tests with pytest
run: pytest

64
.github/workflows/docs.yml vendored
View File

@ -1,64 +0,0 @@
name: Build and Deploy Docs
on:
push:
branches:
- main
permissions:
contents: write
jobs:
build:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- name: Set up Python
uses: actions/setup-python@v3
with:
python-version: "3.x"
- name: Build database
run: |
pip install pyyaml
python thermo_data/combine_data.py thermo_data/combined_data.yaml thermo_data/nasa9*.yaml
python thermo_data/compile_to_bin.py thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin
# echo "Create a dummy file to ensure gaspype does't crash"
# mkdir -p src/gaspype/data
# printf 'gapy\x00\x00\x00\x00' > src/gaspype/data/therm_data.bin
- name: Install gaspype and dependencies
run: |
pip install .[doc_build]
python -m ipykernel install --user --name temp_kernel --display-name "Python (temp_kernel)"
- name: Generate Docs
run: python ./docs/source/generate_class_list.py
- name: Generate Examples
run: python ./docs/source/render_examples.py
- name: Build Docs
run: |
cp LICENSE docs/source/LICENSE.md
cd docs
sphinx-apidoc -o ./source/ ../src/ -M --no-toc
rm ./source/*.rst
make html
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
id: deployment
uses: actions/deploy-pages@v4

44
.github/workflows/publish.yml vendored
View File

@ -1,44 +0,0 @@
name: Publish to PyPI
on:
push:
tags:
- "v*"
jobs:
publish:
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
- name: Ensure this is main branch
if: github.ref == 'refs/heads/main' || startsWith(github.ref, 'refs/tags/v')
run: echo "Proceeding with publish"
- uses: actions/setup-python@v5
with:
python-version: "3.11"
- name: Install build tools
run: python -m pip install --upgrade build twine pyyaml
- name: Prepare data and compile thermo database
run: |
python thermo_data/combine_data.py thermo_data/combined_data.yaml thermo_data/nasa9*.yaml thermo_data/nasa9*.xml
python thermo_data/compile_to_bin.py thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin
- name: Build package
run: python -m build
- name: Publish to PyPI
env:
TWINE_USERNAME: __token__
TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}
run: python -m twine upload dist/*

9
.gitignore vendored
View File

@ -7,11 +7,4 @@ __pycache__
.vscode .vscode
.mypy_cache .mypy_cache
.pytest_cache .pytest_cache
tests/autogenerated_*.py tests/autogenerated_*.py
docs/build/
docs/source/api/
venv/
.venv/
thermo_data/combined_data.yaml
src/gaspype/data/therm_data.bin
*/*/build/

11
CITATION.cff
View File

@ -1,5 +1,4 @@
cff-version: 1.1.0 cff-version: 1.2.0
message: "If you use this software, please cite it as below."
title: Gaspype 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,11 @@ 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: 1.0.1
date-released: "2025-06-24" date-released: "2025-05-09"
#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"

102
README.md
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:
@ -25,7 +24,13 @@ pip install gaspype
Installation with conda: Installation with conda:
``` bash ``` bash
conda install conda-forge::gaspype conda install gaspype
```
Installation for developers with pip:
``` bash
git clone https://github.com/DLR-Institute-of-Future-Fuels/gaspype
pip install -e .[dev]
``` ```
## Getting started ## Getting started
@ -53,7 +58,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
@ -63,8 +68,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}) + \
@ -125,10 +129,7 @@ array([[[0. , 0.5 , 0.5 ],
[0.625 , 0.09375 , 0.28125 ], [0.625 , 0.09375 , 0.28125 ],
[0.71428571, 0.07142857, 0.21428571]]]) [0.71428571, 0.07142857, 0.21428571]]])
``` ```
In some cases not the molecular but the atomic composition is of interest. The ```elements``` class can be used for atom based balances and works similar:
### Elements
In some cases not the molecular but the atomic composition is of interest.
The ```elements``` class can be used for atom based balances and works similar:
``` python ``` python
el = gp.elements({'N': 1, 'Cl': 2}) el = gp.elements({'N': 1, 'Cl': 2})
@ -137,9 +138,7 @@ el.get_mass()
``` ```
np.float64(0.08490700000000001) np.float64(0.08490700000000001)
``` ```
A ```elements``` object can be as well instantiated from a ```fluid``` object. A ```elements``` object can be as well instantiated from a ```fluid``` object. Arithmetic operations between ```elements``` and ```fluid``` result in an ```elements``` object:
Arithmetic operations between ```elements``` and ```fluid``` result in
an ```elements``` object:
``` python ``` python
el2 = gp.elements(fl) + el - 0.3 * fl el2 = gp.elements(fl) + el - 0.3 * fl
el2 el2
@ -151,8 +150,7 @@ N 1.000e+00 mol
O 7.000e-01 mol O 7.000e-01 mol
``` ```
Going from an atomic composition to an molecular composition is possible as well. Going from an atomic composition to an molecular composition is a little bit less straight forward, since there is no universal approach. One way is to calculate the thermodynamic equilibrium for a mixture:
One way is to calculate the thermodynamic equilibrium for a mixture:
``` python ``` python
fs = gp.fluid_system('CH4, H2, CO, CO2, O2') fs = gp.fluid_system('CH4, H2, CO, CO2, O2')
@ -169,13 +167,7 @@ CO2 33.07 %
O2 0.00 % O2 0.00 %
``` ```
The ```equilibrium``` function can be called with a ```fluid``` or ```elements``` object The ```equilibrium``` function can be called with a ```fluid``` or ```elements``` object as first argument. ```fluid``` and ```elements``` referencing a ```fluid_system``` object witch can be be set as shown above during the object instantiation. If not provided, a new one will be created automatically. Providing a ```fluid_system``` gives more control over which molecular species are included in derived ```fluid``` objects. Furthermore arithmetic operations between objects with the same ```fluid_system``` are potentially faster:
as first argument. ```fluid``` and ```elements``` referencing a ```fluid_system``` object
witch can be be set as shown above during the object instantiation. If not provided,
a new one will be created automatically. Providing a ```fluid_system``` gives more
control over which molecular species are included in derived ```fluid``` objects.
Furthermore arithmetic operations between objects with the same ```fluid_system```
are potentially faster:
``` python ``` python
fl3 + gp.fluid({'CH4': 1}, fs) fl3 + gp.fluid({'CH4': 1}, fs)
@ -189,9 +181,7 @@ CO2 18.07 %
O2 0.00 % O2 0.00 %
``` ```
Especially if the ```fluid_system``` of one of the operants has not a subset of Especially if the ```fluid_system``` of one of the operants has not a subset of molecular species of the other ```fluid_system``` a new ```fluid_system``` will be created for the operation which might degrade performance:
molecular species of the other ```fluid_system``` a new ```fluid_system``` will
be created for the operation which might degrade performance:
``` python ``` python
fl3 + gp.fluid({'NH3': 1}) fl3 + gp.fluid({'NH3': 1})
@ -204,46 +194,4 @@ CO2 18.07 %
H2 9.24 % H2 9.24 %
NH3 45.38 % NH3 45.38 %
O2 0.00 % O2 0.00 %
``` ```
## Developer Guide
Contributions are welcome, please open an issue or submit a pull request on GitHub.
To get started with developing the `gaspype` package, follow these steps.
First, clone the repository to your local machine using Git:
```bash
git clone https://github.com/DLR-Institute-of-Future-Fuels/gaspype.git
cd gaspype
```
It's recommended to setup an venv:
```bash
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
in the current directory:
```bash
pip install -e .[dev]
```
Compile binary property database from text based files:
```bash
python thermo_data/combine_data.py thermo_data/combined_data.yaml thermo_data/nasa9*.yaml thermo_data/nasa9*.xml
python thermo_data/compile_to_bin.py thermo_data/combined_data.yaml src/gaspype/data/therm_data.bin
```
Ensure that everything is set up correctly by running the tests:
```bash
pytest
```
## License
This project is licensed under the MIT License - see the [LICENSE](LICENSE) file for details.

20
docs/Makefile
View File

@ -1,20 +0,0 @@
# Minimal makefile for Sphinx documentation
#
# You can set these variables from the command line, and also
# from the environment for the first two.
SPHINXOPTS ?=
SPHINXBUILD ?= sphinx-build
SOURCEDIR = source
BUILDDIR = build
# Put it first so that "make" without argument is like "make help".
help:
@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
.PHONY: help Makefile
# Catch-all target: route all unknown targets to Sphinx using the new
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)

35
docs/make.bat
View File

@ -1,35 +0,0 @@
@ECHO OFF
pushd %~dp0
REM Command file for Sphinx documentation
if "%SPHINXBUILD%" == "" (
set SPHINXBUILD=sphinx-build
)
set SOURCEDIR=source
set BUILDDIR=build
%SPHINXBUILD% >NUL 2>NUL
if errorlevel 9009 (
echo.
echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
echo.installed, then set the SPHINXBUILD environment variable to point
echo.to the full path of the 'sphinx-build' executable. Alternatively you
echo.may add the Sphinx directory to PATH.
echo.
echo.If you don't have Sphinx installed, grab it from
echo.https://www.sphinx-doc.org/
exit /b 1
)
if "%1" == "" goto help
%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
goto end
:help
%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
:end
popd

193
docs/media/soc_inverted.svg
View File

@ -1,193 +0,0 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
width="164.70512mm"
height="86.731506mm"
viewBox="0 0 164.70512 86.731506"
version="1.1"
id="svg1"
inkscape:version="1.3.2 (091e20e, 2023-11-25, custom)"
sodipodi:docname="soc_inverted.svg"
inkscape:export-filename="soc_export.svg"
inkscape:export-xdpi="96"
inkscape:export-ydpi="96"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns="http://www.w3.org/2000/svg"
xmlns:svg="http://www.w3.org/2000/svg">
<sodipodi:namedview
id="namedview1"
pagecolor="#ffffff"
bordercolor="#000000"
borderopacity="0.25"
inkscape:showpageshadow="2"
inkscape:pageopacity="0.0"
inkscape:pagecheckerboard="0"
inkscape:deskcolor="#d1d1d1"
inkscape:document-units="mm"
inkscape:zoom="1.0041055"
inkscape:cx="425.25411"
inkscape:cy="332.13641"
inkscape:window-width="1920"
inkscape:window-height="1009"
inkscape:window-x="1912"
inkscape:window-y="-8"
inkscape:window-maximized="1"
inkscape:current-layer="layer1" />
<defs
id="defs1">
<marker
style="overflow:visible"
id="Triangle"
refX="0"
refY="0"
orient="auto-start-reverse"
inkscape:stockid="Triangle arrow"
markerWidth="2"
markerHeight="1"
viewBox="0 0 1 1"
inkscape:isstock="true"
inkscape:collect="always"
preserveAspectRatio="none">
<path
transform="scale(0.5)"
style="fill:context-stroke;fill-rule:evenodd;stroke:context-stroke;stroke-width:1pt"
d="M 5.77,0 -2.88,5 V -5 Z"
id="path135" />
</marker>
</defs>
<g
inkscape:label="Ebene 1"
inkscape:groupmode="layer"
id="layer1"
transform="translate(-18.027099,-19.242669)">
<path
style="fill:none;stroke:#dadada;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1;marker-start:url(#Triangle);marker-end:url(#Triangle)"
d="M 23.663111,22.545169 V 97.81995 H 179.42973"
id="path1" />
<path
style="fill:none;stroke:#dadada;stroke-width:0.372829;stroke-dasharray:none;stroke-opacity:1"
d="M 170.85884,97.81995 V 95.304579"
id="path2" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3"
cx="40.768127"
cy="46.63377"
r="1.2658331" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3-7"
cx="60.487274"
cy="54.953339"
r="1.2658331" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3-6"
cx="94.775032"
cy="67.060722"
r="1.2658331" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3-8"
cx="133.16524"
cy="76.384407"
r="1.2658331" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3-1"
cx="23.607611"
cy="32.419498"
r="1.2658331" />
<circle
style="fill:#dddddd;fill-opacity:1;stroke:#dadada;stroke-width:0.468334;stroke-dasharray:none;stroke-opacity:1"
id="path3-8-3"
cx="170.49799"
cy="83.681023"
r="1.2658331" />
<path
style="fill:none;fill-opacity:1;stroke:#dadada;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
d="m 23.61653,32.379465 17.141782,14.300346 19.796893,8.151661 34.283565,12.111043 38.38268,9.502507 37.26475,7.266625"
id="path4"
sodipodi:nodetypes="cccccc" />
<text
xml:space="preserve"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
x="86.083488"
y="105.09981"
id="text4"><tspan
sodipodi:role="line"
id="tspan4"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
x="86.083488"
y="105.09981">Relative cell lenghs</tspan></text>
<text
xml:space="preserve"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
x="169.54897"
y="103.00477"
id="text4-1"><tspan
sodipodi:role="line"
id="tspan4-1"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
x="169.54897"
y="103.00477">1</tspan></text>
<text
xml:space="preserve"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
x="22.457224"
y="103.68204"
id="text4-1-1"><tspan
sodipodi:role="line"
id="tspan4-1-2"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
x="22.457224"
y="103.68204">0</tspan></text>
<path
style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1"
d="M 40.710985,49.035087 V 65.438055"
id="path5" />
<path
style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1"
d="m 60.546682,57.469645 v 7.773293"
id="path5-9" />
<path
style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:none;stroke-opacity:1;marker-start:url(#Triangle);marker-end:url(#Triangle)"
d="M 58.396033,62.880678 H 42.878103"
id="path5-0"
sodipodi:nodetypes="cc" />
<text
xml:space="preserve"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
x="48.511959"
y="61.409679"
id="text4-1-1-5"><tspan
sodipodi:role="line"
id="tspan4-1-2-2"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
x="48.511959"
y="61.409679">Δz</tspan></text>
<text
xml:space="preserve"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465;stroke-dasharray:none;stroke-opacity:1"
x="-72.630066"
y="21.160755"
id="text4-9"
transform="rotate(-90)"><tspan
sodipodi:role="line"
id="tspan4-5"
style="font-size:4.23333px;fill:#dedede;fill-opacity:1;stroke:none;stroke-width:0.465"
x="-72.630066"
y="21.160755">Current</tspan></text>
<path
style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:1.2, 0.3;stroke-dashoffset:0;stroke-opacity:1"
d="M 94.728796,67.061137 V 97.890816"
id="path6" />
<path
style="fill:#000000;fill-opacity:1;stroke:#dadada;stroke-width:0.3;stroke-dasharray:1.2, 0.3;stroke-dashoffset:0;stroke-opacity:1"
d="M 133.2181,76.566042 V 97.753468"
id="path6-1"
sodipodi:nodetypes="cc" />
</g>
</svg>
Side by Side

Before

Width:  |  Height:  |  Size: 7.3 KiB

33
docs/source/conf.py
View File

@ -1,33 +0,0 @@
# Configuration file for the Sphinx documentation builder.
#
# For the full list of built-in configuration values, see the documentation:
# https://www.sphinx-doc.org/en/master/usage/configuration.html
# -- Project information -----------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
import os
import sys
sys.path.insert(0, os.path.abspath("../src/"))
project = 'Gaspype'
copyright = '2025, Nicolas Kruse'
author = 'Nicolas Kruse'
# -- General configuration ---------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon", "sphinx_autodoc_typehints", "myst_parser", "sphinx.ext.autosummary"]
templates_path = ['_templates']
exclude_patterns = []
# -- Options for HTML output -------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
# html_theme = 'alabaster'
html_theme = 'pydata_sphinx_theme'
html_static_path = ['_static']
autodoc_inherit_docstrings = True
autoclass_content = 'both'

5
docs/source/floatarray.md
View File

@ -1,5 +0,0 @@
# gaspype.typing.FloatArray
```{eval-rst}
.. autoclass:: gaspype.typing.FloatArray
```

86
docs/source/generate_class_list.py
View File

@ -1,86 +0,0 @@
# This script generates the source md-files for all classes and functions for the docs
import importlib
import inspect
import fnmatch
from io import TextIOWrapper
import os
def write_manual(f: TextIOWrapper, doc_files: list[str], title: str) -> None:
write_dochtree(f, title, doc_files)
def write_classes(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
"""Write the classes to the file."""
module = importlib.import_module(module_name)
classes = [
name for name, obj in inspect.getmembers(module, inspect.isclass)
if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude) and
obj.__doc__ and '(Automatic generated stub)' not in obj.__doc__)
]
if description:
f.write(f'{description}\n\n')
write_dochtree(f, title, classes)
for cls in classes:
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')
f2.write(' :members:\n')
f2.write(' :undoc-members:\n')
f2.write(' :show-inheritance:\n')
f2.write(' :inherited-members:\n')
f2.write('```\n\n')
def write_functions(f: TextIOWrapper, patterns: list[str], module_name: str, title: str, description: str = '', exclude: list[str] = []) -> None:
"""Write the classes to the file."""
module = importlib.import_module(module_name)
functions = [
name for name, _ in inspect.getmembers(module, inspect.isfunction)
if (any(fnmatch.fnmatch(name, pat) for pat in patterns if pat not in exclude))
]
if description:
f.write(f'{description}\n\n')
write_dochtree(f, title, functions)
for func in functions:
if not func.startswith('_'):
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')
f2.write('```\n\n')
def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
f.write('```{toctree}\n')
f.write(':maxdepth: 1\n')
f.write(f':caption: {title}:\n')
#f.write(':hidden:\n')
for text in items:
if not text.startswith('_'):
f.write(f"{text}\n")
f.write('```\n\n')
if __name__ == "__main__":
# Ensure the output directory exists
os.makedirs('docs/source/api', exist_ok=True)
with open('docs/source/api/index.md', 'w') as f:
f.write('# Classes and functions\n\n')
write_classes(f, ['*'], 'gaspype', title='Classes')
write_functions(f, ['*'], 'gaspype', title='Functions')
write_manual(f, ['../ndfloat', '../floatarray'], title='Types')

10
docs/source/index.md
View File

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

5
docs/source/ndfloat.md
View File

@ -1,5 +0,0 @@
# gaspype.typing.NDFloat
```{eval-rst}
.. autoclass:: gaspype.typing.NDFloat
```

66
docs/source/render_examples.py
View File

@ -1,66 +0,0 @@
# This script converts the example md-files as jupyter notebook,
# execute the notebook and convert the notebook back to a md-file
# with outputs included.
import subprocess
from glob import glob
import os
from io import TextIOWrapper
def run_cmd(command: list[str]):
result = subprocess.run(command, capture_output=True, text=True)
print('> ' + ' '.join(command))
print(result.stdout)
assert (not result.stderr or
any('RuntimeWarning: ' in line for line in result.stderr.splitlines()) or
any('[NbConvertApp]' in line or 'Warning' in line for line in result.stderr.splitlines())), 'ERROR: ' + result.stderr
def run_rendering(input_path: str, output_directory: str):
file_name = '.'.join(os.path.basename(input_path).split('.')[:-1])
assert file_name
print(f'- Convert {input_path} ...')
run_cmd(['notedown', input_path, '--to', 'notebook', '--output', f'{output_directory}/{file_name}.ipynb', '--run'])
run_cmd(['jupyter', 'nbconvert', '--to', 'markdown', f'{output_directory}/{file_name}.ipynb', '--output', f'{file_name}.md'])
run_cmd(['python', 'tests/md_to_code.py', 'script', f'{input_path}', f'{output_directory}/{file_name}.py'])
def write_dochtree(f: TextIOWrapper, title: str, items: list[str]):
f.write('```{toctree}\n')
f.write(':maxdepth: 1\n')
#f.write(':hidden:\n')
#f.write(f':caption: {title}:\n')
for text in items:
if not text.startswith('_'):
f.write(f"{text}\n")
f.write('```\n\n')
def render_examples(filter: str, example_file: str):
files = glob(filter)
names = ['.'.join(os.path.basename(path).split('.')[:-1]) for path in files]
with open(example_file, 'w') as f:
f.write('# Gaspype examples\n\n')
write_dochtree(f, '', [n for n in names if n.lower() != 'readme'])
f.write('## Download Jupyter Notebooks\n\n')
for path, name in zip(files, names):
if name.lower() != 'readme':
run_rendering(path, 'docs/source/api')
notebook = name + '.ipynb'
f.write(f'- [{notebook}]({notebook})\n\n')
f.write('## Download plain python files\n\n')
for path, name in zip(files, names):
if name.lower() != 'readme':
script_name = name + '.py'
f.write(f'- [{script_name}]({script_name})\n\n')
if __name__ == "__main__":
render_examples('examples/*.md', 'docs/source/api/examples.md')

3
docs/source/repo.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).

25
examples/README.md
View File

@ -1,25 +0,0 @@
# Example scripts
Examples can be looked-up in the
[documentation](https://dlr-institute-of-future-fuels.github.io/gaspype/)
rendered with results.
The gaspype examples from this directory are available in the documentation as
downloadable Jupyter Notebooks or plain python scripts with comments.
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/api/soec_methane.ipynb --run
# Converting the Jupyter Notebook to Markdown and a folder with image
# 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
the following command:
``` bash
jupyter nbconvert --to markdown new_example.ipynb --NbConvertApp.use_output_suffix=False --ClearOutputPreprocessor.enabled=True --output-dir examples/ --output new_example.md
```

45
examples/hydrogen_production_thermodynamics.md
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.

60
examples/methane_mixtures.md
View File

@ -1,60 +0,0 @@
# Methane Mixtures
This example shows equilibria of methane mixed with steam and CO2
```python
import gaspype as gp
import numpy as np
import matplotlib.pyplot as plt
```
Setting temperature and pressure:
```python
t = 900 + 273.15
p = 1e5
fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4', 'O2'])
```
Equilibrium calculation for methane steam mixtures:
```python
ratio = np.linspace(0.01, 1.5, num=64)
fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
equilibrium_h2o = gp.equilibrium(fl, t, p)
```
```python
fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
ax.set_xlabel("H2O/CH4")
ax.set_ylabel("molar fraction")
ax.set_ylim(0, 1.1)
#ax.set_xlim(0, 100)
ax.plot(ratio, equilibrium_h2o.get_x())
ax.legend(fs.active_species)
```
Equilibrium calculation for methane CO2 mixtures:
```python
fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'CO2': 1}, fs)
equilibrium_co2 = gp.equilibrium(fl, t, p)
```
```python
fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
ax.set_xlabel("CO2/CH4")
ax.set_ylabel("molar fraction")
ax.set_ylim(0, 1.1)
#ax.set_xlim(0, 100)
ax.plot(ratio, equilibrium_co2.get_x())
ax.legend(fs.active_species)
```

124
examples/soec_syngas.md
View File

@ -1,124 +0,0 @@
# SOEC Co-Electrolysis
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 necessarily realistic. For example,
a utilization of 0.95 causes issues with the formation of solid carbon.
```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
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({'H2O': 2, 'CO2': 1}, fs)
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 * utilization * air_dilution
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)
```
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'])
```
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:
```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 = 1.3 # 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, '-')
```

141
examples/sofc_methane.md
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} %")
```

58
examples/sulfur_oxygen_equalibrium.md
View File

@ -1,58 +0,0 @@
# Sulfur Oxygen Equilibrium
This example shows equilibrium calculations for sulfur/oxygen mixtures.
```python
import gaspype as gp
import numpy as np
import matplotlib.pyplot as plt
```
List possible sulfur/oxygen species:
```python
gp.species(element_names = 'S, O')
```
Or more specific by using regular expressions:
```python
gp.species('S?[2-3]?O?[2-5]?', use_regex=True)
```
Calculation of the molar equilibrium fractions for sulfur and oxygen depending on the oxygen to sulfur ratio:
```python
fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2'])
oxygen_ratio = np.linspace(0.5, 3, num=128)
el = gp.elements({'S': 1}, fs) + oxygen_ratio * gp.elements({'O': 1}, fs)
composition = gp.equilibrium(el, 800+273.15, 1e4)
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:
```python
fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2'])
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)
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)
```

90
examples/thermodynamic_carbon_activity.md
View File

@ -1,90 +0,0 @@
# Carbon Activity
This example shows the calculation of the carbon activity for methane mixtures
in thermodynamic equilibrium.
```python
import gaspype as gp
import numpy as np
import matplotlib.pyplot as plt
```
Setting temperatures and pressure:
```python
t_range = np.array([600, 700, 800, 900, 1100, 1500]) # °C
p = 1e5 # Pa
fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4'])
```
Equilibrium calculation for methane steam mixtures:
```python
ratio = np.linspace(0.01, 1.5, num=128)
fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs)
```
gaspype.carbon_activity supports currently only 0D fluids therefore we build this helper function:
```python
def partial_c_activity(fl: gp.fluid, t: float, p: float):
fls = fl.array_composition.shape
eq_fl = gp.equilibrium(fl, t, p)
ret = np.zeros(fls[0])
for i in range(fls[0]):
ret[i] = gp.carbon_activity(gp.fluid(eq_fl.array_composition[i,:], fs), t, p)
return ret
```
Now we use the helper function to calculate the carbon activitx for all
compositions in equilibrium_h2o times all temperatures in t_range:
```python
carbon_activity = np.vstack([partial_c_activity(fl, tc + 273.15, p) for tc in t_range])
```
Plot carbon activities, a activity of > 1 means there is thermodynamically the formation of sold carbon favored.
```python
fig, ax = plt.subplots(figsize=(6, 4), dpi=120)
ax.set_xlabel("H2O/CH4")
ax.set_ylabel("carbon activity")
ax.set_ylim(1e-1, 1e3)
ax.set_yscale('log')
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])
```

25
pyproject.toml
View File

@ -1,6 +1,6 @@
[project] [project]
name = "gaspype" name = "gaspype"
version = "1.1.3" version = "1.0.1"
authors = [ authors = [
{ name="Nicolas Kruse", email="nicolas.kruse@dlr.de" }, { name="Nicolas Kruse", email="nicolas.kruse@dlr.de" },
] ]
@ -13,14 +13,13 @@ classifiers = [
"Operating System :: OS Independent", "Operating System :: OS Independent",
] ]
dependencies = [ dependencies = [
"pyyaml>=6.0.1",
"numpy>2.0.0", "numpy>2.0.0",
"scipy>1.12.0", "scipy>1.12.0",
] ]
[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"
[build-system] [build-system]
@ -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/*.yaml"]
[project.optional-dependencies] [project.optional-dependencies]
dev = [ dev = [
@ -40,22 +39,8 @@ dev = [
"pytest", "pytest",
"pandas", "pandas",
"cantera", "cantera",
"pyyaml>=6.0.1",
"types-PyYAML", "types-PyYAML",
"scipy-stubs", "scipy-stubs"
"matplotlib"
]
doc_build = [
"sphinx",
"pydata_sphinx_theme",
"sphinx-autodoc-typehints",
"myst-parser",
"pandas",
"matplotlib",
"ipykernel",
"jupyter",
"nbconvert",
"notedown"
] ]
[tool.mypy] [tool.mypy]

1073
src/gaspype/__init__.py
View File

File diff suppressed because it is too large Load Diff

820
src/gaspype/_main.py
View File

@ -1,820 +0,0 @@
import numpy as np
from numpy.typing import NDArray
from typing import Sequence, Any, TypeVar, Iterator, overload, Callable
from math import log as ln, ceil
from scipy.linalg import null_space
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, ArrayIndices
T = TypeVar('T', 'fluid', 'elements')
_data = pkgutil.get_data(__name__, 'data/therm_data.bin')
assert _data is not None, 'Could not load thermodynamic data'
_species_db = db_reader(_data)
def species(pattern: str = '*', element_names: str | list[str] = [], use_regex: bool = False) -> list[str]:
"""Returns a alphabetically sorted list of all available species
filtered by a pattern if supplied
Args:
pattern: Optional filter for specific molecules. Placeholder
characters: # A number including non written ones: 'C#H#' matches 'CH4';
$ Arbitrary element name; * Any sequence of characters
element_names: restrict results to species that contain only the specified elements.
The elements can be supplied as list of strings or as comma separated string.
use_regex: using regular expression for the pattern
Returns:
List of species
"""
if isinstance(element_names, str):
elements = {s.strip() for s in element_names.split(',')}
else:
assert isinstance(element_names, list), 'type of element_names must be list or str'
elements = set(element_names)
for el in elements:
assert el in atomic_weights, f'element {el} unknown'
if not use_regex:
el_pattern = '|'.join([el for el in atomic_weights.keys()])
pattern = pattern.replace('*', '.*')
pattern = pattern.replace('#', '\\d*')
pattern = pattern.replace('$', '(' + el_pattern + ')')
pattern = '^' + pattern + '(,.*)?$'
if element_names == []:
return [sn for sn in _species_db.names if re.fullmatch(pattern, sn)]
else:
return [
s.name for s in _species_db
if re.fullmatch(pattern, s.name) and
(len(elements) == 0 or set(s.composition.keys()).issubset(elements))]
class fluid_system:
"""A class to represent a fluid_system defined by a set of selected species.
Attributes:
species_names (list[str]): List of selected species in the fluid_system
array_molar_mass (FloatArray): Array of the molar masses of the species in the fluid_system
array_element_composition (FloatArray): Array of the element composition of the species in the fluid_system.
Dimension is: (number of species, number of elements)
array_atomic_mass (FloatArray): Array of the atomic masses of the elements in the fluid_system
"""
def __init__(self, species: list[str] | str, t_min: int = 250, t_max: int = 2000):
"""Instantiates a fluid_system.
Args:
species: List of species names to be available in the constructed
fluid_system (as list of strings or a comma separated string)
t_min: Lower bound of the required temperature range in Kelvin
t_max: Upper bound of the required temperature range in Kelvin
"""
if isinstance(species, str):
species = [s.strip() for s in species.split(',')]
temperature_base_points = range(int(t_min), ceil(t_max))
data_shape = (len(temperature_base_points), len(species))
self._cp_array = np.zeros(data_shape)
self._h_array = np.zeros(data_shape)
self._s_array = np.zeros(data_shape)
# self._g_array = np.zeros(data_shape)
self._g_rt_array = np.zeros(data_shape)
self._t_offset = int(t_min)
self.species = species
self.active_species = species
element_compositions: list[dict[str, int]] = list()
for i, s in enumerate(species):
species_data = _species_db.read(s)
if not species_data:
raise Exception(f'Species {s} not found')
element_compositions.append(species_data.composition)
assert species_data.model == 9, 'Only NASA9 polynomials are supported'
for t1, t2, a in zip(species_data.t_range[:-1], species_data.t_range[1:], species_data.data):
for j, T in enumerate(temperature_base_points):
if t2 >= T >= t1:
self._cp_array[j, i] = R * (a[0]*T**-2 + a[1]*T**-1 + a[2] + a[3]*T
+ a[4]*T**2 + a[5]*T**3 + a[6]*T**4)
self._h_array[j, i] = R*T * (-a[0]*T**-2 + a[1]*ln(T)/T + a[2]
+ a[3]/2*T + a[4]/3*T**2 + a[5]/4*T**3
+ a[6]/5*T**4 + a[7]/T)
self._s_array[j, i] = R * (-a[0]/2*T**-2 - a[1]*T**-1 + a[2]*ln(T)
+ a[3]*T + a[4]/2*T**2 + a[5]/3*T**3
+ a[6]/4*T**4 + a[8])
#self._g_array[j, i] = self._h_array[j, i] - self._s_array[j, i] * T
self._g_rt_array[j, i] = (self._h_array[j, i] / T - self._s_array[j, i]) / R
# TODO: Check if temperature range is not available
# 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: 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: 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
def get_species_h(self, t: float | FloatArray) -> FloatArray:
"""Get the molar enthalpies for all species in the fluid system
Args:
t: Temperature in Kelvin (can be an array)
Returns:
Array with the enthalpies of each specie in J/mol
"""
return lookup(self._h_array, t, self._t_offset)
def get_species_s(self, t: float | FloatArray) -> FloatArray:
"""Get the molar entropies for all species in the fluid system
Args:
t: Temperature in Kelvin (can be an array)
Returns:
Array with the entropies of each specie in J/mol/K
"""
return lookup(self._s_array, t, self._t_offset)
def get_species_cp(self, t: float | FloatArray) -> FloatArray:
"""Get the isobaric molar heat capacity for all species in the fluid system
Args:
t: Temperature in Kelvin (can be an array)
Returns:
Array with the heat capacities of each specie in J/mol/K
"""
return lookup(self._cp_array, t, self._t_offset)
# def get_species_g(self, t: float | NDArray[_Float]) -> NDArray[_Float]:
# return lookup(self._g_array, t, self._t_offset)
def get_species_g_rt(self, t: float | FloatArray) -> FloatArray:
"""Get specific gibbs free energy divided by RT for all species in the
fluid system (g/R/T == (h/T-s)/R )
Args:
t: Temperature in Kelvin (can be an array)
Returns:
Array of gibbs free energy divided by RT (dimensionless)
"""
return lookup(self._g_rt_array, t, self._t_offset)
def get_species_references(self) -> str:
"""Get a string with the references for all fluids of the fluid system
Returns:
String with the references
"""
return '\n'.join([f'{s:<12}: {_species_db[s].ref_string}' for s in self.species])
def __add__(self, other: 'fluid_system') -> 'fluid_system':
assert isinstance(other, self.__class__)
return self.__class__(self.species + other.species)
def __repr__(self) -> str:
return ('Fluid system\n Species: ' + ', '.join(self.species) +
'\n Elements: ' + ', '.join(self.elements))
class fluid:
"""A class to represent a fluid defined by a composition of
one or more species.
Attributes:
species (list[str]): List of species names in the associated fluid_system
array_composition (FloatArray): Array of the molar amounts of the species in the fluid
array_element_composition (FloatArray): Array of the element composition in the fluid
array_fractions (FloatArray): Array of the molar fractions of the species in the fluid
total (FloatArray | float): Array of the sums of the molar amount of all species
fs (fluid_system): Reference to the fluid_system used for this fluid
shape (_Shape): Shape of the fluid array
elements (list[str]): List of elements in the fluid_system
"""
__array_priority__ = 100
def __init__(self, composition: dict[str, float] | list[float] | FloatArray,
fs: fluid_system | None = None,
shape: Sequence[int] | None = None):
"""Instantiates a fluid.
Args:
composition: A dict of species names with their composition, e.g.
{'O2':0.5,'H2O':0.5} or a list/numpy-array of compositions.
The array can be multidimensional, the size of the last dimension
must match the number of species defined for the fluid_system.
The indices of the last dimension correspond to the indices in
the active_species list of the fluid_system.
fs: Reference to a fluid_system. Is optional if composition is
defined by a dict. If not specified a new fluid_system with
the components from the dict is created.
shape: Tuple or list for the dimensions the fluid array. Can
only be used if composition argument is a dict. Otherwise
the dimensions are specified by the composition argument.
"""
if fs is None:
assert isinstance(composition, dict), 'fluid system must be specified if composition is not a dict'
fs = fluid_system(list(composition.keys()))
if isinstance(composition, list):
composition = np.array(composition)
if isinstance(composition, dict):
missing_species = [s for s in composition if s not in fs.species]
if len(missing_species):
raise Exception(f'Species {missing_species[0]} is not part of the fluid system')
species_composition = [composition[s] if s in composition.keys() else 0 for s in fs.species]
comp_array = np.array(species_composition, dtype=NDFloat)
if shape is not None:
comp_array = comp_array * np.ones(list(shape) + [len(fs.species)], dtype=NDFloat)
else:
assert shape is None, 'specify shape by the shape of the composition array.'
assert composition.shape[-1] == len(fs.species), f'composition.shape[-1] ({composition.shape[-1]}) must be {len(fs.species)}'
comp_array = composition
self.array_composition: FloatArray = comp_array
self.total: FloatArray | float = np.sum(self.array_composition, axis=-1, dtype=NDFloat)
self.array_fractions: FloatArray = self.array_composition / (np.expand_dims(self.total, -1) + epsy)
self.shape: Shape = self.array_composition.shape[:-1]
self.fs = fs
self.array_elemental_composition: FloatArray = np.dot(self.array_composition, fs.array_species_elements)
self.species = fs.species
self.elements = fs.elements
def get_composition_dict(self) -> dict[str, float]:
"""Get a dict of the molar amount of each fluid species
Returns:
Returns a dict of floats with the molar amount of each fluid species in mol
"""
return {s: c for s, c in zip(self.fs.species, self.array_composition)}
def get_fractions_dict(self) -> dict[str, float]:
"""Get a dict of the molar fractions of each fluid species
Returns:
Returns a dict of floats with the molar fractions of each fluid species
"""
return {s: c for s, c in zip(self.fs.species, self.array_fractions)}
def get_h(self, t: float | FloatArray) -> FloatArray | float:
"""Get specific enthalpy of the fluid at the given temperature
Enthalpy is referenced to 25 °C and includes enthalpy of formation.
Therefore the enthalpy of H2 and O2 is 0 at 25 °C, but the enthalpy
of water vapor at 25 °C is −241 kJ/mol (enthalpy of formation).
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Enthalpies in J/mol
"""
return np.sum(self.fs.get_species_h(t) * self.array_fractions, axis=-1, dtype=NDFloat)
def get_H(self, t: float | FloatArray) -> FloatArray | float:
"""Get absolute enthalpy of the fluid at the given temperature
Enthalpy is referenced to 25 °C and includes enthalpy of formation.
Therefore the enthalpy of H2 and O2 is 0 at 25 °C, but the enthalpy
of water vapor at 25 °C is −241 kJ/mol (enthalpy of formation).
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Enthalpies in J
"""
return np.sum(self.fs.get_species_h(t) * self.array_composition, axis=-1, dtype=NDFloat)
def get_s(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get molar entropy of the fluid at the given temperature and pressure
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Entropy in J/mol/K
"""
x = self.array_fractions
s = self.fs.get_species_s(t)
return np.sum(x * (s - R * np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
def get_S(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get absolute entropy of the fluid at the given temperature and pressure
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Entropy in J/K
"""
x = self.array_fractions
n = self.array_composition
s = self.fs.get_species_s(t)
return np.sum(n * (s - R * np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
def get_cp(self, t: float | FloatArray) -> FloatArray | float:
"""Get molar heat capacity at constant pressure
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Heat capacity in J/mol/K
"""
return np.sum(self.fs.get_species_cp(t) * self.array_fractions, axis=-1, dtype=NDFloat)
def get_g(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get molar gibbs free energy (h - Ts)
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Gibbs free energy in J/mol
"""
x = self.array_fractions
grt = self.fs.get_species_g_rt(t)
return R * t * np.sum(x * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
def get_G(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get absolute gibbs free energy (H - TS)
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Gibbs free energy in J
"""
x = self.array_fractions
n = self.array_composition
grt = self.fs.get_species_g_rt(t)
return R * t * np.sum(n * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
def get_g_rt(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get specific gibbs free energy divided by RT: g/R/T == (h/T-s)/R
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressures(s) in Pascal. Fluid shape and shape of the temperature
must be broadcastable
Returns:
Gibbs free energy divided by RT (dimensionless)
"""
x = self.array_fractions
grt = self.fs.get_species_g_rt(t)
return np.sum(x * (grt + np.log(np.expand_dims(p / p0, -1) * x + epsy)), axis=-1, dtype=NDFloat)
def get_v(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get Absolute fluid volume
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressure in Pa. Fluid shape and shape of the pressure
must be broadcastable
Returns:
Volume of the fluid in m³
"""
return R / p * t * self.total
def get_vm(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get molar fluid volume
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressure in Pa. Fluid shape and shape of the pressure
must be broadcastable
Returns:
Molar volume of the fluid in m³/mol
"""
return R / p * t
def get_mass(self) -> FloatArray | float:
"""Get Absolute fluid mass
Returns:
Mass of the fluid in kg
"""
return np.sum(self.array_composition * self.fs.array_molar_mass, axis=-1, dtype=NDFloat)
def get_molar_mass(self) -> FloatArray | float:
"""Get molar fluid mass
Returns:
Mass of the fluid in kg/mol
"""
return np.sum(self.array_fractions * self.fs.array_molar_mass, axis=-1, dtype=NDFloat)
def get_density(self, t: float | FloatArray, p: float | FloatArray) -> FloatArray | float:
"""Get mass based fluid density
Args:
t: Absolute temperature(s) in Kelvin. Fluid shape and shape of the temperature
must be broadcastable
p: Pressure in Pa. Fluid shape and shape of the pressure
must be broadcastable
Returns:
Density of the fluid in kg/m³
"""
return np.sum(self.array_fractions * self.fs.array_molar_mass, axis=-1, dtype=NDFloat) / (R * t) * p
def get_x(self, species: str | list[str] | None = None) -> FloatArray:
"""Get molar fractions of fluid species
Args:
species: A single species name, a list of species names or None for
returning the molar fractions of all species
Returns:
Returns an array of floats with the molar fractions 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_fractions
elif isinstance(species, str):
assert species in self.fs.species, f'Species {species} is not part of the fluid system'
return self.array_fractions[..., 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_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)
def __sub__(self, other: T) -> T:
return array_operation(self, other, np.subtract)
def __truediv__(self, other: int | float | NDArray[Any]) -> 'fluid':
if isinstance(other, np.ndarray):
k = np.expand_dims(other, -1)
else:
k = np.array(other, dtype=NDFloat)
return self.__class__(self.array_composition / k, self.fs)
def __mul__(self, other: int | float | NDArray[Any]) -> 'fluid':
k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
return self.__class__(self.array_composition * k, self.fs)
def __rmul__(self, other: int | float | NDArray[Any]) -> 'fluid':
k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
return self.__class__(self.array_composition * k, self.fs)
def __neg__(self) -> 'fluid':
return self.__class__(-self.array_composition, self.fs)
# def __array__(self) -> FloatArray:
# return self.array_composition
@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)]
else:
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'
aec = self.array_composition.reshape(-1, len(self.fs.species))
return iter({s: c for s, c in zip(self.fs.species, aec[i, :])} for i in range(aec.shape[0]))
def __repr__(self) -> str:
if len(self.array_fractions.shape) == 1:
lines = [f'{s:16} {c * 100:5.2f} %' for s, c in zip(self.fs.species, self.array_fractions)]
return f'{"Total":16} {self.total:8.3e} mol\n' + '\n'.join(lines)
else:
array_disp = self.array_fractions.__repr__()
padding = int(array_disp.find('\n') / (len(self.fs.species) + 1))
return ('Total mol:\n' + self.total.__repr__() +
'\nSpecies:\n' + ' ' * int(padding / 2) +
''.join([(' ' * (padding - len(s))) + s for s in self.fs.species]) +
'\nMolar fractions:\n' + self.array_fractions.__repr__())
class elements:
"""Represent a fluid by composition of elements.
Attributes:
array_element_composition (FloatArray): Array of the element composition
"""
__array_priority__ = 100
def __init__(self, composition: fluid | dict[str, float] | list[str] | list[float] | FloatArray,
fs: fluid_system | None = None, shape: Sequence[int] | None = None):
"""Instantiates an elements object.
Args:
composition: A fluid object, a dict of element names with their
composition, e.g.: {'O':1,'H':2} or a list/numpy-array of compositions.
The array can be multidimensional, the size of the last dimension
must match the number of elements used in the fluid_system.
The indices of the last dimension correspond to the indices in
the active_species list of the fluid_system.
fs: Reference to a fluid_system.
shape: Tuple or list for the dimensions the fluid array. Can
only be used if composition argument is a dict. Otherwise
the dimensions are specified by the composition argument.
"""
if isinstance(composition, list):
composition = np.array(composition)
if isinstance(composition, fluid):
new_composition: FloatArray = np.dot(composition.array_composition, composition.fs.array_species_elements)
if fs:
self.array_elemental_composition = reorder_array(new_composition, composition.fs.elements, fs.elements)
else:
self.array_elemental_composition = new_composition
fs = composition.fs
elif isinstance(composition, dict) and fs is None:
fs = fluid_system(species(element_names=list(composition.keys())))
else:
assert fs, 'fluid system must be specified if composition is not specified by a fluid'
if isinstance(composition, dict):
missing_elements = [s for s in composition if s not in fs.elements]
if len(missing_elements):
raise Exception(f'Element {missing_elements[0]} is not part of the fluid system')
self.array_elemental_composition = np.array([composition[s] if s in composition.keys() else 0 for s in fs.elements])
if shape is not None:
self.array_elemental_composition = self.array_elemental_composition * np.ones(list(shape) + [len(fs.species)])
elif isinstance(composition, np.ndarray):
assert shape is None, 'specify shape by the shape of the composition array.'
assert composition.shape[-1] == len(fs.elements), f'composition.shape[-1] ({composition.shape[-1]}) must be {len(fs.elements)}'
self.array_elemental_composition = composition
self.shape: Shape = self.array_elemental_composition.shape[:-1]
self.fs = fs
self.elements = fs.elements
def get_elemental_composition(self) -> dict[str, float]:
"""Get a dict of the molar amount of each element
Returns:
Returns a dict of floats with the molar amount of each element in mol
"""
return {s: c for s, c in zip(self.fs.elements, self.array_elemental_composition)}
def get_mass(self) -> FloatArray | float:
"""Get absolute mass of elements
Returns:
Mass of the fluid in kg
"""
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)
def __sub__(self, other: 'fluid | elements') -> 'elements':
return array_operation(self, other, np.subtract)
def __truediv__(self, other: int | float | FloatArray) -> 'elements':
k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
ttes = self.array_elemental_composition / k
return self.__class__(self.array_elemental_composition / k + ttes, self.fs)
def __mul__(self, other: int | float | FloatArray) -> 'elements':
k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
return self.__class__(self.array_elemental_composition * k, self.fs)
def __rmul__(self, other: int | float | FloatArray) -> 'elements':
k = np.expand_dims(other, -1) if isinstance(other, np.ndarray) else other
return self.__class__(self.array_elemental_composition * k, self.fs)
def __neg__(self) -> 'elements':
return self.__class__(-self.array_elemental_composition, self.fs)
def __array__(self) -> FloatArray:
return self.array_elemental_composition
@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)]
else:
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'
aec = self.array_elemental_composition.reshape(-1, len(self.fs.elements))
return iter({s: c for s, c in zip(self.fs.elements, aec[i, :])} for i in range(aec.shape[0]))
def __repr__(self) -> str:
if len(self.array_elemental_composition.shape) == 1:
lines = [f'{s:16} {c:5.3e} mol' for s, c in zip(self.fs.elements, self.array_elemental_composition)]
return '\n'.join(lines)
else:
array_disp = self.array_elemental_composition.__repr__()
padding = int(array_disp.find('\n') / (len(self.fs.elements) + 1))
return ('Elements:\n' + ' ' * int(padding / 2) +
''.join([(' ' * (padding - len(s))) + s for s in self.fs.elements]) +
'\nMols:\n' + self.array_elemental_composition.__repr__())
def lookup(prop_array: FloatArray,
temperature: FloatArray | float,
t_offset: float) -> FloatArray:
"""linear interpolates values from the given prop_array
Args:
prop_array: Array of the temperature depended property
temperature: Absolute temperature(s) in Kelvin. Must
be broadcastable to prop_array.
Returns:
Interpolates values based on given temperature
"""
t = np.array(temperature) - t_offset
t_lim = np.minimum(np.maximum(0, t), prop_array.shape[0] - 2)
f = np.expand_dims(t - np.floor(t_lim), axis=-1)
ti1 = t_lim.astype(int)
return f * prop_array[ti1 + 1, :] + (1 - f) * prop_array[ti1, :]
def reorder_array(arr: FloatArray, old_index: list[str], new_index: list[str]) -> FloatArray:
"""Reorder the last dimension of an array according to a provided list of species
names in the old oder and a list in the new order.
Args:
arr: Array to be reordered
old_index: List of species names in the current order
new_index: List of species names in the new order
Returns:
Array with the last dimension reordered
"""
ret_array = np.zeros([*arr.shape[:-1], len(new_index)])
for i, k in enumerate(old_index):
ret_array[..., new_index.index(k)] = arr[..., i]
return ret_array
@overload
def array_operation(self: elements, other: elements | fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> elements:
pass
@overload
def array_operation(self: fluid, other: T, func: Callable[[FloatArray, FloatArray], FloatArray]) -> T:
pass
@overload
def array_operation(self: T, other: fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> T:
pass
def array_operation(self: elements | fluid, other: elements | fluid, func: Callable[[FloatArray, FloatArray], FloatArray]) -> elements | fluid:
"""Perform an array operation on two fluid or elements objects.
The operation is provided by a Callable that takes two arguments.
Args:
self: First fluid or elements object
other: Second fluid or elements object
func: Callable function to perform the operation
Returns:
A new fluid or elements object with the result of the
"""
assert isinstance(other, elements) or isinstance(other, fluid)
if self.fs is other.fs:
if isinstance(self, elements) or isinstance(other, elements):
return elements(func(self.array_elemental_composition, other.array_elemental_composition), self.fs)
else:
return fluid(func(self.array_composition, other.array_composition), self.fs)
elif set(self.fs.species) >= set(other.fs.species):
if isinstance(self, elements) or isinstance(other, elements):
el_array = reorder_array(other.array_elemental_composition, other.fs.elements, self.fs.elements)
return elements(func(self.array_elemental_composition, el_array), self.fs)
else:
el_array = reorder_array(other.array_composition, other.fs.species, self.fs.species)
return fluid(func(self.array_composition, el_array), self.fs)
elif set(self.fs.species) < set(other.fs.species):
if isinstance(self, elements) or isinstance(other, elements):
el_array = reorder_array(self.array_elemental_composition, self.fs.elements, other.fs.elements)
return elements(func(el_array, other.array_elemental_composition), other.fs)
else:
el_array = reorder_array(self.array_composition, self.fs.species, other.fs.species)
return fluid(func(el_array, other.array_composition), other.fs)
else:
new_fs = fluid_system(sorted(list(set(self.fs.species) | set(other.fs.species))))
if isinstance(self, elements) or isinstance(other, elements):
el_array1 = reorder_array(self.array_elemental_composition, self.fs.elements, new_fs.elements)
el_array2 = reorder_array(other.array_elemental_composition, other.fs.elements, new_fs.elements)
return elements(func(el_array1, el_array2), new_fs)
else:
el_array1 = reorder_array(self.array_composition, self.fs.species, new_fs.species)
el_array2 = reorder_array(other.array_composition, other.fs.species, new_fs.species)
return fluid(func(el_array1, el_array2), new_fs)

150
src/gaspype/_operations.py
View File

@ -1,150 +0,0 @@
from math import exp
import numpy as np
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:
"""Stack a list of fluid or elements objects along a new axis
Args:
arrays: List of arrays
axis: Axis to stack the fluid objects along
Returns:
A new array object stacked along the new axis
"""
a0 = arrays[0]
assert all(a.fs == a0.fs for a in arrays), 'All objects must have the same fluid system'
assert axis <= len(a0.shape), f'Axis must be smaller or equal to len(shape) ({len(a0.shape)})'
return a0.__class__(np.stack(
[a.array_elemental_composition if isinstance(a, elements) else a.array_composition for a in arrays],
axis=axis), a0.fs)
def concat(arrays: list[T], axis: int = 0) -> T:
"""Concatenate a list of fluid or elements objects along an existing axis
Args:
arrays: List of arrays
axis: Axis to concatenate the fluid objects along
Returns:
A new array object stacked along the specified axis
"""
a0 = arrays[0]
assert all(f.fs == a0.fs for f in arrays), 'All fluid objects must have the same fluid system'
assert axis < len(a0.shape), f'Axis must be smaller than shape len({a0.shape})'
return a0.__class__(np.concatenate(
[a.array_elemental_composition if isinstance(a, elements) else a.array_composition for a in arrays],
axis=axis), a0.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
additional carbon formation is thermodynamic favored. At a value < 1 a
depletion of solid carbon is favored.
Args:
f: Fluid or elements object
t: Temperature in Kelvin
p: Pressure in Pascal
Returns:
The activity of carbon in the fluid
"""
# Values for solid state carbon (graphite) from NIST-JANAF Tables
# https://janaf.nist.gov/pdf/JANAF-FourthEd-1998-Carbon.pdf
# https://janaf.nist.gov/pdf/JANAF-FourthEd-1998-1Vol1-Intro.pdf
# Polynomial is valid for T from 100 to 2500 K
pgef = np.array([-6.76113852E-02, 2.02187857E+00, -2.38664605E+01,
1.43575462E+02, -4.51375503E+02, 6.06219665E+02])
# Gibbs free energy divided by RT for carbon
g_rtc = -np.sum(pgef * np.log(np.expand_dims(t, -1))**np.array([5, 4, 3, 2, 1, 0])) / R
g_rt = f.fs.get_species_g_rt(t)
x = equilibrium(f, t, p).array_fractions
i_co = f.fs.species.index('CO')
i_co2 = f.fs.species.index('CO2')
i_h2 = f.fs.species.index('H2')
i_h2o = f.fs.species.index('H2O')
i_ch4 = f.fs.species.index('CH4')
if min(x[i_co], x[i_co2]) > min([x[i_ch4], x[i_h2o], x[i_h2]]) and min(x[i_co], x[i_co2]) > 0:
# 2 CO -> CO2 + C(s) (Boudouard reaction)
lnalpha = (2 * g_rt[i_co] - (g_rt[i_co2] + g_rtc)) + np.log(
x[i_co]**2 / x[i_co2] * (p / p0))
elif min([x[i_ch4], x[i_h2o], x[i_co]]) > 1E-4:
# CH4 + 2 CO -> 2 H2O + 3 C(s)
lnalpha = ((g_rt[i_ch4] + 2 * g_rt[i_co] - 3 * g_rtc - 2 * g_rt[i_h2o]) + np.log(
x[i_ch4] * x[i_co]**2 / x[i_h2o]**2 * (p / p0))) / 3
elif min(x[i_h2], x[i_ch4]) > 0:
# if x[i_h2] or x[i_ch4] is small compared to the precision of the
# component concentrations the result can be inaccurate
# CH4 -> 2 H2 + C(s)
# CH4 + CO2 -> 2 H2 + 2 CO
# 2 H2O - O2 -> 2 H2
lnalpha = (g_rt[i_ch4] - (2 * g_rt[i_h2] + g_rtc)) + np.log(
x[i_ch4] / x[i_h2]**2 / (p / p0))
elif x[i_h2] == 0:
# equilibrium on carbon side
lnalpha = 10
else:
# equilibrium on non-carbon side
lnalpha = -10
return exp(lnalpha)
def oxygen_partial_pressure(f: fluid | elements, t: float, p: float) -> FloatArray | float:
_oxygen_data = fluid({'O2': 1})
def get_oxygen(x: FloatArray) -> float:
g_rt = f.fs.get_species_g_rt(t)
g_rt_o2 = _oxygen_data.fs.get_species_g_rt(t)[0]
i_co = f.fs.species.index('CO') if 'C' in f.fs.elements else None
i_co2 = f.fs.species.index('CO2') if 'C' in f.fs.elements else None
i_o2 = f.fs.species.index('O2') if 'O2' in f.fs.species else None
i_h2o = f.fs.species.index('H2O') if 'H' in f.fs.elements else None
i_h2 = f.fs.species.index('H2') if 'H' in f.fs.elements else None
i_ch4 = f.fs.species.index('CH4') if 'CH4' in f.fs.species else None
ox_ref_val = max([float(v) for v in (x[i_h2] * x[i_h2o], x[i_co2] * x[i_co], x[i_ch4]) if v.shape == tuple()] + [0])
# print([i_o2, x[i_o2], ox_ref_val])
if i_o2 is not None and x[i_o2] > ox_ref_val:
# print('o O2')
return float(x[i_o2] * p)
elif i_ch4 is not None and x[i_ch4] > x[i_co2] * 100 and x[i_ch4] > x[i_h2o] * 100:
# print('o ch4')
# 2CH4 + O2 <--> 4H2 + 2CO
lnpo2 = 4 * g_rt[i_h2] + 2 * g_rt[i_co] - 2 * g_rt[i_ch4] - g_rt_o2 + np.log(x[i_h2]**4 * x[i_co]**2 / x[i_ch4]**2) - 2 * np.log(p / p0)
elif (i_co is None and i_h2 is not None) or (i_h2 is not None and i_co is not None and (x[i_h2] * x[i_h2o] > x[i_co2] * x[i_co])):
# print('o h2o')
# 2H2 + O2 <--> 2H2O
lnpo2 = 2 * (g_rt[i_h2o] - g_rt[i_h2] + np.log(x[i_h2o] / x[i_h2])) - g_rt_o2 - np.log(p / p0)
else:
assert i_co is not None
# print('o co2')
# 2CO + O2 <--> 2CO2
lnpo2 = 2 * (g_rt[i_co2] - g_rt[i_co] + np.log(x[i_co2] / x[i_co])) - g_rt_o2 - np.log(p / p0)
return exp(lnpo2) * p
x = equilibrium(f, t, p).array_fractions
if len(x.shape):
return np.apply_along_axis(get_oxygen, -1, x)
else:
return get_oxygen(x)

149
src/gaspype/_phys_data.py
View File

@ -1,149 +0,0 @@
import struct
from typing import Generator, Iterator
from dataclasses import dataclass
def split_on_space(data: bytes, offset: int, end: int) -> Iterator[str]:
"""Splits a byte array into ASCII strings based on spaces.
Args:
data: The byte array to split.
offset: The starting index in the byte array.
end: The ending index in the byte array.
Yields:
str: ASCII strings found in the byte array.
"""
start = offset
for i in range(offset, end):
if data[i] == 0x20: # ASCII space character
if start < i:
yield data[start:i].decode('ascii')
start = i + 1
if start < end:
yield data[start:end].decode('ascii')
@dataclass
class SpeciesData:
"""Class to hold the physical data for a species.
Attributes:
name: Name of the species.
composition: Dictionary of species composition with element symbols as keys and their counts as values.
model: Number of polynomial coefficients used in the model.
ref_string: Reference string for the data source.
t_range: List of temperatures nodes marking intervals.
data: List of lists containing physical data for each temperature interval.
"""
name: str
composition: dict[str, int]
model: int
ref_string: str
t_range: list[float]
data: list[list[float]]
class db_reader():
"""Class to read and parse the binary gas phase species database.
Attributes:
names: An iterator over the names of species in the database.
"""
header_len = 8
def __init__(self, inp_data: bytes):
""" Initializes the database reader with binary data.
Args:
inp_data: The binary data of the gas phase species database.
"""
assert inp_data[:4] == b'gapy', 'Unknown data format'
self._bin_data = inp_data
self._name_lengths = struct.unpack('<I', self._bin_data[4:8])[0]
species_names = split_on_space(self._bin_data, db_reader.header_len, db_reader.header_len + self._name_lengths)
self._index = {s: i for i, s in enumerate(species_names)}
@property
def names(self) -> Iterator[str]:
return iter(self._index.keys())
def __iter__(self) -> Generator[SpeciesData, None, None]:
return (d for d in (self.read(species) for species in self._index.keys()) if d is not None)
def __contains__(self, name: str) -> bool:
""" Checks if a species name is present in the database.
Args:
name: The name of the species to check.
Returns:
bool: True if the species name is present
"""
return name in self._index
def __getitem__(self, name: str) -> SpeciesData:
species_data = self.read(name)
assert species_data, f"Species '{name}' not found in the database."
return species_data
def read(self, name: str) -> SpeciesData | None:
""" Reads the physical data for a given species name from the binary data.
Args:
name (str): The name of the species to read data for.
Returns:
phys_data: An instance of the phys_data class containing the physical data.
"""
if name not in self._index:
return None
head_offset = self._name_lengths + db_reader.header_len + self._index[name] * db_reader.header_len
head = struct.unpack('<I4B', self._bin_data[head_offset:head_offset + db_reader.header_len])
offset = head[0]
composition_count = head[1]
model = head[2]
temperature_count = head[3]
ref_string_len = head[4]
td_data_num = (temperature_count - 1) * model
data_len = composition_count * 3 + (temperature_count + td_data_num) * 4 + ref_string_len
format_string = '<' + '2sB' * composition_count + f'{temperature_count}f{td_data_num}f{ref_string_len}s'
bindat = struct.unpack(format_string, self._bin_data[offset:offset + data_len])
comp = {bindat[i * 2].strip().decode('ASCII'): bindat[i * 2 + 1] for i in range(composition_count)}
noffs = composition_count * 2
t_range = list(bindat[noffs:noffs + temperature_count])
noffs += temperature_count
data = [list(bindat[(noffs + i * model):(noffs + (i + 1) * model)]) for i in range(temperature_count - 1)]
ref = bindat[-1].decode('utf-8')
return SpeciesData(name, comp, model, ref, t_range, data)
"""
Atomic weights values are used from CIAAW.
when a single value is given. Available online at
http://www.ciaaw.org/atomic-weights.htm
When a range of values is given in the CIAAW table, the "conventional
atomic weight" from the IUPAC Periodic Table is used. Available
online at https://iupac.org/wp-content/uploads/2018/12/IUPAC_Periodic_Table-01Dec18.pdf
"""
atomic_weights = {'H': 1.008, 'He': 4.002602, 'Li': 6.94, 'Be': 9.0121831, 'B': 10.81, 'C': 12.011,
'N': 14.007, 'O': 15.999, 'F': 18.998403163, 'Ne': 20.1797, 'Na': 22.98976928,
'Mg': 24.305, 'Al': 26.9815384, 'Si': 28.085, 'P': 30.973761998, 'S': 32.06,
'Cl': 35.45, 'Ar': 39.95, 'K': 39.0983, 'Ca': 40.078, 'Sc': 44.955908,
'Ti': 47.867, 'V': 50.9415, 'Cr': 51.9961, 'Mn': 54.938043, 'Fe': 55.845,
'Co': 58.933194, 'Ni': 58.6934, 'Cu': 63.546, 'Zn': 65.38, 'Ga': 69.723,
'Ge': 72.63, 'As': 74.921595, 'Se': 78.971, 'Br': 79.904, 'Kr': 83.798,
'Rb': 85.4678, 'Sr': 87.62, 'Y': 88.90584, 'Zr': 91.224, 'Nb': 92.90637,
'Mo': 95.95, 'Ru': 101.07, 'Rh': 102.90549, 'Pd': 106.42, 'Ag': 107.8682,
'Cd': 112.414, 'In': 114.818, 'Sn': 118.71, 'Sb': 121.76, 'Te': 127.6,
'I': 126.90447, 'Xe': 131.293, 'Cs': 132.90545196, 'Ba': 137.327, 'La': 138.90547,
'Ce': 140.116, 'Pr': 140.90766, 'Nd': 144.242, 'Sm': 150.36, 'Eu': 151.964,
'Gd': 157.25, 'Tb': 158.925354, 'Dy': 162.5, 'Ho': 164.930328, 'Er': 167.259,
'Tm': 168.934218, 'Yb': 173.045, 'Lu': 174.9668, 'Hf': 178.49, 'Ta': 180.94788,
'W': 183.84, 'Re': 186.207, 'Os': 190.23, 'Ir': 192.217, 'Pt': 195.084,
'Au': 196.96657, 'Hg': 200.592, 'Tl': 204.38, 'Pb': 207.2, 'Bi': 208.9804,
'Th': 232.0377, 'Pa': 231.03588, 'U': 238.02891}

167
src/gaspype/_solver.py
View File

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

24
src/gaspype/constants.py
View File

@ -1,24 +0,0 @@
"""
# Constants
This module contains physical constants used in the gas phase calculations.
kB: Boltzmann constant (J/K)
NA: Avogadro's number (1/mol)
R: Ideal gas constant (J/mol/K)
F: Faraday constant (C/mol)
p0: Standard pressure 1e5 Pa
t0: Standard temperature 298.15 K (25 °C)
p_atm: Standard atmosphere 1 atm = 101325 Pa
epsy: Small value for numerical stability (1e-18)
"""
kB = 1.380649e-23 # J/K
NA = 6.02214076e23 # 1/mol
R = kB * NA # J/mol/K
F = 96485.3321233100 # C/mol
p0 = 1e5 # Pa
t0 = 273.15 + 25 # K
p_atm = 101325 # Pa
epsy = 1e-30

90
src/gaspype/data/atomic_prop.yaml Normal file
View File

@ -0,0 +1,90 @@
atomic_weights:
data:
Ag: 107.8682
Al: 26.9815384
Ar: 39.95
As: 74.921595
Au: 196.96657
B: 10.81
Ba: 137.327
Be: 9.0121831
Bi: 208.9804
Br: 79.904
C: 12.011
Ca: 40.078
Cd: 112.414
Ce: 140.116
Cl: 35.45
Co: 58.933194
Cr: 51.9961
Cs: 132.90545196
Cu: 63.546
Dy: 162.5
Er: 167.259
Eu: 151.964
F: 18.998403163
Fe: 55.845
Ga: 69.723
Gd: 157.25
Ge: 72.63
H: 1.008
He: 4.002602
Hf: 178.49
Hg: 200.592
Ho: 164.930328
I: 126.90447
In: 114.818
Ir: 192.217
K: 39.0983
Kr: 83.798
La: 138.90547
Li: 6.94
Lu: 174.9668
Mg: 24.305
Mn: 54.938043
Mo: 95.95
N: 14.007
Na: 22.98976928
Nb: 92.90637
Nd: 144.242
Ne: 20.1797
Ni: 58.6934
O: 15.999
Os: 190.23
P: 30.973761998
Pa: 231.03588
Pb: 207.2
Pd: 106.42
Pr: 140.90766
Pt: 195.084
Rb: 85.4678
Re: 186.207
Rh: 102.90549
Ru: 101.07
S: 32.06
Sb: 121.76
Sc: 44.955908
Se: 78.971
Si: 28.085
Sm: 150.36
Sn: 118.71
Sr: 87.62
Ta: 180.94788
Tb: 158.925354
Te: 127.6
Th: 232.0377
Ti: 47.867
Tl: 204.38
Tm: 168.934218
U: 238.02891
V: 50.9415
W: 183.84
Xe: 131.293
Y: 88.90584
Yb: 173.045
Zn: 65.38
Zr: 91.224
note: Atomic weights values are used from CIAAW.when a single value is given. Available
online athttp://www.ciaaw.org/atomic-weights.htmWhen a range of values is given
in the CIAAW table, the "conventionalatomic weight" from the IUPAC Periodic Table
is used. Availableonline at https://iupac.org/wp-content/uploads/2018/12/IUPAC_Periodic_Table-01Dec18.pdf

4503
src/gaspype/data/therm_prop.yaml Normal file
View File

File diff suppressed because it is too large Load Diff

562
src/gaspype/data/units.yaml Normal file
View File

@ -0,0 +1,562 @@
prefixes:
- symbol: 'y'
name: yocto
value: 1.0e-24
- symbol: z
name: zepto
value: 1.0e-21
- symbol: a
name: atto
value: 1.0e-18
- symbol: f
name: femto
value: 1.0e-15
- symbol: p
name: pico
value: 1.0e-12
- symbol: 'n'
name: nano
value: 1.0e-9
- symbol: µ
name: micro
value: 1.0e-6
- symbol: m
name: milli
value: 1.0e-3
- symbol: c
name: centi
value: 1.0e-2
- symbol: d
name: deci
value: 1.0e-2
- symbol: da
name: deca
value: 1.0e+1
- symbol: h
name: hecto
value: 1.0e+2
- symbol: k
name: kilo
value: 1.0e+3
- symbol: M
name: mega
value: 1.0e+6
- symbol: G
name: giga
value: 1.0e+9
- symbol: T
name: tera
value: 1.0e+12
- symbol: P
name: peta
value: 1.0e+15
- symbol: E
name: exa
value: 1.0e+18
- symbol: Z
name: zetta
value: 1.0e+21
- symbol: 'Y'
name: yotta
value: 1.0e+24
units:
- symbol: Pa
alt_identifier: [pascal]
si_base_unit:
kg: 1
m: -1
s: -2
typical_si_prefixes: [m, h, k, M, G]
si_factor: 1.0
usage: [pressure, strain]
- symbol: 'N'
alt_identifier: [newton]
si_base_unit:
kg: 1
m: 1
s: -2
typical_si_prefixes: [m, k, M, G]
si_factor: 1.0
usage: [force]
- symbol: lbf
alt_identifier: [pound]
si_base_unit:
kg: 1
m: 1
s: -2
si_factor: 4.4482
usage: [force]
- symbol: dyn
alt_identifier: [dyne]
si_base_unit:
kg: 1
m: 1
s: -2
typical_si_prefixes: []
si_factor: 1.0e-5
usage: [force]
- symbol: atm
alt_identifier: [atmosphere]
si_base_unit:
kg: 1
m: -1
s: -2
typical_si_prefixes: []
si_factor: 101325.0
usage: [pressure]
- symbol: t
alt_identifier: [ton]
si_base_unit:
kg: 1
typical_si_prefixes: []
si_factor: 1000.0
usage: [mass]
- symbol: u
alt_identifier: [atomic mass unit]
si_base_unit:
kg: 1
typical_si_prefixes: []
si_factor: 1.66053906660e-27
usage: [mass]
- symbol: g
alt_identifier: [gram]
si_base_unit:
kg: 1
typical_si_prefixes: [k]
si_factor: 0.001
usage: [mass]
- symbol: bar
alt_identifier: [bar]
si_base_unit:
kg: 1
m: -1
s: -2
typical_si_prefixes: []
si_factor: 100000.0
usage: [pressure]
- symbol: kgf
alt_identifier: [kilogram-force]
si_base_unit:
kg: 1
m: 1
s: -2
typical_si_prefixes: []
si_factor: 9.80665
usage: [force]
- symbol: psi
alt_identifier: [pound-force per square inch]
si_base_unit:
kg: 1
m: -1
s: -2
typical_si_prefixes: []
si_factor: 6894.75729
usage: [pressure]
- symbol: oz
alt_identifier: [ounce]
si_base_unit:
kg: 1
typical_si_prefixes: []
si_factor: 0.0283495
usage: [mass]
- symbol: carat
alt_identifier: [carat]
si_base_unit:
kg: 1
typical_si_prefixes: []
si_factor: 0.0002
usage: [mass]
- symbol: grain
alt_identifier: [grain]
si_base_unit:
kg: 1
typical_si_prefixes: []
si_factor: 0.00006479891
usage: [mass]
- symbol: torr
alt_identifier: [torr, mmHg]
si_base_unit:
kg: 1
m: -1
s: -2
typical_si_prefixes: []
si_factor: 133.322
usage: [pressure]
- symbol: mol
alt_identifier: [mole]
si_base_unit:
mol: 1
typical_si_prefixes: [m, k]
si_factor: 1.0
usage: [amount of substance]
- symbol: m
alt_identifier: [meter]
si_base_unit:
m: 1
typical_si_prefixes: [p, 'n', µ, m, c, d, k]
si_factor: 1.0
usage: [length]
- symbol: mm
alt_identifier: [millimeter]
si_base_unit:
m: 1
typical_si_prefixes: []
si_factor: 0.001
usage: [length]
- symbol: in
alt_identifier: [inch]
si_base_unit:
m: 1
typical_si_prefixes: []
si_factor: 0.0254
usage: [length]
- symbol: ft
alt_identifier: [foot]
si_base_unit:
m: 1
typical_si_prefixes: []
si_factor: 0.3048
usage: [length]
- symbol: yd
alt_identifier: [yard]
si_base_unit:
m: 1
typical_si_prefixes: []
si_factor: 0.9144
usage: [length]
- symbol: mi
alt_identifier: [mile]
si_base_unit:
m: 1
typical_si_prefixes: []
si_factor: 1609.34
usage: [length]
- symbol: m²
alt_identifier: [square meter]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 1.0
usage: [area]
- symbol: cm²
alt_identifier: [square centimeter]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 1.0e+4
usage: [area]
- symbol: mm²
alt_identifier: [square millimeter]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 1.0e+6
usage: [area]
- symbol: km²
alt_identifier: [square kilometer]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 1.0e-6
usage: [area]
- symbol: ac
alt_identifier: [acre]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 4046.86
usage: [area]
- symbol: ha
alt_identifier: [hectare]
si_base_unit:
m: 2
typical_si_prefixes: []
si_factor: 10000.0
usage: [area]
- symbol: l
alt_identifier: [liter]
si_base_unit:
m: 3
typical_si_prefixes: ['n', µ, m, c, d, h]
si_factor: 0.001
usage: [volume]
- symbol: m³
alt_identifier: [cubic meter]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 1.0
usage: [volume]
- symbol: cm³
alt_identifier: [cubic centimeter]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 1.0e+6
usage: [volume]
- symbol: dm³
alt_identifier: [cubic decimeter]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 1.0e+3
usage: [volume]
- symbol: mm³
alt_identifier: [cubic millimeter]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 1.0e+9
usage: [volume]
- symbol: km³
alt_identifier: [cubic kilometer]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 1.0e-9
usage: [volume]
- symbol: gal
alt_identifier: [gallon]
si_base_unit:
m: 3
typical_si_prefixes: []
si_factor: 0.00378541
usage: [volume]
- symbol: K
alt_identifier: [kelvin]
si_base_unit:
K: 3
typical_si_prefixes: [m]
si_factor: 1.0
usage: [temperature]
- symbol: C
alt_identifier: [celsius]
si_base_unit:
K: 3
typical_si_prefixes: [m]
si_factor: 1.0
offset: 273.15
usage: [temperature]
- symbol: F
alt_identifier: [fahrenheit]
si_base_unit:
K: 1
typical_si_prefixes: []
si_factor: 0.5555555555
offset: 255.372222
usage: [temperature]
- symbol: R
alt_identifier: [rankine]
si_base_unit:
K: 1
typical_si_prefixes: []
si_factor: 0.5555555555
offset: 0
usage: [temperature]
- symbol: V
alt_identifier: [volt]
si_base_unit:
kg: 1
m: 2
K: -3
A: -1
typical_si_prefixes: ['n', µ, m, k, M]
si_factor: 0.55555555555555
usage: [voltage]
- symbol: A
alt_identifier: [ampere]
si_base_unit:
A: 1
typical_si_prefixes: [p, 'n', µ, m, k]
si_factor: 0.555555555555555
usage: [current]
- symbol: Barrer
alt_identifier: [barrer]
si_base_unit:
mol: 1
s: 1
kg: -1
typical_si_prefixes: []
si_factor: 3.35e-16
usage: [permeability]
- symbol: C
alt_identifier: [coulomb]
si_base_unit:
s: 1
A: 1
typical_si_prefixes: []
si_factor: 1.0
usage: [charge]
- symbol: Ah
alt_identifier: [ampere-hour, amp-hour]
si_base_unit:
s: 1
A: 1
typical_si_prefixes: [m]
si_factor: 3600.0
usage: [charge]
- symbol: F
alt_identifier: [farad]
si_base_unit:
kg: -1
m: -2
s: 4
A: 2
typical_si_prefixes: [p, 'n', µ, m]
si_factor: 1.0
usage: [capacitance]
- symbol: S
alt_identifier: [siemens, moh]
si_base_unit:
kg: -1
m: -2
s: 3
A: 2
typical_si_prefixes: [µ, m]
si_factor: 1.0
usage: [conductance]
- symbol: Ω
alt_identifier: [ohm]
si_base_unit:
kg: 11
m: 2
s: -3
A: -2
typical_si_prefixes: [m, k, M, G]
si_factor: 1.0
usage: [resistance]
- symbol: s
alt_identifier: [second, sec]
si_base_unit:
s: 1
typical_si_prefixes: ['n', µ, m]
si_factor: 1
usage: [time]
- symbol: min
alt_identifier: [minute]
si_base_unit:
s: 1
typical_si_prefixes: []
si_factor: 60
usage: [time]
- symbol: h
alt_identifier: [hour]
si_base_unit:
s: 1
typical_si_prefixes: []
si_factor: 3600
usage: [time]
- symbol: d
alt_identifier: [day]
si_base_unit:
s: 1
typical_si_prefixes: []
si_factor: 86400
usage: [time]
- symbol: wk
alt_identifier: [week]
si_base_unit:
s: 1
typical_si_prefixes: []
si_factor: 604800
usage: [time]
- symbol: yr
alt_identifier: [year]
si_base_unit:
s: 1
typical_si_prefixes: []
si_factor: 31557600
usage: [time]
- symbol: J
alt_identifier: [joule]
si_base_unit:
kg: 1
m: 2
s: -2
typical_si_prefixes: ['n', µ, m, k, M, G, T, P]
si_factor: 1
usage: [energy]
- symbol: cal
alt_identifier: [calorie]
si_base_unit:
kg: 1
m: 2
s: -2
typical_si_prefixes: [k]
si_factor: 4.184
usage: [energy]
- symbol: Wh
alt_identifier: [watt-hour]
si_base_unit:
kg: 1
m: 2
s: -2
typical_si_prefixes: [m, k, M, G, T]
si_factor: 3600
usage: [energy]
- symbol: eV
alt_identifier: [electronvolt]
si_base_unit:
kg: 1
m: 2
s: -2
typical_si_prefixes: [k, M]
si_factor: 1.602176634e-19
usage: [energy]

5
src/gaspype/membrane.py
View File

@ -1,7 +1,4 @@
from .constants import R, F from . import R, F, oxygen_partial_pressure, fluid, elements, FloatArray
from ._operations import oxygen_partial_pressure
from ._main import fluid, elements
from .typing import FloatArray
import numpy as np import numpy as np
# Each O2 molecule is transported as two ions and each ion has a charged of 2e # Each O2 molecule is transported as two ions and each ion has a charged of 2e

0
src/gaspype/py.typed
View File

10
src/gaspype/typing.py
View File

@ -1,10 +0,0 @@
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, ...]

42
tests/benchmark_cp.py
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])

55
tests/benchmark_cp_comp.py
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])

54
tests/benchmark_equalibrium.py
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])

129
tests/md_to_code.py
View File

@ -1,129 +0,0 @@
import re
from typing import Generator, Iterable
from dataclasses import dataclass
import sys
@dataclass
class markdown_segment:
code_block: bool
language: str
text: str
def convert_to(target_format: str, md_filename: str, out_filename: str, language: str = 'python'):
with open(md_filename, "r") as f_in, open(out_filename, "w") as f_out:
segments = segment_markdown(f_in)
if target_format == 'test':
f_out.write('\n'.join(segments_to_test(segments, language)))
elif target_format == 'script':
f_out.write('\n'.join(segments_to_script(segments, language)))
elif target_format == 'striped_markdown':
f_out.write('\n'.join(segments_to_striped_markdown(segments, language)))
else:
raise ValueError('Unknown target format')
def segment_markdown(markdown_file: Iterable[str]) -> Generator[markdown_segment, None, None]:
regex = re.compile(r"(?:^```\s*(?P<language>(?:\w|-)*)$)", re.MULTILINE)
block_language: str = ''
code_block = False
line_buffer: list[str] = []
for line in markdown_file:
match = regex.match(line)
if match:
if line_buffer:
yield markdown_segment(code_block, block_language, ''.join(line_buffer))
line_buffer.clear()
block_language = match.group('language')
code_block = not code_block
else:
line_buffer.append(line)
if line_buffer:
yield markdown_segment(code_block, block_language, '\n'.join(line_buffer))
def segments_to_script(segments: Iterable[markdown_segment], test_language: str = "python") -> Generator[str, None, None]:
for segment in segments:
if segment.code_block:
if segment.language == test_language:
yield segment.text
else:
for line in segment.text.splitlines():
yield '# | ' + line
yield ''
else:
for line in segment.text.strip(' \n').splitlines():
yield '# ' + line
yield ''
def segments_to_striped_markdown(segments: Iterable[markdown_segment], test_language: str = "python") -> Generator[str, None, None]:
for segment in segments:
if segment.code_block:
if segment.language == test_language:
yield "``` " + test_language
yield segment.text
yield "```"
elif segment.language:
for line in segment.text.splitlines():
yield '# | ' + line
yield ''
else:
for line in segment.text.strip(' \n').splitlines():
yield '# ' + line
yield ''
def segments_to_test(segments: Iterable[markdown_segment], script_language: str = "python") -> Generator[str, None, None]:
ret_block_flag = False
yield 'def run_test():'
for segment in segments:
if segment.code_block:
if segment.language == script_language:
lines = [line for line in segment.text.splitlines() if line.strip()]
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])
yield ''
yield ' print("---------------------------------------------------------")'
yield ''
if ret_block_flag:
yield from [' ' + str(line) for line in segment.text.splitlines()[:-1]]
yield f' print("-- Result (({ret_block_flag})):")'
yield f' print(({ret_block_flag}).__repr__().strip())'
else:
yield from [' ' + str(line) for line in segment.text.splitlines()]
elif ret_block_flag:
yield ' ref_str = r"""'
yield from [str(line) for line in segment.text.splitlines()]
yield '"""'
yield f' print("-- Reference (({ret_block_flag})):")'
yield ' print(ref_str.strip())'
yield f' assert ({ret_block_flag}).__repr__().strip() == ref_str.strip()'
ret_block_flag = False
yield '\nif __name__ == "__main__":'
yield ' run_test()'
if __name__ == "__main__":
format = sys.argv[1]
assert format in ['test', 'script']
convert_to(sys.argv[1], sys.argv[2], sys.argv[3])

20
tests/test_db_reader.py
View File

@ -1,20 +0,0 @@
from gaspype._phys_data import db_reader
def test_db_reader():
with open('src/gaspype/data/therm_data.bin', 'rb') as f:
db = db_reader(f.read())
assert 'HCl' in db
assert 'TEST' not in db
assert db['HCl'].name == 'HCl'
assert db['CH4'].composition == {'C': 1, 'H': 4}
assert db['C2H5OH'].composition == {'C': 2, 'H': 6, 'O': 1}
assert db['H2O'].model == 9
for species in db:
print(species.name)
assert species.model == 9
assert len(species.name) > 0
assert len(species.composition) > 0
assert any(el in species.name for el in species.composition.keys())

70
tests/test_doc_examples.py
View File

@ -1,28 +1,52 @@
import md_to_code import re
from glob import glob from typing import Generator
import importlib
import os
def convert_markdown_file(md_filename: str, out_filename: str):
with open(md_filename, "r") as f_in:
with open(out_filename, "w") as f_out:
f_out.write('def run_test():\n')
for block in markdown_to_code([line for line in f_in]):
f_out.write(block + '\n')
def markdown_to_code(lines: list[str], language: str = "python") -> Generator[str, None, None]:
regex = re.compile(
r"(?P<start>^```\s*(?P<block_language>(\w|-)*)\n)(?P<code>.*?\n)(?P<end>```)",
re.DOTALL | re.MULTILINE,
)
blocks = [
(match.group("block_language"), match.group("code"))
for match in regex.finditer("".join(lines))
]
ret_block_flag = False
for block_language, block in blocks:
if block_language == language:
lines = [line for line in block.splitlines() if line.strip()]
ret_block_flag = lines[-1] if '=' not in lines[-1] else None
yield ''
yield ' print("---------------------------------------------------------")'
yield ''
if ret_block_flag:
yield from [' ' + str(line) for line in block.splitlines()[:-1]]
else:
yield from [' ' + str(line) for line in block.splitlines()]
yield f' print("-- Result (({ret_block_flag})):")'
yield f' print(({ret_block_flag}).__repr__().strip())'
elif ret_block_flag:
yield ' ref_str = """'
yield from [str(line) for line in block.splitlines()]
yield '"""'
yield f' print("-- Reference (({ret_block_flag})):")'
yield ' print(ref_str.strip())'
yield f' assert ({ret_block_flag}).__repr__().strip() == ref_str.strip()'
def test_readme(): def test_readme():
md_to_code.convert_to('test', 'README.md', 'tests/autogenerated_readme.py') convert_markdown_file('README.md', 'tests/autogenerated_readme.py')
import autogenerated_readme import autogenerated_readme
autogenerated_readme.run_test() autogenerated_readme.run_test()
def test_example_code():
filter = 'examples/*.md'
files = glob(filter)
for path in files:
file_name = '.'.join(os.path.basename(path).split('.')[:-1])
if not file_name.lower() == 'readme':
print(f"> Test Example {file_name} ...")
md_to_code.convert_to('test', path, f'tests/autogenerated_{file_name}.py')
mod = importlib.import_module(f'autogenerated_{file_name}')
mod.run_test()
if __name__ == "__main__":
test_readme()
test_example_code()

20
tests/test_edge_case_results.py
View File

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

13
tests/test_fluid_references.py
View File

@ -1,13 +0,0 @@
import gaspype as gp
def test_fluid_references():
fs = gp.fluid_system('Cl, CH4, H2O, C2H6, C3H8')
reference_string = """Cl : Hf:Cox,1989. Moore,1971. Moore,1970a. Gordon,1999. [g 7/97]
CH4 : Gurvich,1991 pt1 p44 pt2 p36. [g 8/99]
H2O : Hf:Cox,1989. Woolley,1987. TRC(10/88) tuv25. [g 8/89]
C2H6 : Ethane. Pamidimukkala,1982. [g 7/00]
C3H8 : Hf:TRC(10/85) w1350. Chao,1973. [g 2/00]"""
assert reference_string == fs.get_species_references(), 'fs.get_species_references() == ' + fs.get_species_references()

17
tests/test_results_cantera.py
View File

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

42
tests/test_slicing.py
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)

11
tests/test_to_pandas.py
View File

@ -21,14 +21,3 @@ def test_elements():
df = pd.DataFrame(list(gp.elements(fl * np.array([1, 2, 3, 4])))) df = pd.DataFrame(list(gp.elements(fl * np.array([1, 2, 3, 4]))))
assert df.shape == (4, 2) assert df.shape == (4, 2)
def test_iter():
fl = gp.fluid({'O2': 1, 'H2': 2, 'H2O': 3})
fl2 = fl * np.array([1, 2, 3, 4])
for i, f in enumerate(fl2):
if i == 1:
assert f == {'O2': np.float64(2.0), 'H2': np.float64(4.0), 'H2O': np.float64(6.0)}
if i == 3:
assert f == {'O2': np.float64(4.0), 'H2': np.float64(8.0), 'H2O': np.float64(12.0)}

66
thermo_data/README.md
View File

@ -1,66 +0,0 @@
# Data preparation
Gaspype uses a binary format for thermodynamic data to minimize memory usage and startup time.
## Combine and prepare data
To prepare the data from YAML and XML sources:
``` bash
python combine_data.py combined_data.yaml nasa9*.yaml nasa9*.xml
```
General Syntax is: ```python combine_data.py OUTPUT_YAML_FILE INPUT_FILE_PATTERN_1 INPUT_FILE_PATTERN_2 ...```
Example output format for a single file entry:
``` yaml
- name: HCl
composition: {Cl: 1, H: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [20625.88287, -309.3368855, 5.27541885, -0.00482887422, 6.1957946e-06, -3.040023782e-09, 4.91679003e-13, -10677.82299, -7.309305408]
- [915774.951, -2770.550211, 5.97353979, -0.000362981006, 4.73552919e-08, 2.810262054e-12, -6.65610422e-16, 5674.95805, -16.42825822]
note: Gurvich,1989 pt1 p186 pt2 p93. [tpis89]
```
Compile ```combined_data.yaml``` to the binary gaspype format:
``` bash
python compile_to_bin.py combined_data.yaml ../src/gaspype/data/therm_data.bin
```
General syntax is: ```python compile_to_bin.py YAML_INPUT_FILE BINARY_OUTPUT_FILE```
The binary format is structured like this, it uses little-endian and IEEE 754 floats:
```
[4 Byte magic number: 'gapy']
[4 Byte: 32 Bit integer for length of all species names (NAMES_LENGTH)]
[NAMES_LENGTH Bytes: ASCII encoded string with all species names separated by space]
[Index
[For each species
[4 Bytes: 32 Bit uint with offset of species data in file]
[1 Byte: 8 Bit uint with number of elements]
[1 Byte: 8 Bit uint for polygon length, value = 9]
[1 Byte: 8 Bit uint for number of temperature supporting points (NUM_TEMPS)]
[1 Byte: 8 Bit uint for length of reference string (REF_LEN)]
]
]
[Data
[For each species
[For each species element
[2 Byte: element name in ASCII, 0x20 padded]
[1 Byte: 8 Bit uint for number of atoms]
]
[For Range(NUM_TEMPS)
[4 Byte: 32 Bit float with temperature supporting point]
]
[For Range(NUM_TEMPS - 1)
[36 Bytes: 9 x 32 Bit float with NASA9-Polynomial for a temperature interval]
]
[REF_LEN Bytes: ASCII string of the data reference]
]
]
```
## Notes
- Original source of the data compilation: https://ntrs.nasa.gov/citations/20020085330
- nasa9_*.yaml files are exported from https://cearun.grc.nasa.gov/ThermoBuild/ and
converted with ck2yaml (https://cantera.org/stable/userguide/thermobuild.html)
- nasa9polynomials.xml is from: https://github.com/guillemborrell/thermopy/tree/master/databases

80
thermo_data/combine_data.py
View File

@ -1,80 +0,0 @@
import yaml
import xml.etree.ElementTree as ET
import glob
import sys
def main():
output_file = sys.argv[1]
input_files = sys.argv[2:]
inp_therm_prop: list[dict] = []
for glop_filter in input_files:
for file_name in glob.glob(glop_filter):
print(f'Processing file: {file_name}...')
if file_name.endswith('.xml'):
inp_therm_prop += get_xml_data(file_name)
else:
with open(file_name, 'r') as f:
inp_therm_prop += yaml.safe_load(f)['species']
therm_prop = []
added_species_names: set[str] = set()
for species in inp_therm_prop:
name = species['name']
assert isinstance(name, str)
for element in species['composition']:
name = name.replace(element.upper(), element)
if name not in added_species_names and 'E' not in species['composition']:
species['name'] = name
therm_prop.append(species)
added_species_names.add(name)
with open(output_file, 'w') as f:
f.write('species:\n')
for species in sorted(therm_prop, key=lambda s: s['name']):
f.write('\n- ' + yaml.dump({'name': species['name']}, default_flow_style=False))
f.write(f' composition: {yaml.dump(species["composition"], default_flow_style=True)}')
f.write(' thermo:\n')
f.write(f' model: {species["thermo"]["model"]}\n')
f.write(f' temperature-ranges: {yaml.dump(species["thermo"]["temperature-ranges"], default_flow_style=True)}')
f.write(' data:\n')
for d in species['thermo']['data']:
f.write(f' - {d}\n')
f.write(' ' + yaml.dump({'note': species['thermo']['note']}, default_flow_style=False))
print('Added: ', ', '.join(sorted(added_species_names)))
def get_xml_data(file_name: str) -> list[dict]:
xml_data_list: list[dict] = []
tree = ET.parse(file_name)
root = tree.getroot()
for i, child in enumerate(root):
elements = {el[0]: int(el[1]) for c in child.find('elements').findall('element') for el in c.items()}
values_temp = [tr for tr in child.findall('T_range')]
t_ranges = [float(v.attrib['Tlow']) for v in values_temp] + [float(values_temp[-1].attrib['Thigh'])]
data = [[float(v.text) for v in tr] for tr in values_temp]
is_gas = child.find('condensed').text == 'False'
if is_gas:
xml_data_list.append({
'name': child.attrib['inp_file_name'],
'composition': elements,
'thermo': {
'model': 'NASA9',
'temperature-ranges': t_ranges,
'data': data,
'note': child.find('comment').text
}
})
return xml_data_list
if __name__ == "__main__":
main()

71
thermo_data/compile_to_bin.py
View File

@ -1,71 +0,0 @@
import yaml
import struct
import sys
import os
def main():
input_file = sys.argv[1]
output_file = sys.argv[2]
assert not output_file.endswith('.yml') and not output_file.endswith('.yaml'), 'Binary output file should not have yaml-extension'
with open(input_file) as f:
data = yaml.safe_load(f)
data = list({v['name']: v for v in data['species']}.values())
species_names = ' '.join(s['name'] for s in data)
header_len = 8
offset = 8 + len(species_names) + len(data) * header_len
header_list = []
body_list = []
added_list = set()
for dat in data:
if dat['name'] not in added_list:
composition_count = len(dat['composition'])
model = {'NASA9': 9}[dat['thermo']['model']]
temperatures = dat['thermo']['temperature-ranges']
ref_string = dat['thermo']['note'].encode('utf-8')
header = struct.pack('<I4B', offset, composition_count, model, len(temperatures), len(ref_string))
assert len(header) == header_len
header_list.append(header)
composition = [el for k, v in dat['composition'].items() for el in [k.ljust(2).encode('ASCII'), v]]
data_vals = [d for darr in dat['thermo']['data'] for d in darr[:model]]
assert len(dat['thermo']['data']) == len(temperatures) - 1, f"Temperature data length mismatch for {dat['name']}."
if any(len(d) != model for d in dat['thermo']['data']):
print(f"Warning: Data length mismatch for {dat['name']}. Expected {model} coefficients, got {len(dat['thermo']['data'][0])}.")
format_string = '<' + '2sB' * composition_count + f'{len(temperatures)}f{len(data_vals)}f{len(ref_string)}s'
body = struct.pack(format_string, *(composition + temperatures + data_vals + [ref_string]))
body_list.append(body)
offset += len(body)
added_list.add(dat['name'])
# Ensure the output directory exists
os.makedirs(os.path.dirname(output_file), exist_ok=True)
with open(output_file, 'wb') as f:
f.write(b'gapy')
f.write(struct.pack('<I', len(species_names)))
f.write(species_names.encode('ASCII'))
for dat in header_list:
f.write(dat)
for dat in body_list:
f.write(dat)
if __name__ == '__main__':
main()

2805
thermo_data/nasa9_hydrocarbons.yaml
View File

File diff suppressed because it is too large Load Diff

293
thermo_data/nasa9_hydrocarbons_si.yaml
View File

@ -1,293 +0,0 @@
generator: ck2yaml
input-files: [mythermo_si.txt]
cantera-version: 3.0.1
date: Sun, 04 Aug 2024 14:04:08 +0000
units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}
species:
- name: H
composition: {H: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, 2.547370801e+04, -0.446682853]
- [60.7877425, -0.1819354417, 2.500211817, -1.226512864e-07, 3.73287633e-11,
-5.68774456e-15, 3.410210197e-19, 2.547486398e+04, -0.448191777]
- [2.173757694e+08, -1.312035403e+05, 33.991742, -3.81399968e-03, 2.432854837e-07,
-7.69427554e-12, 9.64410563e-17, 1.067638086e+06, -274.2301051]
note: D0(H2):Herzberg,1970. Moore,1972. Gordon,1999. [g 6/97]
- name: HO2
composition: {H: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-7.59888254e+04, 1329.383918, -4.67738824, 0.02508308202, -3.006551588e-05,
1.895600056e-08, -4.82856739e-12, -5873.35096, 51.9360214]
- [-1.810669724e+06, 4963.19203, -1.039498992, 4.56014853e-03, -1.061859447e-06,
1.144567878e-10, -4.76306416e-15, -3.20081719e+04, 40.6685092]
note: Hf:Hills,1984 & NASA data. Jacox,1998 p153. [g 4/02]
- name: H2
composition: {H: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [4.07832321e+04, -800.918604, 8.21470201, -0.01269714457, 1.753605076e-05,
-1.20286027e-08, 3.36809349e-12, 2682.484665, -30.43788844]
- [5.60812801e+05, -837.150474, 2.975364532, 1.252249124e-03, -3.74071619e-07,
5.9366252e-11, -3.6069941e-15, 5339.82441, -2.202774769]
- [4.96688412e+08, -3.147547149e+05, 79.8412188, -8.41478921e-03, 4.75324835e-07,
-1.371873492e-11, 1.605461756e-16, 2.488433516e+06, -669.572811]
note: Ref-Elm. Gurvich,1978 pt1 p103 pt2 p31. [tpis78]
- name: H2O
composition: {H: 2, O: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-3.94796083e+04, 575.573102, 0.931782653, 7.22271286e-03, -7.34255737e-06,
4.95504349e-09, -1.336933246e-12, -3.30397431e+04, 17.24205775]
- [1.034972096e+06, -2412.698562, 4.64611078, 2.291998307e-03, -6.83683048e-07,
9.42646893e-11, -4.82238053e-15, -1.384286509e+04, -7.97814851]
note: Hf:Cox,1989. Woolley,1987. TRC(10/88) tuv25. [g 8/89]
- name: H2O2
composition: {H: 2, O: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-9.27953358e+04, 1564.748385, -5.97646014, 0.0327074452, -3.93219326e-05,
2.509255235e-08, -6.46504529e-12, -2.494004728e+04, 58.7717418]
- [1.489428027e+06, -5170.82178, 11.2820497, -8.04239779e-05, -1.818383769e-08,
6.94726559e-12, -4.8278319e-16, 1.418251038e+04, -46.5085566]
note: Hf:Gurvich,1989 pt1 p127. Gurvich,1978 pt1 p121. [g 6/99]
- name: O
composition: {O: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [-7953.6113, 160.7177787, 1.966226438, 1.01367031e-03, -1.110415423e-06,
6.5175075e-10, -1.584779251e-13, 2.840362437e+04, 8.40424182]
- [2.619020262e+05, -729.872203, 3.31717727, -4.28133436e-04, 1.036104594e-07,
-9.43830433e-12, 2.725038297e-16, 3.39242806e+04, -0.667958535]
- [1.779004264e+08, -1.082328257e+05, 28.10778365, -2.975232262e-03,
1.854997534e-07, -5.79623154e-12, 7.191720164e-17, 8.89094263e+05,
-218.1728151]
note: D0(O2):Brix,1954. Moore,1976. Gordon,1999. [g 5/97]
- name: OH
composition: {O: 1, H: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [-1998.85899, 93.0013616, 3.050854229, 1.529529288e-03, -3.157890998e-06,
3.31544618e-09, -1.138762683e-12, 2991.214235, 4.67411079]
- [1.017393379e+06, -2509.957276, 5.11654786, 1.30529993e-04, -8.28432226e-08,
2.006475941e-11, -1.556993656e-15, 2.019640206e+04, -11.01282337]
- [2.847234193e+08, -1.859532612e+05, 50.082409, -5.14237498e-03, 2.875536589e-07,
-8.22881796e-12, 9.56722902e-17, 1.468393908e+06, -402.355558]
note: 'D0(H-OH): Ruscic,2002. Gurvich,1978 pt1 p110 pt2 p37. [g 4/02]'
- name: O2
composition: {O: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [-3.42556342e+04, 484.700097, 1.119010961, 4.29388924e-03, -6.83630052e-07,
-2.0233727e-09, 1.039040018e-12, -3391.45487, 18.4969947]
- [-1.037939022e+06, 2344.830282, 1.819732036, 1.267847582e-03, -2.188067988e-07,
2.053719572e-11, -8.19346705e-16, -1.689010929e+04, 17.38716506]
- [4.9752943e+08, -2.866106874e+05, 66.9035225, -6.16995902e-03, 3.016396027e-07,
-7.4214166e-12, 7.27817577e-17, 2.293554027e+06, -553.062161]
note: Ref-Elm. Gurvich,1989 pt1 p94 pt2 p9. [tpis89]
- name: O3
composition: {O: 3}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-1.282314507e+04, 589.821664, -2.547496763, 0.02690121526, -3.52825834e-05,
2.312290922e-08, -6.04489327e-12, 1.348368701e+04, 38.5221858]
- [-3.86966248e+07, 1.023344994e+05, -89.615516, 0.0370614497, -4.13763874e-06,
-2.725018591e-10, 5.24818811e-14, -6.51791818e+05, 702.910952]
note: Gurvich,1989 pt1 p101 pt2 p15. [g 8/01]
- name: Si
composition: {Si: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04]
data:
- [98.3614081, 154.6544523, 1.87643667, 1.320637995e-03, -1.529720059e-06,
8.95056277e-10, -1.95287349e-13, 5.26351031e+04, 9.69828888]
- [-6.16929885e+05, 2240.683927, -0.444861932, 1.710056321e-03, -4.10771416e-07,
4.55888478e-11, -1.889515353e-15, 3.95355876e+04, 26.79668061]
- [-9.28654894e+08, 5.44398989e+05, -120.6739736, 0.01359662698, -7.60649866e-07,
2.149746065e-11, -2.474116774e-16, -4.29379212e+06, 1086.382839]
note: Hf:Cox,1989. NIST data version1.1 [Online]1997. Gordon,1999. [g
8/97]
- name: SiH
composition: {Si: 1, H: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-6426.6763, 74.1725121, 3.9734916, -4.14940888e-03, 1.022918384e-05,
-8.59238636e-09, 2.567093743e-12, 4.2817458e+04, 2.24693715]
- [4.04208649e+05, -2364.796524, 7.62749914, -2.496591233e-03, 1.10843641e-06,
-1.943991955e-10, 1.136251507e-14, 5.70473768e+04, -24.48054429]
note: Gurvich,1991 pt1 p257 pt2 p234. [tpis91]
- name: SiH2
composition: {Si: 1, H: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-2.063863564e+04, 330.58622, 2.099271145, 3.54253937e-03, 3.37887667e-06,
-5.38384562e-09, 2.081191273e-12, 3.011784298e+04, 12.8233357]
- [4.62403937e+06, -1.14343611e+04, 12.6488087, 9.1148995e-04, -8.76661154e-07,
1.646297357e-10, -9.96509037e-15, 1.072475101e+05, -66.0607807]
note: Gurvich,1991 pt1 p260. Fredin,1985. Dubois,1968. [g 3/01]
- name: SiH3
composition: {Si: 1, H: 3}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [4341.14282, 227.7185085, 0.650825035, 0.01221438558, -4.34760427e-06,
-1.774916828e-09, 1.184191367e-12, 2.259993826e+04, 19.68347482]
- [6.05632122e+05, -4721.25406, 13.29129523, -1.256824868e-03, 2.68828594e-07,
-3.010741582e-11, 1.370945857e-15, 4.97442064e+04, -61.405031]
note: Gurvich,1991 pt1 p261 pt2 p236. [g 3/99]
- name: SiH4
composition: {Si: 1, H: 4}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [7.87299329e+04, -552.608705, 2.498944303, 0.01442118274, -8.46710731e-06,
2.726164641e-09, -5.43675437e-13, 6269.66906, 4.96546183]
- [1.29037874e+06, -7813.39978, 18.28851664, -1.975620946e-03, 4.15650215e-07,
-4.59674561e-11, 2.072777131e-15, 4.76688795e+04, -98.0169746]
note: Silane. Gurvich,1991 pt1 p263 pt 2 p237. [tpis91]
- name: SiO
composition: {Si: 1, O: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-4.72277105e+04, 806.313764, -1.636976133, 0.01454275546, -1.723202046e-05,
1.04239734e-08, -2.559365273e-12, -1.666585903e+04, 33.557957]
- [-1.765134162e+05, -31.9917709, 4.47744193, 4.59176471e-06, 3.55814315e-08,
-1.327012559e-11, 1.613253297e-15, -1.35084236e+04, -0.838695733]
note: Gurvich,1991 pt1 p247 pt2 p227 [tpis91]
- name: SiO2
composition: {Si: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-3.36294878e+04, 473.407892, 0.2309770671, 0.01850230806, -2.242786671e-05,
1.364981554e-08, -3.35193503e-12, -4.22648749e+04, 22.95803206]
- [-1.464031193e+05, -626.144106, 7.96456371, -1.854119096e-04, 4.09521467e-08,
-4.69720676e-12, 2.17805428e-16, -3.79183477e+04, -20.45285414]
note: Gurvich,1991 pt1 p256 pt2 p233. [tpis91]
- name: Si2
composition: {Si: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [1.237596221e+04, -102.4904376, 4.35484852, 1.281063335e-03, -2.531991623e-06,
2.265694244e-09, -7.00129014e-13, 6.90694285e+04, 3.2511252]
- [1.370060657e+06, -4207.06004, 9.33743289, -2.749217168e-03, 9.58634596e-07,
-1.372449748e-10, 6.7650281e-15, 9.51088454e+04, -31.6838519]
note: Gurvich,1991 pt1 p240 pt2 p225. [tpis91]
- name: Si3
composition: {Si: 3}
thermo:
model: NASA9
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [-1.114208177e+04, 157.5785843, 2.486135003, 0.01631637255, -2.208240021e-05,
1.372008287e-08, -3.2623307e-12, 7.32825385e+04, 15.88081347]
- [-1.699395561e+06, 4697.81538, 2.618198124, 1.959082075e-03, -2.581160603e-07,
6.10344486e-12, 6.08630924e-16, 4.27791681e+04, 25.86540384]
note: Gurvich,1991 pt1 p246 pt2 p226. [g 7/95]
- name: H2O(cr)
composition: {H: 2, O: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 273.15]
data:
- [-4.02677748e+05, 2747.887946, 57.3833663, -0.826791524, 4.41308798e-03,
-1.054251164e-05, 9.69449597e-09, -5.53031499e+04, -190.2572063]
note: Ice. Gordon,1982. [g11/99]
- name: H2O(L)
composition: {H: 2, O: 1}
thermo:
model: NASA9
temperature-ranges: [273.15, 373.15, 600.0]
data:
- [1.326371304e+09, -2.448295388e+07, 1.879428776e+05, -767.899505,
1.761556813, -2.151167128e-03, 1.092570813e-06, 1.101760476e+08, -9.77970097e+05]
- [1.263631001e+09, -1.680380249e+07, 9.27823479e+04, -272.237395, 0.447924376,
-3.91939743e-04, 1.425743266e-07, 8.11317688e+07, -5.13441808e+05]
note: Liquid. Cox,1989. Haar,1984. Keenan,1984. Stimson,1969. [g 8/01]
- name: Si(cr)
composition: {Si: 1}
thermo:
model: NASA9
temperature-ranges: [200.0, 298.15, 1690.0]
data:
- [-2.323538208e+04, 0.0, 2.10202168, 1.809220552e-03, 0.0, 0.0, 0.0,
-785.063521, -10.38427318]
- [-5.23255974e+04, 0.0, 2.850169415, 3.97516697e-04, 0.0, 0.0, 0.0,
-1042.947234, -14.38964187]
note: Cubic. Ref-Elm. Gurvich,1991 pt1 p236 pt2 p220. [tpis91]
- name: Si(L)
composition: {Si: 1}
thermo:
model: NASA9
temperature-ranges: [1690.0, 6000.0]
data:
- [0.0, 0.0, 3.271389414, 0.0, 0.0, 0.0, 0.0, 4882.66711, -13.26611073]
note: Liquid. Ref-Elm. Gurvich,1991 pt1 p236 pt2 p220. [tpis91]
- name: SiO2(a-qz)
composition: {Si: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [200.0, 848.0]
data:
- [-5.7768955e+05, 7214.66111, -31.45730294, 0.0741217715, -8.67007782e-06,
-1.080461312e-07, 8.31632491e-11, -1.462398375e+05, 184.2424399]
note: Alpha-quartz,hexagonal. Gurvich,1991 pt1 p250 pt2 p228. [tpis91]
- name: SiO2(b-qz)
composition: {Si: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [848.0, 1200.0]
data:
- [2.317635074e+04, 0.0, 7.026511484, 1.241925261e-03, 0.0, 0.0, 0.0,
-1.117012474e+05, -35.80751356]
note: Beta-quartz,hexagonal. Gurvich,1991 pt1 p250 pt2 p228. [tpis91]
- name: SiO2(b-crt)
composition: {Si: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [1200.0, 1996.0]
data:
- [-5.356419079e+05, 0.0, 9.331036946, -7.306503931e-04, 3.339944266e-07,
0.0, 0.0, -1.134326721e+05, -49.98768383]
note: Beta-cristobalite,cubic. Gurvich,1991 pt1 p250 pt2 p228. [tpis91]
- name: SiO2(L)
composition: {Si: 1, O: 2}
thermo:
model: NASA9
temperature-ranges: [1996.0, 6000.0]
data:
- [0.0, 0.0, 10.04268442, 0.0, 0.0, 0.0, 0.0, -1.140002976e+05, -55.54279592]
note: Liquid. Gurvich,1991 pt1 p250 pt2 p228. [tpis91]

3370
thermo_data/nasa9_hydrocarbons_si_halogens.yaml
View File

File diff suppressed because it is too large Load Diff

1
thermo_data/nasa9polynomials.xml
View File

File diff suppressed because one or more lines are too long