Added YAML support

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.3"
 manifest_format = "2.0"
-project_hash = "3cb89865eaa8c104f3136e377a234921cf7a5336"
+project_hash = "ecf22361c4d51aae12eedd45cc504aeab2544978"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra"]
@@ -261,6 +261,12 @@ version = "1.10.2+0"
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[[deps.Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.16.1+2"
+
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "libblastrampoline_jll"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -429,6 +435,12 @@ version = "1.4.0"
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[[deps.StringEncodings]]
+deps = ["Libiconv_jll"]
+git-tree-sha1 = "33c0da881af3248dafefb939a21694b97cfece76"
+uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
+version = "0.3.6"
+
 [[deps.StructArrays]]
 deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
 git-tree-sha1 = "b03a3b745aa49b566f128977a7dd1be8711c5e71"
@@ -486,6 +498,12 @@ git-tree-sha1 = "d670a70dd3cdbe1c1186f2f17c9a68a7ec24838c"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
 version = "1.12.2"
 
+[[deps.YAML]]
+deps = ["Base64", "Dates", "Printf", "StringEncodings"]
+git-tree-sha1 = "dbc7f1c0012a69486af79c8bcdb31be820670ba2"
+uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+version = "0.4.8"
+
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

Project.toml

@@ -1,12 +1,14 @@
 name = "Apophis"
 uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054"
-authors = ["moataz-sabry <[email protected]>"]
-version = "1.1.0"
+authors = ["moataz-sabry <[email protected]>", "Ahmed Hassan <[email protected]>"]
+version = "1.2.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

README.md

@@ -1,5 +1,5 @@
 # Apophis
-A package for the sensitivity analysis of gas-phase chemical and kinetics
+A package for sensitivity analysis of gas chemical and kinetics
 
 ## Installation
 ```julia
@@ -44,6 +44,12 @@ julia> production_rates(gas) # ω̇ [mol/(cm³⋅s)]
   0.0
 
 ```
+## To-Do
+- Change default units to SI
+- Add routines for other reaction auxillary parameters
+- Correctly implement pressure-dependant reactions
+- Implement diffusion routines
+
 ## References
 - [Chemkin User Manual](https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf)
 - [Chemkin Transport Manual](https://www3.nd.edu/~powers/ame.60636/transport.pdf)

src/Apophis.jl

@@ -1,11 +1,13 @@
 module Apophis
 
-using Base.Iterators: filter, flatten, partition, take
-using Base: OneTo, rest
+using Base.Iterators: filter, flatten, reverse, take, partition
+using Base: Fix1, Fix2, OneTo, rest
 using Combinatorics
+using LinearAlgebra
 using SparseArrays
 using SplitApplyCombine: combinedimsview, mapview
 using Unitful
+using YAML: load_file
 using Zygote
 # using ThreadsX
 

src/calc.jl

@@ -14,4 +14,4 @@ update!(gas::Gas{<:Number}, d::Symbol) = ((update_thermodynamics, update_reactio
 update!(gas::Gas{<:Number}, ds::Vararg{Symbol}) = foreach(d -> update!(gas, d), ds)
 
 _kinetics_sensitivity((; mechanism, state)::Gas{<:Number}, d::Int = 1) = @inbounds map(reaction -> _kinetics_sensitivity(reaction, state)[d], mechanism.reactions)
-kinetics_sensitivity(gas::Gas{N}, d::Int = 1) where {N<:Number} = @inbounds _kinetics_sensitivity(gas, d) |> dqdk -> map(s -> sum(dqdk[i] * ν for ((; i), ν) in s.inreactions; init=zero(N)), species(gas)) |> Diagonal
+kinetics_sensitivity(gas::Gas{N}, d::Int = 1) where {N<:Number} = stoichiometry_matrix(gas) * Diagonal(_kinetics_sensitivity(gas, d))

src/reactions.jl

@@ -80,14 +80,15 @@ end
 Troe() = nothing
 Troe(itr) = Troe(itr...)
 
+Troe(a::N, T₃::N, T₁::N, ::Nothing) where {N<:Number} = Troe(a, T₃, T₁, zero(N))
 Troe(a::N, T₃::N, T₁::N) where {N<:Number} = Troe(a, T₃, T₁, zero(N))
 
 @inline ((; a, T₃, T₁, T₂)::Troe{N})(T::N) where {N<:Number} = (1.0 - a) * exp(-T / T₃) + a * exp(-T / T₁) + (iszero(T₂) ? zero(N) : exp(-T₂ / T))
 @inline ((; a, T₃, T₁, T₂)::Troe{N})(::Val{:dT}, T::N) where {N<:Number} = (a - one(N)) * exp(-T / T₃) / T₃ - a * exp(-T / T₁) / T₁ + (iszero(T₂) ? zero(N) : exp(-T₂ / T) * T₂ / T^2)
 
 (troe::Troe{N})(::Val{:dg}, T::N) where {N<:Number} = gradient(troe -> troe(T), troe) |> only
 
-function troe_function(Fc::N, Pᵣ::N) where N<:Number
+function troe_function(Fc::N, Pᵣ::N) where {N<:Number}
     log10Pᵣ = log10(Pᵣ)
     log10Fc = log10(Fc)
 
@@ -99,7 +100,7 @@ function troe_function(Fc::N, Pᵣ::N) where N<:Number
     return F
 end
 
-function troe_function(::Val{:dT}, Fc::N, Pᵣ::N) where N<:Number
+function troe_function(::Val{:dT}, Fc::N, Pᵣ::N) where {N<:Number}
     log10Pᵣ = log10(Pᵣ)
     dlog10PᵣdPᵣ = log(10.0)Pᵣ |> inv
     log10Fc = log10(Fc)
@@ -131,7 +132,7 @@ function troe_function(::Val{:dT}, Fc::N, Pᵣ::N) where N<:Number
     return F, dFdFc, dFdPᵣ
 end
 
-function troe_function(::Val{:dC}, Fc::N, Pᵣ::N) where N<:Number
+function troe_function(::Val{:dC}, Fc::N, Pᵣ::N) where {N<:Number}
     log10Pᵣ = log10(Pᵣ)
     dlog10PᵣdPᵣ = log(10.0)Pᵣ |> inv
     log10Fc = log10(Fc)
@@ -162,7 +163,7 @@ struct ElementaryReaction{N<:Number} <: AbstractReaction{N}
     reactants::Vector{Pair{Species{N}, N}}
     products::Vector{Pair{Species{N}, N}}
     reaction_order::N
-    forward_rate_parameters::Arrhenius{N}
+    forward_rate_parameters::Maybe{Arrhenius{N}}
     reverse_rate_parameters::Maybe{Arrhenius{N}}
     plog_parameters::Maybe{Plog{N}}
     rates::ReactionRates{N}
@@ -198,7 +199,7 @@ end
 
 const Reaction{N} = Union{ElementaryReaction{N}, ThreeBodyReaction{N}, FallOffReaction{N}} where {N<:Number}
 
-@inline total_molar_concentration(C::Vector{N}, α::Vector{Pair{Species{N},N}}) where N<:Number = @inbounds sum(C[k] * f for ((; k), f) in α)
+@inline total_molar_concentration(C::Vector{N}, α::Vector{Pair{Species{N},N}}) where {N<:Number} = @inbounds sum(C[k] * f for ((; k), f) in α)
 
 forward_rate((; rates, forward_rate_parameters, plog_parameters)::ElementaryReaction{N}, (; T, P)::State{N}) where {N<:Number} = setindex!(rates.kf.val, isnothing(plog_parameters) ? forward_rate_parameters(T) : plog_parameters(T, P), 1)
 forward_rate((; rates, forward_rate_parameters)::ThreeBodyReaction{N}, (; T)::State{N}) where {N<:Number} = setindex!(rates.kf.val, forward_rate_parameters(T), 1)
@@ -284,6 +285,35 @@ function forward_rate(v::Val{:dP}, (; rates, forward_rate_parameters, plog_param
     return nothing
 end
 
+_dkfdA((; forward_rate_parameters)::Union{ElementaryReaction{N}, ThreeBodyReaction{N}}, (; T)::State{N}) where {N<:Real} = forward_rate_parameters(Val(:dg), T) |> first
+
+function _dkfdA((; high_pressure_parameters, low_pressure_parameters, enhancement_factors, troe_parameters)::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Real} ## Add dkfdAₒ later
+    k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T)
+    M = total_molar_concentration(C, enhancement_factors)
+    Pᵣ = reduced_pressure(kₒ, M, k∞)
+    t = inv(one(N) + Pᵣ)
+
+    dk∞dA∞ = high_pressure_parameters(Val(:dg), T) |> first
+    dPᵣdk∞ = reduced_pressure(Val(:dk∞), kₒ, M, k∞)
+
+    if isnothing(troe_parameters)
+        dkfdk∞ = Pᵣ * t
+        dkfdPᵣ = k∞ * t^2
+        dkfdA∞ = dkfdk∞ * dk∞dA∞ + dkfdPᵣ * dPᵣdk∞ * dk∞dA∞
+    else
+        Fc = troe_parameters(T)
+        F, dFdPᵣ = troe_function(Val(:dC), Fc, Pᵣ)
+
+        dkfdF = k∞ * Pᵣ * t 
+        dkfdk∞ = F * Pᵣ * t
+        dkfdPᵣ = F * k∞ * t^2
+
+        dFdA = dFdPᵣ * (dPᵣdk∞ * dk∞dA∞)
+        dkfdA∞ = dkfdF * dFdA + dkfdk∞ * dk∞dA∞ + dkfdPᵣ * dPᵣdk∞ * dk∞dA∞
+    end
+    return dkfdA∞
+end
+
 @inline change_enthalpy((; reactants, products)::Reaction{<:Number}, i::Int = 1) = 
     @inbounds sum(getfield(s.thermo.h, i)[] * ν for (s, ν) in flatten((reactants, products)))
 
@@ -325,7 +355,7 @@ end
 @inline function reverse_rate(reaction::Reaction{N}, (; T)::State{N}) where {N<:Number}
     (; rates, isreversible, reverse_rate_parameters) = reaction
     isreversible || return nothing
-    isnothing(reverse_rate_parameters) || (kr.val[] = reverse_rate_parameters(T);
+    isnothing(reverse_rate_parameters) || (rates.kr.val[] = reverse_rate_parameters(T);
         return nothing
     )
     Kc = equilibrium_constants(reaction, T)
@@ -357,7 +387,7 @@ reverse_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, state::State{N}) where
 @inline function reverse_rate(::Val{:dP}, reaction::Reaction{N}, (; T)::State{N}) where {N<:Number}
     (; rates, isreversible, reverse_rate_parameters) = reaction
     isreversible || return nothing
-    isnothing(reverse_rate_parameters) || (kr.val[] = reverse_rate_parameters(T);
+    isnothing(reverse_rate_parameters) || (rates.kr.val[] = reverse_rate_parameters(T);
         return nothing
     )
     Kc = equilibrium_constants(reaction, T)
@@ -394,7 +424,7 @@ function progress_rate(reaction::Reaction{N}, (; C)::State{N}) where {N<:Number}
     ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
 
     M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
-    q.val[] = M * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
+    @inbounds q.val[] = M * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
     return nothing
 end
 
@@ -404,17 +434,21 @@ function progress_rate(::Val{:dT}, reaction::Reaction{N}, (; C)::State{N}) where
     ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
 
     M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
-    q.dT[] = M * (kf.dT[] * ∏ᴵ - kr.dT[] * ∏ᴵᴵ)
+    @inbounds q.val[] = M * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
+    @inbounds q.dT[] = M * (kf.dT[] * ∏ᴵ - kr.dT[] * ∏ᴵᴵ)
     return nothing
 end
 
 function progress_rate(::Val{:dC}, reaction::ElementaryReaction{N}, (; C)::State{N}) where {N<:Number}
     (; kf, kr, q) = reaction.rates
+    ∏ᴵ = step(reaction.reactants, C)
+    ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
     for k in eachindex(q.dC)
         d∏ᴵdCₖ = step(reaction.reactants, C, k)
         d∏ᴵᴵdCₖ = step(reaction.products, C, k)
         @inbounds q.dC[k] = (kf.val[] * d∏ᴵdCₖ) - (kr.val[] * d∏ᴵᴵdCₖ)
     end
+    @inbounds q.val[] = kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ
     return reaction
 end
 
@@ -432,6 +466,7 @@ function progress_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, (; C)::State{
     for ((; k), dMdCₖ) in reaction.enhancement_factors
         @inbounds q.dC[k] += dMdCₖ * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
     end
+    @inbounds q.val[] = M * (kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ)
     return reaction
 end
 
@@ -445,6 +480,7 @@ function progress_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; C)::State{N}
         d∏ᴵᴵdCₖ = step(reaction.products, C, k) ## get it out of the loop somehow in case of not reversible
         @inbounds q.dC[k] = (kf.val[] * d∏ᴵdCₖ - kr.val[] * d∏ᴵᴵdCₖ) + (kf.dC[k] * ∏ᴵ - kr.dC[k] * ∏ᴵᴵ)
     end
+    @inbounds q.val[] = kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ
     return reaction
 end
 
@@ -454,7 +490,7 @@ function progress_rate(::Val{:dP}, reaction::Reaction{N}, (; C)::State{N}) where
     ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
 
     M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
-    q.dP[] = M * (kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ)
+    @inbounds q.dP[] = M * (kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ)
     return nothing
 end