From 52a72f3c4c9e1b0929a2939a3f066b3671181178 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 3 Feb 2024 16:16:55 -0800 Subject: [PATCH 01/20] include electronchange in reaction reconstruction in phase --- src/Interface.jl | 2 +- src/Phase.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index 466022abc..04896ea4b 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -209,7 +209,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/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, From c8bb7871b9986729b1fa722ba7295ab17ff6a47f Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 3 Feb 2024 16:18:25 -0800 Subject: [PATCH 02/20] add Marcus kinetics type and allow dGrxn dependence of kinetics --- src/Calculators/Rate.jl | 46 +++++++++++++++++++++++++++----------- src/Calculators/Ratevec.jl | 8 +++---- 2 files changed, 37 insertions(+), 17 deletions(-) diff --git a/src/Calculators/Rate.jl b/src/Calculators/Rate.jl index c873c58a0..602b98b6b 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,8 +23,8 @@ 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 @@ -34,10 +34,30 @@ export StickingCoefficient q::B = 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)/(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)/(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 +65,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 +85,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 +99,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 +115,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 +126,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 +146,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 +215,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) 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) From aac5b928be5790583764562ba4a017e9bda8dc4a Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 3 Feb 2024 16:19:20 -0800 Subject: [PATCH 03/20] handle dGrxn and d dependence of reactions in interfaces and phase calculations --- src/Interface.jl | 62 +++++++++--------- src/PhaseState.jl | 156 ++++++++++++++++------------------------------ 2 files changed, 88 insertions(+), 130 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index 04896ea4b..b8dbfc963 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) + dGrxn = -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,20 @@ 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 + 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) 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 From f5f2772e625d97af3494a042e2d98b62e42d3668 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 3 Feb 2024 16:20:20 -0800 Subject: [PATCH 04/20] adapt Domains to dGrxn and d dependence of forward rate coefficient calculations --- src/Domain.jl | 243 ++++++++++++++++++++++++++------------------------ 1 file changed, 128 insertions(+), 115 deletions(-) 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} From e8997a7465f698fdc50f4f132e5a7fedca644389 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 3 Feb 2024 16:20:41 -0800 Subject: [PATCH 05/20] add a Marcus reaction to the ORR.rms mechanism for testing --- src/testing/ORR.rms | 7 +++++++ 1 file changed, 7 insertions(+) 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: {} From d84ef84d33f1c9a52986ef135ea3123d3d848f4f Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Fri, 29 Mar 2024 17:25:16 -0700 Subject: [PATCH 06/20] fix variable naminng --- src/Interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index b8dbfc963..52bdeee81 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -51,7 +51,7 @@ export ReactiveInternalInterface function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot::Array{Q,1}) where {Q} Gpart = ArrayPartition(Gs1,Gs2) - dGrxn = -ri.stoichmatrix*Gpart + 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 From 7de59809434ba5d859f074dc4e6bb50b74e99b9d Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sat, 21 Sep 2024 12:59:38 -0700 Subject: [PATCH 07/20] update documentation build action --- .github/workflows/documentation.yml | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index cc731b92a..7e03c2539 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -13,15 +13,26 @@ jobs: contents: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@latest + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 with: version: '1.10' - - name: Install dependencies + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - name: Build run: | current_path=${{ github.workspace }} export JULIA_CONDAPKG_ENV="$current_path/rms_env" - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.build("ReactionMechanismSimulator");' + julia -e 'using Pkg; Pkg.develop(Pkg.PackageSpec(path="../ReactionMechanismSimulator.jl/")); Pkg.build("ReactionMechanismSimulator");' - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token From de4d6badd1f868e982b866c02b07da7b8d5f4e00 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 22 Sep 2024 18:26:16 -0700 Subject: [PATCH 08/20] loosen solve tolerances --- src/TestReactors.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/TestReactors.jl b/src/TestReactors.jl index 2a089f6ae..f8fd94925 100644 --- a/src/TestReactors.jl +++ b/src/TestReactors.jl @@ -77,7 +77,7 @@ using SciMLSensitivity interfaces = [kLAkHCondensationEvaporationWithReservoir(domain, conds)] react = Reactor(domain, y0, (0.0, 140000.01), interfaces; p=p) - sol1 = solve(react.ode, react.recommendedsolver, abstol=1e-18, reltol=1e-6) + sol1 = solve(react.ode, react.recommendedsolver, abstol=1e-16, reltol=1e-6) phaseDict = readinput("../src/testing/TdependentkLAkH.rms") spcs = phaseDict["phase"]["Species"] @@ -91,7 +91,7 @@ using SciMLSensitivity interfaces = [kLAkHCondensationEvaporationWithReservoir(domain, conds)] react = Reactor(domain, y0, (0.0, 140000.01), interfaces; p=p) #Create the reactor object - sol2 = solve(react.ode, react.recommendedsolver, abstol=1e-18, reltol=1e-6) + sol2 = solve(react.ode, react.recommendedsolver, abstol=1e-16, reltol=1e-6) spcnames = getfield.(liq.species, :name) octaneind = findfirst(isequal("octane"), spcnames) From 571ff08dbeb85daad526de3cce72b9b8ca54b338 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 22 Sep 2024 21:03:37 -0700 Subject: [PATCH 09/20] fix getKc call in ReactiveInternalInterface --- src/Interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 52bdeee81..e56c6f39f 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -53,7 +53,7 @@ function getkfskrevs(ri::ReactiveInternalInterface,T1,T2,phi1,phi2,Gs1,Gs2,cstot 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) + Kc = getKcs(ri.domain1.phase,ri.domain2.phase,T1,ri.Nrp1,ri.Nrp2,dGrxns) krevs = kfs./Kc return kfs,krevs end From 8616dfc5169af78e06f6fc6f00a0c1e403a6567a Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 22 Sep 2024 21:05:30 -0700 Subject: [PATCH 10/20] add phi and d to newer domains --- src/Domain.jl | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/Domain.jl b/src/Domain.jl index a304b055b..2a8619253 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -740,6 +740,8 @@ mutable struct FragmentBasedConstantTrhoDomain{N<:AbstractPhase,S<:Integer,W<:Re constantspeciesinds::Array{S,1} T::Float64 rho::Float64 + phi::W + d::W A::Float64 kfs::Array{W,1} krevs::Array{W,1} @@ -766,6 +768,8 @@ function FragmentBasedConstantTrhoDomain(; phase::Z, initialconds::Dict{X,E}, co mass = 0.0 P = 1.0e8 A = 0.0 + phi = 0.0 #default 0 + d = 0.0 fragmentnames = getfield.(getphasespecies(phase), :name) @@ -781,6 +785,10 @@ function FragmentBasedConstantTrhoDomain(; phase::Z, initialconds::Dict{X,E}, co elseif key == "mass" mass = val y0[end] = val + elseif key == "Phi" + phi = val + elseif key == "d" + d = val else ind = findfirst(isequal(key), fragmentnames) @assert typeof(ind) <: Integer "$key not found in fragment list: $fragmentnames" @@ -829,7 +837,7 @@ function FragmentBasedConstantTrhoDomain(; phase::Z, initialconds::Dict{X,E}, co end rxnarray = getreactionindices(phase) return FragmentBasedConstantTrhoDomain(phase, [1, length(fragmentnames), length(fragmentnames) + 1], [1, length(phase.species) + length(phase.reactions)], constspcinds, - T, rho, A, kfs, krevs, kfsnondiff, efficiencyinds, Gs, rxnarray, mu, diffs, jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict("mass" => length(fragmentnames) + 1)), y0, p + T, rho, phi, d, A, kfs, krevs, kfsnondiff, efficiencyinds, Gs, rxnarray, mu, diffs, jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict("mass" => length(fragmentnames) + 1)), y0, p end export FragmentBasedConstantTrhoDomain @@ -842,6 +850,7 @@ mutable struct ConstantTLiqFilmDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Re T::W epsilon::W phi::W + d::W kfs::Array{W,1} krevs::Array{W,1} kfsnondiff::Array{W,1} @@ -866,6 +875,7 @@ function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspe V = 0.0 P = 1.0e8 phi = 0.0 + d = 0.0 epsilon = 0.0 y0 = zeros(length(phase.species) + 1) spnames = [x.name for x in phase.species] @@ -879,6 +889,10 @@ function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspe y0[end] = val elseif key == "epsilon" epsilon = 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" @@ -922,7 +936,7 @@ function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspe end rxnarray = getreactionindices(phase) return ConstantTLiqFilmDomain(phase, [1, length(phase.species), length(phase.species) + 1], [1, length(phase.species) + length(phase.reactions)], constspcinds, - T, epsilon, phi, kfs, krevs, kfsnondiff, efficiencyinds, Gs, rxnarray, mu, diffs, jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}(["V" => length(phase.species) + 1])), y0, p + T, epsilon, phi, d, kfs, krevs, kfsnondiff, efficiencyinds, Gs, rxnarray, mu, diffs, jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}(["V" => length(phase.species) + 1])), y0, p end export ConstantTLiqFilmDomain From ba271543c5f8973809f1302b26e9a4782db4ed38 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Sun, 22 Sep 2024 21:10:05 -0700 Subject: [PATCH 11/20] doc build tweak --- .github/workflows/documentation.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 7e03c2539..673daaabf 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -33,6 +33,7 @@ jobs: current_path=${{ github.workspace }} export JULIA_CONDAPKG_ENV="$current_path/rms_env" julia -e 'using Pkg; Pkg.develop(Pkg.PackageSpec(path="../ReactionMechanismSimulator.jl/")); Pkg.build("ReactionMechanismSimulator");' + julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate();' - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token From d3119bd6a2b3c92fc69b2a3775acfd62052c6099 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 18:21:06 -0700 Subject: [PATCH 12/20] phase electron change fixes --- src/Phase.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/Phase.jl b/src/Phase.jl index 1ef17426a..f4c955bcd 100644 --- a/src/Phase.jl +++ b/src/Phase.jl @@ -125,7 +125,7 @@ function IdealSurface(species,reactions,sitedensity;name="",diffusionlimited=fal therm = getvecthermo(species) rxnarray = getreactionindices(species,rxns) M,Nrp = getstoichmatrix(rxnarray,species) - echangevec = getfield.(rxns,:electronchange).*F + echangevec = getfield.(rxns,:electronchange) electronchange = convert(Array{Float64,1},echangevec) reversibility = getfield.(rxns,:reversible) forwardability = getfield.(rxns,:forwardable) @@ -170,11 +170,7 @@ function FragmentBasedIdealFilm(species, reactions; name="", diffusionlimited=fa rxnarray = getreactionindices(species,rxns) M,Nrp = getstoichmatrix(rxnarray,species) echangevec = getfield.(rxns,:electronchange) - if all(echangevec .== 0) - electronchange = nothing - else - electronchange = convert(echangevec,Array{Float64,1}) - end + electronchange = convert(Array{Float64,1},echangevec) reversibility = getfield.(rxns,:reversible) forwardability = getfield.(rxns,:forwardable) From 1b6a31be696992ce55536546f83c69a4710e6c1b Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 18:22:13 -0700 Subject: [PATCH 13/20] missing domain tweaks for phi/d --- src/Domain.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Domain.jl b/src/Domain.jl index 2a8619253..5ae98c1de 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -826,9 +826,9 @@ function FragmentBasedConstantTrhoDomain(; phase::Z, initialconds::Dict{X,E}, co end 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) + 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(getphasespecies(phase)), length(getphasespecies(phase))) @@ -926,8 +926,8 @@ function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspe 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) + kfs, krevs = getkfkrevs(phase, T, P, C, N, ns, Gs, diffs, V, phi, d) + kfsnondiff = getkfs(phase, T, P, C, ns, V, phi, d) p = vcat(Gs, kfsnondiff) if sparse jacobian = zeros(typeof(T), length(phase.species), length(phase.species)) From 1a8ddf43d219333e85789a92e7bb28727f9118bd Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 18:22:49 -0700 Subject: [PATCH 14/20] fix typo --- src/PhaseState.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PhaseState.jl b/src/PhaseState.jl index 6abe48b73..3fbe03963 100644 --- a/src/PhaseState.jl +++ b/src/PhaseState.jl @@ -233,7 +233,7 @@ end return kfsdiff, krev else len = length(phase.reactions) - kfs = zeros(typeof(Gs[1]),len)ss + 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,dGrxns[i],diffs,V,phi,d) From 12a65ce09a282b717c66a9404ef9a1c51178aa89 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 18:23:01 -0700 Subject: [PATCH 15/20] tweak docs build --- .github/workflows/documentation.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 673daaabf..35f0ce791 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -32,7 +32,6 @@ jobs: run: | current_path=${{ github.workspace }} export JULIA_CONDAPKG_ENV="$current_path/rms_env" - julia -e 'using Pkg; Pkg.develop(Pkg.PackageSpec(path="../ReactionMechanismSimulator.jl/")); Pkg.build("ReactionMechanismSimulator");' julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate();' - name: Build and deploy env: From a871c47f67a6e5041c46cf190853ce24f047c784 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 19:01:06 -0700 Subject: [PATCH 16/20] tweak docs --- .github/workflows/documentation.yml | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 35f0ce791..3f974185b 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -26,18 +26,8 @@ jobs: restore-keys: | ${{ runner.os }}-test-${{ env.cache-name }}- ${{ runner.os }}-test- - ${{ runner.os }}- - - uses: julia-actions/julia-buildpkg@v1 - - name: Build - run: | - current_path=${{ github.workspace }} - export JULIA_CONDAPKG_ENV="$current_path/rms_env" - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate();' - - name: Build and deploy - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - DOCUMENTER_DEBUG: true - run: | - current_path=${{ github.workspace }} - export JULIA_CONDAPKG_ENV="$current_path/rms_env" - julia --color=yes --project=docs docs/make.jl + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-docdeploy@latest + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} From a80848b572ad43c91c7f7c15f114a7e32196f59b Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 19:04:41 -0700 Subject: [PATCH 17/20] fix FragmentBasedReactiveFilmGrowthInterfaceConstantT interface --- src/Interface.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index e56c6f39f..563eedceb 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -456,9 +456,12 @@ function FragmentBasedReactiveFilmGrowthInterfaceConstantT(domainfilm, domain2, rxnarray, fragmentbasedrxnarray = getfragmentbasedinterfacereactioninds(domainfilm, domain2, reactions) - kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), 0.0, 0.0) - Kc = getKc.(reactions, domainfilm.phase, domain2.phase, Ref(domainfilm.Gs), Ref(domain2.Gs), T, 0.0) - krevs = kfs ./ Kc + M,Nrp1,Nrp2 = getstoichmatrix(domainfilm,domain2,reactions) + Gpart = ArrayPartition(domainfilm.Gs,domain2.Gs) + dGrxns = -M*Gpart + kfs = getkf.(reactions,nothing,T,0.0,0.0,Ref([]),0.0,0.0,dGrxns,0.0) + Kc = getKcs(domain1.phase,domain2.phase,T,Nrp1,Nrp2,dGrxns) + krevs = kfs./Kc M, Nrp1, Nrp2 = getstoichmatrix(domainfilm, domain2, reactions) reversibility = Array{Bool,1}(getfield.(reactions, :reversible)) From 8c1357d26dd85a64be055bd79480bcc07337ea17 Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 19:53:10 -0700 Subject: [PATCH 18/20] fix --- src/Interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 563eedceb..504231b73 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -460,7 +460,7 @@ function FragmentBasedReactiveFilmGrowthInterfaceConstantT(domainfilm, domain2, Gpart = ArrayPartition(domainfilm.Gs,domain2.Gs) dGrxns = -M*Gpart kfs = getkf.(reactions,nothing,T,0.0,0.0,Ref([]),0.0,0.0,dGrxns,0.0) - Kc = getKcs(domain1.phase,domain2.phase,T,Nrp1,Nrp2,dGrxns) + Kc = getKcs(domainfilm.phase,domain2.phase,T,Nrp1,Nrp2,dGrxns) krevs = kfs./Kc M, Nrp1, Nrp2 = getstoichmatrix(domainfilm, domain2, reactions) From dcb23c094db9b5eedff4ec6b8e823f11d9124f1b Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 19:53:30 -0700 Subject: [PATCH 19/20] revert documentation --- .github/workflows/documentation.yml | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 3f974185b..cc731b92a 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -13,21 +13,20 @@ jobs: contents: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@latest with: version: '1.10' - - uses: actions/cache@v1 + - name: Install dependencies + run: | + current_path=${{ github.workspace }} + export JULIA_CONDAPKG_ENV="$current_path/rms_env" + julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.build("ReactionMechanismSimulator");' + - name: Build and deploy env: - cache-name: cache-artifacts - with: - path: ~/.julia/artifacts - key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} - restore-keys: | - ${{ runner.os }}-test-${{ env.cache-name }}- - ${{ runner.os }}-test- - - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-docdeploy@latest - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_DEBUG: true + run: | + current_path=${{ github.workspace }} + export JULIA_CONDAPKG_ENV="$current_path/rms_env" + julia --color=yes --project=docs docs/make.jl From dd7856f60e92ea0c1c3aad94b1a58ffeebafc64b Mon Sep 17 00:00:00 2001 From: Matt Johnson Date: Tue, 24 Sep 2024 20:23:11 -0700 Subject: [PATCH 20/20] add documentation workflow --- .github/workflows/documentation.yml | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index cc731b92a..3ec7ad619 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -12,11 +12,21 @@ jobs: permissions: contents: write runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@latest + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 with: version: '1.10' + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 - name: Install dependencies run: | current_path=${{ github.workspace }}