deng-mit / arrhenius.jl Goto Github PK

View Code? Open in Web Editor NEW

55.0 10.0 20.0 3.18 MB

Differentiable Reacting Flow Modeling Software

Home Page: https://deng-mit.github.io/Arrhenius.jl/dev/

License: MIT License

Julia 16.44% Python 2.79% Jupyter Notebook 80.77%

combustion sciml julialang machine-learning deeplearning auto-differentiation cantera reaction flame energy

arrhenius.jl's People

Contributors

Stargazers

Watchers

Forkers

minhbau stjordanis suxy15 tjp-karpowski rsuryanarayan deepkinetics ourchinaindata jatropj undi123 lkampoli naikless 13269252368 avinashresearch1 jiweiqi myx98 sezerh24 younghwacho 66my lsyh366

arrhenius.jl's Issues

Implement a multiplier to rate constants

Pass a vector that contains factors to create the solution, and multiply it in wdot_func

Handle falloff reaction without troe coefficients

Arrhenius.jl/src/Kinetics.jl

Line 23 in cb0b2ee

if reaction.Troe_A[j] > 1.e-12

Need better approach to handle such situations. Currently, I have to manually go to the yaml file and put zero Troe coefficients.

Implement an equilibrium solver

This is very essential for a lot of the cases

Wrong definition of X2C

The current equation for X2C(gas, X, ρ_mass)

C = X * ρ_mass ./ gas.MW

seems incorrect to me. It is the same as for Y2C, which can't be the case.
Instead it should read

C = X * ρ_mass / sum(X.*gas.MW)

Set up documentation

https://sparkbox.com/foundry/github_wiki_tutorial_for_technical_wiki_documentation

Discretization for one-dimensional flame

Goal:

implement the discretization for one-dimensional flame, then one can compute the residuals and the sensitivity for a given solution. Further, one can also solve the equations.

Equations:

see Chemkin manual Sec. 12.6

Notes:

similar to the one did for sensBVP for ignition, but more complex with the central difference.

Investigate if it is beneficial to use sparse matrix for stoichiometric coeffcients

Arrhenius.jl/src/Solution.jl

Line 56 in 51a2f3d

product_stoich_coeffs = sparse(npz["product_stoich_coeffs"])

Perhaps testing one small mechanism of GRI3.0, and a large one for n-heptane.

Divide the interpretation of kinetic model into three parts

Arrhenius.jl/src/Solution.jl

Line 9 in 51a2f3d

function CreateSolution(mech)

Currently, the interpretation of kinetic, thermodynamic, and transport are coded together. It will be great to spill them into three parts, to make the code more modulized.

Determine if we need those @inbounds

If those @inbounds are not necessary, we shall remove them to improve the readiness of the source code.

Arrhenius.jl/src/Kinetics.jl

Lines 2 to 53 in 03e31ea

 function wdot_func(reaction, T, C, S0, h_mole; get_qdot=false) 

 @inbounds _kf = @. @view(reaction.Arrhenius_coeffs[:, 1]) * exp( 

 @view(reaction.Arrhenius_coeffs[:, 2]) * log(T) - 

 @view(reaction.Arrhenius_coeffs[:, 3]) * (4184.0 / R / T), 

 ) 

 for i in reaction.index_three_body 

 @inbounds _kf[i] *= dot(@view(reaction.efficiencies_coeffs[:, i]), C) 

 end 

 for (j, i) in enumerate(reaction.index_falloff) 

 @inbounds A0, b0, Ea0 = reaction.Arrhenius_0[j, :] 

 @inbounds k0 = A0 * exp(b0 * log(T) - Ea0 * 4184.0 / R / T) 

 @inbounds Pr = 

 k0 * dot(@view(reaction.efficiencies_coeffs[:, i]), C) / _kf[i] 

 lPr = log10(Pr) 

 _kf[i] *= (Pr / (1 + Pr)) 

 # reference: 

 # http://web.mit.edu/2.62/cantera/doc/html/classCantera_1_1Troe4.html#a38aa787421d426dfd0a587fd6fc8108e 

 if reaction.index_falloff_Troe[j] > 0 

 k = reaction.index_falloff_Troe[j] 

 @inbounds F_cent = 

 (1 - reaction.Troe_[k, 1]) * exp(-T / reaction.Troe_[k, 4]) + 

 reaction.Troe_[k, 1] * exp(-T / reaction.Troe_[k, 2]) + 

 exp(-reaction.Troe_[k, 3] / T) 

 lF_cent = log10(F_cent) 

 _C = -0.4 - 0.67 * lF_cent 

 N = 0.75 - 1.27 * lF_cent 

 @inbounds f1 = (lPr + _C) / (N - 0.14 * (lPr + _C)) 

 @inbounds _kf[i] *= exp(log(10.0) * lF_cent / (1 + f1^2)) 

 end 

 end 

 @inbounds ΔS_R = reaction.vk' * S0 / R 

 @inbounds ΔH_RT = reaction.vk' * h_mole / (R * T) 

 @inbounds Keq = 

 @. exp(ΔS_R - ΔH_RT + log(one_atm / R / T) * reaction.vk_sum) 

 @inbounds _kr = @. _kf / Keq * reaction.is_reversible 

 for i = 1:reaction.n_reactions 

 @inbounds for j in reaction.i_reactant[i] 

 @inbounds _kf[i] *= C[j]^reaction.reactant_orders[j, i] 

 end 

 if reaction.is_reversible[i] 

 @inbounds for j in reaction.i_product[i] 

 @inbounds _kr[i] *= C[j]^reaction.product_stoich_coeffs[j, i] 

 end 

 end 

 end

