Advanced Thermodynamics

Robust density solver for Helmholtz Equations of State

30 October 2021

In the series about helium-containing mixtures, I created a simple density solver and a simple flash procedure compatible with multiparameter Helmholtz Equations of State like GERG-2008 [1]. What I found was that the simple solvers lacked robustness: sometimes they would fail to converge or they may have converged to the wrong solution. In this post, I'd like to investigate more robust methods of solving for the density. If I find an acceptable method, I may then be able to pursue a robust method for solving the more-complicated flash problem.

In Part 9 of the helium series, I showed the beginnings of a robust solver based on interval root finding. It looked like the following:

function ρ_robust(m::FluidMixture,p_target,T,x)
    # find all density roots between 1.e-4 mol/m³ and 1.e5 mol/m³
    return IntervalRootFinding.roots(ρ->p(m,ρ,T,x)-p_target,1.e-4..1.e5)
end

• it returns all roots, even unstable/unphysical ones
• it's extremely slow — 10 or 100 thousand times slower in some cases than my quick-and-dirty density solver

methane = AdvancedThermo.get_fluid_component("methane")
m = AdvancedThermo.create_mixture([methane])
T = 180.0
x = [1.0]
@time AdvancedThermo.ρ_robust(m,1.e6,T,x)
    2.209085 seconds (6.28 M allocations: 244.578 MiB, 5.04% gc time)
    1-element Vector{Root{Interval{Float64}}}:
    Root([736.006, 736.007], :unique)

The basic benchmarking above indicates that 244 MiB of allocations are happening when the solver runs. Allocating memory (and the subsequent cleanup) is very time-consuming in Julia, so the first thing to investigate is where the allocations are happening. I'll post the script here but more information can be found online (e.g. here).

julia --project --track-allocation=user
using AdvancedThermo
methane = AdvancedThermo.get_fluid_component("methane")
m = AdvancedThermo.create_mixture([methane])
T = 180.0
x = [1.0]
AdvancedThermo.ρ_robust(m,1.e6,T,x)
using Profile
Profile.clear_malloc_data()
AdvancedThermo.ρ_robust(m,1.e6,T,x)
[Ctrl-D]

    - function ρ_robust(m::FluidMixture,p_target,T,x)
    -     # find all density roots between 1.e-4 mol/m³ and 1.e5 mol/m³
100384     return IntervalRootFinding.roots(ρ->p(m,ρ,T,x)-p_target,1.e-4..1.e5)
    - end
    - function a(m::FluidMixture,v,T,x)
587936     return R*T*α(m,inv(v),T,x)
    - end
    - function ∂a∂v(m::FluidMixture,v,T,x)
    -     f = v -> a(m,v,T,x)
69520     return ForwardDiff.derivative(f,v)
    - end
    - function p(m::FluidMixture,ρ,T,x)
239728     return -∂a∂v(m,inv(ρ),T,x)
    - end
    - function α(m::FluidMixture,ρ,T,x)
216352     return αᵒ(m,ρ,T,x) + αʳ(m,ρ,T,x)
    - end
    - function αᵒ(fluid::FluidComponent,ρ,T)
    0     δᵢ = ρ/fluid.ρc
    0     τᵢ = fluid.Tc/T
277776     s = [
    -         1.0, 
    -         τᵢ, 
    -         log(τᵢ), 
    -         log(abs(sinh(fluid.vi[4]*τᵢ))), 
    -        -log(abs(cosh(fluid.vi[5]*τᵢ))), 
    -         log(abs(sinh(fluid.vi[6]*τᵢ))), 
    -        -log(abs(cosh(fluid.vi[7]*τᵢ)))]
526816     return log(δᵢ) + fluid.ni's
    - end
    - function αᵒ(m::FluidMixture,ρ,T,x)
    -     result = 0.0
    0     for i = 1:length(x)
    0         if x[i]>0.0
463872             result += x[i]*(αᵒ(m.components[i],ρ,T)+log(x[i]))
    -         end
    -     end
    0     return result
    - end
    - function αʳ(fluid::FluidComponent,ρ,T)
    -     δ = ρ/fluid.ρc
    -     τ = fluid.Tc/T
    -     return αʳ(δ,τ,fluid)
    - end
    - function αʳ(m::FluidMixture,ρ,T,x)
    0     δ = ρ/ρᵣ(x,m.ρc,m.βv,m.γv)
    0     τ = Tᵣ(x,m.Tc,m.βT,m.γT)/T
    0     return αʳ(δ,τ,x,m)
    - end
    - function nαʳ(m::FluidMixture,VTn)
    -     ntotal = sum(VTn[3:end])
    -     return ntotal*αʳ(m,ntotal/VTn[1],VTn[2],VTn[3:end]./ntotal)
    - end
    - function αʳ(δ,τ,x,m::FluidMixture)
278080     return sum(x.*αʳ.(δ,τ,m.components))
    - end
    - function αʳ(δ,τ,f::FluidComponent)
    0     return αʳ(δ,τ,f.np,f.dp,f.tp) + αʳ(δ,τ,f.ne,f.de,f.te,f.γe,f.ce)
    - end
    - function αʳ(δ,τ,n,d,t)
