Multiple bug fixes and adjustments

Project.toml

@@ -1,23 +1,18 @@
 name = "Apophis"
 uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054"
 authors = ["moataz-sabry <[email protected]>", "Ahmed Hassan <[email protected]>"]
-version = "1.3.3"
+version = "1.3.4"
 
 [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"
 
 [compat]
-Combinatorics = "1.0.2"
 PyCall = "1.95.1"
 SplitApplyCombine = "1.2.2"
-ThreadsX = "0.1.11"
 Unitful = "1.12.3"
 YAML = "0.4.8"
 Zygote = "0.6.55"

src/Apophis.jl

@@ -2,14 +2,11 @@ module Apophis
 
 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
 
 abstract type AbstractSpecies{N<:Number} end
 abstract type AbstractReaction{N<:Number} end
@@ -22,7 +19,6 @@ const Tᵣ = 300.0 # K
 const d = 0.14
 
 const Maybe{T} = Union{T, Nothing}
-const Diffusions = (:binary_diffusions, :mean_diffusions)
 const Energies{N<:Number} = NamedTuple{(:val, :dT), Tuple{Vector{N}, Vector{N}}}
 const Thermodynamics{N<:Number} = NamedTuple{(:cₚ, :h, :s), NTuple{3, Energies{N}}}
 const Rates{N<:Number} = NamedTuple{(:val, :dT, :dC), NTuple{3, Vector{N}}}

src/gas.jl

@@ -51,15 +51,15 @@ end
 """
 Main function to create a `Gas` object.
 """
-function Gas(name::Union{Nothing, Symbol,String} = nothing; 
-    kinetics_path = "", 
-    thermo_path = "", 
-    transport_path = "",
-    as::Type{<:Number} = Float64, init...)
-
-    mechanism = read_mechanism(name; kinetics_path = kinetics_path, thermo_path = thermo_path, transport_path = transport_path, as = as)
+function Gas(name::Union{Nothing, Symbol, String} = nothing;
+    kinetics_path::Maybe{String} = nothing,
+    thermo_path::Maybe{String} = nothing,
+    transport_path::Maybe{String} = nothing,
+    as::Type{<:Number} = Float64,
+    init...)
+
+    mechanism = read_mechanism(name; kinetics_path, thermo_path, transport_path, as)
     gas = Gas(mechanism, State(mechanism))
-
     set!(gas; init...)
     return gas
 end

src/reactions.jl

@@ -208,8 +208,9 @@ function forward_rate(reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N
 
     k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T)
     M = total_molar_concentration(C, enhancement_factors)
-    Pᵣ = reduced_pressure(kₒ, M, k∞)
+    iszero(M) && (@inbounds rates.kf.val[] = zero(N); return nothing)
 
+    Pᵣ = reduced_pressure(kₒ, M, k∞)
     F = isnothing(troe_parameters) ? one(N) : troe_parameters(T) |> Fix2(troe_function, Pᵣ)
     @inbounds rates.kf.val[] = F * k∞ * Pᵣ * inv(one(N) + Pᵣ)
     return nothing
@@ -230,6 +231,7 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State
 
     (k∞, dk∞dT), (kₒ, dkₒdT) = (high_pressure_parameters, low_pressure_parameters)(v, T)
     M = total_molar_concentration(C, enhancement_factors)
+    iszero(M) && (@inbounds rates.kf.val[] = rates.kf.dT[] = zero(N); return nothing)
 
     Pᵣ = reduced_pressure(kₒ, M, k∞)
     t = inv(one(N) + Pᵣ)
@@ -250,7 +252,7 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State
         dFdT = dFdFc * dFcdT + dFdPᵣ * dPᵣdT
 
         @inbounds rates.kf.val[] = F * k∞ * Pᵣ * t
-        @inbounds rates.kf.dT[] = t * ((dk∞dT * Pᵣ + dPᵣdT * k∞ * t)F + dFdT * k∞ * Pᵣ)
+        @inbounds rates.kf.dT[] = t * (F * (dk∞dT * Pᵣ + dPᵣdT * k∞ * t) + dFdT * k∞ * Pᵣ)
     end
     return nothing
 end
@@ -263,6 +265,7 @@ function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State
 
     k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T)
     M = total_molar_concentration(C, enhancement_factors)
