Added routines for unit parsing

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.0"
+version = "1.3.1"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/Apophis.jl

@@ -14,8 +14,8 @@ using Zygote
 abstract type AbstractSpecies{N<:Number} end
 abstract type AbstractReaction{N<:Number} end
 
-const R = 8314.46261815324 # J K⁻¹ kmol⁻¹
-const Rc = 1.987261815324 # cal K⁻¹ kmol⁻¹
+const R = 8.31446261815324e3 # J K⁻¹ kmol⁻¹
+const Rc = 1.987261815324 # cal K⁻¹ mol⁻¹
 const kB = 1.380649e-23 # J K⁻¹
 const Pa = 101325.0 # Pa
 const Tᵣ = 300.0 # K

src/reactions.jl

@@ -4,10 +4,11 @@ struct Arrhenius{N<:Number} ## [A] := M^(∑νᵣ - 1) ⋅ s^-1; where M = cm^3
     E::N
 end
 
-function Arrhenius(itr; order=1)
+function Arrhenius(itr; order=1, mechunits = chemkin_default_units)
+    activation_energy_unit, length_unit, quantity_unit = (mechunits[k] for k in ("activation-energy", "length", "quantity"))
     isempty(itr) && return nothing
     A, β, E = itr
-    return Arrhenius(A * rate_units(order), β, E)
+    return Arrhenius(unitfy_rate(A, order, length_unit, quantity_unit), β, unitfy_activation_energy(E, activation_energy_unit))
 end
 
 ((; A, β, E)::Arrhenius{N})(T::N) where {N<:Number} = A * T^β * exp(-E * inv(Rc * T))

src/reader.jl

@@ -2,6 +2,7 @@
 Contains regular expressions for finding data. To-Do: brief explanation.
 =#
 const regextionary = Dict{Symbol, Union{Regex, Function}}(
+    :units => r"(?:\G(?!^)|reactions|reac)\s+\K(?:(?!reac|reactions|\n)CAL\/MOLE|KCAL\/MOLE|JOULES\/MOLE|KELVINS|EVOLTS|MOLES|MOLECULES)"i,
     :elements => r"(?:\G(?!^)|elements|elem)\s+\K(?:(?!end|elements|elem|thermo|species|spec)\S){1,2}"i,
     :species => r"(?:\G(?!^)|species|spec)\s+\K(?:(?!end|species|spec|thermo|reactions|reac)\S){1,16}"i,
     :reactions => r"(\s*(.+<?\s*=\s*>?.+?)(?<=\S)(?!\S)(?s:.+?))(?=\g'2'|end|\z)"i,
@@ -33,6 +34,13 @@ const weights = Dict{Symbol, Float64}( ## Reference: Google Search
     :He => 4.002602
 ) ## contains molecular weight of common elements in combustion.
 
+const chemkin_default_units = Dict{String, String}(
+    "length" => "cm",
+    "time" => "s",
+    "activation-energy" => "cal/mol",
+    "quantity" => "mol"
+)
+
 without_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.
@@ -118,6 +126,12 @@ Finds the properties of the species found in the given kinetics `data`. Return t
 find_species(kinetics::T, thermo::T, transport::Maybe{T}, N::Type{<:Number}) where {T<:String} = Species{N}[_find_species(k, Symbol(m.match), thermo, transport, N) for (k, m) in enumerate(eachmatch(:species, kinetics))]
 find_species(species::Vector{<:Dict}, N::Type{<:Number}) = Species{N}[_find_species(k, s, N) for (k, s) in enumerate(species)]
 
+function find_mechanism_units(data::String) 
+    units_found = Tuple(s for s in eachmatch(:units, data))
+    mechanism_units = Dict{String, String}(contains(u, r"^MOLES|MOLECULES$") ? "activation-energy" => u : "quantity" => u for u in units_found)
+    return mechanism_units
+end
+
 #=
 Finds index of the given `species` in a species list.
 =#
@@ -144,7 +158,7 @@ end
 Parses the number of moles for the given reactant or product pair. A pair means ...
 =#
 function parse_moles(s::Int, pair::AbstractString, species_list::Vector{Species{N}}) where {N<:Number}
-    moles_string, species_string = match(regextionary[:moles], pair)
+    moles_string, species_string = match(:moles, pair)
     species = assign(species_list, species_string)
     moles = (-1)^s * (isnothing(moles_string) ? one(N) : parse(N, moles_string))
     return species, moles
@@ -171,19 +185,19 @@ function decompose(reaction::AbstractString, species_list::Vector{Species{N}}) w
     return species_moles
 end
 
-get_arrhenius(arrhenius::Dict, order::Int) = Arrhenius(arrhenius[p] for p in ("A", "b", "Ea"); order)
+get_arrhenius(arrhenius::Dict, order::Int, mechunits::Dict) = Arrhenius(arrhenius[p] for p in ("A", "b", "Ea"); order, mechunits)
 
 #=
 Finds the auxillary parameters in the given reaction `data`. Auxillary parameters could be `LOW`, `TROE` or `REV` parameters.
 Returns the found parameters as object of type `Arrhenius` or `Troe`
 =#
-find_auxillaries(type::Symbol, data::AbstractString, N::Type{<:Number}) = (isequal(type, :troe) ? Troe : Arrhenius)(parse(N, p.match) for p in eachmatch(regextionary[type], data))
+find_auxillaries(type::Symbol, data::AbstractString, N::Type{<:Number}) = (isequal(type, :troe) ? Troe : Arrhenius)(parse(N, p.match) for p in eachmatch(type, data))
 
 #=
 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(regextionary[: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"dyn/cm^2") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)])