Ideally, we could come up with a guideline on when we need explicitly to use the keyword, and when the compiler will automatically do the code optimization.

Prepare a list of resources for learning Julia

As Julia is not widely adopted in the Combustion community, it will be great to have instructions for beginners from different backgrounds. For example,

without any programming experience
knowing Fortran/Chemkin well
knowing Cantera's python/Matlab interface well

Fail to load CH4_Kazakov_without_NOx.yaml in PSR.jl

hi,

I try just to compile the PSR.jl resulting in the error below. Can you recompile the DRM22 mech (CH4_Kazakov_without_NOx.yaml) / or fix it ?

Thanks, Dmitry

julia> gas = CreateSolution("CH4_Kazakov_without_NOx.yaml")
ERROR: MethodError: no method matching Array{Float64,2}(::Array{Float64,3})
Closest candidates are:
Array{Float64,2}(::AbstractArray{S,N}) where {T, N, S} at array.jl:562
Array{Float64,2}(::UndefInitializer, ::Int64, ::Int64) where T at boot.jl:408
Array{Float64,2}(::UndefInitializer, ::Int64...) where {T, N} at boot.jl:412
...
Stacktrace:
[1] convert(::Type{Array{Float64,2}}, ::Array{Float64,3}) at .\array.jl:554
[2] Arrhenius.Transport(::Int64, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,3}) at C:\Users\lysenkod.julia\packages\Arrhenius\BOKjU\src\DataStructure.jl:23
[3] CreateSolution(::String) at C:\Users\lysenkod.julia\packages\Arrhenius\BOKjU\src\Solution.jl:142
[4] top-level scope at REPL[31]:1

Implement a PSR in the example folder

For example
https://github.com/DENG-MIT/NN-PSR/blob/main/node_pyrolysis.jl

Arrhenius.jl for Rapid Compression Machine

Code is almost ready, research can come online soon

Determine the choice of sparse matrix based on the matrix size

Need a careful benchmarking to optimize the choice of whether using sparse matrix according to the size of the matrix and its sparsity

Determine the tolerance for combustion simulations

References

ReactionMechanismGenerator/ReactionMechanismSimulator.jl#87

Interpret kinetic data from yaml directly and bypass .npz

Arrhenius.jl/src/Solution.jl

Line 56 in 51a2f3d

product_stoich_coeffs = sparse(npz["product_stoich_coeffs"])

Basically, following the code in ReacTorch. So far, it works well, and this task is not urgent at all.

Use NN to solve PSR

Use a NN to produce S-Curve

Anaconda Environment for interpreter.py

Hello,
I am trying to convert my own mechanisms to the required npz file.
With my current setup i can't manage to to this.
python interpreter.py -i h2o2.yaml

C:\Users\neuma\anaconda3\envs\reac-env\Lib\site-packages\torch\__init__.py:749: UserWarning: torch.set_default_tensor_type() is deprecated as of PyTorch 2.1, please use torch.set_default_dtype() and torch.set_default_device() as alternatives. (Triggered internally at ..\torch\csrc\tensor\python_tensor.cpp:433.)
  _C._set_default_tensor_type(t)
