Small adjustments

Project.toml

@@ -1,7 +1,7 @@
 name = "Apophis"
 uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054"
 authors = ["moataz-sabry <[email protected]>", "Ahmed Hassan <[email protected]>"]
-version = "1.3.5"
+version = "1.3.6"
 
 [deps]
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

README.md

@@ -9,8 +9,9 @@ pkg> add Apophis
 ### Reading a Mechanism
 ```julia
 julia> using Apophis
+julia> using Unitful: K, Pa
 
-julia> gas = Gas(:GRI3; T = 1000, P = 101325, Y = "CH4: 0.05, O2: 0.20, N2: 0.75")
+julia> gas = Gas(:GRI3; T = 1000K, P = 101325Pa, Y = "CH4: 0.05, O2: 0.20, N2: 0.75")
 # Creates a Gas based on the given mechanism data files
 
               T: 1000 K  //  P: 101325 Pa  //  ρ: 0.3372 kg/m³
@@ -42,34 +43,33 @@ julia> pressure(gas) # P [Pa]
 julia> update(gas);
 # Updates internal variables based on the current gas state
 
-julia> production_rates(gas) # ω̇ [kmol/(m³⋅s)]
+julia> production_rates(gas; view=false) # ω̇ [kmol/(m³⋅s)]
 53-element Vector{Float64}:
   0.0
-  3.7072899145013144e-11
-  1.3302904784766728e-16
- -8.845602918023088e-8
+  8.034465079230524e-7
+  1.6467013518389382e-12
+ -2.8148788209137565e-5
   0.0
   0.0
-  8.845602904821038e-8
+  2.8148786697190487e-5
   0.0
   0.0
   0.0
   ⋮
   0.0
   0.0
   0.0
- -1.434894871257002e-16
+ -2.6060314085092194e-12
   0.0
   0.0
   0.0
   0.0
   0.0
-
 ```
 ## To-Do
 - Add routines for other reaction auxillary parameters
 - 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)