+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)])
 
 #=
 Finds the enhancement factors in the given reaction `data`. Returns vector of pairs with species index as keys and enhancement factors as values.
@@ -216,21 +230,22 @@ function find_enhancements(data::Union{AbstractString, Dict}, species_list::Vect
 end
 
 order(part::Vector{<:Pair}) = sum(Int ∘ abs ∘ last, part)
-rate_units(order::Int) = uparse("(m^3/kmol)^($order-1)/s") |> Fix2(ustrip, 1uparse("(cm^3/mol)^($order-1)/s"))
+unitfy_rate(A::Number, order::Int, l::String, q::String) = uparse("(m^3/kmol)^($order-1)/s") |> Fix2(ustrip, A * uparse("($l^3/$q)^($order-1)/s"))
+unitfy_activation_energy(E::Number, e::String) = ustrip(uparse("cal/mol"), E * uparse("$e"))
 #=
 Finds the kinetic parameters of an `ElementaryReaction` in the given reaction `data`.
 =#
 function find_kinetics(::Type{ElementaryReaction}, data::AbstractString, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
     forward_order, reverse_order = Tuple(order(p) for p in components)
-    forward_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:high], data); order = forward_order)
-    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:rev], data); order = reverse_order)
+    forward_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(:high, data); order = forward_order)
+    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(:rev, data); order = reverse_order)
     plog_parameters = find_plog(data, forward_order, N)
     return forward_parameters, reverse_paramters, plog_parameters
 end
 
