diff --git a/src/Domain.jl b/src/Domain.jl index 82d2e05e..30b47f2a 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -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 @@ -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]] @@ -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]] diff --git a/src/Interface.jl b/src/Interface.jl index 466022ab..c2d7f173 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -389,11 +389,65 @@ function VaporLiquidMassTransferInternalInterfaceConstantT(domaingas, domainliq, @assert isa(domaingas.phase, IdealGas) @assert isa(domainliq.phase, IdealDiluteSolution) @assert getfield.(domaingas.phase.species, :name) == getfield.(domainliq.phase.species, :name) + T = domainliq.T phase = domainliq.phase ignoremasstransferspcinds = getinterfaceignoremasstransferspcinds(domaingas, domainliq, ignoremasstransferspcnames) kLAs = [kLA(T=T) for kLA in getfield.(phase.species, :liquidvolumetricmasstransfercoefficient)] kHs = [kH(T=T) for kH in getfield.(phase.species, :henrylawconstant)] + + non_empty_indexes_kH = findall(x -> typeof(x)!=EmptyHenryLawConstant,[y.henrylawconstant for y in phase.species]) + Tmax_kH = [x.henrylawconstant.kH.smspl.Xorig[end] for x in phase.species[non_empty_indexes_kH]] + Tmin_kH = [x.henrylawconstant.kH.smspl.Xorig[1] for x in phase.species[non_empty_indexes_kH]] + largeTs_kH = non_empty_indexes_kH[findall(x->xx>T, Tmin_kH)] + + if !isempty(largeTs_kH) + species_names = getfield.(phase.species[largeTs_kH], :name) + @warn "The current temperature ($(T) K) is above the given range. kH(T) will be extrapolated for the following species:\n$(join(species_names, "\n"))" + negative_kH_indexes = findall(x->x<0, kHs) + if !isempty(negative_kH_indexes) + @warn "The extrapolated kH(T) is negative for the following species. The highest temperature in the given range will be used to calculate kH instead. \n$(join(species_names, "\n"))" + kHs[negative_kH_indexes].=[x.henrylawconstant.kH(x.henrylawconstant.kH.smspl.Xorig[end]) for x in phase.species[negative_kH_indexes]] + end + end + + if !isempty(smallTs_kH) + species_names = getfield.(phase.species[smallTs_kH], :name) + @warn "The current temperature ($(T) K) is below the given range. kH(T) will be extrapolated for the following species:\n$(join(species_names, "\n"))" + negative_kH_indexes = findall(x->x<0, kHs) + if !isempty(negative_kH_indexes) + @warn "The extrapolated kH(T) is negative for the following species. The lowest temperature in the given range will be used to calculate kH instead. \n$(join(species_names, "\n"))" + kHs[negative_kH_indexes].=[x.henrylawconstant.kH(x.henrylawconstant.kH.smspl.Xorig[1]) for x in phase.species[negative_kH_indexes]] + end + end + + non_empty_indexes_kLA = findall(x -> typeof(x)!=EmptyLiquidVolumetricMassTransferCoefficient,[y.liquidvolumetricmasstransfercoefficient for y in phase.species]) + Tmax_kLA = [x.liquidvolumetricmasstransfercoefficient.kLA.smspl.Xorig[end] for x in phase.species[non_empty_indexes_kLA]] + Tmin_kLA = [x.liquidvolumetricmasstransfercoefficient.kLA.smspl.Xorig[1] for x in phase.species[non_empty_indexes_kLA]] + largeTs_kLA = non_empty_indexes_kLA[findall(x->xx>T, Tmin_kLA)] + + if !isempty(largeTs_kLA) + species_names = getfield.(phase.species[largeTs_kLA], :name) + @warn "The current temperature ($(T) K) is above the given range. kLA(T) will be extrapolated for the following species:\n$(join(species_names, "\n"))" + negative_kLA_indexes = findall(x->x<0, kLAs) + if !isempty(negative_kLA_indexes) + @warn "The extrapolated kLA(T) is negative for the following species. The highest temperature in the given range will be used to calculate kLA instead. \n$(join(species_names, "\n"))" + kLAs[negative_kLA_indexes].=[x.liquidvolumetricmasstransfercoefficient.kLA(x.liquidvolumetricmasstransfercoefficient.kLA.smspl.Xorig[end]) for x in phase.species[negative_kLA_indexes]] + end + end + + if !isempty(smallTs_kLA) + species_names = getfield.(phase.species[smallTs_kLA], :name) + @warn "The current temperature ($(T) K) is below the given range. kLA(T) will be extrapolated for the following species:\n$(join(species_names, "\n"))" + negative_kLA_indexes = findall(x->x<0, kLAs) + if !isempty(negative_kLA_indexes) + @warn "The extrapolated kLA(T) is negative for the following species. The lowest temperature in the given range will be used to calculate kLA instead. \n$(join(species_names, "\n"))" + kLAs[negative_kLA_indexes].=[x.liquidvolumetricmasstransfercoefficient.kLA(x.liquidvolumetricmasstransfercoefficient.kLA.smspl.Xorig[1]) for x in phase.species[negative_LA_indexes]] + end + end + return VaporLiquidMassTransferInternalInterfaceConstantT(domaingas, domainliq, ignoremasstransferspcnames, ignoremasstransferspcinds, kLAs, kHs, [1, length(domainliq.phase.species)], [0, 0], ones(length(domainliq.phase.species))), ones(length(domainliq.phase.species)) end export VaporLiquidMassTransferInternalInterfaceConstantT @@ -403,15 +457,15 @@ function getkLAkHs(vl::VaporLiquidMassTransferInternalInterfaceConstantT, Tgas, end function evaluate(vl::VaporLiquidMassTransferInternalInterfaceConstantT, dydt, Vgas, Vliq, Tgas, Tliq, N1, N2, P1, P2, Cvave1, Cvave2, ns1, ns2, Us1, Us2, cstot, p) - kLAs, kHs = getkLAkHs(vl, Tgas, Tliq) - @views @inbounds @fastmath evap = kLAs * Vliq .* cstot[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] #evap_i = kLA_i * Vliq * cliq_i - @views @inbounds @fastmath cond = kLAs * Vliq .* cstot[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] * R * Tgas ./ kHs #cond_i = kLA_i * Vliq * Pgas_i / kH_i, Pgas_i = cgas_i * R * Tgas - netevap = (evap .- cond) + kLAs, kHs = getkLAkHs(vl, vl.domaingas.T, vl.domainliq.T) + @views @inbounds @fastmath evap = kLAs * vl.domainliq.V .* cstot[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] #evap_i = kLA_i * Vliq * cliq_i + @views @inbounds @fastmath cond = kLAs * vl.domainliq.V .*cstot[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] * R * vl.domaingas.T ./ kHs #cond_i = kLA_i * Vliq * Pgas_i / kH_i, Pgas_i = cgas_i * R * Tgas + net_evap = (evap .- cond) @simd for ind in vl.ignoremasstransferspcinds - @inbounds netevap[ind] = 0.0 + net_evap[ind] = 0.0 end - @views @inbounds @fastmath dydt[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] .+= netevap - @views @inbounds @fastmath dydt[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] .-= netevap + @views @inbounds @fastmath dydt[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] .+= net_evap + @views @inbounds @fastmath dydt[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] .-= net_evap end export evaluate @@ -538,4 +592,4 @@ function evaluate(ri::FragmentBasedReactiveFilmGrowthInterfaceConstantT, dydt, V end end -export evaluate \ No newline at end of file +export evaluate diff --git a/src/Simulation.jl b/src/Simulation.jl index f80f614d..884ea3b8 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -203,11 +203,11 @@ end export concentrations -getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantTVDomain,FragmentBasedConstantTrhoDomain},K<:Real,Q,G,L} = bsol.domain.T +getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantTVDomain,ConstantTVGasDomain,FragmentBasedConstantTrhoDomain},K<:Real,Q,G,L} = bsol.domain.T getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ParametrizedVDomain,ConstantPDomain,ParametrizedPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]] getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedTConstantVDomain,ParametrizedTPDomain},K<:Real,Q,G,L} = bsol.domain.T(t) export getT -getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.domain.V +getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ConstantTVGasDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.domain.V getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.domain.V(t) getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ParametrizedTPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]] getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedPDomain,ConstantPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[4]] @@ -215,11 +215,12 @@ getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:FragmentBasedConstantTrhoDomain, export getV getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantPDomain},K<:Real,Q,G,L} = bsol.domain.P getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTVDomain,ParametrizedTConstantVDomain,FragmentBasedConstantTrhoDomain},K<:Real,Q,G,L} = 1.0e6 +getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTVGasDomain,K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]] getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedTPDomain,ParametrizedPDomain},K<:Real,Q,G,L} = bsol.domain.P(t) getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ParametrizedVDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[4]] export getP getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L} = bsol.domain.P / (R * bsol.domain.T) -getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.N(t) / bsol.domain.V +getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ConstantTVGasDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.N(t) / bsol.domain.V getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:FragmentBasedConstantTrhoDomain,K<:Real,Q,G,L} = bsol.N(t) / getV(bsol, t) getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.N(t) / bsol.domain.V(t) getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedTPDomain,K<:Real,Q,G,L} = bsol.domain.P(t) / (R * bsol.domain.T(t)) @@ -253,7 +254,11 @@ 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)]) + Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims]) + Nrxns_interface = [length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter, :reactions)] + if Nrxns_interface != Union{}[] + Nxrns += sum(Nrxns_interface) + end Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims]) cstot = zeros(Nspcs) vns = Array{Any,1}(undef, length(domains))