Skip to content

Commit

Permalink
enable 4 party reaction for edge analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
@hwpang
hwpang committed Apr 19, 2022
1 parent 52d4b46 commit 1674cca
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 32 deletions.
24 changes: 12 additions & 12 deletions src/EdgeAnalysis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters)
rxnindexedge += length(coreedgedomains[i].phase.reactions)

indend = length(domains[i].phase.reactions)
reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:3,:]
productindices[:,ind:ind+indend-1] = domains[i].rxnarray[4:6,:]
reactantindices[:,ind:ind+indend-1] = domains[i].rxnarray[1:4,:]
productindices[:,ind:ind+indend-1] = domains[i].rxnarray[5:8,:]
ind += indend
end

Expand All @@ -220,8 +220,8 @@ function getkeyselectioninds(coreedgedomains,coreedgeinters,domains,inters)
index += length(coreedgeinters[i].phase.reactions)

indend = length(inters[i].reactions)
reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:3,:]
productindices[:,ind:ind+indend] = inters[i].rxnarray[4:6,:]
reactantindices[:,ind:ind+indend] = inters[i].rxnarray[1:4,:]
productindices[:,ind:ind+indend] = inters[i].rxnarray[5:8,:]
ind += indend
end
end
Expand Down Expand Up @@ -251,8 +251,8 @@ function getkeyselectioninds(coreedgedomain::AbstractDomain,coreedgeinters,domai
end
end

@views @inbounds reactantindices = coreedgedomain.rxnarray[1:3,:]
@views @inbounds productindices = coreedgedomain.rxnarray[4:6,:]
@views @inbounds reactantindices = coreedgedomain.rxnarray[1:4,:]
@views @inbounds productindices = coreedgedomain.rxnarray[5:8,:]
coretoedgespcmap = Dict([i=>findfirst(isequal(spc),coreedgedomain.phase.species) for (i,spc) in enumerate(domain.phase.species)])
@simd for j = 3:length(domain.indexes)
coretoedgespcmap[domain.indexes[j]] = coreedgedomain.indexes[j]
Expand Down Expand Up @@ -286,15 +286,15 @@ function processfluxes(sim::SystemSimulation,
if any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:3
for j = 1:4
if d.rxnarray[j,i] != 0
corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index]
corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index]
else
break
end
end
for j = 4:6
for j = 5:8
if d.rxnarray[j,i] != 0
corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index]
corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index]
Expand All @@ -311,15 +311,15 @@ function processfluxes(sim::SystemSimulation,
if any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:3
for j = 1:4
if d.rxnarray[j,i] != 0
corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i+index]
corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i+index]
else
break
end
end
for j = 4:6
for j = 5:8
if d.rxnarray[j,i] != 0
corespeciesproductionrates[d.rxnarray[j,i]] += frts[i+index]
corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i+index]
Expand Down Expand Up @@ -359,15 +359,15 @@ function processfluxes(sim::Simulation,corespcsinds,corerxninds,edgespcsinds,edg
if any(d.rxnarray[:,i].>length(corespeciesconcentrations))
continue
end
for j = 1:3
for j = 1:4
if d.rxnarray[j,i] != 0
corespeciesconsumptionrates[d.rxnarray[j,i]] += frts[i]
corespeciesproductionrates[d.rxnarray[j,i]] += rrts[i]
else
break
end
end
for j = 4:6
for j = 5:8
if d.rxnarray[j,i] != 0
corespeciesproductionrates[d.rxnarray[j,i]] += frts[i]
corespeciesconsumptionrates[d.rxnarray[j,i]] += rrts[i]
Expand Down
4 changes: 2 additions & 2 deletions src/Phase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@ Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
export broadcastable

function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
arr = zeros(Int64,(6,length(rxns)))
arr = zeros(Int64,(8,length(rxns)))
names = [spc.name for spc in spcs]
for (i,rxn) in enumerate(rxns)
inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants]
Expand All @@ -268,7 +268,7 @@ function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
end
for (j,spc) in enumerate(rxn.products)
ind = findfirst(isequal(spc),spcs)
arr[j+3,i] = ind
arr[j+4,i] = ind
rxn.productinds[j] = ind
end
if hasproperty(rxn.kinetics,:efficiencies) && length(rxn.kinetics.nameefficiencies) > 0
Expand Down
50 changes: 32 additions & 18 deletions src/Simulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,15 +292,19 @@ function rops!(ropmat,rarray,cs,kfs,krevs,V,start)
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
elseif @inbounds rarray[3,i] == 0
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
else
elseif @inbounds rarray[4,i] == 0
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
else
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]]
end
if @inbounds rarray[5,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
elseif @inbounds rarray[6,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
if @inbounds rarray[6,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]
elseif @inbounds rarray[7,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]
elseif rarray[8,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]
else
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]]
end
@fastmath R = (fR - rR)*V

Expand All @@ -309,35 +313,45 @@ function rops!(ropmat,rarray,cs,kfs,krevs,V,start)
@inbounds @fastmath ropmat[i+start,rarray[2,i]] -= R
if @inbounds rarray[3,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[3,i]] -= R
if @inbounds rarray[4,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[4,i]] -= R
end
end
end
@inbounds @fastmath ropmat[i+start,rarray[4,i]] += R
if @inbounds rarray[5,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[5,i]] += R
if @inbounds rarray[6,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[6,i]] += R
@inbounds @fastmath ropmat[i+start,rarray[5,i]] += R
if @inbounds rarray[6,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[6,i]] += R
if @inbounds rarray[7,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[7,i]] += R
if @inbounds rarray[8,i] != 0
@inbounds @fastmath ropmat[i+start,rarray[8,i]] += R
end
end
end
end
end

function rops!(ropvec,rarray,cs,kfs,krevs,V,start,ind)
for i = 1:length(kfs)
c = count(isequal(ind),rarray[4:6,i])-count(isequal(ind),rarray[1:3,i])
c = count(isequal(ind),rarray[5:8,i])-count(isequal(ind),rarray[1:4,i])
if c != 0.0
if @inbounds rarray[2,i] == 0
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
elseif @inbounds rarray[3,i] == 0
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
else
elseif @inbounds rarray[4,i] == 0
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
else
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]]
end
if @inbounds rarray[5,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
elseif @inbounds rarray[6,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
if @inbounds rarray[6,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]
elseif @inbounds rarray[7,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]
elseif rarray[8,i] == 0
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]
else
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]]
end
@fastmath R = (fR - rR)*V
@fastmath @inbounds ropvec[i+start] = c*R
Expand Down

0 comments on commit 1674cca

Please sign in to comment.