Compare commits
No commits in common. "372ad7fe5a5f6c8d1e2959a4d96a7b64fabe777a" and "d6c6bcac0ffd033fb2b1dc9d79603a96e41cfc7e" have entirely different histories.
372ad7fe5a
...
d6c6bcac0f
|
|
@ -12,7 +12,7 @@ jobs:
|
||||||
|
|
||||||
strategy:
|
strategy:
|
||||||
matrix:
|
matrix:
|
||||||
python-version: ["3.10", 3.11, 3.12, 3.13, 3.14]
|
python-version: ["3.10", 3.11, 3.12, 3.13]
|
||||||
|
|
||||||
steps:
|
steps:
|
||||||
- name: Check out code
|
- name: Check out code
|
||||||
|
|
|
||||||
|
|
@ -14,5 +14,4 @@ venv/
|
||||||
.venv/
|
.venv/
|
||||||
thermo_data/combined_data.yaml
|
thermo_data/combined_data.yaml
|
||||||
src/gaspype/data/therm_data.bin
|
src/gaspype/data/therm_data.bin
|
||||||
*/*/build/
|
*/*/build/
|
||||||
examples/*.ipynb
|
|
||||||
|
|
@ -11,11 +11,11 @@ The conversion is done like the following automated by the
|
||||||
[docs/source/render_examples.py](../docs/source/render_examples.py) script:
|
[docs/source/render_examples.py](../docs/source/render_examples.py) script:
|
||||||
``` bash
|
``` bash
|
||||||
# Converting markdown with code sections to Jupyter Notebook and run it:
|
# Converting markdown with code sections to Jupyter Notebook and run it:
|
||||||
notedown examples/sofc_methane.md --to notebook --output docs/source/api/sofc_methane.ipynb --run
|
notedown examples/soec_methane.md --to notebook --output docs/source/api/soec_methane.ipynb --run
|
||||||
|
|
||||||
# Converting the Jupyter Notebook to Markdown and a folder with image
|
# Converting the Jupyter Notebook to Markdown and a folder with image
|
||||||
# files placed in docs/source/api/:
|
# files placed in docs/source/api/:
|
||||||
jupyter nbconvert --to markdown docs/source/api/sofc_methane.ipynb --output sofc_methane.md
|
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
|
A new example Markdown file can be created from a Jupyter Notebook running
|
||||||
|
|
|
||||||
|
|
@ -50,10 +50,10 @@ doc_build = [
|
||||||
"myst-parser",
|
"myst-parser",
|
||||||
"pandas",
|
"pandas",
|
||||||
"matplotlib",
|
"matplotlib",
|
||||||
"ipykernel==6.29",
|
"ipykernel",
|
||||||
"jupyter==1.1",
|
"jupyter",
|
||||||
"nbconvert==7.16",
|
"nbconvert",
|
||||||
"notedown==1.5"
|
"notedown"
|
||||||
]
|
]
|
||||||
|
|
||||||
[tool.mypy]
|
[tool.mypy]
|
||||||
|
|
|
||||||
|
|
@ -3,25 +3,22 @@ import numpy as np
|
||||||
import time
|
import time
|
||||||
import gaspype as gp
|
import gaspype as gp
|
||||||
|
|
||||||
try:
|
|
||||||
import cea
|
|
||||||
CEA_AVAILABLE = True
|
|
||||||
except ImportError:
|
|
||||||
CEA_AVAILABLE = False
|
|
||||||
|
|
||||||
gas = ct.Solution("gri30.yaml")
|
gas = ct.Solution("gri30.yaml")
|
||||||
composition = {"H2": 0.3, "H2O": 0.3, "N2": 0.4}
|
composition = {"H2": 0.3, "H2O": 0.3, "N2": 0.4}
|
||||||
|
|
||||||
n_species = gas.n_species
|
n_species = gas.n_species
|
||||||
n_states = 1_000_000
|
n_states = 1_000_000
|
||||||
|
|
||||||
|
# Random temperatures and pressures
|
||||||
temperatures = np.linspace(300.0, 2500.0, n_states)
|
temperatures = np.linspace(300.0, 2500.0, n_states)
|
||||||
pressures = np.full(n_states, ct.one_atm)
|
pressures = np.full(n_states, ct.one_atm)
|
||||||
|
|
||||||
|
# Create a SolutionArray with many states at once
|
||||||
states = ct.SolutionArray(gas, len(temperatures))
|
states = ct.SolutionArray(gas, len(temperatures))
|
||||||
|
|
||||||
time.sleep(0.5)
|
time.sleep(0.5)
|
||||||
|
|
||||||
|
# Vectorized assignment
|
||||||
t0 = time.perf_counter()
|
t0 = time.perf_counter()
|
||||||
states.TPX = temperatures, pressures, composition
|
states.TPX = temperatures, pressures, composition
|
||||||
cp_values = states.cp_mole
|
cp_values = states.cp_mole
|
||||||
|
|
@ -30,36 +27,16 @@ elapsed = time.perf_counter() - t0
|
||||||
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized cantera)")
|
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)
|
print("First 5 Cp values (J/mol-K):", cp_values[:5] / 1000)
|
||||||
|
|
||||||
|
|
||||||
|
# Vectorized fluid creation
|
||||||
fluid = gp.fluid(composition)
|
fluid = gp.fluid(composition)
|
||||||
|
|
||||||
time.sleep(0.5)
|
time.sleep(0.5)
|
||||||
|
|
||||||
|
# Benchmark: calculate Cp for all states at once
|
||||||
t0 = time.perf_counter()
|
t0 = time.perf_counter()
|
||||||
cp_values = fluid.get_cp(t=temperatures)
|
cp_values = fluid.get_cp(t=temperatures)
|
||||||
elapsed = time.perf_counter() - t0
|
elapsed = time.perf_counter() - t0
|
||||||
|
|
||||||
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
|
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
|
||||||
print("First 5 Cp values (J/mol·K):", cp_values[:5])
|
print("First 5 Cp values (J/mol·K):", cp_values[:5])
|
||||||
|
|
||||||
if CEA_AVAILABLE:
|
|
||||||
cea_mix = cea.Mixture(['H2', 'H2O', 'N2'])
|
|
||||||
mole_fracs = np.array([0.3, 0.3, 0.4])
|
|
||||||
MW = np.array([2.016, 18.015, 28.014])
|
|
||||||
mass_weights = mole_fracs * MW / (mole_fracs * MW).sum()
|
|
||||||
avg_MW = np.sum(mole_fracs * MW)
|
|
||||||
p_bar = cea.units.atm_to_bar(1.0)
|
|
||||||
|
|
||||||
time.sleep(0.5)
|
|
||||||
|
|
||||||
# the current NASA CEA Python API does not provide a NumPy-style
|
|
||||||
# vectorized interface for thermodynamic property lookups
|
|
||||||
t0 = time.perf_counter()
|
|
||||||
cea_cp = np.zeros(n_states)
|
|
||||||
for i in range(n_states):
|
|
||||||
cea_cp[i] = cea_mix.calc_property(cea.FROZEN_CP, mass_weights, temperatures[i], p_bar)
|
|
||||||
elapsed = time.perf_counter() - t0
|
|
||||||
|
|
||||||
cea_cp_molar = cea_cp * avg_MW / 1000
|
|
||||||
|
|
||||||
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (CEA)")
|
|
||||||
print("First 5 Cp values (J/mol-K):", cea_cp_molar[:5])
|
|
||||||
|
|
|
||||||
|
|
@ -3,12 +3,6 @@ import numpy as np
|
||||||
import time
|
import time
|
||||||
import gaspype as gp
|
import gaspype as gp
|
||||||
|
|
||||||
try:
|
|
||||||
import cea
|
|
||||||
CEA_AVAILABLE = True
|
|
||||||
except ImportError:
|
|
||||||
CEA_AVAILABLE = False
|
|
||||||
|
|
||||||
gas = ct.Solution("gri30.yaml")
|
gas = ct.Solution("gri30.yaml")
|
||||||
n_species = gas.n_species
|
n_species = gas.n_species
|
||||||
n_states = 1_000_000
|
n_states = 1_000_000
|
||||||
|
|
@ -59,26 +53,3 @@ elapsed = time.perf_counter() - t0
|
||||||
|
|
||||||
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
|
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (vectorized Gaspype)")
|
||||||
print("First 5 Cp values (J/mol·K):", cp_values[:5])
|
print("First 5 Cp values (J/mol·K):", cp_values[:5])
|
||||||
|
|
||||||
|
|
||||||
if CEA_AVAILABLE:
|
|
||||||
MW = np.array([2.016, 18.015, 28.014])
|
|
||||||
mass_weights = fractions * MW / (fractions * MW).sum(axis=1)[:, None]
|
|
||||||
avg_MW = np.sum(fractions * MW, axis=1)
|
|
||||||
p_bar = cea.units.atm_to_bar(1.0)
|
|
||||||
cea_mix = cea.Mixture(['H2', 'H2O', 'N2'])
|
|
||||||
|
|
||||||
time.sleep(0.5)
|
|
||||||
|
|
||||||
# the current NASA CEA Python API does not provide a NumPy-style
|
|
||||||
# vectorized interface for thermodynamic property lookups
|
|
||||||
t0 = time.perf_counter()
|
|
||||||
cea_cp = np.zeros(n_states)
|
|
||||||
for i in range(n_states):
|
|
||||||
cea_cp[i] = cea_mix.calc_property(cea.FROZEN_CP, mass_weights[i], temperatures[i], p_bar)
|
|
||||||
elapsed = time.perf_counter() - t0
|
|
||||||
|
|
||||||
cea_cp_molar = cea_cp * avg_MW / 1000
|
|
||||||
|
|
||||||
print(f"Computed {n_states} Cp values in {elapsed:.4f} seconds (CEA)")
|
|
||||||
print("First 5 Cp values (J/mol-K):", cea_cp_molar[:5])
|
|
||||||
|
|
|
||||||
|
|
@ -4,16 +4,9 @@ import numpy as np
|
||||||
import time
|
import time
|
||||||
from gaspype import fluid_system
|
from gaspype import fluid_system
|
||||||
|
|
||||||
try:
|
|
||||||
import cea
|
|
||||||
CEA_AVAILABLE = True
|
|
||||||
except ImportError:
|
|
||||||
CEA_AVAILABLE = False
|
|
||||||
|
|
||||||
# -----------------------
|
# -----------------------
|
||||||
# Settings
|
# Settings
|
||||||
# -----------------------
|
# -----------------------
|
||||||
|
|
||||||
n_temps = 1000
|
n_temps = 1000
|
||||||
temps_C = np.linspace(300, 1000, n_temps) # °C
|
temps_C = np.linspace(300, 1000, n_temps) # °C
|
||||||
temperatures = temps_C + 273.15 # K
|
temperatures = temps_C + 273.15 # K
|
||||||
|
|
@ -55,40 +48,11 @@ print(f"Gaspype: {elapsed_gaspype:.4f} s")
|
||||||
el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2)
|
el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2)
|
||||||
assert np.all(el_err < 1e-20)
|
assert np.all(el_err < 1e-20)
|
||||||
|
|
||||||
if CEA_AVAILABLE:
|
|
||||||
reac_names = ["CH4", "H2O"]
|
|
||||||
prod_names = ["CH4", "H2O", "CO", "CO2", "H2", "O2", "H", "O", "OH"]
|
|
||||||
|
|
||||||
reac = cea.Mixture(reac_names)
|
|
||||||
prod = cea.Mixture(prod_names)
|
|
||||||
|
|
||||||
solver = cea.EqSolver(prod, reactants=reac)
|
|
||||||
solution = cea.EqSolution(solver)
|
|
||||||
|
|
||||||
input_moles = np.array([8.0, 2.0])
|
|
||||||
input_weights = reac.moles_to_weights(input_moles)
|
|
||||||
|
|
||||||
p_bar = pressure * 1e-5 # Convert to bar
|
|
||||||
|
|
||||||
eq_cea = np.zeros((n_temps, len(species_to_track)))
|
|
||||||
|
|
||||||
time.sleep(0.5)
|
|
||||||
t0 = time.perf_counter()
|
|
||||||
for i, T in enumerate(temperatures):
|
|
||||||
solver.solve(solution, cea.TP, T, p_bar, input_weights)
|
|
||||||
if solution.converged:
|
|
||||||
for j, s in enumerate(species_to_track):
|
|
||||||
eq_cea[i, j] = solution.mole_fractions.get(s, 0.0)
|
|
||||||
elapsed_cea = time.perf_counter() - t0
|
|
||||||
print(f"CEA: {elapsed_cea:.4f} s")
|
|
||||||
|
|
||||||
# -----------------------
|
# -----------------------
|
||||||
# Compare first 5 results
|
# Compare first 5 results
|
||||||
# -----------------------
|
# -----------------------
|
||||||
print("\nFirst 5 equilibrium compositions:")
|
print("First 5 equilibrium compositions (mole fractions):")
|
||||||
for i in range(5):
|
for i in range(5):
|
||||||
print(f"T = {temperatures[i]:.1f} K")
|
print(f"T = {temperatures[i]:.1f} K")
|
||||||
print(" Cantera:", eq_cantera[i])
|
print(" Cantera:", eq_cantera[i])
|
||||||
print(" Gaspype :", eq_gaspype.array_composition[i])
|
print(" Gaspype :", eq_gaspype.array_composition[i])
|
||||||
if CEA_AVAILABLE:
|
|
||||||
print(" CEA :", eq_cea[i])
|
|
||||||
|
|
|
||||||
|
|
@ -1,54 +0,0 @@
|
||||||
"""
|
|
||||||
Testing the interpolation error of the lookup table for Cp, G, and H for multiple species.
|
|
||||||
"""
|
|
||||||
|
|
||||||
import gaspype as gp
|
|
||||||
import numpy as np
|
|
||||||
from gaspype.constants import R
|
|
||||||
|
|
||||||
fl = gp.fluid({'CH4': 1, 'CO2': 1, 'H2O': 1, 'SO2': 1, 'NH3': 1, 'NO2': 1, 'H2S': 1, 'C2H6': 1, 'C3H8': 1, 'H2': 1, 'O2': 1})
|
|
||||||
t_values = np.arange(200, 2000, step=10, dtype=int)
|
|
||||||
|
|
||||||
|
|
||||||
def test_cp_continuity_multi_species():
|
|
||||||
cp_actual = fl.fs.get_species_cp(t_values)
|
|
||||||
cp_neighbor_mean = (fl.fs.get_species_cp(t_values - 1) + fl.fs.get_species_cp(t_values + 1)) / 2
|
|
||||||
diff = (cp_actual - cp_neighbor_mean) / cp_actual / 2
|
|
||||||
avg_diff = np.mean(np.abs(diff))
|
|
||||||
max_diff = np.max(np.abs(diff))
|
|
||||||
print(f'Average difference of Cp: {avg_diff} J/mol/K')
|
|
||||||
print(f'Max difference of Cp: {max_diff} J/mol/K')
|
|
||||||
assert avg_diff < 1e-5
|
|
||||||
assert max_diff < 1e-5
|
|
||||||
|
|
||||||
|
|
||||||
def test_g_continuity_multi_species():
|
|
||||||
g_actual = fl.fs.get_species_g_rt(t_values) * R * np.reshape(t_values, (-1, 1))
|
|
||||||
g_neighbor_mean = (fl.fs.get_species_g_rt(t_values - 1) + fl.fs.get_species_g_rt(t_values + 1)) / 2 * R * np.reshape(t_values, (-1, 1))
|
|
||||||
|
|
||||||
diff = (g_actual - g_neighbor_mean) / g_actual / 2
|
|
||||||
avg_diff = np.mean(np.abs(diff))
|
|
||||||
max_diff = np.max(np.abs(diff))
|
|
||||||
print(f'Average difference of G: {avg_diff} J/mol')
|
|
||||||
print(f'Max difference of G: {max_diff} J/mol')
|
|
||||||
assert avg_diff < 1e-4
|
|
||||||
assert max_diff < 1e-4
|
|
||||||
|
|
||||||
|
|
||||||
def test_h_continuity_multi_species():
|
|
||||||
h_actual = fl.fs.get_species_h(t_values)
|
|
||||||
h_neighbor_mean = (fl.fs.get_species_h(t_values - 1) + fl.fs.get_species_h(t_values + 1)) / 2
|
|
||||||
|
|
||||||
rel_diff = (h_actual - h_neighbor_mean) / h_actual / 2
|
|
||||||
avg_diff = np.mean(np.abs(rel_diff))
|
|
||||||
max_diff = np.max(np.abs(rel_diff))
|
|
||||||
print(f'Average difference of H: {avg_diff} J/mol')
|
|
||||||
print(f'Max difference of H: {max_diff} J/mol')
|
|
||||||
assert avg_diff < 1e-4
|
|
||||||
assert max_diff < 1e-4
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
|
||||||
test_cp_continuity_multi_species()
|
|
||||||
test_g_continuity_multi_species()
|
|
||||||
test_h_continuity_multi_species()
|
|
||||||
Loading…
Reference in New Issue