Merge pull request #181 from ReactionMechanismGenerator/fix_edge_analysis

Minor fixes for edge analysis

Author: mjohnson541

.DS_Store

@@ -234,13 +234,27 @@ Precalculate important indices and maps for use in edge analysis
 function getkeyselectioninds(coreedgedomain::AbstractDomain,coreedgeinters,domain,inters)
     corespcsinds = 1:length(domain.phase.species)
     edgespcsinds = length(domain.phase.species)+1:length(coreedgedomain.phase.species)
-    corerxninds = 1:length(domain.phase.reactions)
-    edgerxninds = length(domain.phase.reactions)+1:length(coreedgedomain.phase.reactions)
-    reactantindices = coreedgedomain.rxnarray[1:3,:]
-    productindices = coreedgedomain.rxnarray[4:6,:]
+    corerxninds = zeros(Int64,length(domain.phase.reactions))
+    corerxncount = 1
+    edgerxninds = zeros(Int64,length(coreedgedomain.phase.reactions))
+    edgerxncount = 1
+    coretoedgerxnmap = Dict{Int64,Int64}()
+    for (j,rxn) in enumerate(coreedgedomain.phase.reactions)
+        coreind = findfirst(x->rxn.reactants==x.reactants && rxn.products==x.products && rxn.kinetics==x.kinetics,domain.phase.reactions)
+        if coreind === nothing
+            edgerxninds[edgerxncount] = j
+            edgerxncount += 1
+        else
+            coretoedgerxnmap[coreind] = j
+            corerxninds[corerxncount] = j
+            corerxncount += 1
+        end
+    end
+
+    @views @inbounds reactantindices = coreedgedomain.rxnarray[1:3,:]
+    @views @inbounds productindices = coreedgedomain.rxnarray[4:6,:]
     coretoedgespcmap = Dict([i=>findfirst(isequal(spc),coreedgedomain.phase.species) for (i,spc) in enumerate(domain.phase.species)])
-    coretoedgerxnmap = Dict([i=>findfirst(isequal(rxn),coreedgedomain.phase.reactions) for (i,rxn) in enumerate(domain.phase.reactions)])
-    for j = 3:length(domain.indexes)
+    @simd for j = 3:length(domain.indexes)
         coretoedgespcmap[domain.indexes[j]] = coreedgedomain.indexes[j]
     end
     return corespcsinds,corerxninds,edgespcsinds,edgerxninds,reactantindices,productindices,coretoedgespcmap,coretoedgerxnmap
@@ -499,10 +513,11 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds,
     if charrate == 0 && length(edgereactionrates) > 0
         maxspeciesindex = argmax(edgespeciesrates)
         maxspeciesrate = edgespeciesrates[maxspeciesindex]
-        name = sim.names[maxspeciesindex]
+        index = edgespcsinds[maxspeciesindex]
+        name = sim.names[index]
         @info "at time $t s, species $name was added to model core to avoid singularity"
-        push!(invalidobjects,sim.species[maxspeciesindex])
-        return (false,true)
+        push!(invalidobjects,sim.species[index])
+        return (false,true,0.0)
     end
 
     if branchfactor != 0.0 && !firsttime