+- [Chemkin Transport Manual](https://www3.nd.edu/~powers/ame.60636/transport.pdf)

src/Apophis.jl

@@ -4,7 +4,7 @@ using Base.Iterators: filter, flatten, partition, reverse, take
 using Base: Fix1, Fix2, OneTo, rest
 using PrettyTables
 using SparseArrays
-using SplitApplyCombine: combinedimsview, mapview
+using SplitApplyCombine: combinedims, combinedimsview, mapview
 using Unitful
 using YAML: load_file
 using Zygote

src/calc.jl

@@ -6,7 +6,7 @@ update_thermodynamics(v::Val, (; mechanism, state)::Gas{<:Number}) = foreach(spe
 update_reaction_rates((; mechanism, state)::Gas{<:Number}) = foreach(reaction -> _update_reaction_rates(reaction, state), mechanism.reactions)
 update_reaction_rates(v::Val, (; mechanism, state)::Gas{<:Number}) = foreach(reaction -> _update_reaction_rates(v, reaction, state), mechanism.reactions)
 
-update_production_rates((; mechanism)::Gas{<:Number}) = foreach(species -> _update_production_rates(species), mechanism.species)
+update_production_rates((; mechanism)::Gas{<:Number}) = foreach(_update_production_rates, mechanism.species)
 update_production_rates(v::Val, (; mechanism)::Gas{<:Number}) = foreach(species -> _update_production_rates(v, species), mechanism.species)
 
 update(gas::Gas{<:Number}) = ((update_thermodynamics, update_reaction_rates, update_production_rates)(gas); return gas)

src/reactions.jl

@@ -162,7 +162,7 @@ struct ElementaryReaction{N<:Number, S<:AbstractSpecies{N}} <: AbstractReaction{
     reactants::Vector{Pair{S, N}}
     products::Vector{Pair{S, N}}
     reaction_order::N
-    forward_rate_parameters::Maybe{Arrhenius{N}}
+    forward_rate_parameters::Arrhenius{N}
     reverse_rate_parameters::Maybe{Arrhenius{N}}
     plog_parameters::Maybe{Plog{N}}
     rates::ReactionRates{N}
@@ -198,9 +198,9 @@ end
 
 const Reaction{N<:Number, S<:AbstractSpecies{N}} = Union{ElementaryReaction{N, S}, ThreeBodyReaction{N, S}, FallOffReaction{N, S}}
 
-total_molar_concentration(C::Vector{N}, α::Vector{Pair{Species{N},N}}) where {N<:Number} = @inbounds sum(C[k] * f for ((; k), f) in α)
+total_molar_concentration(C::Vector{N}, α::Vector{Pair{Species{N}, N}}) where {N<:Number} = @inbounds sum(C[k] * f for ((; k), f) in α; init=zero(N))
 
-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, plog_parameters)::ElementaryReaction{N}, (; T, P)::State{N}) where {N<:Number} = setindex!(rates.kf.val, forward_rate_parameters(T), 1) # 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)
 
 function forward_rate(reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Number}
@@ -217,7 +217,7 @@ function forward_rate(reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N
 end
 
 function forward_rate(v::Val{:dT}, (; rates, forward_rate_parameters, plog_parameters)::ElementaryReaction{N}, (; T, P)::State{N}) where {N<:Number}
-    @inbounds rates.kf.val[], rates.kf.dT[] = isnothing(plog_parameters) ? forward_rate_parameters(v, T) : plog_parameters(v, T, P)
+    @inbounds rates.kf.val[], rates.kf.dT[] = forward_rate_parameters(v, T) # isnothing(plog_parameters) ? forward_rate_parameters(v, T) : plog_parameters(v, T, P)
     return nothing
 end
 
@@ -282,16 +282,16 @@ function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State
     return nothing
 end
 
-function forward_rate(v::Val{:dP}, (; rates, forward_rate_parameters, plog_parameters)::ElementaryReaction{N}, (; T, P)::State{N}) where {N<:Number}
-    rates.kf.val[], rates.kf.dT[] = isnothing(plog_parameters) ? (forward_rate_parameters(T), zero(N)) : plog_parameters(v, T, P)
-    return nothing
-end
+# function forward_rate(v::Val{:dP}, (; rates, forward_rate_parameters, plog_parameters)::ElementaryReaction{N}, (; T, P)::State{N}) where {N<:Number}
+#     rates.kf.val[], rates.kf.dT[] = isnothing(plog_parameters) ? (forward_rate_parameters(T), zero(N)) : plog_parameters(v, T, P)
+#     return nothing
+# end
 
-change_enthalpy((; reactants, products)::AbstractReaction{<:Number}, i::Int=1) =
-    @inbounds sum(getfield(s.thermo.h, i)[] * ν for (s, ν) in flatten((reactants, products)))
+change_enthalpy((; reactants, products)::AbstractReaction{N}, i::Int=1) where {N<:Number} =
+    @inbounds sum(getfield(s.thermo.h, i)[] * ν for (s, ν) in flatten((reactants, products)); init=zero(N))
 
-change_entropy((; reactants, products)::AbstractReaction{<:Number}, i::Int=1) =
-    @inbounds sum(getfield(s.thermo.s, i)[] * ν for (s, ν) in flatten((reactants, products)))
+change_entropy((; reactants, products)::AbstractReaction{N}, i::Int=1) where {N<:Number} =
+    @inbounds sum(getfield(s.thermo.s, i)[] * ν for (s, ν) in flatten((reactants, products)); init=zero(N))
 
 function equilibrium_constants(reaction::AbstractReaction{N}, T::N) where {N<:Number}
     ∑ν = reaction.reaction_order
@@ -359,17 +359,17 @@ end
 reverse_rate(::Val{:dC}, reaction::ElementaryReaction{N}, state::State{N}) where {N<:Number} = reverse_rate(reaction, state)
 reverse_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, state::State{N}) where {N<:Number} = reverse_rate(reaction, state)
 
-function reverse_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; T)::State{N}) where {N<:Number}
-    (; rates, isreversible, reverse_rate_parameters) = reaction
-    isreversible || return nothing
-    isnothing(reverse_rate_parameters) || (@inbounds rates.kr.val[] = reverse_rate_parameters(T);
-        return nothing
-    )
+# function reverse_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; T)::State{N}) where {N<:Number}
+#     (; rates, isreversible, reverse_rate_parameters) = reaction
+#     isreversible || return nothing
+#     isnothing(reverse_rate_parameters) || (@inbounds rates.kr.val[] = reverse_rate_parameters(T);
+#         return nothing
+#     )
 