+    iszero(M) && (@inbounds rates.kf.dC .= rates.kf.val[] = zero(N); return nothing)
 
     Pᵣ = reduced_pressure(kₒ, M, k∞)
     t = inv(one(N) + Pᵣ)
@@ -272,8 +275,8 @@ function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State
 
     @inbounds rates.kf.val[] = F * k∞ * Pᵣ * t
     dkfdPᵣ = F * k∞ * t^2
-    dkfdF = k∞ * Pᵣ * t 
-    for ((; k), dMdCₖ) in enhancement_factors ## SparseArrays?
+    dkfdF = k∞ * Pᵣ * t
+    for ((; k), dMdCₖ) in enhancement_factors
         @inbounds rates.kf.dC[k] = (dkfdPᵣ + dkfdF * dFdPᵣ) * dPᵣdM * dMdCₖ
     end
     return nothing
@@ -424,9 +427,9 @@ end
 
 function progress_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, (; C)::State{N}) where {N<:Number}
     (; kf, kr, q) = reaction.rates
+    M = total_molar_concentration(C, reaction.enhancement_factors)
     ∏ᴵ = step(reaction.reactants, C)
     ∏ᴵᴵ = reaction.isreversible ? step(reaction.products, C) : zero(N)
-    M = total_molar_concentration(C, reaction.enhancement_factors)
 
     for k in eachindex(q.dC)
         d∏ᴵdCₖ = step(reaction.reactants, C, k)

src/reader.jl

@@ -41,7 +41,7 @@ const chemkin_default_units = Dict{String, String}(
     "quantity" => "mol"
 )
 
-without_spaces(text::AbstractString) = replace(text, r"\s+" => "")
+remove_spaces(text::AbstractString) = replace(text, r"\s+" => "")
 #=
 Reads a file data located at `file_path` as String. Filters out any exclamation marks `!` and the text following it.
 =#
