Fixed plogs

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.6"
+version = "1.3.7"
 
 [deps]
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

README.md

@@ -1,5 +1,5 @@
 # Apophis
-Welcome to Apophis! This package is designed to provide sensitivity analysis of combustion processes using efficient and accurate derivatives. It is built to be fast, reliable, and easy to use, with a range of features and tools that can help to get the most out of combustion simulations.
+Welcome to Apophis! This package is designed to provide sensitivity analysis of combustion processes using efficient and accurate derivatives. It is built to be fast, reliable, and easy to use; but still early in development!
 
 ## Installation
 ```julia

src/Apophis.jl

@@ -22,7 +22,7 @@ const d = 0.14
 const Maybe{T} = Union{T, Nothing}
 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}}}
+const Rates{N<:Number} = NamedTuple{(:val, :dT, :dP, :dC), NTuple{4, Vector{N}}}
 const SpeciesRates{N<:Number} = NamedTuple{(:ω̇,), Tuple{Rates{N}}}
 const ReactionRates{N<:Number} = NamedTuple{(:kf, :kr, :q), NTuple{3, Rates{N}}}
 

src/reactions.jl

@@ -31,46 +31,49 @@ struct Plog{N<:Number}
     arrhenius::Vector{Arrhenius{N}}
 end
 
-Plog(itr...; order=1) = isempty(first(itr)) ? nothing : Plog(first(itr), Arrhenius{Float64}[Arrhenius(p; order) for p in zip(rest(itr, 2)...)])
+Plog(itr...; order=1) = isempty(first(itr)) ? nothing : Plog(first(itr) * ustrip(u"Pa", 1u"atm"), Arrhenius{Float64}[Arrhenius(p; order) for p in zip(rest(itr, 2)...)])
+interpolate(kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = log(kᵢ) + (log(kᵢ₊₁) - log(kᵢ)) * (log(P) - log(Pᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ)) |> exp
 
-function isin(p::N, x::Vector{N}) where {N<:Number}
-    for i in eachindex(x)
-        x[i] ≤ p && x[i+1] > p && return i
-    end
+function interpolate(::Val{:dT}, kᵢ::N, kᵢ₊₁::N, dkᵢdT::N, dkᵢ₊₁dT::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number}
+    k = interpolate(kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
+    dkdT = k * (dkᵢdT / kᵢ + (dkᵢ₊₁dT / kᵢ₊₁ - dkᵢdT / kᵢ) * (log(P) - log(Pᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ)))
+    return k, dkdT
 end
 
-isout(p::N, x::Vector{N}) where {N<:Number} = p > last(x) ? lastindex(x) : p < first(x) ? firstindex(x) : nothing ## on the assumption that plogs are sorted, To-Do: sort in advance
-interpolate(kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = log(kᵢ) + (log(kᵢ₊₁) - log(kᵢ)) * (log(P) - log(Pᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ))
-interpolate(::Val{:dP}, kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number} = (log(kᵢ₊₁) - log(kᵢ)) / (log(Pᵢ₊₁) - log(Pᵢ))P
+function interpolate(::Val{:dP}, kᵢ::N, kᵢ₊₁::N, Pᵢ::N, Pᵢ₊₁::N, P::N) where {N<:Number}
+    k = interpolate(kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
+    dkdP = k * (log(kᵢ₊₁) - log(kᵢ)) / (P * (log(Pᵢ₊₁) - log(Pᵢ)))
+    return k, dkdP
+end
 
 function ((; pressures, arrhenius)::Plog{N})(T::N, P::N) where {N<:Number}
-    i = isout(P, pressures)
-    isnothing(i) || return arrhenius[i](T)
+    i = searchsortedlast(pressures, P)
+    iszero(i) && return first(arrhenius)(T)
+    i == lastindex(pressures) && return last(arrhenius)(T)
 
-    i = isin(P, pressures)
-    Pᵢ, Pᵢ₊₁ = P[i], P[i+1]
-    kᵢ, kᵢ₊₁ = arrhenius[i](T), arrhenius[i+1](T)
+    Pᵢ, Pᵢ₊₁ = pressures[i], pressures[i+1] # maybe log(P[i]) and log(P[i+1])
+    kᵢ, kᵢ₊₁ = (arrhenius[i], arrhenius[i+1])(T)
     return interpolate(kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
 end
 
 function ((; pressures, arrhenius)::Plog{N})(v::Val{:dT}, T::N, P::N) where {N<:Number}
-    i = isout(P, pressures)
-    isnothing(i) || return arrhenius[i](v, T)
+    i = searchsortedlast(pressures, P)
+    iszero(i) && return first(arrhenius)(v, T)
+    i == lastindex(pressures) && return last(arrhenius)(v, T)
 
-    i = isin(P, pressures)
-    Pᵢ, Pᵢ₊₁ = P[i], P[i+1]
-    kᵢ, dkᵢdT, kᵢ₊₁, dkᵢ₊₁dT = arrhenius[i](v, T), arrhenius[i+1](v, T)
-    return interpolate(kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P), interpolate(dkᵢdT / kᵢ, dkᵢ₊₁dT / kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
+    Pᵢ, Pᵢ₊₁ = pressures[i], pressures[i+1]
+    (kᵢ, dkᵢdT), (kᵢ₊₁, dkᵢ₊₁dT) = (arrhenius[i], arrhenius[i+1])(v, T)
+    return interpolate(v, kᵢ, kᵢ₊₁, dkᵢdT, dkᵢ₊₁dT, Pᵢ, Pᵢ₊₁, P)
 end
 
 function ((; pressures, arrhenius)::Plog{N})(v::Val{:dP}, T::N, P::N) where {N<:Number}
-    i = isout(P, pressures)
-    isnothing(i) || return arrhenius[i](T), zero(N)
+    i = searchsortedlast(pressures, P)
+    iszero(i) && return first(arrhenius)(T), zero(N)
+    i == lastindex(pressures) && return last(arrhenius)(T), zero(N)
 
-    i = isin(P, pressures)
-    Pᵢ, Pᵢ₊₁ = P[i], P[i+1]
-    kᵢ, kᵢ₊₁ = arrhenius[i](T), arrhenius[i+1](T)
-    return interpolate(kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P), interpolate(v, kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
+    Pᵢ, Pᵢ₊₁ = pressures[i], pressures[i+1]
+    kᵢ, kᵢ₊₁ = (arrhenius[i], arrhenius[i+1])(T)
+    return interpolate(v, kᵢ, kᵢ₊₁, Pᵢ, Pᵢ₊₁, P)
 end
 
 struct Troe{N<:Number}
@@ -200,7 +203,7 @@ const Reaction{N<:Number, S<:AbstractSpecies{N}} = Union{ElementaryReaction{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 α; init=zero(N))
 
-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, 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)
 
 function forward_rate(reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Number}
@@ -217,7 +220,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[] = forward_rate_parameters(v, T) # isnothing(plog_parameters) ? forward_rate_parameters(v, T) : plog_parameters(v, T, P)
+    @inbounds rates.kf.val[], rates.kf.dT[] = isnothing(plog_parameters) ? forward_rate_parameters(v, T) : plog_parameters(v, T, P)
     return nothing
 end
 
@@ -257,8 +260,14 @@ function forward_rate(v::Val{:dT}, reaction::FallOffReaction{N}, (; T, C)::State
     return nothing
 end
 
-forward_rate(::Val{:dC}, reaction::ElementaryReaction{N}, state::State{N}) where {N<:Number} = forward_rate(reaction, state)
-forward_rate(::Val{:dC}, reaction::ThreeBodyReaction{N}, state::State{N}) where {N<:Number} = forward_rate(reaction, state)
+forward_rate(::Val{:dP}, reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = forward_rate(reaction, state)
+
+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.dP[] = isnothing(plog_parameters) ? (forward_rate_parameters(T), zero(N)) : plog_parameters(v, T, P)
+    return nothing
+end
+
+forward_rate(::Val{:dC}, reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = forward_rate(reaction, state)
 
 function forward_rate(v::Val{:dC}, reaction::FallOffReaction{N}, (; T, C)::State{N}) where {N<:Number}
     (; rates, high_pressure_parameters, low_pressure_parameters, troe_parameters, enhancement_factors) = reaction
@@ -282,11 +291,6 @@ 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
-
 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))
 
@@ -356,20 +360,21 @@ function reverse_rate(v::Val{:dT}, reaction::AbstractReaction{N}, (; T)::State{N
     return nothing
 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)
+reverse_rate(::Val{:dP}, reaction::AbstractReaction{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::ElementaryReaction{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
+
+reverse_rate(::Val{:dC}, reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = reverse_rate(reaction, state)
 
 function reverse_rate(::Val{:dC}, reaction::FallOffReaction{N}, (; T)::State{N}) where {N<:Number}
     (; rates, isreversible, reverse_rate_parameters) = reaction
@@ -392,8 +397,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; 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)
+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
@@ -405,6 +410,18 @@ function progress_rate(reaction::AbstractReaction{N}, (; C)::State{N}) where {N<
     return nothing
 end
 
+progress_rate(::Val{:dP}, reaction::AbstractReaction{N}, state::State{N}) where {N<:Number} = progress_rate(reaction, state)
+
+function progress_rate(::Val{:dP}, 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)
+
+    @inbounds q.val[] = kf.val[] * ∏ᴵ - kr.val[] * ∏ᴵᴵ
+    @inbounds q.dP[] = kf.dP[] * ∏ᴵ - kr.dP[] * ∏ᴵᴵ
+    return nothing
+end
+
 function progress_rate(::Val{:dT}, reaction::AbstractReaction{N}, (; C)::State{N}) where {N<:Number}
     (; kf, kr, q) = reaction.rates
     ∏ᴵ = step(reaction.reactants, C)
@@ -462,15 +479,5 @@ 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)
-
-#     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)

src/reader.jl

@@ -107,7 +107,8 @@ function _find_species(k::Int, species::Symbol, thermo::T, transport::Maybe{T},
     thermodynamics_parameters = find_thermodynamic_parameters(species, thermo, 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)))...)))...))
+        (; zip((:ω̇,), tuple((; zip((:val, :dT, :dP, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...)
+    )
 end
 
 function _find_species(k::Int, species::Dict, N::Type{<:Number})
@@ -118,7 +119,8 @@ function _find_species(k::Int, species::Dict, N::Type{<:Number})
     thermodynamic_parameters = find_thermodynamic_parameters(species["thermo"], 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)))...)))...))
+        (; zip((:ω̇,), tuple((; zip((:val, :dT, :dP, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 1), zeros(N, 0)))...)))...)
+    )
 end
 
 #=
@@ -178,14 +180,14 @@ function _combine_moles!(pairs::Vector{<:Pair}, species::Species)
 end
 
 combine_moles!(pairs::Vector{<:Pair}) = foreach(pair -> count(otherpair -> isequal(otherpair.first, pair.first), pairs) > 1 && _combine_moles!(pairs, pair.first), pairs)
-
+remove_rbracket(s::AbstractString) = replace(s, r"\(\+\w+\K\)" => "") # temporary
 #=
 Decomposes the given `reaction` into pairs of species and number of moles, for both of reactants and products.
 =#
 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(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))
+    sides_components = Tuple((remove_rbracket ∘ remove_spaces)(side) |> side -> split(side, r"\+|\(\+"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
 end
@@ -201,8 +203,8 @@ 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"] |> remove_spaces |> uparse |> Fix1(ustrip, u"Pa") for p in eachindex(plogs)], [plogs[p] |> Fix2(get_arrhenius, order) for p in eachindex(plogs)])
+find_plog(data::AbstractString, order::Real, 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::Real, 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.
@@ -234,8 +236,8 @@ function find_enhancements(data::Union{AbstractString, Dict}, species_list::Vect
     return enhancements
 end
 
-order(part::Vector{<:Pair}) = sum(Int ∘ abs ∘ last, part)
-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"))
+order(part::Vector{<:Pair}) = sum(abs ∘ last, part)
+unitfy_rate(A::Number, order::Real, 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`.
@@ -299,7 +301,7 @@ Find the reactions of the given kinetics `data`. Return the reactions found as a
 =#
 function _find_reaction(i::Int, data::RegexMatch, species::Vector{Species{N}}) 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)))...)
+    duals = (; zip((:kf, :kr, :q), ((; zip((:val, :dT, :dP, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 1), zeros(N, K)))...) for _ in OneTo(3)))...)
 
     reaction_data, reaction_string = data
     equation = Symbol(reaction_string)
@@ -318,14 +320,14 @@ end
 
 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)))...)
+    duals = (; zip((:kf, :kr, :q), ((; zip((:val, :dT, :dP, :dC), (zeros(N, 1), zeros(N, 1), zeros(N, 1), zeros(N, K)))...) for _ in OneTo(3)))...)
 
     reaction_string = reaction["equation"]
     equation = Symbol(reaction_string)
 
     isreversible = check_reversibility(reaction_string)
     components = decompose(reaction_string, species)
-    order = sum(last, flatten(components))
+    order = sum(last, flatten(components)) 
 
     Type = find_type(reaction_string, components)