Skip to content
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

Jump to bottom

A new domain ConstantTVGasDomain #258

Open
wants to merge 10 commits into
base: main
Choose a base branch
from
Open

A new domain ConstantTVGasDomain #258

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
5bc511a
remove constantTVGasDomain
donerancl Apr 4, 2024
467c909
add constantTVGasDomain
donerancl Apr 29, 2024
d5b47f0
still does not work
donerancl May 7, 2024
88f9d42
clean up
donerancl May 7, 2024
c22d723
more cleanup
donerancl May 7, 2024
a32ad1d
added checks and warnings to VaporLiquidMassTransferInternalInterface…
donerancl May 13, 2024
683c51a
cleanup
donerancl May 13, 2024
39d5bdc
fixed rops bug
donerancl May 13, 2024
39319ca
fixed rops bug
donerancl May 13, 2024
512174c
Update Interface.jl
donerancl Jul 1, 2024
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
279 changes: 279 additions & 0 deletions src/Domain.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,95 @@ function ConstantTPDomain(; phase::E2, initialconds::Dict{X,X2}, constantspecies
end
export ConstantTPDomain


mutable struct ConstantTVGasDomain{N<:AbstractPhase,S<:Integer,W<:Real, W2<:Real, I<:Integer, Q<:AbstractArray} <: AbstractConstantKDomain
phase::N
indexes::Q #assumed to be in ascending order
parameterindexes::Q
constantspeciesinds::Array{S,1}
T::W
V::W
kfs::Array{W,1}
krevs::Array{W,1}
efficiencyinds::Array{I,1}
Gs::Array{W,1}
rxnarray::Array{Int64,2}
mu::W
diffusivity::Array{W,1}
jacobian::Array{W,2}
sensitivity::Bool
alternativepformat::Bool
jacuptodate::MArray{Tuple{1},Bool,1,1}
t::MArray{Tuple{1},W2,1,1}
p::Array{W,1}
thermovariabledict::Dict{String,Int64}
end

function ConstantTVGasDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies::Array{X3,1}=Array{String,1}(),
sparse::Bool=false,sensitivity::Bool=false) where {E<:Real,E2<:AbstractPhase,Q<:AbstractInterface,W<:Real,X,X2,X3}
#set conditions and initialconditions
T = 0.0
V = 0.0
y0 = zeros(length(phase.species)+1)
spnames = [x.name for x in phase.species]
for (key,val) in initialconds
if key == "T"
T = val
elseif key == "V"
V = val
else
ind = findfirst(isequal(key),spnames)
@assert typeof(ind)<: Integer "$key not found in species list: $spnames"
y0[ind] = val
end
end
@assert T != 0.0
@assert V != 0.0
ns = y0[1:end-1]
N = sum(ns)

if length(constantspecies) > 0
spcnames = getfield.(phase.species,:name)
constspcinds = [findfirst(isequal(k),spcnames) for k in constantspecies]
else
constspcinds = Array{Int64,1}()
end
efficiencyinds = [rxn.index for rxn in phase.reactions if typeof(rxn.kinetics)<:AbstractFalloffRate && length(rxn.kinetics.efficiencies) > 0]
Gs = calcgibbs(phase,T)
if :solvent in fieldnames(typeof(phase)) && typeof(phase.solvent) != EmptySolvent
mu = phase.solvent.mu(T)
else
mu = 0.0
end
if phase.diffusionlimited
diffs = [x(T=T,mu=mu,P=P) for x in getfield.(phase.species,:diffusion)]
else
diffs = Array{typeof(T),1}()
end
P = N*R*T/V
C = P/(R*T)