-    Kc = equilibrium_constants(reaction, T)
-    @inbounds rates.kr.dP[] = rates.kf.dP[] / Kc
-    return nothing
-end
+#     Kc = equilibrium_constants(reaction, T)
+#     @inbounds rates.kr.dP[] = rates.kf.dP[] / Kc
+#     return nothing
+# end
 
 function reverse_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; T)::State{N}) where {N<:Number}
     (; rates, isreversible, reverse_rate_parameters) = reaction
@@ -392,8 +392,8 @@ function reverse_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; T)::State{N})
     return nothing
 end
 
-step(part::Vector{Pair{Species{N},N}}, C::Vector{N}) where {N<:Number} = @inbounds prod(C[k]^abs(ν) for ((; k), ν) in part) ## C in small values like in units of mol/cm^3 may lead to slightly different results due to this operation
-step(part::Vector{Pair{Species{N},N}}, C::Vector{N}, j::Int) where {N<:Number} = @inbounds in(j, k for ((; k), ν) in part) ? prod(k ≠ j ? C[k]^abs(ν) : abs(ν) * C[k]^(abs(ν) - 1) for ((; k), ν) in part) : zero(N)
+step(part::Vector{Pair{Species{N},N}}, C::Vector{N}) where {N<:Number} = @inbounds prod(C[k]^abs(ν) for ((; k), ν) in part; init=one(N)) ## init=one(N) as there is no reaction without reactants or products => "part" is never empty, has at least one species
+step(part::Vector{Pair{Species{N},N}}, C::Vector{N}, j::Int) where {N<:Number} = @inbounds in(j, k for ((; k), ν) in part) ? prod(k ≠ j ? C[k]^abs(ν) : abs(ν) * C[k]^(abs(ν) - 1) for ((; k), ν) in part; init=one(N)) : zero(N)
 
 function progress_rate(reaction::AbstractReaction{N}, (; C)::State{N}) where {N<:Number}
     (; kf, kr, q) = reaction.rates
@@ -462,15 +462,15 @@ function progress_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; C)::State{N}
     return reaction
 end
 
-function progress_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; C)::State{N}) where {N<:Number}
-    (; kf, kr, q) = reaction.rates
-    ∏ᴵ = step(reaction.reactants, C)
-    ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
+# function progress_rate(::Val{:dP}, reaction::AbstractReaction{N}, (; C)::State{N}) where {N<:Number}
+#     (; kf, kr, q) = reaction.rates
+#     ∏ᴵ = step(reaction.reactants, C)
+#     ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
 
-    M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
-    @inbounds q.dP[] = M * (kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ)
-    return nothing
-end
+#     M = reaction isa ThreeBodyReaction ? total_molar_concentration(C, reaction.enhancement_factors) : one(N)
+#     @inbounds q.dP[] = M * (kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ)
+#     return nothing
+# end
 
 _update_reaction_rates(reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = (forward_rate, reverse_rate, progress_rate)(reaction, state)
 _update_reaction_rates(v::Val, reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = (forward_rate, reverse_rate, progress_rate)(v, reaction, state)