-function find_kinetics(::Type{ElementaryReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
+function find_kinetics(::Type{ElementaryReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}, mechunits::Dict) where {N<:Number}
     forward_order = first(components) |> order
-    forward_parameters, plog_parameters = haskey(reaction, "rate-constants") ? (nothing, find_plog(reaction["rate-constants"], forward_order, N)) : (get_arrhenius(reaction["rate-constant"], forward_order), nothing)
+    forward_parameters, plog_parameters = haskey(reaction, "rate-constants") ? (nothing, find_plog(reaction["rate-constants"], forward_order, N)) : (get_arrhenius(reaction["rate-constant"], forward_order, mechunits), nothing)
     return forward_parameters, nothing, plog_parameters
 end
 
@@ -240,15 +255,15 @@ Finds the kinetic parameters of a `ThreeBodyReaction` in the given reaction `dat
 function find_kinetics(::Type{ThreeBodyReaction}, data::AbstractString, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
     enhancements = find_enhancements(data, species, components)
     forward_order, reverse_order = Tuple(order(p) for p in components)
-    forward_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:high], data); order = forward_order + 1)
-    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:rev], data); order = reverse_order + 1)
+    forward_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(:high, data); order = forward_order + 1)
+    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(:rev, data); order = reverse_order + 1)
     return forward_parameters, reverse_paramters, enhancements
 end
 
-function find_kinetics(::Type{ThreeBodyReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
+function find_kinetics(::Type{ThreeBodyReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}, mechunits::Dict) where {N<:Number}
     enhancements = find_enhancements(reaction, species, components)
     forward_order = first(components) |> order
-    forward_parameters = reaction["rate-constant"] |> Fix2(get_arrhenius, forward_order + 1)
+    forward_parameters = reaction["rate-constant"] |> params -> get_arrhenius(params, forward_order + 1, mechunits)
     return forward_parameters, nothing, enhancements
 end
 
@@ -258,18 +273,18 @@ Finds the kinetic parameters of a `FallOffReaction` in the given reaction `data`
 function find_kinetics(::Type{FallOffReaction}, data::AbstractString, species_list::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
     enhancements = find_enhancements(data, species_list, components)
     forward_order, reverse_order = Tuple(order(p) for p in components)
-    high_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:high], data); order = forward_order)
-    low_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:low], data); order = forward_order + 1)
-    troe_parameters = Troe(parse(N, p.match) for p in eachmatch(regextionary[:troe], data))
-    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(regextionary[:rev], data); order = reverse_order)
+    high_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(:high, data); order = forward_order)
+    low_parameters = Arrhenius(parse(N, p.match) for p in eachmatch(:low, data); order = forward_order + 1)
+    troe_parameters = Troe(parse(N, p.match) for p in eachmatch(:troe, data))
+    reverse_paramters = Arrhenius(parse(N, p.match) for p in eachmatch(:rev, data); order = reverse_order)
     return high_parameters, low_parameters, troe_parameters, reverse_paramters, enhancements
 end
 
-function find_kinetics(::Type{FallOffReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}) where {N<:Number}
+function find_kinetics(::Type{FallOffReaction}, reaction::Dict, species::Vector{Species{N}}, components::NTuple{2, Vector{<:Pair}}, mechunits::Dict) where {N<:Number}
     enhancements = find_enhancements(reaction, species, components)
     forward_order = first(components) |> order
-    high_parameters = reaction["high-P-rate-constant"] |> Fix2(get_arrhenius, forward_order)
-    low_parameters = reaction["low-P-rate-constant"] |> Fix2(get_arrhenius, forward_order + 1)
+    high_parameters = reaction["high-P-rate-constant"] |> params -> get_arrhenius(params, forward_order, mechunits)
+    low_parameters = reaction["low-P-rate-constant"] |> params -> get_arrhenius(params, forward_order + 1, mechunits)
     troe_parameters = get(reaction, "Troe", nothing) |> troe -> isnothing(troe) ? Troe() : Troe(get(troe, p, nothing) for p in ("A", "T3", "T1", "T2"))
     return high_parameters, low_parameters, troe_parameters, nothing, enhancements
 end
@@ -296,7 +311,7 @@ function _find_reaction(i::Int, data::RegexMatch, species::Vector{Species{N}}) w
     return reaction
 end
 
-function _find_reaction(i::Int, reaction::Dict, species::Vector{Species{N}}) where {N<:Number}
+function _find_reaction(i::Int, reaction::Dict, species::Vector{Species{N}}, mechunits::Dict) where {N<:Number}
     K = length(species)
     duals = (; zip((:kf, :kr, :q), ((; zip((:val, :dT, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, K)))...) for _ in OneTo(3)))...)
 
@@ -309,7 +324,7 @@ function _find_reaction(i::Int, reaction::Dict, species::Vector{Species{N}}) whe
 
     Type = find_type(reaction_string, components)
 
-    kinetics = find_kinetics(Type, reaction, species, components)
+    kinetics = find_kinetics(Type, reaction, species, components, mechunits)
     reaction = Type(i, equation, isreversible, components..., order, kinetics..., duals)
     foreach(i -> push!(i.first.inreactions, Pair{Reaction{N, Species{N}}, N}(reaction, i.second)), flatten(components))
     return reaction
@@ -319,7 +334,7 @@ end
 Finds the reactions of the given kinetics `data`. Returns a vector of the reactions.
 =#
 find_reaction(data::String, species::Vector{Species{N}}) where {N<:Number} = Reaction{N, Species{N}}[_find_reaction(i, r, species) for (i, r) in enumerate(eachmatch(:reactions, data))]
-find_reaction(reactions::Vector{<:Dict}, species::Vector{Species{N}}) where {N<:Number} = Reaction{N, Species{N}}[_find_reaction(i, r, species) for (i, r) in enumerate(reactions)]
+find_reaction(reactions::Vector{<:Dict}, species::Vector{Species{N}}, mechunits::Dict) where {N<:Number} = Reaction{N, Species{N}}[_find_reaction(i, r, species, mechunits) for (i, r) in enumerate(reactions)]
 
 function get_stoichiometry_matrix(species::Vector{<:Species}) ## fail in case a species doesn't appear in any reaction
     K = [s.k for s in species for _ in s.inreactions]
@@ -345,7 +360,9 @@ function _read_mechanism(mechanism::Dict, N::Type{<:Number})
     species = find_species(mechanism["species"], N)
     K = length(species)
     foreach(s -> append!(s.rates.ω̇.dC, zeros(N, K)), species)
-    reactions = find_reaction(mechanism["reactions"], species)
+
+    mechunits = mechanism["units"]
+    reactions = find_reaction(mechanism["reactions"], species, mechunits)
     stoichiometry_matrix = get_stoichiometry_matrix(species)
     return (; species, reactions, stoichiometry_matrix)
 end

test/Manifest.toml

@@ -18,9 +18,9 @@ version = "3.5.0"
 
 [[deps.Apophis]]
 deps = ["Combinatorics", "LinearAlgebra", "SparseArrays", "SplitApplyCombine", "ThreadsX", "Unitful", "YAML", "Zygote"]
-git-tree-sha1 = "a55a15b9c7bb06f790702f5fe7bdbe4e9a2bfced"
-repo-rev = "tests-clean-up"
-repo-url = "/Users/sabry/Developer/ApophisOrg/Apophis.jl"
+git-tree-sha1 = "bae642067496fe63974edb0118a362dbde9a4de2"
+repo-rev = "default-units"
+repo-url = "https://github.com/ApophisOrg/Apophis.jl.git"
 uuid = "cf72b054-1e85-4718-a3f4-02cfdc10a054"
 version = "1.2.3"