Processing mechanism h2o2.yaml
C:\Users\neuma\anaconda3\envs\reac-env\Lib\site-packages\reactorch\solution.py:76: UserWarning: Creating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with numpy.array() before converting to a tensor. (Triggered internally at ..\torch\csrc\utils\tensor_new.cpp:277.)
  self.molecular_weights = torch.Tensor([self.gas.molecular_weights]).T.to(self.device)
Traceback (most recent call last):
  File "C:\Users\neuma\Documents\Arrhenius\interpreter.py", line 39, in <module>
    sol = rt.Solution(mech_yaml=mech_yaml, device=torch.device('cpu'),
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\neuma\anaconda3\envs\reac-env\Lib\site-packages\reactorch\solution.py", line 88, in __init__
    self.set_reactions()
  File "C:\Users\neuma\anaconda3\envs\reac-env\Lib\site-packages\reactorch\import_kinetics.py", line 77, in set_reactions
    self.reaction[i] = {'equation': self.gas.reaction_equations(i)}
                                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "build\\python\\cantera\\kinetics.pyx", line 240, in cantera.kinetics.Kinetics.reaction_equations
TypeError: 'int' object is not iterable

My guess is that my cantera/pytorch/reactorch version is uncompatible. I would want to avoid modifying the interpreter.py script. If you could provide me with a requirement.txt to duplicate your environments, that would be great.

conda list of my current environment yields:

packages in environment at C:\Users\neuma\anaconda3\envs\reac-env:

Name                    Version                   Build  Channel
blas                      1.0                         mkl
bzip2                     1.0.8                h2bbff1b_6
ca-certificates           2024.3.11            haa95532_0
cantera                   3.0.0           py311h859585b_1    cantera
contourpy                 1.2.1                    pypi_0    pypi
cycler                    0.12.1                   pypi_0    pypi
filelock                  3.15.1                   pypi_0    pypi
fonttools                 4.53.0                   pypi_0    pypi
fsspec                    2024.6.0                 pypi_0    pypi
hdf5                      1.12.1               h51c971a_3
icc_rt                    2022.1.0             h6049295_2
intel-openmp              2021.4.0                 pypi_0    pypi
jinja2                    3.1.4                    pypi_0    pypi
kiwisolver                1.4.5                    pypi_0    pypi
libcantera                3.0.0                h82bb817_1    cantera
libffi                    3.4.4                hd77b12b_1
markupsafe                2.1.5                    pypi_0    pypi
matplotlib                3.9.0                    pypi_0    pypi
mkl                       2021.4.0                 pypi_0    pypi
mkl-service               2.4.0           py311h2bbff1b_1
mkl_fft                   1.3.8           py311h2bbff1b_0
mkl_random                1.2.4           py311h59b6b97_0
mpmath                    1.3.0                    pypi_0    pypi
networkx                  3.3                      pypi_0    pypi
numpy                     1.26.4          py311hdab7c0b_0
numpy-base                1.26.4          py311hd01c5d8_0
openssl                   3.0.14               h827c3e9_0
packaging                 24.1                     pypi_0    pypi
pillow                    10.3.0                   pypi_0    pypi
pip                       24.0            py311haa95532_0
pyparsing                 3.1.2                    pypi_0    pypi
python                    3.11.9               he1021f5_0
python-dateutil           2.9.0.post0              pypi_0    pypi
reactorch                 0.1.1                    pypi_0    pypi
ruamel.yaml               0.17.21         py311h2bbff1b_0
scipy                     1.13.1                   pypi_0    pypi
setuptools                69.5.1          py311haa95532_0
six                       1.16.0                   pypi_0    pypi
sqlite                    3.45.3               h2bbff1b_0
sympy                     1.12.1                   pypi_0    pypi
tbb                       2021.12.0                pypi_0    pypi
tk                        8.6.14               h0416ee5_0
torch                     2.3.1+cu118              pypi_0    pypi
torchaudio                2.3.1+cu118              pypi_0    pypi
torchvision               0.18.1+cu118             pypi_0    pypi
typing-extensions         4.12.2                   pypi_0    pypi
tzdata                    2024a                h04d1e81_0
vc                        14.2                 h2eaa2aa_1
vs2015_runtime            14.29.30133          h43f2093_3
wheel                     0.43.0          py311haa95532_0
xz                        5.4.6                h8cc25b3_1
zlib                      1.2.13               h8cc25b3_1

Add Unit Testing for development

As the package becomes more and more complex, it is important to adopt unit testing. We can following the example of RMS.jl at https://github.com/ReactionMechanismGenerator/ReactionMechanismSimulator.jl/blob/ac2a3bd88504b86d5bc60dedc09533839a9c9cbb/src/Calculators/TestThermo.jl#L1-L60.

Some thoughts:

We could put the test files inside the src folder, with a consistent prefix, so that one can easily spot the actual source file.
We can test the functions against Cantera via PyCall. This should work well for most of the functionality.
@test Cps ≈ Cpexplist rtol=1e-3 where 1e-3 is a good default value for tolerance.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

add elemental matrix into struct gas

something like

E_C = [10, 0, 0, 1, 1, 3, 2, 3, 5, 6, 7, 0, 0, 0, 0];
E_H = [16, 1, 2, 3, 4, 5, 4, 6, 6, 6, 8, 0, 0, 1, 2];
E_N = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0];
E_O = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2];

