Use getphasespecies to generalize post processing codes

Author: hwpang

src/AutomaticMechanismAnalysis.jl

@@ -130,7 +130,7 @@ function follow(sim,rxnind,spcind,ropp,ropl,rts,forward;steptol=1e-2,branchtol=5
             return rpout
         else
             rpout = newrps[argmax([getsimilarity(sim.species[rpout.spcsinds[end]],
-                            sim.domain.phase.species[nrp.spcsinds[end]]) for nrp in newrps])]
+                            getphasespecies(sim.domain.phase)[nrp.spcsinds[end]]) for nrp in newrps])]
         end
     end
 end
@@ -361,7 +361,7 @@ function gettimescaleconnectivity(ssys::SystemSimulation,t)
     tsc = ones(length(ssys.species),length(ssys.species))*Inf
     domains = getfield.(ssys.sims,:domain)
     Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
-    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
+    Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
     cstot = zeros(Nspcs)
     vns = Array{Any,1}(undef,length(domains))
     vcs = Array{Any,1}(undef,length(domains))

src/Debugging.jl

@@ -129,7 +129,7 @@ function analyzecrash(ssys::SystemSimulation;tol=1e6)
             end
         end
         rdindex += length(sim.domain.phase.reactions)
-        for (i,dy) in enumerate(dydt[sdindex:sdindex-1+length(sim.domain.phase.species)])
+        for (i,dy) in enumerate(dydt[sdindex:sdindex-1+length(getphasespecies(sim.domain.phase))])
             if i > length(sim.species)
                 continue
             end
@@ -140,7 +140,7 @@ function analyzecrash(ssys::SystemSimulation;tol=1e6)
                 push!(spcs,getdebugspecies(sim.species[i],T;dy=dy,ratio=ratio,tol=tol,index=i))
             end
         end
-        sdindex += length(sim.domain.phase.species)
+        sdindex += length(getphasespecies(sim.domain.phase))
     end
     return DebugMech(spcs,rxns)
 end

src/Domain.jl

@@ -3304,7 +3304,7 @@ export getreactionindices
     senstooriginspcind = Array{Int64,1}(undef,length(sensspcinds))
     senstooriginrxnind = Array{Int64,1}(undef,length(sensrxninds))
     for (i,spcind) in enumerate(sensspcinds)
