Merge pull request #149 from ReactionMechanismGenerator/edge_analysis_improvements

Edge analysis improvements

Author: mjohnson541

src/Domain.jl

@@ -96,8 +96,8 @@ function ConstantTPDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies::A
         jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
     end
     rxnarray = getreactionindices(phase)
-    return ConstantTPDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-        T,P,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,false,MVector(false),MVector(0.0),p, Dict("V"=>phase.species[end].index+1)), y0, p
+    return ConstantTPDomain(phase,[1,length(phase.species),length(phase.species)+1],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+        T,P,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,false,MVector(false),MVector(0.0),p, Dict("V"=>length(phase.species)+1)), y0, p
 end
 export ConstantTPDomain
 
@@ -164,8 +164,8 @@ function ConstantVDomain(;phase::Z,initialconds::Dict{X,E},constantspecies::Arra
         jacobian=zeros(typeof(T),length(phase.species)+2,length(phase.species)+2)
     end
     rxnarray = getreactionindices(phase)
-    return ConstantVDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1,phase.species[end].index+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-    V,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>phase.species[end].index+1,"P"=>phase.species[end].index+2)), y0, p
+    return ConstantVDomain(phase,[1,length(phase.species),length(phase.species)+1,length(phase.species)+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    V,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>length(phase.species)+1,"P"=>length(phase.species)+2)), y0, p
 end
 export ConstantVDomain
 
@@ -232,8 +232,8 @@ function ConstantPDomain(;phase::Z,initialconds::Dict{X,E},constantspecies::Arra
         jacobian=zeros(typeof(T),length(phase.species)+2,length(phase.species)+2)
     end
     rxnarray = getreactionindices(phase)
-    return ConstantPDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1,phase.species[end].index+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-    P,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>phase.species[end].index+1,"V"=>phase.species[end].index+2)), y0, p
+    return ConstantPDomain(phase,[1,length(phase.species),length(phase.species)+1,length(phase.species)+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    P,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>length(phase.species)+1,"V"=>length(phase.species)+2)), y0, p
 end
 export ConstantPDomain
 
@@ -297,8 +297,8 @@ function ParametrizedTPDomain(;phase::Z,initialconds::Dict{X,Any},constantspecie
     N = sum(ns)
     V = N*R*Tfcn(0.0)/Pfcn(0.0)
     y0 = zeros(length(phase.species)+1)
-    y0[phase.species[1].index:phase.species[end].index] = ns
-    y0[phase.species[end].index+1] = V
+    y0[1:length(phase.species)] = ns
+    y0[length(phase.species)+1] = V
     p = vcat(zeros(length(phase.species)),ones(length(phase.reactions)))
     if length(constantspecies) > 0
         spcnames = getfield.(phase.species,:name)
@@ -313,8 +313,8 @@ function ParametrizedTPDomain(;phase::Z,initialconds::Dict{X,Any},constantspecie
         jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
     end
     rxnarray = getreactionindices(phase)
-    return ParametrizedTPDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-    Tfcn,Pfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("V"=>phase.species[end].index+1)), y0, p
+    return ParametrizedTPDomain(phase,[1,length(phase.species),length(phase.species)+1],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    Tfcn,Pfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("V"=>length(phase.species)+1)), y0, p
 end
 export ParametrizedTPDomain
 
@@ -390,8 +390,8 @@ function ParametrizedVDomain(;phase::Z,initialconds::Dict{X,Any},constantspecies
         jacobian=zeros(typeof(T),length(phase.species)+2,length(phase.species)+2)
     end
     rxnarray = getreactionindices(phase)
-    return ParametrizedVDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1,phase.species[end].index+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-    Vfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>phase.species[end].index+1,"P"=>phase.species[end].index+2)), y0, p
+    return ParametrizedVDomain(phase,[1,length(phase.species),length(phase.species)+1,length(phase.species)+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    Vfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>length(phase.species)+1,"P"=>length(phase.species)+2)), y0, p
 end
 export ParametrizedVDomain
 
