Skip to content

Commit

Permalink
Split into two functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@hwpang
hwpang committed Jan 23, 2024
1 parent 7819ee5 commit 8f262ec
Showing 1 changed file with 56 additions and 5 deletions.
61 changes: 56 additions & 5 deletions src/Simulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,51 @@ function rops(bsol::Q, t::X) where {Q<:Simulation,X<:Real}
end

function rops(ssys::SystemSimulation, t)
domains = getfield.(ssys.sims, :domain)
Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims]) + sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter, :reactions)])
Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
cstot = zeros(Nspcs)
vns = Array{Any,1}(undef, length(domains))
vcs = Array{Any,1}(undef, length(domains))
vT = Array{Any,1}(undef, length(domains))
vP = Array{Any,1}(undef, length(domains))
vV = Array{Any,1}(undef, length(domains))
vC = Array{Any,1}(undef, length(domains))
vN = Array{Any,1}(undef, length(domains))
vmu = Array{Any,1}(undef, length(domains))
vkfs = Array{Any,1}(undef, length(domains))
vkrevs = Array{Any,1}(undef, length(domains))
vHs = Array{Any,1}(undef, length(domains))
vUs = Array{Any,1}(undef, length(domains))
vGs = Array{Any,1}(undef, length(domains))
vdiffs = Array{Any,1}(undef, length(domains))
vCvave = Array{Any,1}(undef, length(domains))
vphi = Array{Any,1}(undef, length(domains))
ropmat = spzeros(Nrxns, Nspcs)
start = 0
for (k, sim) in enumerate(ssys.sims)
vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
rops!(ropmat, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start)
start += length(vkfs[k])
end
for inter in ssys.interfaces
if inter isa FragmentBasedReactiveFilmGrowthInterfaceConstantT
kfs, krevs = getkfskrevs(inter)
rops!(ropmat, nothing, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], inter.domainfilm.indexes[3], start)
start += length(kfs)
elseif hasproperty(inter, :reactions)
kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
start += length(kfs)
end
end
return ropmat
end

function massrops(ssys::SystemSimulation, t)
@assert any([domain isa FragmentBasedConstantTrhoDomain for domain in domains]), "massrops only works for FragmentBasedConstantTrhoDomain"

domains = getfield.(ssys.sims, :domain)
Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims]) + sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter, :reactions)])
Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims])
Expand All @@ -276,6 +321,8 @@ function rops(ssys::SystemSimulation, t)
else
ropmat = spzeros(Nrxns, Nspcs)
end
massropvec = spzeros(Nrxns, 1)
ropmat = spzeros(Nrxns, Nspcs)
start = 0
for (k, sim) in enumerate(ssys.sims)
vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
Expand All @@ -286,15 +333,15 @@ function rops(ssys::SystemSimulation, t)
for inter in ssys.interfaces
if inter isa FragmentBasedReactiveFilmGrowthInterfaceConstantT
kfs, krevs = getkfskrevs(inter)
rops!(ropmat, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], inter.domainfilm.indexes[3], start)
rops!(ropmat, massropvec, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], inter.domainfilm.indexes[3], start)
start += length(kfs)
elseif hasproperty(inter, :reactions)
kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
start += length(kfs)
end
end
return ropmat
return ropmat, massropvec
end

"""
Expand Down Expand Up @@ -423,7 +470,7 @@ function rops!(ropvec, rarray::Array{Int64,2}, cs, kfs, krevs, V, start, ind)
end


function rops!(ropmat, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int64,2}, cs, kfs, krevs, V, Mws, fragmentindexes, massindex, start)
function rops!(ropmat, massropvec, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int64,2}, cs, kfs, krevs, V, Mws, fragmentindexes, massindex, start)
numfragmentbasedreacprod, numrxns = size(fragmentbasedrxnarray)
half = Int(numfragmentbasedreacprod / 2)
for i = 1:length(kfs)
Expand All @@ -433,7 +480,9 @@ function rops!(ropmat, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int6
if fragmentbasedrxnarray[j, i] != 0
@fastmath ropmat[i+start, fragmentbasedrxnarray[j, i]] -= R
if !(fragmentbasedrxnarray[j, i] in fragmentindexes)
ropmat[i+start, massindex] += R * Mws[fragmentbasedrxnarray[j, i]]
if !(massropvec == nothing)
@fastmath @inbounds massropvec[i+start] -= R * Mws[fragmentbasedrxnarray[j, i]]
end
end
end
end
Expand All @@ -442,7 +491,9 @@ function rops!(ropmat, rarray::Array{Int64,2}, fragmentbasedrxnarray::Array{Int6
if fragmentbasedrxnarray[j, i] != 0
@fastmath ropmat[i+start, fragmentbasedrxnarray[j, i]] += R
if !(fragmentbasedrxnarray[j, i] in fragmentindexes)
ropmat[i+start, massindex] -= R * Mws[fragmentbasedrxnarray[j, i]]
if !(massropvec == nothing)
@fastmath @inbounds massropvec[i+start] += R * Mws[fragmentbasedrxnarray[j, i]]
end
end
end
end
Expand Down

0 comments on commit 8f262ec

Please sign in to comment.