-        spc = domain.phase.species[spcind]
+        spc = getphasespecies(domain.phase)[spcind]
         sensspcnames[i] = spc.name
         @inbounds sensspcs[i] = Species(
             name=spc.name,

src/Plotting.jl

@@ -12,9 +12,9 @@ function plotmolefractions(bsol::Q, tf::V; t0::Z=1e-15,N::Z2=1000,tol::Z3=0.01,e
     maxes = maximum(xs,dims=2)
     spnames = []
     for i = 1:length(maxes)
-        if maxes[i] > tol && !(bsol.domain.phase.species[i].name in exclude)
+        if maxes[i] > tol && !(getphasespecies(bsol.domain.phase)[i].name in exclude)
             plot(ts,xs[i,:])
-            push!(spnames,bsol.domain.phase.species[i].name)
+            push!(spnames,getphasespecies(bsol.domain.phase)[i].name)
         end
     end
     legend(spnames)
@@ -32,9 +32,9 @@ function plotmolefractions(bsol::Q; tol::V=0.01, exclude::M=Array{String,1}()) w
     maxes = maximum(xs,dims=2)
     spnames = []
     for i = 1:length(maxes)
-        if maxes[i] > tol && !(bsol.domain.phase.species[i].name in exclude)
+        if maxes[i] > tol && !(getphasespecies(bsol.domain.phase)[i].name in exclude)
             plot(bsol.sol.t,xs[i,:])
-            push!(spnames,bsol.domain.phase.species[i].name)
+            push!(spnames,getphasespecies(bsol.domain.phase)[i].name)
         end
     end
     legend(spnames)
@@ -63,7 +63,7 @@ function plotradicalmolefraction(bsol;t=0.0,ts=[])
         ts = range(0.0,stop=t, length=100);
     end
     radmf = zeros(length(ts))
-    rads = [x.radicalelectrons for x in bsol.domain.phase.species]
+    rads = [x.radicalelectrons for x in getphasespecies(bsol.domain.phase)]
     for (i,t) in enumerate(ts)
         radmf[i] = dot(rads,molefractions(bsol,t))
     end
@@ -75,7 +75,7 @@ end
 export plotradicalmolefraction
 
 function plotmaxthermoforwardsensitivity(bsol, spcname; N=0, tol= 1e-2)
-    spnames = getfield.(bsol.domain.phase.species,:name)
+    spnames = getfield.(getphasespecies(bsol.domain.phase),:name)
     values = Array{Float64,1}()
     outnames = Array{String,1}()
     for spn in spnames
@@ -132,7 +132,7 @@ associated with each reaction
 N reactions are included all of which must have absolute value greater than abs(maximum prod or loss rate)*tol
 """
 function plotrops(bsol::Y,name::X,t::Z;N=0,tol=0.01) where {Y<:Simulation, X<:AbstractString, Z<:Real}
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     end
     rop = rops(bsol,name,t)
@@ -195,7 +195,7 @@ reactions with maximum (over time) production value greater than max production*
 maximum (over time) loss value greater than maximum loss*tol are included
 """
 function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     end
     if length(rxnrates) == 0 || length(ts) == 0
@@ -243,7 +243,7 @@ reactions with maximum (over time) production value greater than max production*
 maximum (over time) loss value greater than maximum loss*tol are included
 """
 function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     end
     if length(rxnrates) == 0 || length(ts) == 0
@@ -291,7 +291,7 @@ reactions with maximum (over time) production value greater than max production*
 maximum (over time) loss value greater than maximum loss*tol are included
 """
 function plotrops(bsol::Y,name::X;rxnrates=Array{Float64,1}(),ts=Array{Float64,1}(),tol=0.05) where {Y<:Simulation, X<:AbstractString}
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     end
     if length(rxnrates) == 0 || length(ts) == 0
@@ -431,7 +431,7 @@ function plotthermoadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01) whe
     inds = inds[1:k]
     xs = Array{Float64,1}(1:length(inds))
     barh(xs,reverse(dpvals[inds]))
-    yticks(xs,reverse(getfield.(bsol.domain.phase.species[inds],:name)))
+    yticks(xs,reverse(getfield.(getphasespecies(bsol.domain.phase)[inds],:name)))
     if name in ["T","V"]
         xlabel("d$name/dG mol/kcal")
     else
@@ -445,7 +445,7 @@ function plotrateadjointsensitivities(bsol::Y,name::X,dps::Z;N=0,tol=0.01) where
     if !(name in ["T","V"] || name in bsol.names)
         error("Name $name not in domain")
     end
-    dpvals = dps[length(bsol.domain.phase.species)+1:end]
+    dpvals = dps[length(getphasespecies(bsol.domain.phase))+1:end]
     inds = reverse(sortperm(abs.(dpvals)))
     if N == 0
         N = length(inds)
@@ -492,7 +492,7 @@ end
 export plottimescales
 
 function plotrxntransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0,rxntol=1e-6)
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     elseif !isnothing(dSdt) && (sum(dim > 1 for dim in size(dSdt)) > 1 || maximum(size(dSdt)) != length(bsol.p))
         error("dSdt must be a vector of length number of parameters")
@@ -546,7 +546,7 @@ end
 export plotrxntransitorysensitivities
 
 function plotthermotransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,tol=1e-3,N=0)
-    if !(name in getfield.(bsol.domain.phase.species,:name))
+    if !(name in getfield.(getphasespecies(bsol.domain.phase),:name))
         error("Species $name not in domain")
     end
     ind = findfirst(isequal(name),bsol.names)
