From eb855f97b2e2b7ba5f95c0dad655f15a234f6c73 Mon Sep 17 00:00:00 2001 From: Hao-Wei Pang Date: Mon, 3 Apr 2023 19:15:57 -0400 Subject: [PATCH 1/6] Calculate core species net consumption rates & core radical net termination rates --- src/EdgeAnalysis.jl | 40 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index a919cfeb..f640db64 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -330,6 +330,10 @@ function processfluxes(sim::SystemSimulation, corespeciesconcentrations = cs[corespcsinds] corespeciesconsumptionrates = zeros(length(corespeciesconcentrations)) corespeciesproductionrates = zeros(length(corespeciesconcentrations)) + # record the net consumption rates for core adicals, defined as sum_rxn max(L_rxn-P_rxn,0.0) + corespeciesnetconsumptionrates = zeros(length(corespeciesconcentrations)) + # record the net termination rates for core radicals, defined as sum_rxn max(L_rxn-P_rxn,0.0) * min(radicalchange_rxn,0.0) + coreradicalnetterminationrates = zeros(length(corespeciesconcentrations)) #process core species consumption and production rates index = 1 @@ -338,18 +342,28 @@ function processfluxes(sim::SystemSimulation, if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end + net_forward_rate = max(frts[i+index]-rrts[i+index], 0.0) for j = 1:4 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index] @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_forward_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_forward_rate * abs(min(d.phase.reactions[i].radicalchange, 0.0)) + end else break end end + net_reverse_rate = max(rrts[i+index]-frts[i+index], 0.0) for j = 5:8 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index] @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_reverse_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_reverse_rate * abs(min(-d.phase.reactions[i].radicalchange, 0.0)) + end else break end @@ -363,18 +377,28 @@ function processfluxes(sim::SystemSimulation, if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end + net_forward_rate = max(frts[i+index]-rrts[i+index], 0.0) for j = 1:4 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index] @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_forward_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_forward_rate * abs(min(d.phase.reactions[i].radicalchange, 0.0)) + end else break end end + net_reverse_rate = max(rrts[i+index]-frts[i+index], 0.0) for j = 5:8 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index] @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_reverse_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_reverse_rate * abs(min(-d.phase.reactions[i].radicalchange, 0.0)) + end else break end @@ -384,7 +408,7 @@ function processfluxes(sim::SystemSimulation, end end - return dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates,edgereactionrates,edgerxnradrateratios,corespeciesrateratios,edgespeciesrateratios,corereactionrates,corespeciesconcentrations,corespeciesproductionrates,corespeciesconsumptionrates + return dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates,edgereactionrates,edgerxnradrateratios,corespeciesrateratios,edgespeciesrateratios,corereactionrates,corespeciesconcentrations,corespeciesproductionrates,corespeciesconsumptionrates,corespeciesnetconsumptionrates,coreradicalnetterminationrates end """ @@ -416,6 +440,8 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg @inbounds corespeciesconcentrations = cs[corespcsinds] corespeciesconsumptionrates = zeros(length(corespeciesconcentrations)) corespeciesproductionrates = zeros(length(corespeciesconcentrations)) + corespeciesnetconsumptionrates = zeros(length(corespeciesconcentrations)) + coreradicalnetterminationrates = zeros(length(corespeciesconcentrations)) #process core species consumption and production rates d = sim.domain @@ -423,25 +449,35 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg if @inbounds any(d.rxnarray[:,i].>length(corespeciesconcentrations)) continue end + net_forward_rate = max(frts[i]-rrts[i], 0.0) for j = 1:4 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i] @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_forward_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_forward_rate * abs(min(d.phase.reactions[i].radicalchange, 0.0)) + end else break end end + net_reverse_rate = max(rrts[i]-frts[i], 0.0) for j = 5:8 if @inbounds d.rxnarray[j,i] != 0 @inbounds corespeciesproductionrates[d.rxnarray[j,i]] += frts[i] @inbounds corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i] + corespeciesnetconsumptionrates[d.rxnarray[j,i]] += net_reverse_rate + if d.phase.species[d.rxnarray[j,i]].radicalelectrons == 1 + coreradicalnetterminationrates[d.rxnarray[j,i]] += net_reverse_rate * abs(min(-d.phase.reactions[i].radicalchange, 0.0)) + end else break end end end - return dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates,edgereactionrates,edgerxnradrateratios,corespeciesrateratios,edgespeciesrateratios,corereactionrates,corespeciesconcentrations,corespeciesproductionrates,corespeciesconsumptionrates + return dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates,edgereactionrates,edgerxnradrateratios,corespeciesrateratios,edgespeciesrateratios,corereactionrates,corespeciesconcentrations,corespeciesproductionrates,corespeciesconsumptionrates,corespeciesnetconsumptionrates,coreradicalnetterminationrates end export processfluxes From 74ab4abbc831b63184eb43d5f986ef2a04ef6467 Mon Sep 17 00:00:00 2001 From: Hao-Wei Pang Date: Tue, 18 Apr 2023 18:53:27 -0400 Subject: [PATCH 2/6] Calculate the ratio of net consumption/termination of core radicals contributed by each edge reaction compared to all core reactions --- src/EdgeAnalysis.jl | 68 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index f640db64..d9d79bf9 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -535,6 +535,74 @@ end export calcbranchingnumbers +""" +Calculate the ratio of net consumption/termination of core radicals contributed by each edge reaction compared to all core reactions +""" +function calcradconstermratios(sim,reactantinds,productinds,corespcsinds,corerxninds,edgerxninds,edgereactionrates,corespeciesnetconsumptionrates,coreradicalnetterminationrates) + radconstermratios = zeros(length(edgereactionrates)) + for ind in 1:length(edgereactionrates) + chain_termination = false + chain_transfer = false + + index = edgerxninds[ind] + reactionrate = edgereactionrates[ind] + + if reactionrate > 0 + reactantside = reactantinds[:,index] + productside = productinds[:,index] + else + reactantside = productinds[:,index] + productside = reactantinds[:,index] + end + + productrade = [sim.species[i].radicalelectrons for i in productside if i != 0] + reactantrade = [sim.species[i].radicalelectrons for i in reactantside if i != 0] + + if length(productrade) == 1 && length(reactantrade) == 2 + if productrade[1] == 0 && reactantrade[1] == 1 && reactantrade[2] == 1 + chain_termination = true + end + elseif length(reactantrade) == length(productrade) && length(productrade) == 2 + if (0 in reactantrade && 1 in reactantrade) && (0 in productrade && 1 in productrade) + chain_transfer = true + end + end + + if !(chain_transfer || chain_termination) + continue + end + + for spcindex in reactantside + if spcindex == 0 + continue + elseif spcindex in corespcsinds + if sim.species[spcindex].radicalelectrons != 1 + continue + end + + if chain_transfer + consumption = corespeciesnetconsumptionrates[spcindex] + consumption = max(consumption, 1e-40) # avoid divide by zeros + radconstermratio = abs(reactionrate) / consumption + end + + if chain_termination + termination = coreradicalnetterminationrates[spcindex] + termination = max(termination, 1e-40) # avoid divide by zeros + radconstermratio = abs(reactionrate) / termination + end + + if radconstermratio > radconstermratios[ind] + radconstermratios[ind] = radconstermratio + end + end + end + end + return radconstermratios +end + +export calcradconstermratios + """ determine species pairings that are concentrated enough that they should be reacted """ From debeca45387dca70cefd25d06f1ccd962647a120 Mon Sep 17 00:00:00 2001 From: Hao-Wei Pang Date: Mon, 3 Apr 2023 13:20:40 -0400 Subject: [PATCH 3/6] feed in deadendradicaltolerance add deadendradicaltolerance to select object --- src/EdgeAnalysis.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index d9d79bf9..b4c85fdb 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -655,7 +655,7 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds, trimolecularthreshold,maxedgespeciesrateratios,tolmovetocore,tolradmovetocore,tolinterruptsimulation, ignoreoverallfluxcriterion,filterreactions,maxnumobjsperiter,branchfactor,branchingratiomax, branchingindex,terminateatmaxobjects,termination,y0,invalidobjects,firsttime, - filterthreshold,transitorydict,checktransitory) + filterthreshold,transitorydict,checktransitory,deadendradicaltolerance) rxnarray = vcat() t = sim.sol.t[end] @@ -935,7 +935,7 @@ function selectobjects(react,edgereact,coreedgedomains,coreedgeinters,domains,in corep,coreedgep,tolmovetocore,tolradmovetocore,tolinterruptsimulation,ignoreoverallfluxcriterion,filterreactions, maxnumobjsperiter,tolbranchrxntocore,branchingratiomax, branchingindex,terminateatmaxobjects,termination, - filterthreshold,transitorydict,transitorystepperiod; + filterthreshold,transitorydict,transitorystepperiod,deadendradicaltolerance; atol=1e-20,rtol=1e-6,solver=CVODE_BDF()) (corespcsinds,corerxninds,edgespcsinds,edgerxninds,reactantindices, @@ -980,7 +980,7 @@ function selectobjects(react,edgereact,coreedgedomains,coreedgeinters,domains,in tolinterruptsimulation,ignoreoverallfluxcriterion,filterreactions, maxnumobjsperiter,branchfactor,branchingratiomax, branchingindex,terminateatmaxobjects,termination,y0,invalidobjects,firsttime, - filterthreshold,transitorydict,checktransitory) + filterthreshold,transitorydict,checktransitory,deadendradicaltolerance) if firsttime firsttime = false end From 1a761f08644bd0acb7d88c29f1edbe0af7481c51 Mon Sep 17 00:00:00 2001 From: Hao-Wei Pang Date: Mon, 3 Apr 2023 13:18:09 -0400 Subject: [PATCH 4/6] output corespeciesnetconsumptionrates, coreradicalnetterminationrates --- src/EdgeAnalysis.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index b4c85fdb..9c070a9c 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -670,7 +670,8 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds, (dydt,rts,frts,rrts,cs,corespeciesrates,charrate,edgespeciesrates, edgereactionrates,edgerxnradrateratios,corespeciesrateratios,edgespeciesrateratios, corereactionrates,corespeciesconcentrations, - corespeciesproductionrates,corespeciesconsumptionrates) = processfluxes(sim,corespcsinds,corerxninds,edgespcsinds,edgerxninds) + corespeciesproductionrates,corespeciesconsumptionrates, + corespeciesnetconsumptionrates, coreradicalnetterminationrates) = processfluxes(sim,corespcsinds,corerxninds,edgespcsinds,edgerxninds) for i = 1:length(edgespeciesrateratios) if @inbounds edgespeciesrateratios[i] > maxedgespeciesrateratios[i] From c58037a92bb9c10ce9ad91547f3589e134e16165 Mon Sep 17 00:00:00 2001 From: Hao-Wei Pang Date: Thu, 13 Apr 2023 20:22:57 -0400 Subject: [PATCH 5/6] use deadend radical algorithm to identify objects --- src/EdgeAnalysis.jl | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/EdgeAnalysis.jl b/src/EdgeAnalysis.jl index 9c070a9c..67b6f9e7 100644 --- a/src/EdgeAnalysis.jl +++ b/src/EdgeAnalysis.jl @@ -689,6 +689,10 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds, return (false,true,0.0) end + if deadendradicaltolerance != 0.0 && !firsttime + radconstermratios = calcradconstermratios(sim,reactantinds,productinds,corespcsinds,corerxninds,edgerxninds,edgereactionrates,corespeciesnetconsumptionrates,coreradicalnetterminationrates) + end + if branchfactor != 0.0 && !firsttime branchingnums = calcbranchingnumbers(sim,reactantinds,productinds,corespcsinds,corerxninds,edgerxninds,edgereactionrates, corespeciesrateratios,corespeciesconsumptionrates,branchfactor,branchingratiomax,branchingindex) @@ -771,6 +775,38 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds, tempnewobjecttype = [] end + # dead-end radical species selection algorithm should be evaluated before branching algorithm as + # branching algorithm tends to pick up propagation reactions first, and if these propagation reactions are added + # into the core for many iterations in a row, it can create dead-end radicals + if deadendradicaltolerance != 0.0 && !firsttime + for (i,ind) in enumerate(edgerxninds) + drr = radconstermratios[i] + if drr > deadendradicaltolerance + obj = sim.reactions[ind] + if !(obj in newobjects || obj in invalidobjects) + push!(tempnewobjects,obj) + push!(tempnewobjectinds,ind) + push!(tempnewobjectvals,drr) + push!(tempnewobjecttype,"deadendradical") + end + end + end + + sortedinds = reverse(sortperm(tempnewobjectvals)) + + for q in sortedinds + push!(newobjects,tempnewobjects[q]) + push!(newobjectinds,tempnewobjectinds[q]) + push!(newobjectvals,tempnewobjectvals[q]) + push!(newobjecttype,tempnewobjecttype[q]) + end + + tempnewobjects = [] + tempnewobjectinds = Array{Int64,1}() + tempnewobjectvals = Array{Float64,1}() + tempnewobjecttype = [] + end + if !ignoreoverallfluxcriterion #radical rate ratios for (i,ind) in enumerate(edgerxninds) rr = edgerxnradrateratios[i] @@ -902,6 +938,9 @@ function identifyobjects!(sim,corespcsinds,corerxninds,edgespcsinds, @inbounds spcname = transitoryoutdict[ind] @inbounds tol = transitorydict[spcname] @info "at time $t sec, reaction $rstr at a normalized transitory sensitivity from $spcname of $sens exceeded the threshold of $tol for moving to model core" + elseif newobjecttype[i] == "deadendradical" + tol = deadendradicaltolerance + @info "at time $t sec, reaction $rstr at a net radical consumption/termination ratio of $val exceeded the threshold of $tol for moving to model core" end end end From df1224e20e7d1933f5e24ca000366a5cc202ed00 Mon Sep 17 00:00:00 2001 From: hwpang Date: Mon, 3 Apr 2023 13:20:54 -0400 Subject: [PATCH 6/6] Add deadendradicaltolerance to test --- src/TestEdgeAnalysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TestEdgeAnalysis.jl b/src/TestEdgeAnalysis.jl index 0b8ca210..b02bde86 100644 --- a/src/TestEdgeAnalysis.jl +++ b/src/TestEdgeAnalysis.jl @@ -55,7 +55,7 @@ using Sundials reactedge = Reactor(coreedgedomain,coreedgey0,(0.0,1e6);p=coreedgep); (terminated,resurrected,invalidobjects,unimolecularthreshold,bimolecularthreshold, trimolecularthreshold,maxedgespeciesrateratios) = selectobjects(react,reactedge,coreedgedomain,[],coredomain, - [],corep,coreedgep,0.03,Inf,0.03,false,true,5,0.005,1.0,1.0,true,termination,1.0e8,Dict(),20) + [],corep,coreedgep,0.03,Inf,0.03,false,true,5,0.005,1.0,1.0,true,termination,1.0e8,Dict(),20,Inf) @test terminated == false @test invalidobjects[1].name == "[CH2]CCC" @test unimolecularthreshold[5] == true