From 5bc511a00f1c8159a3468726505c2bdbbf8311e1 Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Thu, 4 Apr 2024 10:41:55 -0400 Subject: [PATCH 01/10] remove constantTVGasDomain --- src/Domain.jl | 261 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 261 insertions(+) diff --git a/src/Domain.jl b/src/Domain.jl index 82d2e05e..0e193b6d 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -102,6 +102,94 @@ 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,V,C,N,ns,Gs,diffs,P,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 ConstantTPDomain(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 ConstantTPDomain + 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 @@ -820,6 +908,7 @@ end export FragmentBasedConstantTrhoDomain +<<<<<<< HEAD mutable struct ConstantTLiqFilmDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:Integer,Q<:AbstractArray} <: AbstractConstantKDomain phase::N indexes::Q #assumed to be in ascending order @@ -843,6 +932,172 @@ mutable struct ConstantTLiqFilmDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Re t::MArray{Tuple{1},W2,1,1} p::Array{W,1} thermovariabledict::Dict{String,Int64} +======= + +@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,d.P,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,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::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::ConstantTPDomain{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]] + V = y[d.indexes[3]] + N = d.P*V/(R*d.T) + cs = ns./V + C = N/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,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0) + end + return ns,cs,d.T,d.P,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 +>>>>>>> 02018d3 (add ConstantTVGasDomain) end function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspecies::Array{X2,1}=Array{String,1}(), @@ -985,7 +1240,13 @@ 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 +<<<<<<< HEAD @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p +======= + + +@inline function calcthermo(d::ConstantTPDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p +>>>>>>> 02018d3 (add ConstantTVGasDomain) ns = y[d.indexes[1]:d.indexes[2]] V = y[d.indexes[3]] N = d.P * V / (R * d.T) From 467c9091b7df927b575528a9f3ea7e26aff6f9f1 Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Mon, 29 Apr 2024 17:02:07 -0400 Subject: [PATCH 02/10] add constantTVGasDomain --- src/Domain.jl | 329 +++++++++++++++++++++++--------------------------- 1 file changed, 154 insertions(+), 175 deletions(-) diff --git a/src/Domain.jl b/src/Domain.jl index 0e193b6d..c193e5c9 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -188,7 +188,19 @@ function ConstantTVGasDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies return ConstantTPDomain(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 ConstantTPDomain +export ConstantTVGasDomain +@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 + struct ConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:Integer,Q<:AbstractArray} <: AbstractVariableKDomain phase::N @@ -908,7 +920,6 @@ end export FragmentBasedConstantTrhoDomain -<<<<<<< HEAD mutable struct ConstantTLiqFilmDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:Integer,Q<:AbstractArray} <: AbstractConstantKDomain phase::N indexes::Q #assumed to be in ascending order @@ -932,172 +943,6 @@ mutable struct ConstantTLiqFilmDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Re t::MArray{Tuple{1},W2,1,1} p::Array{W,1} thermovariabledict::Dict{String,Int64} -======= - -@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,d.P,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,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::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::ConstantTPDomain{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]] - V = y[d.indexes[3]] - N = d.P*V/(R*d.T) - cs = ns./V - C = N/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,d.P,C,N,ns,d.Gs,d.diffusivity,V,0.0) - end - return ns,cs,d.T,d.P,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 ->>>>>>> 02018d3 (add ConstantTVGasDomain) end function ConstantTLiqFilmDomain(; phase::Z, initialconds::Dict{X,E}, constantspecies::Array{X2,1}=Array{String,1}(), @@ -1240,13 +1085,7 @@ 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 -<<<<<<< HEAD @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p -======= - - -@inline function calcthermo(d::ConstantTPDomain{W,Y},y::J,t::Q,p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p ->>>>>>> 02018d3 (add ConstantTVGasDomain) ns = y[d.indexes[1]:d.indexes[2]] V = y[d.indexes[3]] N = d.P * V / (R * d.T) @@ -1322,6 +1161,146 @@ 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::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,d.P,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,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::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]] @@ -3853,4 +3832,4 @@ end else return eval(d)(phase=IdealDiluteSolution(sensspcs, sensrxns, domain.phase.solvent, name="phase"), initialconds=initialconds)[1], sensspcnames, senstooriginspcind, senstooriginrxnind end -end +end \ No newline at end of file From d5b47f0176e2cadd769ff60f36925a90e0f54a3b Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Tue, 7 May 2024 13:53:13 -0400 Subject: [PATCH 03/10] still does not work --- src/Domain.jl | 95 ++++++++++++++++++++++++++++++++++------------- src/Interface.jl | 14 +++---- src/Reactor.jl | 48 ++++++++++++++++++++++++ src/Simulation.jl | 7 ++-- 4 files changed, 129 insertions(+), 35 deletions(-) diff --git a/src/Domain.jl b/src/Domain.jl index c193e5c9..a6442a9c 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -144,8 +144,6 @@ function ConstantTVGasDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies y0[ind] = val end end - - @assert T != 0.0 @assert V != 0.0 ns = y0[1:end-1] @@ -173,7 +171,7 @@ function ConstantTVGasDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies C = P/(R*T) y0[end] = P - kfs,krevs = getkfkrevs(phase,T,V,C,N,ns,Gs,diffs,P,0.0) + 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 @@ -185,21 +183,12 @@ function ConstantTVGasDomain(;phase::E2,initialconds::Dict{X,X2},constantspecies jacobian=zeros(typeof(T),length(phase.species),length(phase.species)) end rxnarray = getreactionindices(phase) - return ConstantTPDomain(phase,[1,length(phase.species),length(phase.species)+1],[1,length(phase.species)+length(phase.reactions)],constspcinds, + 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 -@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 + + struct ConstantVDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,I<:Integer,Q<:AbstractArray} <: AbstractVariableKDomain @@ -1162,10 +1151,21 @@ end 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) + N = P * d.V / (R * d.T) cs = ns./d.V C = N/d.V if !d.alternativepformat @@ -1180,13 +1180,13 @@ end 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,d.P,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 + 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 + 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 @@ -1198,14 +1198,15 @@ end 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 + 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) + N = P * d.V / (R * d.T) + cs = ns./d.V C = N/d.V if !d.alternativepformat @@ -1225,7 +1226,7 @@ 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) + N = P * d.V / (R * d.T) cs = ns./d.V C = N/d.V if !d.alternativepformat @@ -1239,13 +1240,13 @@ end 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,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 + 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) + N = P * d.V / (R * d.T) cs = ns./d.V C = N/d.V kfs = similar(y,length(d.phase.reactions)) @@ -1267,7 +1268,7 @@ 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) + N = P * d.V / (R * d.T) cs = ns./d.V C = N/d.V if !d.alternativepformat @@ -1284,7 +1285,7 @@ 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) + N = P * d.V / (R * d.T) cs = ns./d.V C = N/d.V if !d.alternativepformat @@ -2138,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]] @@ -2441,6 +2480,8 @@ end end end + + @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantTPDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) @@ -2534,6 +2575,7 @@ end end end + @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantVDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) @@ -2693,6 +2735,9 @@ end return jac end + + + @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantPDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) diff --git a/src/Interface.jl b/src/Interface.jl index 466022ab..313aa7ed 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -403,15 +403,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 @fastmath evap = kLAs * vl.domainliq.V .* cstot[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] #evap_i = kLA_i * Vliq * cliq_i + @views @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 @fastmath dydt[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] .+= net_evap + @views @fastmath dydt[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] .-= net_evap end export evaluate diff --git a/src/Reactor.jl b/src/Reactor.jl index 549e9ead..5816c677 100644 --- a/src/Reactor.jl +++ b/src/Reactor.jl @@ -1395,6 +1395,54 @@ end end end + +@inline function _jacobianynswrtP!(jac::S, Pind::Int64, rxnarray::Array{Int64,2}, rxnind::Int64, cs::Array{Float64,1}, kf::Float64, krev::Float64) where {S<:AbstractArray} + k = kf + if rxnarray[2, rxnind] == 0 + nothing + elseif rxnarray[3, rxnind] == 0 + @inbounds @fastmath deriv = -k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] + @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv + _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) + elseif rxnarray[4, rxnind] == 0 + @inbounds @fastmath deriv = -2.0 * k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] * cs[rxnarray[3, rxnind]] + @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[3, rxnind], Pind] -= deriv + _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) + else + @inbounds @fastmath deriv = -3.0 * k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] * cs[rxnarray[3, rxnind]] * cs[rxnarray[4, rxnind]] + @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[3, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[4, rxnind], Pind] -= deriv + _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) + end + k = krev + if rxnarray[6, rxnind] == 0 + nothing + elseif rxnarray[7, rxnind] == 0 + @inbounds @fastmath deriv = -k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] + @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv + _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) + elseif rxnarray[8, rxnind] == 0 + @inbounds @fastmath deriv = -2.0 * k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] * cs[rxnarray[7, rxnind]] + @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[7, rxnind], Pind] -= deriv + _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) + else + @inbounds @fastmath deriv = -3.0 * k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] * cs[rxnarray[7, rxnind]] * cs[rxnarray[8, rxnind]] + @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[7, rxnind], Pind] -= deriv + @inbounds jac[rxnarray[8, rxnind], Pind] -= deriv + _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) + end +end + """ This function calculates the ns partials in jacobiany involving k derivatives. dkdx is either dkdni and dkdV. x is either ni or V. """ diff --git a/src/Simulation.jl b/src/Simulation.jl index f80f614d..3460e956 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)) From 88f9d42de6cd8e780bb3310e009608ccfceab8dd Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Tue, 7 May 2024 13:57:29 -0400 Subject: [PATCH 04/10] clean up --- src/Interface.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Interface.jl b/src/Interface.jl index 313aa7ed..49fb0266 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -404,14 +404,14 @@ 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, vl.domaingas.T, vl.domainliq.T) - @views @fastmath evap = kLAs * vl.domainliq.V .* cstot[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] #evap_i = kLA_i * Vliq * cliq_i - @views @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 + @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 net_evap[ind] = 0.0 end - @views @fastmath dydt[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] .+= net_evap - @views @fastmath dydt[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] .-= net_evap + @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 +538,4 @@ function evaluate(ri::FragmentBasedReactiveFilmGrowthInterfaceConstantT, dydt, V end end -export evaluate \ No newline at end of file +export evaluate From c22d723c1adda8c0f3b81f177f6ab248471fad38 Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Tue, 7 May 2024 13:59:35 -0400 Subject: [PATCH 05/10] more cleanup --- src/Reactor.jl | 48 ------------------------------------------------ 1 file changed, 48 deletions(-) diff --git a/src/Reactor.jl b/src/Reactor.jl index 5816c677..549e9ead 100644 --- a/src/Reactor.jl +++ b/src/Reactor.jl @@ -1395,54 +1395,6 @@ end end end - -@inline function _jacobianynswrtP!(jac::S, Pind::Int64, rxnarray::Array{Int64,2}, rxnind::Int64, cs::Array{Float64,1}, kf::Float64, krev::Float64) where {S<:AbstractArray} - k = kf - if rxnarray[2, rxnind] == 0 - nothing - elseif rxnarray[3, rxnind] == 0 - @inbounds @fastmath deriv = -k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] - @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv - _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) - elseif rxnarray[4, rxnind] == 0 - @inbounds @fastmath deriv = -2.0 * k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] * cs[rxnarray[3, rxnind]] - @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[3, rxnind], Pind] -= deriv - _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) - else - @inbounds @fastmath deriv = -3.0 * k * cs[rxnarray[1, rxnind]] * cs[rxnarray[2, rxnind]] * cs[rxnarray[3, rxnind]] * cs[rxnarray[4, rxnind]] - @inbounds jac[rxnarray[1, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[2, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[3, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[4, rxnind], Pind] -= deriv - _spreadreactantpartials!(jac, deriv, rxnarray, rxnind, Pind) - end - k = krev - if rxnarray[6, rxnind] == 0 - nothing - elseif rxnarray[7, rxnind] == 0 - @inbounds @fastmath deriv = -k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] - @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv - _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) - elseif rxnarray[8, rxnind] == 0 - @inbounds @fastmath deriv = -2.0 * k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] * cs[rxnarray[7, rxnind]] - @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[7, rxnind], Pind] -= deriv - _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) - else - @inbounds @fastmath deriv = -3.0 * k * cs[rxnarray[5, rxnind]] * cs[rxnarray[6, rxnind]] * cs[rxnarray[7, rxnind]] * cs[rxnarray[8, rxnind]] - @inbounds jac[rxnarray[5, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[6, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[7, rxnind], Pind] -= deriv - @inbounds jac[rxnarray[8, rxnind], Pind] -= deriv - _spreadproductpartials!(jac, deriv, rxnarray, rxnind, Pind) - end -end - """ This function calculates the ns partials in jacobiany involving k derivatives. dkdx is either dkdni and dkdV. x is either ni or V. """ From a32ad1dd0fd584bd1315d80f67d4d3ae78e15e76 Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Mon, 13 May 2024 12:34:04 -0400 Subject: [PATCH 06/10] added checks and warnings to VaporLiquidMassTransferInternalInterfaceConstantT to prevent kH and kLA being extrapolated to negative values --- src/Interface.jl | 57 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 49fb0266..3c4c803d 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 @@ -405,7 +459,8 @@ 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, 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 + @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 + cond[findall(x->x<0,kHs)].=0 net_evap = (evap .- cond) @simd for ind in vl.ignoremasstransferspcinds net_evap[ind] = 0.0 From 683c51a68cc8ca584e15d7124dc75cc69c14956e Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Mon, 13 May 2024 12:40:59 -0400 Subject: [PATCH 07/10] cleanup --- src/Domain.jl | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/Domain.jl b/src/Domain.jl index a6442a9c..30b47f2a 100644 --- a/src/Domain.jl +++ b/src/Domain.jl @@ -2480,8 +2480,6 @@ end end end - - @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantTPDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) @@ -2575,7 +2573,6 @@ end end end - @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantVDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) @@ -2735,9 +2732,6 @@ end return jac end - - - @inline function jacobiany!(jac::Q, y::U, p::W, t::Z, domain::D, interfaces::Q3, colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:ConstantPDomain} ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR = calcthermo(domain, y, t, p) jacobianynsderiv!(jac, domain, domain.rxnarray, domain.efficiencyinds, cs, kfs, krevs, T, V, C) @@ -3877,4 +3871,4 @@ end else return eval(d)(phase=IdealDiluteSolution(sensspcs, sensrxns, domain.phase.solvent, name="phase"), initialconds=initialconds)[1], sensspcnames, senstooriginspcind, senstooriginrxnind end -end \ No newline at end of file +end From 39d5bdc706b990f9ce09eccf2a70034c9d5d4bca Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Mon, 13 May 2024 13:28:41 -0400 Subject: [PATCH 08/10] fixed rops bug --- src/Simulation.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Simulation.jl b/src/Simulation.jl index 3460e956..dd1312f3 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -254,7 +254,10 @@ 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) Nspcs = sum([length(getphasespecies(sim.domain.phase)) for sim in ssys.sims]) cstot = zeros(Nspcs) vns = Array{Any,1}(undef, length(domains)) From 39319caf919e05ce18faed51628ba3c870a58cbd Mon Sep 17 00:00:00 2001 From: Anna Doner Date: Mon, 13 May 2024 13:29:16 -0400 Subject: [PATCH 09/10] fixed rops bug --- src/Simulation.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Simulation.jl b/src/Simulation.jl index dd1312f3..884ea3b8 100644 --- a/src/Simulation.jl +++ b/src/Simulation.jl @@ -258,6 +258,7 @@ function rops(ssys::SystemSimulation, t) 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)) From 512174c3303437dfd59d1eb0a08d5feca38d87a8 Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Mon, 1 Jul 2024 13:47:56 -0400 Subject: [PATCH 10/10] Update Interface.jl --- src/Interface.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Interface.jl b/src/Interface.jl index 3c4c803d..c2d7f173 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -460,7 +460,6 @@ function evaluate(vl::VaporLiquidMassTransferInternalInterfaceConstantT, dydt, V 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 - cond[findall(x->x<0,kHs)].=0 net_evap = (evap .- cond) @simd for ind in vl.ignoremasstransferspcinds net_evap[ind] = 0.0