@@ -582,7 +582,7 @@ function plotthermotransitorysensitivities(bsol,name,t;dSdt=nothing,tau=nothing,
     inds = inds[1:k]
     xs = Array{Float64,1}(1:length(inds))
     barh(xs,reverse(dSdt[inds]))
-    yticks(xs,reverse(getfield.(bsol.domain.phase.species[inds],:name)))
+    yticks(xs,reverse(getfield.(getphasespecies(bsol.domain.phase)[inds],:name)))
     xlabel("Normalized Transitory Sensitivities for $name")
     return
 end

src/Reactor.jl

@@ -100,7 +100,7 @@ function Reactor(domains::T,y0s::W1,tspan::W2,interfaces::Z=Tuple(),ps::X=SciMLB
     n = Nvars
     k = 1
     for (j,domain) in enumerate(domains)
-        Nspcs = length(domain.phase.species)
+        Nspcs = length(getphasespecies(domain.phase))
         Ntherm = length(domain.indexes) - 2
         for i = 1:8, j = 1:size(domain.rxnarray)[2]
             if domain.rxnarray[i,j] != 0
@@ -827,12 +827,12 @@ end
 export jacobianpforwarddiff!
 
 function jacobianpforwarddiff(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
-    J = zeros(length(y),length(domain.phase.species)+length(domain.phase.reactions))
+    J = zeros(length(y),length(getphasespecies(domain.phase))+length(domain.phase.reactions))
     jacobianpforwarddiff!(J,y,p,t,domain,interfaces,colorvec)
 end
 
 function jacobianpforwarddiff(y::U,p::W,t::Z,domains::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:Tuple}
-    J = zeros(length(y),length(domains.phase.species)+length(domains.phase.reactions))
+    J = zeros(length(y),length(getphasespecies(domains.phase))+length(domains.phase.reactions))
     jacobianpforwarddiff!(J,y,p,t,domains,interfaces,colorvec)
 end
 
@@ -858,7 +858,7 @@ end
 export jacobiany
 
 function jacobianp(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
-    J = zeros(length(y),length(domain.phase.species)+length(domain.phase.reactions))
+    J = zeros(length(y),length(getphasespecies(domain.phase))+length(domain.phase.reactions))
     jacobianp!(J,y,p,t,domain,interfaces,colorvec)
     return J
 end

src/Simulation.jl

@@ -20,7 +20,7 @@ struct Simulation{Q<:AbstractODESolution,W<:AbstractDomain,M,L<:AbstractArray,G<
 end
 
 function Simulation(sol::Q,domain::W,interfaces=[],p=nothing) where {Q<:AbstractODESolution,W<:AbstractDomain}
-    names = getfield.(domain.phase.species,:name)
+    names = getfield.(getphasespecies(domain.phase),:name)
     Ns = sum(hcat(sol.interp.u...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
     if hasproperty(sol.interp,:du)
         Nderivs = sum(hcat(sol.interp.du...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
@@ -32,7 +32,7 @@ function Simulation(sol::Q,domain::W,interfaces=[],p=nothing) where {Q<:Abstract
     if p === nothing
         p = domain.p
     end
-    return Simulation(sol,domain,interfaces,names,F,Ns,domain.phase.species,domain.phase.reactions,p)
+    return Simulation(sol,domain,interfaces,names,F,Ns,getphasespecies(domain.phase),domain.phase.reactions,p)
 end
 
 function Simulation(sol::Q,domain::W,reducedmodelmappings::ReducedModelMappings,interfaces=[],p=nothing) where {Q<:AbstractODESolution,W<:AbstractDomain} 
@@ -68,8 +68,8 @@ function Simulation(sol::Q,domain::W,reducedmodelmappings::ReducedModelMappings,
     t = sol.t
     sol = SciMLBase.build_solution(sol.prob, sol.alg, t, u, retcode = :Success)
 
-    species_range = 1:length(domain.phase.species)
-    names = getfield.(domain.phase.species,:name)
+    species_range = 1:length(getphasespecies(domain.phase))
+    names = getfield.(getphasespecies(domain.phase),:name)
     Ns = sum(hcat(sol.interp.u...)[species_range,:],dims=1)
     if hasproperty(sol.interp,:du)
         Nderivs = sum(hcat(sol.interp.du...)[species_range,:],dims=1)
@@ -81,7 +81,7 @@ function Simulation(sol::Q,domain::W,reducedmodelmappings::ReducedModelMappings,
     if p === nothing
         p = domain.p
     end
-    return Simulation(sol,domain,interfaces,names,F,Ns,domain.phase.species,domain.phase.reactions,p)
+    return Simulation(sol,domain,interfaces,names,F,Ns,getphasespecies(domain.phase),domain.phase.reactions,p)
 end
 
 export Simulation
@@ -119,7 +119,7 @@ function SystemSimulation(sol,domains,interfaces,p)
     species = Array{Species,1}()
     for sim in sims
         append!(names,sim.names)
-        append!(species,sim.domain.phase.species)
+        append!(species,getphasespecies(sim.domain.phase))
         append!(reactions,sim.domain.phase.reactions)
     end
     for inter in interfaces
@@ -140,7 +140,7 @@ export iterate
 Broadcast.broadcastable(p::T) where {T<:AbstractSimulation} = Ref(p)
 export broadcastable
 
-spcindex(bsol::Z,name::Q) where {Z<:Simulation,Q<:AbstractString} = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
+spcindex(bsol::Z,name::Q) where {Z<:Simulation,Q<:AbstractString} = findfirst(isequal(name),getfield.(getphasespecies(bsol.domain.phase),:name))
 export spcindex
 
 function molefractions(bsol::Q,name::W,t::E) where {Q<:AbstractSimulation, W<:String, E<:Real}
@@ -234,7 +234,7 @@ this outputs a sparse matrix of  num reactions xnum species containing the produ
 rate of that species associated with that reaction
 """
 function rops(bsol::Q,t::X) where {Q<:Simulation,X<:Real}
-    ropmat = spzeros(length(bsol.domain.phase.reactions),length(bsol.domain.phase.species))
+    ropmat = spzeros(length(bsol.domain.phase.reactions),length(getphasespecies(bsol.domain.phase)))
     cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
     V = getdomainsize(bsol,t)
     @simd for i in 1:length(bsol.domain.phase.reactions)
@@ -253,7 +253,7 @@ end
 function rops(ssys::SystemSimulation,t)
     domains = getfield.(ssys.sims,:domain)
     Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
-    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
+    Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
     cstot = zeros(Nspcs)
     vns = Array{Any,1}(undef,length(domains))
     vcs = Array{Any,1}(undef,length(domains))
@@ -302,7 +302,7 @@ function rops(bsol::Y,name::X,t::Z) where {Y<:Simulation, X<:AbstractString, Z<:
     rop = spzeros(length(bsol.domain.phase.reactions))
     cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
     V = getdomainsize(bsol,t)
-    ind = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
+    ind = findfirst(isequal(name),getfield.(getphasespecies(bsol.domain.phase),:name))
     @assert !isa(ind,Nothing) "species $name not in species array"
     for (i,rxn) in enumerate(bsol.domain.phase.reactions)
         c = 0
@@ -320,7 +320,7 @@ function rops(ssys::SystemSimulation,name,t)
     domains = getfield.(ssys.sims,:domain)
     ind = findfirst(isequal(name),ssys.names)
     Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
-    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
+    Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
     cstot = zeros(Nspcs)
     vns = Array{Any,1}(undef,length(domains))
     vcs = Array{Any,1}(undef,length(domains))
@@ -609,7 +609,7 @@ function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::Strin
     @assert denominator in bsol.names
     indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
     inddeno = findfirst(isequal(denominator),bsol.names)+bsol.domain.parameterindexes[1]-1
-    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2 
+    Nvars = length(getphasespecies(bsol.domain.phase))+length(bsol.domain.indexes)-2 
     Nrxns = length(bsol.domain.phase.reactions)
     x,dp = extract_local_sensitivities(bsol.sol,t)
     s = dp[inddeno][indnum]
@@ -626,7 +626,7 @@ function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::Strin
     @assert denominator in bsol.names
     indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
     inddeno = findfirst(isequal(denominator),bsol.names)+bsol.domain.parameterindexes[1]-1
-    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2 
+    Nvars = length(getphasespecies(bsol.domain.phase))+length(bsol.domain.indexes)-2 
     Nrxns = length(bsol.domain.phase.reactions)
     x,dp = extract_local_sensitivities(bsol.sol,t)
     svals = dp[inddeno][bsol.domain.indexes[1]:bsol.domain.indexes[2]]
@@ -645,7 +645,7 @@ function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::Strin
     @assert numerator in bsol.names
     indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
     inddeno = denominator+bsol.domain.parameterindexes[1]-1
-    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2 
+    Nvars = length(getphasespecies(bsol.domain.phase))+length(bsol.domain.indexes)-2 
     Nrxns = length(bsol.domain.phase.reactions)
     x,dp = extract_local_sensitivities(bsol.sol,t)
     s = dp[inddeno+length(bsol.domain.phase.species)][indnum]
@@ -665,7 +665,7 @@ function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::Strin
     @assert numerator in bsol.names
     indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
     inddeno = denominator+bsol.domain.parameterindexes[1]-1
-    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2 
+    Nvars = length(getphasespecies(bsol.domain.phase))+length(bsol.domain.indexes)-2 
     Nrxns = length(bsol.domain.phase.reactions)
     x,dp = extract_local_sensitivities(bsol.sol,t)
     svals = dp[inddeno+length(bsol.domain.phase.species)][bsol.domain.indexes[1]:bsol.domain.indexes[2]]
@@ -713,7 +713,7 @@ function rates(ssys::Q,t::X) where {Q<:SystemSimulation,X<:Real}
     rts = zeros(length(ssys.reactions))
     domains = getfield.(ssys.sims,:domain)
     Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([hasproperty(inter,:reactions) ? length(inter.reactions) : 0 for inter in ssys.interfaces])
-    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
+    Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
     cstot = zeros(Nspcs)
     vns = Array{Any,1}(undef,length(domains))
     vcs = Array{Any,1}(undef,length(domains))

src/fluxdiagrams.jl

@@ -133,7 +133,7 @@ function makefluxdiagrams(bsol,ts;centralspecieslist=Array{String,1}(),superimpo
     colorscheme="viridis",removeunconnectednodes=false)
 
     if hasproperty(bsol,:domain)
-        specieslist = bsol.domain.phase.species
+        specieslist = getphasespecies(bsol.domain.phase)
         reactionlist = bsol.domain.phase.reactions
     else 
         specieslist = bsol.species