y0[end] = P
kfs,krevs = getkfkrevs(phase,T,P,C,N,ns,Gs,diffs,V,0.0)
kfsp = deepcopy(kfs)
for ind in efficiencyinds
kfsp[ind] = 1.0
end
p = vcat(deepcopy(Gs),kfsp)
if sparse
jacobian=spzeros(typeof(T),length(phase.species),length(phase.species))
else
jacobian=zeros(typeof(T),length(phase.species),length(phase.species))
end
rxnarray = getreactionindices(phase)
return ConstantTVGasDomain(phase,[1,length(phase.species),length(phase.species)+1],[1,length(phase.species)+length(phase.reactions)],constspcinds,
T,V,kfs,krevs,efficiencyinds,Gs,rxnarray,mu,diffs,jacobian,sensitivity,false,MVector(false),MVector(0.0),p, Dict("P"=>length(phase.species)+1)), y0, p
end
export ConstantTVGasDomain




struct ConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:Integer,Q<:AbstractArray} <: AbstractVariableKDomain
phase::N
indexes::Q #assumed to be in ascending order
Expand Down Expand Up @@ -1061,6 +1150,158 @@ end
return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
end


@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W3=SciMLBase.NullParameters()) where {W3<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:Array{Float64,1},Q} #no parameter input
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
for ind in d.efficiencyinds #efficiency related rates may have changed
d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, 0.0)
end
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:Array{Float64,1},Q<:Float64} #uses parameter input
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
if !d.alternativepformat
@views kfps = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
@views nothermochg= d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
else
kfps = d.p[length(d.phase.species)+1:end].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
nothermochg= d.Gs == d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
@views nokfchg = count(d.kfs .!= kfps) <= length(d.efficiencyinds) && all(kfps[d.efficiencyinds] .== 1.0)
if nothermochg && nokfchg
for ind in d.efficiencyinds #efficiency related rates may have changed
d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,0.0;f=kfps[ind])
end
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
elseif nothermochg
d.kfs = kfps
for ind in d.efficiencyinds #efficiency related rates may have changed
d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,0.0;f=kfps[ind])
end
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
else #need to handle thermo changes
d.kfs .= kfps
if !d.alternativepformat
d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
else
d.Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,0.0;kfs=d.kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
d.kfs[ind],d.krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,d.Gs,d.diffusivity,d.V,0.0;f=kfps[ind])
end
return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end
end

@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::Array{W3,1},t::Q,p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,W3<:ForwardDiff.Dual,Q} #Autodiff y
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)

cs = ns./d.V
C = N/d.V
if !d.alternativepformat
kfs = convert(typeof(y),p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
else
kfs = convert(typeof(y),d.p[length(d.phase.species)+1:end].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
krevs = convert(typeof(y),getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;kfs=kfs)[2])
for ind in d.efficiencyinds #efficiency related rates may have changed
kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;f=kfs[ind])
end
return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end

@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
if !d.alternativepformat
kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
else
kfs = d.p[length(d.phase.species)+1:end].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;kfs=kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;f=kfs[ind])
end
return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end

@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Autodiff p
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
kfs = similar(y,length(d.phase.reactions))
krevs = similar(y,length(d.phase.reactions))
if !d.alternativepformat
kfs .= p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
else
kfs .= d.p[length(d.phase.species)+1:end].*p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
end
krevs .= getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;kfs=kfs)[2]
for ind in d.efficiencyinds #efficiency related rates may have changed
kfs[ind],krevs[ind] = getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;f=kfs[ind])
end
return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end

@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J,Q} #Tracker/reversediff
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
if !d.alternativepformat
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0)[1]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
else
Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0)[1]*d.p[length(d.phase.species)+ind]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
end
krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;kfs=kfs)[2]
return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end

@inline function calcthermo(d::ConstantTVGasDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Tracker/reversediff
ns = y[d.indexes[1]:d.indexes[2]]
P = y[d.indexes[3]]
N = P * d.V / (R * d.T)
cs = ns./d.V
C = N/d.V
if !d.alternativepformat
Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0)[1]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
else
Gs = d.p[1:length(d.phase.species)].+p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind],d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0)[1]*d.p[length(d.phase.species)+ind]*p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
end
krevs = getkfkrevs(d.phase,d.T,P,C,N,ns,Gs,d.diffusivity,d.V,0.0;kfs=kfs)[2]
return ns,cs,d.T,P,d.V,C,N,d.mu,kfs,krevs,Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),Array{Float64,1}(),0.0,Array{Float64,1}(),0.0
end




@inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
ns = y[d.indexes[1]:d.indexes[2]]
T = y[d.indexes[3]]
Expand Down Expand Up @@ -1898,6 +2139,44 @@ end
end
end

@inline function calcdomainderivatives!(d::Q, dydt::Z7, interfaces::Z12; t::Z10, T::Z4, P::Z9, Us::Array{Z,1}, Hs::Array{Z11,1}, V::Z2, C::Z3, ns::Z5, N::Z6, Cvave::Z8) where {Q<:ConstantTVGasDomain,Z12,Z11,Z10,Z9,Z8<:Real,Z7,W<:IdealGas,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5}
@views @fastmath @inbounds dydt[d.indexes[3]] = sum(dydt[d.indexes[1]:d.indexes[2]]) * R * T / V
for ind in d.constantspeciesinds #make dydt zero for constant species
@inbounds dydt[ind] = 0.0
end
for inter in interfaces
if isa(inter, Inlet) && d == inter.domain
dydt[d.indexes[1]:d.indexes[2]] .+= inter.y .* inter.F(t)
dydt[d.indexes[3]] += inter.F(t) * R * T / P
elseif isa(inter, Outlet) && d == inter.domain
dydt[d.indexes[1]:d.indexes[2]] .-= inter.F(t) .* ns ./ N
dydt[d.indexes[3]] -= inter.F(t) * R * T / P
elseif isa(inter, kLAkHCondensationEvaporationWithReservoir) && d == inter.domain
kLAs = map.(inter.kLAs, inter.T)
kHs = map.(inter.kHs, inter.T)
evap = kLAs .* inter.V .* inter.cs
cond = kLAs .* inter.V .* cs * R * T ./ kHs
net_evap = evap .- cond
dydt[d.indexes[1]:d.indexes[2]] .+= net_evap
dydt[d.indexes[3]] += sum(net_evap) * R * T / P
elseif isa(inter, VolumetricFlowRateInlet) && d == inter.domain
dydt[d.indexes[1]:d.indexes[2]] .+= inter.Vin(t) * inter.cs
dydt[d.indexes[3]] += inter.Vin(t)
elseif isa(inter, VolumetricFlowRateOutlet) && d == inter.domain
dydt[d.indexes[1]:d.indexes[2]] .-= inter.Vout(t) * ns / V
dydt[d.indexes[3]] -= inter.Vout(t)
end
end
for inter in interfaces
if isa(inter, VolumeMaintainingOutlet) && d == inter.domain #VolumeMaintainingOutlet has to be evaluated after dVdt has been modified by everything else
@inbounds dVdt = dydt[d.indexes[3]]
@inbounds flow = P * dVdt / (R * T)
@views @inbounds dydt[d.indexes[1]:d.indexes[2]] .-= flow * ns / N
@inbounds dydt[d.indexes[3]] -= dVdt
end
end
end

@inline function calcdomainderivatives!(d::ConstantVDomain{W,Y}, dydt::K, interfaces::Z12; t::Z10, T::Z4, P::Z9, Us::Z, Hs::Z11, V::Z2, C::Z3, ns::Z5, N::Z6, Cvave::Z7) where {Z12,Z11,Z10,Z9,W<:IdealGas,Z7,K,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5}
@views @fastmath @inbounds dydt[d.indexes[3]] = -dot(Us, dydt[d.indexes[1]:d.indexes[2]]) / (N * Cvave) #divide by V to cancel ωV to ω
@views @fastmath @inbounds dydt[d.indexes[4]] = sum(dydt[d.indexes[1]:d.indexes[2]]) * R * T / V + P / T * dydt[d.indexes[3]]
Expand Down
Loading
Loading