This is especially useful for computing elemental conversation subspace.

	function wdot_func(reaction, T, C, S0, h_mole; get_qdot=false)

	@inbounds _kf = @. @view(reaction.Arrhenius_coeffs[:, 1]) * exp(
	@view(reaction.Arrhenius_coeffs[:, 2]) * log(T) -
	@view(reaction.Arrhenius_coeffs[:, 3]) * (4184.0 / R / T),
	)

	for i in reaction.index_three_body
	@inbounds _kf[i] *= dot(@view(reaction.efficiencies_coeffs[:, i]), C)
	end

	for (j, i) in enumerate(reaction.index_falloff)
	@inbounds A0, b0, Ea0 = reaction.Arrhenius_0[j, :]
	@inbounds k0 = A0 * exp(b0 * log(T) - Ea0 * 4184.0 / R / T)
	@inbounds Pr =
	k0 * dot(@view(reaction.efficiencies_coeffs[:, i]), C) / _kf[i]
	lPr = log10(Pr)
	_kf[i] *= (Pr / (1 + Pr))

	# reference:
	# http://web.mit.edu/2.62/cantera/doc/html/classCantera_1_1Troe4.html#a38aa787421d426dfd0a587fd6fc8108e
	if reaction.index_falloff_Troe[j] > 0
	k = reaction.index_falloff_Troe[j]
	@inbounds F_cent =
	(1 - reaction.Troe_[k, 1]) * exp(-T / reaction.Troe_[k, 4]) +
	reaction.Troe_[k, 1] * exp(-T / reaction.Troe_[k, 2]) +
	exp(-reaction.Troe_[k, 3] / T)

	lF_cent = log10(F_cent)
	_C = -0.4 - 0.67 * lF_cent
	N = 0.75 - 1.27 * lF_cent
	@inbounds f1 = (lPr + _C) / (N - 0.14 * (lPr + _C))
	@inbounds _kf[i] = exp(log(10.0) lF_cent / (1 + f1^2))
	end
	end

	@inbounds ΔS_R = reaction.vk' * S0 / R
	@inbounds ΔH_RT = reaction.vk' * h_mole / (R * T)
	@inbounds Keq =
	@. exp(ΔS_R - ΔH_RT + log(one_atm / R / T) * reaction.vk_sum)
	@inbounds _kr = @. _kf / Keq * reaction.is_reversible

	for i = 1:reaction.n_reactions
	@inbounds for j in reaction.i_reactant[i]
	@inbounds _kf[i] *= C[j]^reaction.reactant_orders[j, i]
	end
	if reaction.is_reversible[i]
	@inbounds for j in reaction.i_product[i]
	@inbounds _kr[i] *= C[j]^reaction.product_stoich_coeffs[j, i]
	end
	end
	end

deng-mit / arrhenius.jl Goto Github PK

arrhenius.jl's People

Contributors

Stargazers

Watchers

Forkers

arrhenius.jl's Issues

Goal:

Equations:

Notes:

Recommend Projects

Recommend Topics

Recommend Org