diff --git a/README.md b/README.md
index e29dcf0f4..008cf8d31 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,7 @@ RMS has been used in many applications:
## How to cite
Please include the following citations for ReactionMechanismSimulator.jl in general and for transitory sensitivities and the automatic mechanism analysis toolkit respectively.
-- Johnson, M. S., Pang, H.-W., Payne, A. M., & Green, W. H. (2023). ReactionMechanismSimulator.jl: A Modern Approach to Chemical Kinetic Mechanism Simulation and Analysis. https://doi.org/10.26434/CHEMRXIV-2023-TJ34T
+- Johnson, M. S., Pang, H.-W., Payne, A. M., & Green, W. H. (2023). ReactionMechanismSimulator.jl: A Modern Approach to Chemical Kinetic Mechanism Simulation and Analysis. Int. J. Chem. Kinet. 2024;1-16 https://doi.org/10.1002/kin.21753
- Johnson, M. S., McGill, C. J., & Green, W. H. (2022). Transitory Sensitivity in Automatic Chemical Kinetic Mechanism Analysis. https://doi.org/10.26434/CHEMRXIV-2022-ZSFJC
## Installation
diff --git a/src/Calculators/Rate.jl b/src/Calculators/Rate.jl
index c873c58a0..37b1979f8 100644
--- a/src/Calculators/Rate.jl
+++ b/src/Calculators/Rate.jl
@@ -13,8 +13,8 @@ export AbstractFalloffRate
Ea::Q
unc::P = EmptyRateUncertainty()
end
-@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
-@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
+@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
+@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
export Arrhenius
@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
@@ -23,21 +23,42 @@ export Arrhenius
Ea::Q
unc::P = EmptyRateUncertainty()
end
-@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
-@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
+@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
+@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
export StickingCoefficient
-@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,B} <: AbstractRate
+@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,R<:Real,P<:AbstractRateUncertainty,B} <: AbstractRate
A::N
n::K
Ea::Q
q::B = 0.0
+ V0::R = 0.0
unc::P = EmptyRateUncertainty()
end
-@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
-@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
+@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*(phi-arr.V0))/(R*T))
+@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*(phi-arr.V0))/(R*T))
export Arrheniusq
+@with_kw struct Marcus{N<:Real,K<:Real,Q,P<:AbstractRateUncertainty,B} <: AbstractRate
+ A::N
+ n::K
+ lmbd_i_coefs::Q
+ lmbd_o::K
+ wr::K
+ wp::K
+ beta::B
+ unc::P = EmptyRateUncertainty()
+end
+@inline function (arr::Marcus)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real}
+ @fastmath lmbd = arr.lmbd_o + evalpoly(T,arr.lmbd_i_coefs)
+ @fastmath arr.A*T^arr.n*exp(-lmbd/4.0*(1.0+dGrxn/lmbd)^2/(R*T)-arr.beta*d)
+end
+@inline function (arr::Marcus)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,N<:Real,S<:Real}
+ @fastmath lmbd = arr.lmbd_o + evalpoly(T,arr.lmbd_i_coefs)
+ @fastmath return arr.A*T^arr.n*exp(-lmbd/4.0*(1.0+dGrxn/lmbd)^2/(R*T)-arr.beta*d)
+end
+export Marcus
+
@with_kw struct PdepArrhenius{T<:Real,Q<:AbstractRateUncertainty,Z<:AbstractRate} <: AbstractRate
Ps::Array{T,1}
arrs::Array{Z,1}
@@ -45,7 +66,7 @@ export Arrheniusq
end
PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} = PdepArrhenius(sort(Ps),arrs)
-@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,V<:Real,S<:Real}
inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
if inds[2] == -1
@@ -65,7 +86,7 @@ export PdepArrhenius
unc::Q = EmptyRateUncertainty()
end
-@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
out = 0.0
for arr in marr.arrs
@fastmath out += arr(T)
@@ -79,7 +100,7 @@ export MultiArrhenius
unc::Q = EmptyRateUncertainty()
end
-@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
out = 0.0
for pdar in parr.parrs
@fastmath out += pdar(T=T,P=P)
@@ -95,7 +116,7 @@ export MultiPdepArrhenius
unc::Q = EmptyRateUncertainty()
end
-(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
+(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
export ThirdBody
@with_kw struct Lindemann{N<:Integer,K<:AbstractFloat,Q<:AbstractRateUncertainty} <: AbstractFalloffRate
@@ -106,7 +127,7 @@ export ThirdBody
unc::Q = EmptyRateUncertainty()
end
-@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = lnd.arrlow(T=T)
kinf = lnd.arrhigh(T=T)
@fastmath Pr = k0*C/kinf
@@ -126,7 +147,7 @@ export Lindemann
unc::R = EmptyRateUncertainty()
end
-@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0,dGrxn=0.0,d=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = tr.arrlow(T=T)
kinf = tr.arrhigh(T=T)
@fastmath Pr = k0*C/kinf
@@ -195,7 +216,7 @@ export getredtemp
end
export getredpress
-@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,dGrxn=0.0,d=0.0) where {N<:Real,B<:Real,Q<:Real}
k = 0.0
Tred = getredtemp(ch,T)
Pred = getredpress(ch,P)
@@ -228,7 +249,7 @@ export getkineticstype
@inline extracttypename(typ::Symbol) = string(typ)
@inline extracttypename(typ) = string(typ.name)
-
+
@inline function _calcdkdCeff(tbarr::ThirdBody,T::Float64,Ceff::Float64)
return @fastmath tbarr.arr(T)
end
diff --git a/src/Calculators/Ratevec.jl b/src/Calculators/Ratevec.jl
index 7b5045743..7b9dcacbb 100644
--- a/src/Calculators/Ratevec.jl
+++ b/src/Calculators/Ratevec.jl
@@ -21,7 +21,7 @@ function Arrheniusvec(arrs::T) where {T<:AbstractArray}
end
return Arrheniusvec(A=A,n=n,Ea=Ea)
end
-@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,dGrxns=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
export Arrheniusvec
@@ -86,7 +86,7 @@ export getredtemp
end
export getredpress
-@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,dGrxns=0.0) where {N<:Real,B<:Real,Q<:Real}
k = zeros(N,size(ch.coefs)[1])
Tred = getredtemp(ch,T)
Pred = getredpress(ch,P)
@@ -170,7 +170,7 @@ function Troevec(troes::T) where {T<:AbstractArray}
end
export Troevec
-@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,dGrxns=0.0) where {Q<:Real,R<:Real,S<:Real}
k0 = tr.arrlow(T=T)
kinf = tr.arrhigh(T=T)
@fastmath Pr = k0.*C./kinf
@@ -201,7 +201,7 @@ function PdepArrheniusvec(pdeparrs::T) where {T<:AbstractArray}
end
export PdepArrheniusvec
-@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,dGrxns=0.0) where {Q<:Real,V<:Real,S<:Real}
inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
if inds[2] == -1
return @inbounds parr.arrvecs[inds[1]](T=T)
diff --git a/src/Domain.jl b/src/Domain.jl
index 82d2e05ea..a304b055b 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -6,6 +6,7 @@ using SciMLBase
using ForwardDiff
using Tracker
using ReverseDiff
+using Logging
abstract type AbstractDomain end
export AbstractDomain
@@ -85,7 +86,7 @@ function ConstantTPDomain(; phase::E2, initialconds::Dict{X,X2}, constantspecies
C = P / (R * T)
V = N * R * T / P
y0[end] = V
- kfs, krevs = getkfkrevs(phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
+ kfs, krevs = getkfkrevs(phase, T, P, C, N, ns, Gs, diffs, V, 0.0, 0.0)
kfsp = deepcopy(kfs)
for ind in efficiencyinds
kfsp[ind] = 1.0
@@ -481,6 +482,7 @@ mutable struct ConstantTVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:
T::W
V::W
phi::W
+ d::W
kfs::Array{W,1}
krevs::Array{W,1}
kfsnondiff::Array{W,1}
@@ -504,6 +506,7 @@ function ConstantTVDomain(; phase::Z, initialconds::Dict{X,E}, constantspecies::
V = 0.0
P = 1.0e8
phi = 0.0
+ d = 0.0
y0 = zeros(length(phase.species))
spnames = [x.name for x in phase.species]
for (key, val) in initialconds
@@ -515,6 +518,8 @@ function ConstantTVDomain(; phase::Z, initialconds::Dict{X,E}, constantspecies::
V = val
elseif key == "Phi"
phi = val
+ elseif key == "d"
+ d = val
else
ind = findfirst(isequal(key), spnames)
@assert typeof(ind) <: Integer "$key not found in species list: $spnames"
@@ -545,18 +550,19 @@ function ConstantTVDomain(; phase::Z, initialconds::Dict{X,E}, constantspecies::
diffs = Array{Float64,1}()
end
P = 1.0e8 #essentiallly assuming this is a liquid
- C = N / V
- kfs, krevs = getkfkrevs(phase, T, P, C, N, ns, Gs, diffs, V, phi)
- kfsnondiff = getkfs(phase, T, P, C, ns, V, phi)
- p = vcat(Gs, kfsnondiff)
+ C = N/V
+ dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
+ kfs,krevs = getkfkrevs(phase,T,P,C,N,ns,Gs,diffs,V,phi,d)
+ kfsnondiff = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
+ p = vcat(Gs,kfsnondiff)
if sparse
jacobian = zeros(typeof(T), length(phase.species), length(phase.species))
else
jacobian = zeros(typeof(T), length(phase.species), length(phase.species))
end
rxnarray = getreactionindices(phase)
- return ConstantTVDomain(phase, [1, length(phase.species)], [1, length(phase.species) + length(phase.reactions)], constspcinds,
- T, V, phi, kfs, krevs, kfsnondiff, efficiencyinds, Gs, rxnarray, mu, diffs, jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}()), y0, p
+ return ConstantTVDomain(phase,[1,length(phase.species)],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+ T,V,phi,d,kfs,krevs,kfsnondiff,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,false,MVector(false),MVector(0.0),p,Dict{String,Int64}()), y0, p
end
export ConstantTVDomain
@@ -568,6 +574,7 @@ struct ParametrizedTConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real
T::FT
V::W
phi::W
+ d::W
efficiencyinds::Array{I,1}
rxnarray::Array{Int64,2}
jacobian::Array{W,2}
@@ -584,6 +591,7 @@ function ParametrizedTConstantVDomain(; phase::IdealDiluteSolution, initialconds
P = 1.0e8 #essentiallly assuming this is a liquid
V = 0.0
phi = 0.0
+ d = 0.0
ts = Array{Float64,1}()
ns = zeros(length(phase.species))
spnames = [x.name for x in phase.species]
@@ -599,6 +607,8 @@ function ParametrizedTConstantVDomain(; phase::IdealDiluteSolution, initialconds
ts = val
elseif key == "Phi"
phi = val
+ elseif key == "d"
+ d = val
else
ind = findfirst(isequal(key), spnames)
@assert typeof(ind) <: Integer "$key not found in species list: $spnames"
@@ -629,8 +639,8 @@ function ParametrizedTConstantVDomain(; phase::IdealDiluteSolution, initialconds
jacobian = zeros(typeof(V), length(phase.species) + 1, length(phase.species) + 1)
end
rxnarray = getreactionindices(phase)
- return ParametrizedTConstantVDomain(phase, [1, length(phase.species)], [1, length(phase.species) + length(phase.reactions)], constspcinds,
- Tfcn, V, phi, efficiencyinds, rxnarray, jacobian, sensitivity, MVector(false), MVector(0.0), p, Dict{String,Int64}()), y0, p
+ return ParametrizedTConstantVDomain(phase,[1,length(phase.species)],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+ Tfcn,V,phi,d,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict{String,Int64}()), y0, p
end
export ParametrizedTConstantVDomain
@@ -642,6 +652,7 @@ mutable struct ConstantTAPhiDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,
T::W
A::W
phi::W
+ d::W
kfs::Array{W,1}
krevs::Array{W,1}
efficiencyinds::Array{I,1}
@@ -663,6 +674,7 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
T = 0.0
A = 0.0
phi = 0.0 #default 0.0
+ d = 0.0
y0 = zeros(length(phase.species))
spnames = [x.name for x in phase.species]
for (key, val) in initialconds
@@ -672,6 +684,8 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
A = val
elseif key == "Phi"
phi = val
+ elseif key == "d"
+ d = val
else
ind = findfirst(isequal(key), spnames)
@assert typeof(ind) <: Integer "$key not found in species list: $spnames"
@@ -698,16 +712,16 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
mu = 0.0
end
C = 0.0 #this currently shouldn't matter here, on a surface you shouldn't have pdep
- kfs, krevs = getkfkrevs(phase, T, 0.0, C, N, ns, Gs, [], A, phi)
- p = vcat(Gs, kfs)
+ kfs,krevs = getkfkrevs(phase,T,0.0,C,N,ns,Gs,[],A,phi,d)
+ p = vcat(Gs,kfs)
if sparse
jacobian = spzeros(typeof(T), length(phase.species), length(phase.species))
else
jacobian = zeros(typeof(T), length(phase.species), length(phase.species))
end
rxnarray = getreactionindices(phase)
- return ConstantTAPhiDomain(phase, [1, length(phase.species)], [1, length(phase.species) + length(phase.reactions)], constspcinds,
- T, A, phi, kfs, krevs, efficiencyinds, Gs, rxnarray, mu, Array{Float64,1}(), jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}()), y0, p
+ return ConstantTAPhiDomain(phase,[1,length(phase.species)],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+ T,A,phi,d,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,Array{Float64,1}(),jacobian,sensitivity,false,MVector(false),MVector(0.0),p,Dict{String,Int64}()), y0, p
end
export ConstantTAPhiDomain
@@ -803,11 +817,11 @@ function FragmentBasedConstantTrhoDomain(; phase::Z, initialconds::Dict{X,E}, co
diffs = Array{Float64,1}()
end
- C = N / V
- kfs, krevs = getkfkrevs(phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- kfsnondiff = getkfs(phase, T, P, C, ns, V, 0.0)
-
- p = vcat(Gs, kfsnondiff)
+ C = N/V
+ dGrxns = -(phase.stoichmatrix*Gs)
+ kfs,krevs = getkfkrevs(phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ kfsnondiff = getkfs(phase,T,P,C,ns,V,0.0,dGrxns,0.0)
+ p = vcat(Gs,kfsnondiff)
if sparse
jacobian = zeros(typeof(T), length(getphasespecies(phase)), length(getphasespecies(phase)))
else
@@ -920,7 +934,7 @@ export ConstantTLiqFilmDomain
cs = ns ./ V
C = N / V
for ind in d.efficiencyinds #efficiency related rates may have changed
- d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0)
+ d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0)
end
return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end
@@ -941,13 +955,13 @@ end
@views nokfchg = count(d.kfs .!= kfps) <= length(d.efficiencyinds) && all(kfps[d.efficiencyinds] .== 1.0)
if nothermochg && nokfchg
for ind in d.efficiencyinds #efficiency related rates may have changed
- d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
+ d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;f=kfps[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
elseif nothermochg
d.kfs = kfps
for ind in d.efficiencyinds #efficiency related rates may have changed
- d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
+ d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;f=kfps[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
else #need to handle thermo changes
@@ -957,9 +971,9 @@ end
else
d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
- krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfs)[2]
+ krevs = getkfkrevs(d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;kfs=d.kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
- d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
+ d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;f=kfps[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end
@@ -978,9 +992,9 @@ end
kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
- krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2])
+ krevs = convert(typeof(y),getkfkrevs(d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfs)[2])
for ind in d.efficiencyinds #efficiency related rates may have changed
- kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
+ kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;f=kfs[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end
@@ -998,9 +1012,9 @@ end
kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
- krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
+ krevs = getkfkrevs(d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
- kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
+ kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;f=kfs[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end
@@ -1020,9 +1034,9 @@ end
kfs .= d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
- krevs .= getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
+ krevs .= getkfkrevs(d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
- kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
+ kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;f=kfs[ind])
end
return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end
@@ -1035,13 +1049,13 @@ end
C = N / V
if !d.alternativepformat
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
+ kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0)[1]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
else
- Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
+ Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+ kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0)[1]*d.p[length(d.phase.species)+ind]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
end
- krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
- return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ krevs = getkfkrevs(d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfs)[2]
+ return ns,cs,d.T,d.P,V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Tracker/reversediff
@@ -1052,13 +1066,13 @@ end
C = N / V
if !d.alternativepformat
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
+ kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0)[1]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
else
- Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
+ Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+ kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0)[1]*d.p[length(d.phase.species)+ind]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
end
- krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
- return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ krevs = getkfkrevs(d.phase,d.T,d.P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfs)[2]
+ return ns,cs,d.T,d.P,V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1082,8 +1096,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
- return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,d.V,0.0,0.0)
+ return ns,cs,T,P,d.V,C,N,0.0,kfs,krevs,Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1108,8 +1122,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
- return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,d.V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,d.V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1134,8 +1148,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
- return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,d.V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,d.V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1158,8 +1172,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
- return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,d.P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return ns,cs,T,d.P,V,C,N,0.0,kfs,krevs,Hs,Array{Float64,1}(),Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1183,7 +1197,7 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
+ kfs,krevs = getkfkrevs(d.phase,T,d.P,C,N,ns,Gs,diffs,V,0.0,0.0)
if p != SciMLBase.NullParameters()
return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
else
@@ -1212,8 +1226,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,d.P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,d.P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Hs,Array{Float64,1}(),Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1237,8 +1251,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return ns,cs,T,P,V,C,N,0.0,kfs,krevs,Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1263,8 +1277,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1289,8 +1303,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Us,Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1314,8 +1328,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return ns,cs,T,P,V,C,N,0.0,kfs,krevs,Hs,Array{Float64,1}(),Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1340,8 +1354,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Hs,Array{Float64,1}(),Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
P = d.P(t)
@@ -1365,8 +1379,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Hs,Array{Float64,1}(),Gs,diffs,Cvave,cpdivR,0.0
end
@inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1386,8 +1400,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
- return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,d.phi,d.d)
+ return ns,cs,T,P,V,C,N,mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1408,8 +1422,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
- return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,d.phi,d.d)
+ return @views @fastmath ns,cs,T,P,V,C,N,mu,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),d.phi
end
@inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1430,8 +1444,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
- return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,d.phi,d.d)
+ return @views @fastmath ns,cs,T,P,V,C,N,mu,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),d.phi
end
@inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1455,8 +1469,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return ns,cs,T,P,V,C,N,0.0,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1481,8 +1495,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1507,8 +1521,8 @@ end
else
diffs = Array{Float64,1}()
end
- kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
- return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+ kfs,krevs = getkfkrevs(d.phase,T,P,C,N,ns,Gs,diffs,V,0.0,0.0)
+ return @views @fastmath ns,cs,T,P,V,C,N,0.0,kfs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],krevs.*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)],Array{Float64,1}(),Array{Float64,1}(),Gs,diffs,0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1533,13 +1547,13 @@ end
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
elseif nothermochg
d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
else
d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
else
@views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
@@ -1547,14 +1561,14 @@ end
if nothermochg && nokfchg
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
elseif nothermochg
- d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
else
- d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+ d.Gs .= d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
end
@@ -1573,8 +1587,8 @@ end
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
end
- kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ kfs,krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,d.phi,d.d;kfs=kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
@inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1590,8 +1604,8 @@ end
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
end
- kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
- return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ kfs,krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,d.phi,d.d;kfs=kfsnondiff)
+ return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
@inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q}
@@ -1616,18 +1630,18 @@ end
else
d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
- return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfs)[2]
+ return ns,cs,d.T,P,d.A,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),d.Gs,Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
else
@views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
if nothermochg
return ns, cs, d.T, P, d.A, C, N, d.mu, d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, 0.0, Array{Float64,1}(), d.phi
else
- d.kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
- return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ d.kfs = d.p[length(d.phase.species)+1:end].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+ d.Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+ d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,d.phi,d.d;kfs=d.kfs)[2]
+ return ns,cs,d.T,P,d.A,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),d.Gs,Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
end
end
@@ -1645,8 +1659,8 @@ end
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
end
- krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2])
- return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ krevs = convert(typeof(y),getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.A,d.phi,d.d;kfs=kfs)[2])
+ return ns,cs,d.T,P,d.A,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Gs,Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
@inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1662,8 +1676,8 @@ end
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
end
- krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2]
- return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.A,d.phi,d.d;kfs=kfs)[2]
+ return ns,cs,d.T,P,d.A,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Gs,Array{Float64,1}(),0.0,Array{Float64,1}(),d.phi
end
function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q}
@@ -1692,13 +1706,13 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
elseif nothermochg
d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
else
d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
else
@views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
@@ -1706,14 +1720,14 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
if nothermochg && nokfchg
return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
elseif nothermochg
- d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
else
- d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
- d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
- d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+ d.Gs .= d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+ d.kfs,d.krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,V,0.0,0.0;kfs=d.kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,d.kfs,d.krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
end
@@ -1734,8 +1748,8 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::Array{W2,1}, t::
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
end
- kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+ kfs,krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1830,10 +1844,9 @@ end
Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
end
- kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
- return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+ kfs,krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,V,0.0,0.0;kfs=kfsnondiff)
+ return ns,cs,d.T,P,V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
-
export calcthermo
@inline function calcdomainderivatives!(d::Q, dydt::Z7, interfaces::Z12; t::Z10, T::Z4, P::Z9, Us::Array{Z,1}, Hs::Array{Z11,1}, V::Z2, C::Z3, ns::Z5, N::Z6, Cvave::Z8) where {Q<:AbstractDomain,Z12,Z11,Z10,Z9,Z8<:Real,Z7,W<:IdealGas,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5}
diff --git a/src/Interface.jl b/src/Interface.jl
index 466022abc..e1a9d1225 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -34,26 +34,28 @@ struct ReactiveInternalInterface{T,B,C,C2,N,Q<:AbstractReaction,X} <: AbstractRe
reversibililty::Array{Bool,1}
forwardability::Array{Bool,1}
end
-function ReactiveInternalInterface(domain1, domain2, reactions, A)
- vectuple, vecinds, otherrxns, otherrxninds, posinds = getveckinetics(reactions)
- rxns = vcat(reactions[vecinds], reactions[otherrxninds])
- rxns = [ElementaryReaction(index=i, reactants=rxn.reactants, reactantinds=rxn.reactantinds, products=rxn.products,
- productinds=rxn.productinds, kinetics=rxn.kinetics, radicalchange=rxn.radicalchange, reversible=rxn.reversible, forwardable=rxn.forwardable, pairs=rxn.pairs) for (i, rxn) in enumerate(rxns)]
- rxnarray = getinterfacereactioninds(domain1, domain2, rxns)
- M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
- reversibility = Array{Bool,1}(getfield.(rxns, :reversible))
- forwardability = Array{Bool,1}(getfield.(rxns, :forwardable))
- return ReactiveInternalInterface(domain1, domain2,
- rxns, vectuple, posinds, rxnarray, M, Nrp1, Nrp2, A, [1, length(reactions)],
- [0, 1], ones(length(rxns)), reversibility, forwardability), ones(length(rxns))
+function ReactiveInternalInterface(domain1,domain2,reactions,A)
+ vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
+ rxns = vcat(reactions[vecinds],reactions[otherrxninds])
+ rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
+ productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
+ rxnarray = getinterfacereactioninds(domain1,domain2,rxns)
+ M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
+ reversibility = Array{Bool,1}(getfield.(rxns,:reversible))
+ forwardability = Array{Bool,1}(getfield.(rxns,:forwardable))
+ return ReactiveInternalInterface(domain1,domain2,
+ rxns,vectuple,posinds,rxnarray,M,Nrp1,Nrp2,A,[1,length(reactions)],
+ [0,1],ones(length(rxns)),reversibility,forwardability),ones(length(rxns))
end
export ReactiveInternalInterface
-function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, cstot::Array{Q,1}) where {Q}
- kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
- Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), T1, phi1)
- krevs = kfs ./ Kc
- return kfs, krevs
+function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot::Array{Q,1}) where {Q}
+ Gpart = ArrayPartition(Gs1,Gs2)
+ dGrxns = -ri.stoichmatrix*Gpart
+ kfs = getkfs(ri,T1,0.0,0.0,Array{Q,1}(),ri.A,phi1,dGrxns,0.0)
+ Kc = getKc.(ri.reactions,ri.domain1.phase,ri.domain2.phase,Ref(Gs1),Ref(Gs2),dGrxns,T1,phi1)
+ krevs = kfs./Kc
+ return kfs,krevs
end
function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
@@ -86,18 +88,28 @@ struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReactio
reversibility::Array{Bool,1}
forwardability::Array{Bool,1}
end
-function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A, phi=0.0)
- @assert domain1.T == domain2.T
- reactions = upgradekinetics(reactions, domain1, domain2)
- rxnarray = getinterfacereactioninds(domain1, domain2, reactions)
- kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi)
- Kc = getKc.(reactions, domain1.phase, domain2.phase, Ref(domain1.Gs), Ref(domain2.Gs), T, phi)
- krevs = kfs ./ Kc
- M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
- reversibility = Array{Bool,1}(getfield.(reactions, :reversible))
- forwardability = Array{Bool,1}(getfield.(reactions, :forwardable))
- if isa(reactions, Vector{Any})
- reactions = convert(Vector{ElementaryReaction}, reactions)
+function ReactiveInternalInterfaceConstantTPhi(domain1,domain2,reactions,T,A,phi=0.0)
+ @assert domain1.T == domain2.T
+ reactions = upgradekinetics(reactions,domain1,domain2)
+ rxnarray = getinterfacereactioninds(domain1,domain2,reactions)
+ M,Nrp1,Nrp2 = getstoichmatrix(domain1,domain2,reactions)
+ Gpart = ArrayPartition(domain1.Gs,domain2.Gs)
+ dGrxns = -M*Gpart
+ electronchanges = [hasproperty(reaction, :electronchange) ? reaction.electronchange : 0.0 for reaction in reactions]
+ referencepotentials = [hasproperty(reaction.kinetics, :V0) ? reaction.kinetics.V0 : 0.0 for reaction in reactions]
+ if isa(domain1.phase, IdealSurface)
+ phi = domain1.phi !== nothing ? domain1.phi : phi
+ elseif isa(domain2.phase, IdealSurface)
+ phi = domain2.phi !== nothing ? domain2.phi : phi
+ end
+ dGrxns .+= electronchanges.*(phi.-referencepotentials).*F
+ kfs = getkf.(reactions,nothing,T,0.0,0.0,Ref([]),A,phi,dGrxns,0.0)
+ Kc = getKcs(domain1.phase,domain2.phase,T,Nrp1,Nrp2,dGrxns)
+ krevs = kfs./Kc
+ reversibility = Array{Bool,1}(getfield.(reactions,:reversible))
+ forwardability = Array{Bool,1}(getfield.(reactions,:forwardable))
+ if isa(reactions,Vector{Any})
+ reactions = convert(Vector{ElementaryReaction},reactions)
end
if isa(kfs, Vector{Any})
kfs = convert(Vector{Float64}, kfs)
@@ -209,7 +221,7 @@ function upgradekinetics(rxns, domain1, domain2)
@assert length(spc) == 1
kin = stickingcoefficient2arrhenius(rxn.kinetics, surfdomain.phase.sitedensity, length(rxn.reactants) - 1, spc[1].molecularweight)
newrxns[i] = ElementaryReaction(index=rxn.index, reactants=rxn.reactants, reactantinds=rxn.reactantinds, products=rxn.products,
- productinds=rxn.productinds, kinetics=kin, radicalchange=rxn.radicalchange, reversible=rxn.reversible, forwardable=rxn.forwardable, pairs=rxn.pairs)
+ productinds=rxn.productinds, kinetics=kin, electronchange=rxn.electronchange, radicalchange=rxn.radicalchange, reversible=rxn.reversible, forwardable=rxn.forwardable, pairs=rxn.pairs)
else
newrxns[i] = rxn
end
diff --git a/src/Parse.jl b/src/Parse.jl
index 3698081f3..375851712 100644
--- a/src/Parse.jl
+++ b/src/Parse.jl
@@ -168,7 +168,12 @@ function getatomdictsmiles(smiles)
mol.assign_representative_molecule()
getatomdictfromrmg(mol.mol_repr)
else
- getatomdictfromrdkit(Chem.AddHs(Chem.MolFromSmiles(smiles)))
+ try
+ return getatomdictfromrdkit(Chem.AddHs(Chem.MolFromSmiles(smiles)))
+ catch e
+ println("RDKit parsing failed, using RMG instead", e)
+ return getatomdictfromrmg(molecule.Molecule().from_smiles(smiles))
+ end
end
end
diff --git a/src/Phase.jl b/src/Phase.jl
index 784ed5c80..1ef17426a 100644
--- a/src/Phase.jl
+++ b/src/Phase.jl
@@ -41,7 +41,7 @@ function IdealGas(species,reactions; name="",diffusionlimited=false)
vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
rxns = vcat(reactions[vecinds],reactions[otherrxninds])
rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,
- products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,
+ products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,
radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs,comment=rxn.comment) for (i,rxn) in enumerate(rxns)]
therm = getvecthermo(species)
rxnarray = getreactionindices(species,rxns)
@@ -80,7 +80,7 @@ function IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited
vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
rxns = vcat(reactions[vecinds],reactions[otherrxninds])
rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,
- products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,
+ products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,
radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs,comment=rxn.comment) for (i,rxn) in enumerate(rxns)]
therm = getvecthermo(species)
rxnarray = getreactionindices(species,rxns)
@@ -120,7 +120,7 @@ function IdealSurface(species,reactions,sitedensity;name="",diffusionlimited=fal
vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
rxns = vcat(reactions[vecinds],reactions[otherrxninds])
rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,
- products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,
+ products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,
radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs,comment=rxn.comment) for (i,rxn) in enumerate(rxns)]
therm = getvecthermo(species)
rxnarray = getreactionindices(species,rxns)
@@ -161,7 +161,7 @@ function FragmentBasedIdealFilm(species, reactions; name="", diffusionlimited=fa
vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
rxns = vcat(reactions[vecinds],reactions[otherrxninds])
rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,
- products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,
+ products=rxn.products,productinds=rxn.productinds,kinetics=rxn.kinetics,electronchange=rxn.electronchange,
radicalchange=rxn.radicalchange,reversible=rxn.reversible,forwardable=rxn.forwardable,pairs=rxn.pairs,
fragmentbasedreactants=rxn.fragmentbasedreactants,fragmentbasedreactantinds=rxn.fragmentbasedreactantinds,
fragmentbasedproducts=rxn.fragmentbasedproducts,fragmentbasedproductinds=rxn.fragmentbasedproductinds,
diff --git a/src/PhaseState.jl b/src/PhaseState.jl
index ff2683e8f..6abe48b73 100644
--- a/src/PhaseState.jl
+++ b/src/PhaseState.jl
@@ -4,6 +4,7 @@ using LinearAlgebra
using Tracker
using ReverseDiff
using RecursiveArrayTools
+using Logging
@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
return getGibbs.(getfield.(ph.species,:thermo),T)
@@ -50,15 +51,15 @@ end
export makespcsvector
-@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi)
+@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi,dGrxn,d)
if isdefined(rxn.kinetics,:efficiencies) && length(rxn.kinetics.efficiencies) > 0
@views @inbounds @fastmath C += sum([ns[i]*val for (i,val) in rxn.kinetics.efficiencies])/V
end
- return rxn.kinetics(T=T,P=P,C=C,phi=phi)
+ return rxn.kinetics(T=T,P=P,C=C,phi=phi,dGrxn=dGrxn,d=d)
end
export getkf
-@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
+@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi,dGrxns,d) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
kfs = similar(ns,length(ph.reactions))
i = 1
oldind = 1
@@ -70,7 +71,7 @@ export getkf
i += 1
end
@simd for i in ind+1:length(ph.reactions)
- @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi)
+ @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi,dGrxns[i],d)
end
return kfs
end
@@ -84,13 +85,13 @@ for 2 spc calculates using the Smolchowski equation
for >2 spc calculates using the Generalized Smolchowski equation
Equations from Flegg 2016
"""
-@inline function getDiffusiveRate(spcs::Q,diffs::Array{W,1}) where {Q<:AbstractArray,W<:Real}
+@inline function getDiffusiveRate(spcs::Q,spcsinds::Q2,diffs::Array{W,1}) where {Q<:AbstractArray,W<:Real,Q2<:AbstractArray}
if length(spcs) == 1
return Inf
elseif length(spcs) == 2
- @fastmath @inbounds kf = 4.0*Base.pi*(diffs[spcs[1].index]+diffs[spcs[2].index])*(spcs[1].radius+spcs[2].radius)*Na
+ @fastmath @inbounds kf = 4.0*Base.pi*(diffs[spcsinds[1]]+diffs[spcsinds[2]])*(spcs[1].radius+spcs[2].radius)*Na
else
- @views @inbounds diffusivity = diffs[getfield.(spcs,:index)]
+ @views @inbounds diffusivity = diffs[spcsinds]
N = length(spcs)
@fastmath a = (3.0*length(spcs)-5.0)/2.0
@fastmath Dinv = 1.0./diffusivity
@@ -103,76 +104,26 @@ Equations from Flegg 2016
end
export getDiffusiveRate
-@inline function getKc(rxn::ElementaryReaction,ph::U,T::Z,Gs::Q,phi::V=0.0) where {U<:AbstractPhase,V,Q,Z<:Real}
+@inline function getKc(rxn::ElementaryReaction,ph::U,T::Z,dGrxn::Q,phi::V=0.0) where {U<:AbstractPhase,V,Q,Z<:Real}
Nreact = length(rxn.reactantinds)
Nprod = length(rxn.productinds)
- dGrxn = 0.0
- if Nreact == 1
- @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]
- elseif Nreact == 2
- @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]
- elseif Nreact == 3
- @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]
- elseif Nreact == 4
- @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]+Gs[rxn.reactantinds[4]]
- end
- if Nprod == 1
- @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]
- elseif Nprod == 2
- @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]
- elseif Nprod == 3
- @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]
- elseif Nprod == 4
- @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]+Gs[rxn.productinds[4]]
- end
- return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*(getC0(ph,T))^(Nprod-Nreact)
+ return @inbounds @fastmath exp((dGrxn+rxn.electronchange*(phi*F))/(-R*T))*(getC0(ph,T))^(Nprod-Nreact)
end
-@inline function getKc(rxn::ElementaryReaction,phase1,phase2,Gs1,Gs2,T,phi=0.0) #for constant k interfaces
- dGrxn = 0.0
- dN1 = 0
- dN2 = 0
- for r in rxn.reactants
- isfirst = true
- ind = findfirst(isequal(r),phase1.species)
- if ind === nothing
- isfirst = false
- ind = findfirst(isequal(r),phase2.species)
- dGrxn -= Gs2[ind]
- dN2 -= 1
- else
- dGrxn -= Gs1[ind]
- dN1 -= 1
- end
- end
- for r in rxn.products
- isfirst = true
- ind = findfirst(isequal(r),phase1.species)
- if ind === nothing
- isfirst = false
- ind = findfirst(isequal(r),phase2.species)
- dGrxn += Gs2[ind]
- dN2 += 1
- else
- dGrxn += Gs1[ind]
- dN1 += 1
- end
- end
- return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*getC0(phase1,T)^dN1*getC0(phase2,T)^dN2
+@inline function getKcs(ph::U,T::Z,dGrxns::Q) where {U<:AbstractPhase,Q,Z<:Real}
+ return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp.*log(getC0(ph,T)));
end
-export getKc
-@inline function getKcs(ph::U,T::Z,Gs::Q) where {U<:AbstractPhase,Q,Z<:Real}
- return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
+@inline function getKcs(ph::U,T::Z,dGrxns::Q,phi::V) where {U<:AbstractPhase,Q,Z<:Real,V<:Real}
+ return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp.*log(getC0(ph,T)));
end
-@inline function getKcs(ph::U,T::Z,Gs::Q,phi::V) where {U<:AbstractPhase,Q,Z<:Real,V<:Real}
- return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)).+ph.electronchange.*(phi/(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
+@inline function getKcs(ph,T,Gs1,Gs2,dGrxns)
+ return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ ph.Nrp1.*log(getC0(ph.domain1.phase,T)) .+ ph.Nrp2.*log(getC0(ph.domain2.phase,T)));
end
-@inline function getKcs(ph,T,Gs1,Gs2)
- Gpart = ArrayPartition(Gs1,Gs2)
- return @fastmath @inbounds exp.(ph.stoichmatrix*(Gpart./(R*T)) .+ ph.Nrp1.*log(getC0(ph.domain1.phase,T)) .+ ph.Nrp2.*log(getC0(ph.domain2.phase,T)));
+@inline function getKcs(ph1,ph2,T,Nrp1,Nrp2,dGrxns)
+ return @fastmath @inbounds exp.(dGrxns./(-R*T) .+ Nrp1.*log(getC0(ph1,T)) .+ Nrp2.*log(getC0(ph2,T)));
end
export getKcs
@@ -181,30 +132,30 @@ export getKcs
Calculates the forward and reverse rate coefficients for a given reaction, phase and state
Maintains diffusion limitations if the phase has diffusionlimited=true
"""
-@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
- if signbit(kf)
+@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,dGrxn::Q2,diffs::Q3,V::W5,phi::W8,d;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
+ if signbit(kf)
if signbit(f)
- kf = getkf(rxn,ph,T,P,C,ns,V,phi)
+ kf = getkf(rxn,ph,T,P,C,ns,V,phi,dGrxn,d)
else
- kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
+ kf = getkf(rxn,ph,T,P,C,ns,V,phi,dGrxn,d)*f
end
end
- Kc = getKc(rxn,ph,T,Gs,phi)
+ Kc = getKc(rxn,ph,T,dGrxn,phi)
@fastmath krev = kf/Kc
if ph.diffusionlimited
if length(rxn.reactants) == 1
if length(rxn.products) > 1
- krevdiff = getDiffusiveRate(rxn.products,diffs)
+ krevdiff = getDiffusiveRate(rxn.products,rxn.productinds,diffs)
@fastmath krev = krev*krevdiff/(krev+krevdiff)
@fastmath kf = Kc*krev
end
elseif length(rxn.products) == 1
- kfdiff = getDiffusiveRate(rxn.reactants,diffs)
+ kfdiff = getDiffusiveRate(rxn.reactants,rxn.reactantinds,diffs)
@fastmath kf = kf*kfdiff/(kf+kfdiff)
@fastmath krev = kf/Kc
elseif length(rxn.products) == length(rxn.reactants)
- kfdiff = getDiffusiveRate(rxn.reactants,diffs)
- krevdiff = getDiffusiveRate(rxn.products,diffs)
+ kfdiff = getDiffusiveRate(rxn.reactants,rxn.reactantinds,diffs)
+ krevdiff = getDiffusiveRate(rxn.products,rxn.productinds,diffs)
@fastmath kff = kf*kfdiff/(kf+kfdiff)
@fastmath krevr = krev*krevdiff/(krev+krevdiff)
@fastmath kfr = Kc*krevr
@@ -223,32 +174,33 @@ Maintains diffusion limitations if the phase has diffusionlimited=true
end
export getkfkrev
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7,d;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
+ dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
if !phase.diffusionlimited && kfs === nothing
- kfs = getkfs(phase,T,P,C,ns,V,phi)
+ kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif !phase.diffusionlimited && !(kfs === nothing)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif phase.diffusionlimited && !(kfs === nothing)
len = length(phase.reactions)
krev = zeros(typeof(N),len)
@simd for i = 1:len
- @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+ @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
end
else
len = length(phase.reactions)
kfs = zeros(typeof(N),len)
krev = zeros(typeof(N),len)
@simd for i = 1:len
- @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+ @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
end
end
kfs .*= phase.forwardability
@@ -256,65 +208,67 @@ export getkfkrev
return kfs,krev
end
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7,d;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
+ dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
if !phase.diffusionlimited && kfs === nothing
- kfs = getkfs(phase,T,P,C,ns,V,phi)
+ kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif !phase.diffusionlimited && !(kfs === nothing)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif phase.diffusionlimited && !(kfs === nothing)
len = length(phase.reactions)
krev = similar(kfs)
kfsdiff = similar(kfs)
@simd for i = 1:len
- @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+ @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
end
return kfsdiff, krev
else
len = length(phase.reactions)
- kfs = zeros(typeof(Gs[1]),len)
+ kfs = zeros(typeof(Gs[1]),len)ss
krev = zeros(typeof(Gs[1]),len)
@simd for i = 1:len
- @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+ @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
end
end
return kfs,krev
end
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7,d;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
+ dGrxns = -(phase.stoichmatrix*Gs).+phase.electronchange.*(phi*F)
if !phase.diffusionlimited && kfs === nothing
- kfs = getkfs(phase,T,P,C,ns,V,phi)
+ kfs = getkfs(phase,T,P,C,ns,V,phi,dGrxns,d)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif !phase.diffusionlimited && !(kfs === nothing)
if phi == 0.0
- krev = @fastmath kfs./getKcs(phase,T,Gs)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns)
else
- krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
+ krev = @fastmath kfs./getKcs(phase,T,dGrxns,phi)
end
elseif phase.diffusionlimited && !(kfs === nothing)
len = length(phase.reactions)
krev = similar(kfs)
@simd for i = 1:len
- @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+ @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d;kf=kfs[i])
end
else
len = length(phase.reactions)
kfs = zeros(typeof(Gs[1]),len)
krev = zeros(typeof(Gs[1]),len)
@simd for i = 1:len
- @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+ @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,dGrxns[i],diffs,V,phi,d)
end
end
kfs .*= phase.forwardability
diff --git a/src/testing/ORR.rms b/src/testing/ORR.rms
index a402cb514..566fe3d59 100644
--- a/src/testing/ORR.rms
+++ b/src/testing/ORR.rms
@@ -218,4 +218,11 @@ Reactions:
radicalchange: 0
reactants: [H2O2A]
type: ElementaryReaction
+- kinetics: {A: 6.0e12, lmbd_o: 100000.0, lmbd_i_coefs: [ 2.63844404e+03 9.64948798e+00 7.41268237e-03 -6.10277107e-06], beta: 1.2e10 ,n: 1.0, d: 0.0, type: Marcus}
+ products: [HOOA]
+ radicalchange: 0
+ electronchange: -1
+ reversible: false
+ reactants: [O2A, H+]
+ type: ElementaryReaction
Units: {}