633776     return sum(@.(n*δ^d*τ^t))
    - end
    - function αʳ(δ,τ,n,d,t,γ,c)
1576800     return sum(@.(n*δ^d*τ^t*exp(-γ*δ^c)))
    - end
    - function ρᵣ(x,ρc,βv,γv)
185184     return inv(vᵣ(x,inv.(ρc),βv,γv))
    - end
    - function vᵣ(x,vc,βv,γv)
    -     result = 0.0
    0     for j = 1:length(x)
    0         if x[j]>0.0
    0             for i = 1:length(x)
    0                 if x[i]>0.0
    0                     vᵢⱼ = (vc[i]^(1/3)+vc[j]^(1/3))^3/8.0
    0                     result += ( x[i]*x[j]*βv[i,j]*γv[i,j]*
    -                                 (x[i]+x[j])/(βv[i,j]^2*x[i]+x[j])*vᵢⱼ )
    -                 end
    -             end
    -         end
    -     end
    0     return result
    - end
    - function Tᵣ(x,Tc,βT,γT)
    -     result = 0.0
    0     for j = 1:length(x)
    0         if x[j]>0.0
    0             for i = 1:length(x)
    0                 if x[i]>0.0
    0                     Tᵢⱼ = sqrt(Tc[i]*Tc[j])
    0                     result += ( x[i]*x[j]*βT[i,j]*γT[i,j]*
    -                                 (x[i]+x[j])/(βT[i,j]^2*x[i]+x[j])*Tᵢⱼ )
    -                 end
    -             end
    -         end
    -     end
    0     return result
    - end

    - function setprecision(f::Function, ::Type{Interval}, prec::Integer)
    - 
6431264     old_precision = precision(Interval)
    0     setprecision(Interval, prec)
    - 
    -     try
32156320         return f()
    -     finally
    0         setprecision(Interval, old_precision)
    -     end
    - end

    - @static if Sys.iswindows()  # Windows cannot round properly
    -     function atomic(::Type{Interval{T}}, x::S) where {T<:AbstractFloat, S<:Real}
    -         isinf(x) && return wideinterval(T(x))
    - 
26308416         Interval{T}( parse(T, string(x), RoundDown),
    -                      parse(T, string(x), RoundUp) )
    -     end
    - 
    -     function atomic(::Type{Interval{T}}, x::Union{Irrational,Rational}) where {T<:AbstractFloat}
    0         isinf(x) && return wideinterval(T(x))
    - 
77175168         Interval{T}( T(x, RoundDown), T(x, RoundUp) )
    -         # the rounding up could be done as nextfloat of the rounded down one?
    -         # use @round_up and @round_down here?
    -     end
    - 
    - else
    -     function atomic(::Type{Interval{T}}, x::S) where {T<:AbstractFloat, S<:Real}
    -         isinf(x) && return wideinterval(T(x))
    - 
    -         Interval{T}( T(x, RoundDown), T(x, RoundUp) )
    -         # the rounding up could be done as nextfloat of the rounded down one?
    -         # use @round_up and @round_down here?
    -     end
    - end

    - function ^(a::Interval{BigFloat}, n::Integer)
    0     isempty(a) && return a
4779792     n == 0 && return one(a)
    0     n == 1 && return a
    -     # n == 2 && return sqr(a)
    0     n < 0 && a == zero(a) && return emptyinterval(a)
    - 
    -     T = BigFloat
    - 
    0     if isodd(n) # odd power
    0         isentire(a) && return a
    0         if n > 0
    0             a.lo == 0 && return @round(0, a.hi^n)
    0             a.hi == 0 && return @round(a.lo^n, 0)
14387408             return @round(a.lo^n, a.hi^n)
    -         else
    0             if a.lo ≥ 0
    0                 a.lo == 0 && return @round(a.hi^n, Inf)
    0                 return @round(a.hi^n, a.lo^n)
    - 
    0             elseif a.hi ≤ 0
    0                 a.hi == 0 && return @round(-Inf, a.lo^n)
    0                 return @round(a.hi^n, a.lo^n)
    -             else
    0                 return entireinterval(a)
    -             end
    -         end
    - 
    -     else # even power
    0         if n > 0
    0             if a.lo ≥ 0
22318464                 return @round(a.lo^n, a.hi^n)
    0             elseif a.hi ≤ 0
    0                 return @round(a.hi^n, a.lo^n)
    -             else
    0                 return @round(mig(a)^n, mag(a)^n)
    -             end
    - 
    -         else
    0             if a.lo ≥ 0
    0                 return @round(a.hi^n, a.lo^n)
    0             elseif a.hi ≤ 0
    0                 return @round(a.lo^n, a.hi^n)
    -             else
    0                 return @round(mag(a)^n, mig(a)^n)
    -             end
    -         end
    -     end
    - end

• most of the allocations are coming from within the library code
• the library doesn't seem to offer any features that would allow me to pre-allocate memory
• other julia packages for intervals exist but they don't offer root-finding nor seem to be active

References

[1] Kunz and Wagner, J. Chem. Eng. Data, 2012, 57, 3032 (link to publisher)

[2] Aursand et al., Fluid Phase Equilib., 2017, 436, 98 (link to publisher)