Skip to content

Deadend radical algorithm RMG-RMS twin PR #203

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
May 16, 2023
Merged

Deadend radical algorithm RMG-RMS twin PR #203

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
eb855f9
Calculate core species net consumption rates & core radical net termi…
hwpang Apr 3, 2023
74ab4ab
Calculate the ratio of net consumption/termination of core radicals c…
hwpang Apr 18, 2023
debeca4
feed in deadendradicaltolerance
hwpang Apr 3, 2023
1a761f0
output corespeciesnetconsumptionrates, coreradicalnetterminationrates
hwpang Apr 3, 2023
c58037a
use deadend radical algorithm to identify objects
hwpang Apr 14, 2023
df1224e
Add deadendradicaltolerance to test
hwpang Apr 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
156 changes: 150 additions & 6 deletions src/EdgeAnalysis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

"""
Expand Down Expand Up @@ -416,32 +440,44 @@ 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
@inbounds for i = 1:size(d.rxnarray)[2]
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
Expand Down Expand Up @@ -499,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
"""
Expand Down Expand Up @@ -551,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]
Expand All @@ -566,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]
Expand All @@ -584,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)
Expand Down Expand Up @@ -666,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]
Expand Down Expand Up @@ -797,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
Expand Down Expand Up @@ -831,7 +975,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,
Expand Down Expand Up @@ -876,7 +1020,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
Expand Down
2 changes: 1 addition & 1 deletion src/TestEdgeAnalysis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down