@@ -467,8 +467,8 @@ function ParametrizedPDomain(;phase::Z,initialconds::Dict{X,Any},constantspecies
         jacobian=zeros(typeof(T),length(phase.species)+2,length(phase.species)+2)
     end
     rxnarray = getreactionindices(phase)
-    return ParametrizedPDomain(phase,[phase.species[1].index,phase.species[end].index,phase.species[end].index+1,phase.species[end].index+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
-    Pfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>phase.species[end].index+1,"V"=>phase.species[end].index+2)), y0, p
+    return ParametrizedPDomain(phase,[1,length(phase.species),length(phase.species)+1,length(phase.species)+2],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    Pfcn,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict("T"=>length(phase.species)+1,"V"=>length(phase.species)+2)), y0, p
 end
 export ParametrizedPDomain
 
@@ -550,7 +550,7 @@ function ConstantTVDomain(;phase::Z,initialconds::Dict{X,E},constantspecies::Arr
         jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
     end
     rxnarray = getreactionindices(phase)
-    return ConstantTVDomain(phase,[phase.species[1].index,phase.species[end].index],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    return ConstantTVDomain(phase,[1,length(phase.species)],[1,length(phase.species)+length(phase.reactions)],constspcinds,
         T,V,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
@@ -605,7 +605,7 @@ function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,initialconds::
     end
     N = sum(ns)
     y0 = zeros(length(phase.species))
-    y0[phase.species[1].index:phase.species[end].index] = ns
+    y0[1:length(phase.species)] = ns
     p = vcat(zeros(length(phase.species)),ones(length(phase.reactions)))
     if length(constantspecies) > 0
         spcnames = getfield.(phase.species,:name)
@@ -620,7 +620,7 @@ function ParametrizedTConstantVDomain(;phase::IdealDiluteSolution,initialconds::
         jacobian=zeros(typeof(V),length(phase.species)+1,length(phase.species)+1)
     end
     rxnarray = getreactionindices(phase)
-    return ParametrizedTConstantVDomain(phase,[phase.species[1].index,phase.species[end].index],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    return ParametrizedTConstantVDomain(phase,[1,length(phase.species)],[1,length(phase.species)+length(phase.reactions)],constspcinds,
     Tfcn,V,efficiencyinds,rxnarray,jacobian,sensitivity,MVector(false),MVector(0.0),p,Dict{String,Int64}()), y0, p
 end
 export ParametrizedTConstantVDomain
@@ -697,7 +697,7 @@ function ConstantTAPhiDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies
         jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
     end
     rxnarray = getreactionindices(phase)
-    return ConstantTAPhiDomain(phase,[phase.species[1].index,phase.species[end].index],[1,length(phase.species)+length(phase.reactions)],constspcinds,
+    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
 end
 export ConstantTAPhiDomain

src/EdgeAnalysis.jl

@@ -482,6 +482,7 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
     numcorerxns = length(corerxninds)
     invalidobjectsprintboolean = true
     terminated = false
+    conversion = 0.0
 
     (dydt,rts,frts,rrts,cs,corespeciesratse,charrate,edgespeciesrates,
     edgereactionrates,corespeciesrateratios,edgespeciesrateratios,
@@ -596,10 +597,10 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
             invalidobjectsprintboolean = false
         end
         num = maxnumobjsperiter - length(invalidobjects)
-        newobjects = newobjects[:num]
-        newobjectinds = newobjectinds[:num]
-        newobjectvals = newobjectvals[:num]
-        newobjecttype = newobjecttype[:num]
+        newobjects = newobjects[1:num]
+        newobjectinds = newobjectinds[1:num]
+        newobjectvals = newobjectvals[1:num]
+        newobjecttype = newobjecttype[1:num]
     end
 
     if terminateatmaxobjects && length(invalidobjects) + length(newobjects) >= maxnumobjsperiter
@@ -640,7 +641,7 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
                 @info "At time $t sec, reached target termination rate ratio $ratio"
             end
         else isa(term, TerminationConversion)
-            index = findfirst(isequal(term.species),sim.species)
+            index = findfirst(isequal(term.species.name),sim.names)
             conversion = 1 - (y[index] / y0[index])
             name = sim.species[index].name
             if conversion >= term.conversion
@@ -653,15 +654,15 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
     if terminated
         for term in termination
             if isa(term, TerminationConversion)
-                index = findfirst(isequal(term.species),sim.species)
+                index = findfirst(isequal(term.species.name),sim.names)
                 conversion = 1 - (y[index] / y0[index])
                 name = sim.species[index].name
                 @info "$name conversion: $conversion"
             end
         end
     end
 
-    return (terminated,interrupt) 
+    return (terminated,interrupt,conversion) 
 end
 
 export identifyobjects!
@@ -670,7 +671,7 @@ export identifyobjects!
 run edge analysis to determine objects (species/reactions) that should be added to model core
 """
 function selectobjects(react,coreedgedomains,coreedgeinters,domains,inters,
-                p,tolmovetocore,tolinterruptsimulation,ignoreoverallfluxcriterion,filterreactions,
+                corep,coreedgep,tolmovetocore,tolinterruptsimulation,ignoreoverallfluxcriterion,filterreactions,
                 maxnumobjsperiter,tolbranchrxntocore,branchingratiomax,
                 branchingindex,terminateatmaxobjects,termination,
                 filterthreshold;
@@ -685,6 +686,8 @@ function selectobjects(react,coreedgedomains,coreedgeinters,domains,inters,
     maxedgespeciesrateratios = zeros(length(edgespcsinds))
     invalidobjects = []
     terminated = false
+    conversion = 0.0
+    code = :Success
 
     if tolbranchrxntocore != 0.0
         branchfactor = 1.0/tolbranchrxntocore
@@ -696,20 +699,21 @@ function selectobjects(react,coreedgedomains,coreedgeinters,domains,inters,
     inte = init(react.ode,solver,abstol=atol,reltol=rtol);
 
     t = inte.t
-    sim = getsim(inte,react,coreedgedomains,inters,p,coretoedgespcmap)
+    sim = getsim(inte,react,coreedgedomains,inters,corep,coretoedgespcmap)
 
     y0 = sim.sol[end]
     spcsaddindices = Array{Int64,1}()
     firsttime = true
 
     n = 1
-    while t < tf
+    while t < tf && code == :Success
         for i = 1:n
             step!(inte)
         end
+        code = check_error(inte)
         t = inte.t
-        sim = getsim(inte,react,coreedgedomains,inters,p,coretoedgespcmap)
-        terminated,interrupt = identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
+        sim = getsim(inte,react,coreedgedomains,inters,coreedgep,coretoedgespcmap)
+        terminated,interrupt,conversion = identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
             edgerxninds,reactantindices,productindices,unimolecularthreshold,bimolecularthreshold,
                 trimolecularthreshold,maxedgespeciesrateratios,tolmovetocore,tolinterruptsimulation,ignoreoverallfluxcriterion,filterreactions,
                 maxnumobjsperiter,branchfactor,branchingratiomax,
@@ -722,11 +726,24 @@ function selectobjects(react,coreedgedomains,coreedgeinters,domains,inters,
             break
         end
     end
-    return (terminated,invalidobjects,unimolecularthreshold,
-        bimolecularthreshold,trimolecularthreshold,maxedgespeciesrateratios)
+
+    if code == :Success
+        return (terminated,false,invalidobjects,unimolecularthreshold,
+            bimolecularthreshold,trimolecularthreshold,maxedgespeciesrateratios,t,conversion)
+    else
+        @error "Solver failed with code $code resurrecting job"
+        dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates,edgereactionrates,
+        corespeciesrateratios,edgespeciesrateratios,corereactionrates,corespeciesconcentrations,
+        corespeciesproductionrates,corespeciesconsumptionrates = processfluxes(sim,corespcsinds,corerxninds,edgespcsinds,edgerxninds)
+        ind = edgespcsinds[argmax(edgespeciesrates)]
+        invalidobjects = [sim.species[ind]]
+        return (terminated,true,invalidobjects,unimolecularthreshold,
+            bimolecularthreshold,trimolecularthreshold,maxedgespeciesrateratios,t,conversion)
+    end
+        
 end
 
-export selectspecies
+export selectobjects
 
 """
 calculate threshold rate constants for different numbers of reactants
@@ -746,4 +763,4 @@ function getthresholdrateconstants(sim::Simulation,phase::IdealDiluteSolution,fi
     return 2.08366122e10*T,22.2*T/mu,0.11*T/mu
 end
 
-export getthresholdrateeconstant
+export getthresholdrateconstants

src/Parse.jl

@@ -160,7 +160,7 @@ export getatomdictinchi
 getatomdictadjlist(adjlist) = getatomdictfromrmg(molecule.Molecule().from_adjacency_list(adjlist))
 export getatomdictadjlist
 
-function getspeciesradius(atomdict::Dict{String,Int64},nbonds::Int64)
+function getspeciesradius(atomdict,nbonds::Int64)
     """
     estimates the McGowan radius by calculating the McGowan Volume
     from the number of each type of atom and the number of bonds

src/Phase.jl

@@ -39,17 +39,17 @@ function IdealGas(species,reactions; name="",diffusionlimited=false)
     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,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
     therm = getvecthermo(species)
-    M,Nrp = getstoichmatrix(species,rxns)
+    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
     reversibility = getfield.(rxns,:reversible)
-    rxnarray = getreactionindices(species,rxns)
     return IdealGas(species=species,reactions=rxns,name=name,
-        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple, 
+        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple, 
         veckineticsinds=posinds, vecthermo=therm, otherreactions=otherrxns, electronchange=electronchange, 
         reversibility=reversibility,diffusionlimited=diffusionlimited,)
 end
@@ -79,17 +79,18 @@ function IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited
     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,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
     therm = getvecthermo(species)
-    M,Nrp = getstoichmatrix(species,rxns)
+    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
     reversibility = getfield.(rxns,:reversible)
-    rxnarray = getreactionindices(species,rxns)
+
     return IdealDiluteSolution(species=species,reactions=rxns,solvent=solvent,name=name,
-        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
+        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
         veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,electronchange=electronchange,
         reversibility=reversibility,diffusionlimited=diffusionlimited)
 end
@@ -119,17 +120,17 @@ function IdealSurface(species,reactions,sitedensity;name="",diffusionlimited=fal
     rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
         productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,electronchange=rxn.electronchange,reversible=rxn.reversible,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
     therm = getvecthermo(species)
-    M,Nrp = getstoichmatrix(species,rxns)
+    rxnarray = getreactionindices(species,rxns)
+    M,Nrp = getstoichmatrix(rxnarray,species)
     echangevec = getfield.(rxns,:electronchange).*F
     if all(echangevec .== 0)
         electronchange = nothing
     else 
         electronchange = convert(typeof(Nrp),echangevec)
     end
     reversibility = getfield.(rxns,:reversible)
-    rxnarray = getreactionindices(species,rxns)
     return IdealSurface(species=species,reactions=rxns,name=name,
-        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
+        spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple,
         veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,electronchange=electronchange,
         reversibility=reversibility,sitedensity=sitedensity,diffusionlimited=diffusionlimited)
 end
@@ -138,17 +139,23 @@ export IdealSurface
 """
 construct the stochiometric matrix for the reactions and the reaction molecule # change
 """
-function getstoichmatrix(spcs,rxns)
-    M = spzeros(length(rxns),length(spcs))
-    Nrp = zeros(length(rxns))
-    for (i,rxn) in enumerate(rxns)
-        Nrp[i] = Float64(length(rxn.productinds) - length(rxn.reactantinds))
-        for ind in rxn.reactantinds
-            M[i,ind] += 1.0
-        end
-        for ind in rxn.productinds
-            M[i,ind] -= 1.0
+function getstoichmatrix(rxnarray,spcs)
+    M = spzeros(size(rxnarray)[2],length(spcs))
+    Nrp = zeros(size(rxnarray)[2])
+    for i in 1:size(rxnarray)[2]
+        n = 0.0
+        for (j,ind) in enumerate(rxnarray[:,i])
+            if ind == 0
+                continue
+            elseif j > 3
+                M[i,ind] -= 1.0
+                n += 1.0
+            else
+                M[i,ind] += 1.0
+                n -= 1.0
+            end
         end
+        Nrp[i] = n
     end
     return M,Nrp
 end