Initial commit

This commit is contained in:
Nicolas Kruse 2025-05-25 23:23:02 +02:00
commit 9143bfef5b
10 changed files with 840 additions and 0 deletions

24
.flake8 Normal file
View File

@ -0,0 +1,24 @@
[flake8]
# Specify the maximum allowed line length
max-line-length = 88
# Ignore specific rules
# For example, E501: Line too long, W503: Line break before binary operator
ignore = E501, W503, W504, E226, E265
# Exclude specific files or directories
exclude =
.git,
__pycache__,
build,
dist,
.conda,
tests/autogenerated_*
# Enable specific plugins or options
# Example: Enabling flake8-docstrings
select = C,E,F,W,D
# Specify custom error codes to ignore or enable
per-file-ignores =
tests/*: D

10
.gitignore vendored Normal file
View File

@ -0,0 +1,10 @@
backend/index/*
__pycache__
*.code-workspace
*.egg.info
/src/*.egg-info/*
/dist/*
.vscode
.mypy_cache
.pytest_cache
tests/autogenerated_*.py

18
CITATION.cff Normal file
View File

@ -0,0 +1,18 @@
cff-version: 1.2.0
title: Gaspype
abstract: Gaspype is a performant library for thermodynamic calculations with ideal gases
authors:
- family-names: Kruse
given-names: Nicolas
orcid: "https://orcid.org/0000-0001-6758-2269"
affiliation: "German Aerospace Center (DLR)"
address: "Linder Höhe"
city: Köln
version: 1.0.1
date-released: "2025-05-09"
#identifiers:
# - description: This is the collection of archived snapshots of all versions of Gaspype
# type: doi
# value: ""
license: MIT License
repository-code: "https://github.com/DLR-Institute-of-Future-Fuels/gaspype"

21
LICENSE Normal file
View File

@ -0,0 +1,21 @@
MIT License
Copyright (c) 2025 Nicolas Kruse, German Aerospace Center (DLR)
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

228
README.md Normal file
View File

@ -0,0 +1,228 @@
# Gaspype
The python package provides an performant library for thermodynamic calculations
like equilibrium reactions for several hundred gas species and their mixtures -
written in Python/Numpy.
Species are treated as ideal gases. Therefore the application is limited to moderate
pressures or high temperature applications.
Its designed with goal to be portable to Numpy-style GPU frameworks like JAX and PyTorch.
## Key features
- Tensor operations to prevent bottlenecks by the python interpreter
- Immutable types
- Elegant pythonic interface
- Great readable and compact syntax when using this package
- Good usability in Jupyter Notebook
- High performance for multidimensional fluid arrays
## Installation
Installation with pip:
``` bash
pip install gaspype
```
Installation with conda:
``` bash
conda install conda-forge::gaspype
```
## Getting started
Gaspype provides two main classes: ```fluid``` and ```elements```.
### Fluid
A fluid class describes a mixture of molecular species and their individual molar amounts.
``` python
import gaspype as gp
fl = gp.fluid({'H2O': 1, 'H2': 2})
fl
```
```
Total 3.000e+00 mol
H2O 33.33 %
H2 66.67 %
```
Its' functions provides thermodynamic, mass balance and ideal gas properties of the mixture.
``` python
cp = fl.get_cp(t=800+273.15)
mass = fl.get_mass()
gas_volume = fl.get_v(t=800+273.15, p=1e5)
```
The arguments can be provided as numpy-arrays:
``` python
import numpy as np
t_range = np.linspace(600, 800, 5) + 273.15
fl.get_density(t=t_range, p=1e5)
```
```
array([0.10122906, 0.09574625, 0.09082685, 0.08638827, 0.08236328])
```
A ```fluid``` object can have multiple compositions. A multidimensional ```fluid``` object can be created for example by multiplication with a numpy array:
``` python
fl2 = gp.fluid({'H2O': 1, 'N2': 2}) + \
np.linspace(0, 10, 4) * gp.fluid({'H2': 1})
fl2
```
```
Total mol:
array([ 3. , 6.33333333, 9.66666667, 13. ])
Species:
H2 H2O N2
Molar fractions:
array([[0. , 0.33333333, 0.66666667],
[0.52631579, 0.15789474, 0.31578947],
[0.68965517, 0.10344828, 0.20689655],
[0.76923077, 0.07692308, 0.15384615]])
```
A fluid object can be converted to a pandas dataframe:
``` python
import pandas as pd
pd.DataFrame(list(fl2))
```
| | H2O | N2 | H2
|----|-----|-----|-------
|0 | 1.0 | 2.0 | 0.000000
|1 | 1.0 | 2.0 | 3.333333
|2 | 1.0 | 2.0 | 6.666667
|3 | 1.0 | 2.0 | 10.000000
The broadcasting behavior is not limited to 1D-arrays:
``` python
fl3 = gp.fluid({'H2O': 1}) + \
np.linspace(0, 10, 4) * gp.fluid({'H2': 1}) + \
np.expand_dims(np.linspace(1, 3, 3), axis=1) * gp.fluid({'N2': 1})
fl3
```
```
Total mol:
array([[ 2. , 5.33333333, 8.66666667, 12. ],
[ 3. , 6.33333333, 9.66666667, 13. ],
[ 4. , 7.33333333, 10.66666667, 14. ]])
Species:
H2 H2O N2
Molar fractions:
array([[[0. , 0.5 , 0.5 ],
[0.625 , 0.1875 , 0.1875 ],
[0.76923077, 0.11538462, 0.11538462],
[0.83333333, 0.08333333, 0.08333333]],
[[0. , 0.33333333, 0.66666667],
[0.52631579, 0.15789474, 0.31578947],
[0.68965517, 0.10344828, 0.20689655],
[0.76923077, 0.07692308, 0.15384615]],
[[0. , 0.25 , 0.75 ],
[0.45454545, 0.13636364, 0.40909091],
[0.625 , 0.09375 , 0.28125 ],
[0.71428571, 0.07142857, 0.21428571]]])
```
### 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
el = gp.elements({'N': 1, 'Cl': 2})
el.get_mass()
```
```
np.float64(0.08490700000000001)
```
A ```elements``` object can be as well instantiated from a ```fluid``` object. Arithmetic operations between ```elements``` and ```fluid``` result in an ```elements``` object:
``` python
el2 = gp.elements(fl) + el - 0.3 * fl
el2
```
```
Cl 2.000e+00 mol
H 4.200e+00 mol
N 1.000e+00 mol
O 7.000e-01 mol
```
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:
``` python
fs = gp.fluid_system('CH4, H2, CO, CO2, O2')
el3 = gp.elements({'C': 1, 'H': 2, 'O':1}, fs)
fl3 = gp.equilibrium(el3, t=800)
fl3
```
```
Total 1.204e+00 mol
CH4 33.07 %
H2 16.93 %
CO 16.93 %
CO2 33.07 %
O2 0.00 %
```
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:
``` python
fl3 + gp.fluid({'CH4': 1}, fs)
```
```
Total 2.204e+00 mol
CH4 63.44 %
H2 9.24 %
CO 9.24 %
CO2 18.07 %
O2 0.00 %
```
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:
``` python
fl3 + gp.fluid({'NH3': 1})
```
```
Total 2.204e+00 mol
CH4 18.07 %
CO 9.24 %
CO2 18.07 %
H2 9.24 %
NH3 45.38 %
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]
```
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.

57
pyproject.toml Normal file
View File

@ -0,0 +1,57 @@
[project]
name = "copapy"
version = "1.0.0"
authors = [
{ name="Nicolas Kruse", email="nicolas.kruse@nonan.net" },
]
description = "Copy-Patch Compiler"
readme = "README.md"
requires-python = ">=3.10"
license = "MIT"
classifiers = [
"Programming Language :: Python :: 3",
"Operating System :: OS Independent",
]
dependencies = [
"pelfy>=1.0.2",
]
[project.urls]
Homepage = "https://github.com/nonannet/copapy"
Issues = "https://github.com/nonannet/copapy/issues"
[build-system]
requires = ["setuptools>=61.0", "wheel"]
build-backend = "setuptools.build_meta"
[tool.setuptools.packages.find]
where = ["src"]
[tool.setuptools.package-data]
copapy = ["*.c"]
[project.optional-dependencies]
dev = [
"flake8",
"mypy",
"pytest",
"pandas",
"cantera",
"types-PyYAML",
"scipy-stubs"
]
[tool.mypy]
files = ["src"]
strict = true
warn_return_any = true
warn_unused_configs = true
check_untyped_defs = true
no_implicit_optional = true
show_error_codes = true
[tool.pytest.ini_options]
minversion = "6.0"
addopts = "-ra -q"
testpaths = ["tests"]
pythonpath = ["src"]

230
src/copapy/__init__.py Normal file
View File

@ -0,0 +1,230 @@
import re
import pkgutil
def get_var_name(var, scope=globals()):
return [name for name, value in scope.items() if value is var]
def _get_c_function_definitions(code: str):
ret = re.findall(r".*?void\s+([a-z_1-9]*)\s*\([^\)]*?\)[^\}]*?\{[^\}]*?result_([a-z_]*)\(.*?", code, flags=re.S)
return {r[0]: r[1] for r in ret}
_function_definitions = _get_c_function_definitions(pkgutil.get_data(__name__, 'ops.c').decode('utf-8'))
class Node:
def __repr__(self):
#return f"Node:{self.name}({', '.join(str(a) for a in self.args) if self.args else self.value})"
return f"Node:{self.name}({', '.join(str(a) for a in self.args) if self.args else (self.value if isinstance(self, Const) else '')})"
class Device():
pass
class Net:
def __init__(self, dtype: str, source: Node):
self.dtype = dtype
self.source = source
def __mul__(self, other):
return _add_op('mul', [self, other], True)
def __rmul__(self, other):
return _add_op('mul', [self, other], True)
def __add__(self, other):
return _add_op('add', [self, other], True)
def __radd__(self, other):
return _add_op('add', [self, other], True)
def __sub__ (self, other):
return _add_op('sub', [self, other])
def __rsub__ (self, other):
return _add_op('sub', [other, self])
def __truediv__ (self, other):
return _add_op('div', [self, other])
def __rtruediv__ (self, other):
return _add_op('div', [other, self])
def __repr__(self):
names = get_var_name(self)
return f"{'name:' + names[0] if names else 'id:' + str(id(self))[-5:]}"
class Const(Node):
def __init__(self, value: float | int | bool):
self.dtype, self.value = _get_data_and_dtype(value)
self.name = 'const_' + self.dtype
#if self.name not in _function_definitions:
# raise ValueError(f"Unsupported operand type for a const: {self.dtype}")
self.args = []
class Write(Node):
def __init__(self, net: Net):
self.name = 'write_' + net.dtype
self.args = [net]
#if self.name not in _function_definitions:
# raise ValueError(f"Unsupported operand type for write: {net.dtype}")
class Op(Node):
def __init__(self, typed_op_name: str, args: list[Net]):
assert not args or any(isinstance(t, Net) for t in args), 'args parameter must be of type list[Net]'
self.name = typed_op_name
self.args = args
def _add_op(op: str, args: list[Net], commutative = False):
args = [a if isinstance(a, Net) else const(a) for a in args]
if commutative:
args = sorted(args, key=lambda a: a.dtype)
typed_op = '_'.join([op] + [a.dtype for a in args])
if typed_op not in _function_definitions:
raise ValueError(f"Unsupported operand type(s) for {op}: {' and '.join([a.dtype for a in args])}")
result_type = _function_definitions[typed_op]
result_net = Net(result_type, Op(typed_op, args))
return result_net
#def read_input(hw: Device, test_value: float):
# return Net(type(value))
def const(value: float | int | bool):
new_const = Const(value)
return Net(new_const.dtype, new_const)
def _get_data_and_dtype(value):
if isinstance(value, int):
return ('int', int(value))
elif isinstance(value, float):
return ('float', float(value))
elif isinstance(value, bool):
return ('bool', int(value))
else:
raise ValueError(f'Non supported data type: {type(value).__name__}')
def const_vector3d(x: float, y: float, z: float):
return vec3d((const(x), const(y), const(z)))
class vec3d:
def __init__(self, value: tuple[float, float, float]):
self.value = value
def __add__(self, other):
return vec3d(tuple(a+b for a,b in zip(self.value, other.value)))
def get_multiuse_nets(root: list[Op]):
"""
Finds all nets that get accessed more than one time. Therefore
storage on the heap might be better.
"""
known_nets: set[Net] = set()
def recursiv_node_search(net_list: list[Net]):
for net in net_list:
#print(net)
if net in known_nets:
yield net
else:
known_nets.add(net)
yield from recursiv_node_search(net.source.args)
return set(recursiv_node_search(op.args[0] for op in root))
def get_path_segments(root: list[Op]) -> list[list[Op]]:
"""
List of all possible paths. Ops in order of execution (output at the end)
"""
def recursiv_node_search(op_list: list[Op], path: list[Op]) -> list[Op]:
for op in op_list:
new_path = [op] + path
if op.args:
yield from recursiv_node_search([net.source for net in op.args], new_path)
else:
yield new_path
known_ops: set[Op] = set()
sorted_path_list = sorted(recursiv_node_search(root, []), key=lambda x: -len(x))
for path in sorted_path_list:
sflag = False
for i, net in enumerate(path):
if net in known_ops or i == len(path) - 1:
if sflag:
if i > 0:
yield path[:i+1]
break
else:
sflag = True
known_ops.add(net)
def get_ordered_ops(path_segments: list[list[Op]]) -> list[Op]:
finished_paths: set[int] = set()
for i, path in enumerate(path_segments):
#print(i)
if i not in finished_paths:
for op in path:
for j in range(i + 1, len(path_segments)):
path_stub = path_segments[j]
if op == path_stub[-1]:
for insert_op in path_stub[:-1]:
#print('->', insert_op)
yield insert_op
finished_paths.add(j)
#print('- ', op)
yield op
finished_paths.add(i)
def get_consts(op_list: list[Node]):
net_lookup = {net.source: net for op in op_list for net in op.args}
return [(n.name, net_lookup[n], n.value) for n in op_list if isinstance(n, Const)]
def add_read_ops(op_list):
"""Add read operation before each op where arguments are not allredy possitioned
correctly in the registers
Returns:
Yields a tuples of a net and a operation. The net is the result net
from the retuned operation"""
registers = [None] * 16
net_lookup = {net.source: net for op in op_list for net in op.args}
for op in op_list:
if not op.name.startswith('const_'):
for i, net in enumerate(op.args):
if net != registers[i]:
#if net in registers:
# print('x swap registers')
new_op = Op('read_reg' + str(i) + '_' + net.dtype, [])
yield net, new_op
registers[i] = net
yield net_lookup.get(op), op
if op in net_lookup:
registers[0] = net_lookup[op]
def add_write_ops(op_list, const_list):
"""Add write operation for each new defined net if a read operation is later folowed"""
stored_nets = {c[1] for c in const_list}
read_back_nets = {net for net, op in op_list if op.name.startswith('read_reg')}
for net, op in op_list:
yield net, op
if net in read_back_nets and net not in stored_nets:
yield (net, Op('write_' + net.dtype, [net]))
stored_nets.add(net)

98
src/copapy/binwrite.py Normal file
View File

@ -0,0 +1,98 @@
from enum import Enum
from pelfy import elf_symbol
Command = Enum('Command', [('ALLOCATE_DATA', 1), ('COPY_DATA', 2),
('ALLOCATE_CODE', 3), ('COPY_CODE', 4),
('RELOCATE_FUNC', 5), ('RELOCATE_OBJECT', 6),
('SET_ENTR_POINT', 64), ('END_PROG', 255)])
RelocationType = Enum('RelocationType', [('RELOC_RELATIVE_32', 0)])
def translate_relocation(new_sym_addr: int, new_patch_addr: int, reloc_type: str, r_addend: int) -> int:
if reloc_type in ('R_AMD64_PLT32', 'R_AMD64_PC32'):
# S + A - P
value = new_sym_addr + r_addend - new_patch_addr
return RelocationType.RELOC_RELATIVE_32, value
else:
raise Exception(f"Unknown: {reloc_type}")
def get_variable_data(symbols: list[elf_symbol]) -> tuple[list[tuple[elf_symbol, int, int]], int]:
object_list = []
out_offs = 0
for sym in symbols:
assert sym.info == 'STT_OBJECT'
lengths = sym.fields['st_size']
object_list.append((sym, out_offs, lengths))
out_offs += (lengths + 3) // 4 * 4
return object_list, out_offs
def get_function_data(symbols: list[elf_symbol]) -> tuple[list[tuple[elf_symbol, int, int, int]], int]:
code_list = []
out_offs = 0
for sym in symbols:
assert sym.info == 'STT_FUNC'
lengths = sym.fields['st_size']
#if strip_function:
# assert False, 'Not implemente'
# TODO: Strip functions
# Symbol, start out_offset in symbol, offset in output file, output lengths
# Set in_sym_out_offs and lengths
in_sym_offs = 0
code_list.append((sym, in_sym_offs, out_offs, lengths))
# out_offs += (lengths + 3) // 4 * 4
out_offs += lengths # should be aligned by default?
return code_list, out_offs
def get_function_data_blob(symbols: list[elf_symbol]) -> tuple[list[tuple[elf_symbol, int, int, int]], int]:
code_list = []
out_offs = 0
for sym in symbols:
assert sym.info == 'STT_FUNC'
lengths = sym.fields['st_size']
#if strip_function:
# assert False, 'Not implemente'
# TODO: Strip functions
# Symbol, start out_offset in symbol, offset in output file, output lengths
# Set in_sym_out_offs and lengths
in_sym_offs = 0
code_list.append((sym, in_sym_offs, out_offs, lengths))
out_offs += (lengths + 3) // 4 * 4
return code_list, out_offs
class data_writer():
def __init__(self, byteorder: str):
self._data: list[(str, bytes)] = list()
self.byteorder = byteorder
def write_int(self, value: int, num_bytes: int = 4, signed: bool = False):
self._data.append((f"INT {value}", value.to_bytes(length=num_bytes, byteorder=self.byteorder, signed=signed), 0))
def write_com(self, value: Enum, num_bytes: int = 4):
self._data.append((value.name, value.value.to_bytes(length=num_bytes, byteorder=self.byteorder, signed=False), 1))
def write_byte(self, value: int):
self._data.append((f"BYTE {value}", bytes([value]), 0))
def write_bytes(self, value: bytes):
self._data.append((f"BYTES {len(value)}", value, 0))
def print(self) -> str:
for name, dat, flag in self._data:
if flag:
print('')
print(f"{name:18}" + ' '.join(f'{b:02X}' for b in dat))
def get_data(self) -> bytes:
return b''.join(dat for _, dat, _ in self._data)
def to_file(self, path: str):
with open(path, 'wb') as f:
f.write(self.get_data())

95
src/copapy/ops.c Normal file
View File

@ -0,0 +1,95 @@
//notes:
//Helper functions
void result_int(int ret1){
asm ("");
}
void result_float(float ret1){
asm ("");
}
void result_float_float(float ret1, float ret2){
asm ("");
}
//Operations
void add_int_int(int arg1, int arg2){
result_int(arg1 + arg2);
}
void add_float_float(float arg1, float arg2){
result_float(arg1 + arg2);
}
void add_float_int(float arg1, float arg2){
result_float(arg1 + arg2);
}
void sub_int_int(int arg1, int arg2){
result_int(arg1 - arg2);
}
void sub_float_float(float arg1, float arg2){
result_float(arg1 - arg2);
}
void sub_float_int(float arg1, int arg2){
result_float(arg1 - arg2);
}
void sub_int_float(int arg1, float arg2){
result_float(arg1 - arg2);
}
void mul_int_int(int arg1, int arg2){
result_int(arg1 * arg2);
}
void mul_float_float(float arg1, float arg2){
result_float(arg1 * arg2);
}
void mul_float_int(float arg1, int arg2){
result_float(arg1 + arg2);
}
void div_int_int(int arg1, int arg2){
result_int(arg1 / arg2);
}
void div_float_float(float arg1, float arg2){
result_float(arg1 / arg2);
}
void div_float_int(float arg1, int arg2){
result_float(arg1 / arg2);
}
void div_int_float(int arg1, float arg2){
result_float(arg1 / arg2);
}
//Read global variables from heap
int read_int_ret = 1337;
void read_int(){
result_int(read_int_ret);
}
float read_float_ret = 1337;
void read_float(){
result_float(read_float_ret);
}
void read_float_2(float arg1){
result_float_float(arg1, read_float_ret);
}

59
tests/test_ast_gen.py Normal file
View File

@ -0,0 +1,59 @@
#from diss_helper import dmprint
import pelfy
#from IPython.display import Markdown
from typing import Literal
from copapy import Op, Write, const
import copapy as rc
from copapy.binwrite import data_writer, Command, RelocationType, get_variable_data, get_function_data, translate_relocation
def test_ast_generation():
c1 = const(1.11)
c2 = const(2.22)
#c3 = const(3.33)
#i1 = c1 + c2
#i2 = c2 * i1
#i3 = i2 + 4
#r1 = i1 + i3
#r2 = i3 * i2
i1 = c1 * 2
i2 = i1 + 3
r1 = i1 + i2
r2 = c2 + 4 + c1
out = [Write(r1), Write(r2)]
print(out)
print('--')
path_segments = list(rc.get_path_segments(out))
for p in path_segments:
print(p)
print('--')
ordered_ops = list(rc.get_ordered_ops(path_segments))
for p in path_segments:
print(p)
print('--')
const_list = rc.get_consts(ordered_ops)
for p in const_list:
print(p)
print('--')
output_ops = list(rc.add_read_ops(ordered_ops))
for p in output_ops:
print(p)
print('--')
extended_output_ops = list(rc.add_write_ops(output_ops, const_list))
for p in extended_output_ops:
print(p)
print('--')
if __name__ == "__main__":
test_ast_generation()