@@ -66,7 +66,7 @@ find_composition(composition::Dict, N::Type{<:Number}) = Pair{Symbol, N}[element
 #=
 Finds molecular weight of the given species.
 =#
-find_weight(composition::Vector{Pair{Symbol, N}}) where {N<:Number} = sum(weights[element] * count for (element, count) in composition)
+find_weight(composition::Vector{Pair{Symbol, N}}) where {N<:Number} = sum(N(weights[element]) * count for (element, count) in composition)
 
 #=
 Finds temperature ranges of the given species.
@@ -104,7 +104,7 @@ end
 function _find_species(k::Int, species::Symbol, thermo::T, transport::Maybe{T}, N::Type{<:Number}) where {T<:String}
     transport_parameters = isnothing(transport) ? nothing : find_transport_parameters(species, transport, N)
     thermodynamics_parameters = find_thermodynamic_parameters(species, thermo, N)
-    return Species(k, species, Pair{Reaction{N, Species{N}}, N}[], thermodynamics_parameters..., transport_parameters, (; zip(Diffusions, (zeros(N, 9), zero(N)))...),
+    return Species(k, species, Pair{Reaction{N, Species{N}}, N}[], thermodynamics_parameters..., transport_parameters,
         (; zip((:cₚ, :h, :s), Tuple((; zip((:val, :dT), (zeros(N, 1) for _ in OneTo(2)))...) for _ in OneTo(3)))...), 
         (; zip((:ω̇,), tuple((; zip((:val, :dT, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...))
 end
@@ -115,7 +115,7 @@ function _find_species(k::Int, species::Dict, N::Type{<:Number})
 
     weight = find_weight(composition)
     thermodynamic_parameters = find_thermodynamic_parameters(species["thermo"], N)
-    return Species(k, name, Pair{Reaction{N, Species{N}}, N}[], composition, weight, thermodynamic_parameters, nothing, (; zip(Diffusions, (zeros(N, 9), zero(N)))...),
+    return Species(k, name, Pair{Reaction{N, Species{N}}, N}[], composition, weight, thermodynamic_parameters, nothing,
         (; zip((:cₚ, :h, :s), Tuple((; zip((:val, :dT), (zeros(N, 1) for _ in OneTo(2)))...) for _ in OneTo(3)))...), 
         (; zip((:ω̇,), tuple((; zip((:val, :dT, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...))
 end
@@ -135,7 +135,11 @@ end
 #=
 Finds index of the given `species` in a species list.
 =#
-assign(list::Vector, item::Union{AbstractString, Symbol}) = list[findfirst(i -> isequal(getfield(i, 2), Symbol(item)), list)]
+function assign(list::Vector, item::Union{AbstractString, Symbol})
+    index = findfirst(i -> isequal(getfield(i, 2), Symbol(item)), list)
+    @assert !isnothing(index) "$item is not present in the specified mechanism. Please check the mechanism files and ensure that $item is included"
+    return list[index]
+end
 
 #=
 Checks if the given `reaction` is reversible. Returns true if reversible, false otherwise.
@@ -179,7 +183,7 @@ Decomposes the given `reaction` into pairs of species and number of moles, for b
 =#
 function decompose(reaction::AbstractString, species_list::Vector{Species{N}}) where {N<:Number}
     sides = split(reaction, r"\s*<?=>?\s*", keepempty=false) ## left and right side of the equation
-    sides_components = Tuple(without_spaces(side) |> side -> split(side, r"\+|\(\+|(?<=M)\)"i; keepempty=false) for side in sides) ## products or reactants
+    sides_components = Tuple(remove_spaces(side) |> side -> split(side, r"\+|\(\+|(?<![S|L|G])\)"i; keepempty=false) for side in sides) ## products or reactants
     species_moles = Tuple(Pair{Species{N},N}[Pair(parse_moles(s, pair, species_list)...) for pair in components if pair ≠ "M"] for (s, components) in enumerate(sides_components))
     foreach(combine_moles!, species_moles)
     return species_moles
@@ -197,7 +201,7 @@ find_auxillaries(type::Symbol, data::AbstractString, N::Type{<:Number}) = (isequ
 Finds the `PLOG` parameters of an `ElementaryReaction` in the given reaction `data`. Returns a `Plog` object with the found parameters.
 =#
 find_plog(data::AbstractString, order::Int, N::Type{<:Number}) = Plog(([parse(N, m[p].match) for m in partition(eachmatch(:plog, data), 4)] for p in OneTo(4))...; order)
-find_plog(plogs::Vector{<:Dict}, order::Int, N::Type{<:Number}) = Plog{N}([plogs[p]["P"] |> without_spaces |> uparse |> Fix1(ustrip, u"Pa") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)])
+find_plog(plogs::Vector{<:Dict}, order::Int, N::Type{<:Number}) = Plog{N}([plogs[p]["P"] |> remove_spaces |> uparse |> Fix1(ustrip, u"Pa") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)])
 
 #=
 Finds the enhancement factors in the given reaction `data`. Returns vector of pairs with species index as keys and enhancement factors as values.
@@ -343,10 +347,10 @@ function get_stoichiometry_matrix(species::Vector{<:Species}) ## fail in case a
     return sparse(K, I, V)
 end
 
-function _read_mechanism(kinetics::T, thermo::T, transport::T, N::Type{<:Number}) where {T<:String}
+function _read_mechanism(kinetics::T, thermo::T, transport::Maybe{T}, N::Type{<:Number}) where {T<:String}
     kinetics_data = read_data(kinetics)
     thermo_data = read_data(thermo)
-    trans_data = isfile(transport) ? read_data(transport) : nothing
+    trans_data = isnothing(transport) ? nothing : isfile(transport) ? read_data(transport) : nothing
 
     species = find_species(kinetics_data, thermo_data, trans_data, N)
     K = length(species)
@@ -367,18 +371,25 @@ function _read_mechanism(mechanism::Dict, N::Type{<:Number})
     return (; species, reactions, stoichiometry_matrix)
 end
 
-function read_mechanism(name::Union{Nothing,Symbol,String} = nothing; kinetics_path = "", thermo_path = "", transport_path = "", as::Type{<:Number} = Float64)
-    name isa String && last(name, 5) == ".yaml" && return _read_mechanism(load_file(name), as)
+function read_mechanism(name::Union{Nothing, Symbol, String} = nothing;
+    kinetics_path::Maybe{String} = nothing,
+    thermo_path::Maybe{String} = nothing,
+    transport_path::Maybe{String} = nothing,
+    as::Type{<:Number} = Float64)
 
-    if name isa Symbol
-        mechanism_path = pkgdir(Apophis, "test/mechanisms/$(name)")
-        kinetics = mechanism_path * "/kinetics.dat"
-        thermo = mechanism_path * "/thermo.dat"
-        transport = mechanism_path * "/transport.dat"
-    else
+    name isa String && endswith(name, r"\.(:?yaml|yml)") && return _read_mechanism(load_file(name), as)
+
+    if name isa Symbol || name isa String
+        mechanism_path = pkgdir(Apophis, "test/mechanisms/$name")
+        kinetics = joinpath(mechanism_path, "kinetics.dat")
+        thermo = joinpath(mechanism_path, "thermo.dat")
+        transport = joinpath(mechanism_path, "transport.dat")
+    elseif !isnothing(kinetics_path) && !isnothing(thermo_path)
         kinetics = kinetics_path
         thermo = thermo_path
-        transport = transport_path
+        transport = isnothing(transport_path) ? nothing : transport_path
+    else
+        throw(ArgumentError("Invalid arguments"))
     end
     return _read_mechanism(kinetics, thermo, transport, as)
 end

src/sensitivity.jl

@@ -20,6 +20,8 @@ _dkfdA((; forward_rate_parameters)::Union{ElementaryReaction{N}, ThreeBodyReacti
 function _dkfdA((; high_pressure_parameters, low_pressure_parameters, enhancement_factors, troe_parameters)::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Number} ## Add dkfdAₒ later
     k∞, kₒ = (high_pressure_parameters, low_pressure_parameters)(T)
     M = total_molar_concentration(C, enhancement_factors)
+    iszero(M) && return zero(N)
+
     Pᵣ = reduced_pressure(kₒ, M, k∞)
     t = inv(one(N) + Pᵣ)
 

src/species.jl

@@ -54,7 +54,6 @@ struct Species{N<:Number, R<:AbstractReaction{N}} <: AbstractSpecies{N}
     weight::N
     nasa_polynomial::NasaPolynomial{N}
     transport_parameters::Maybe{TransportParameters{N}}
-    diffusions::NamedTuple{Diffusions, Tuple{Vector{N}, N}}
     thermo::Thermodynamics{N}
     rates::SpeciesRates{N}
 end

src/utils.jl

@@ -14,7 +14,7 @@ Base.show(io::IO, (; state)::Gas) = foreach(f -> println(io, "$f: $(getfield(sta
 
 function find_init((; mechanism, state)::Gas{N}, I) where {N<:Number}
     S = zero(state.X)
-    for (s, f) in eachmatch(regextionary[:inputs], I)
+    for (s, f) in eachmatch(:inputs, I)
         species = assign(mechanism.species, s)
         S[species.k] = parse(N, f)
     end
@@ -37,7 +37,7 @@ function set!(gas::Gas{N}; init...) where {N<:Number}
         Xᵢ = init[:X]
         X = Xᵢ isa AbstractVector ? Xᵢ : find_init(gas, init[:X])
         P = init[:P] |> (P -> P isa Quantity ? ustrip(u"Pa", P) : P) |> N
-        TPX!(gas, T, P, X) |> N
+        TPX!(gas, T, P, X)
     elseif haskey(init, :Y) && haskey(init, :ρ)
         Xᵢ = init[:X]
         X = Xᵢ isa AbstractVector ? Xᵢ : find_init(gas, init[:X])