Further performance enhancement

[Sush1090]

Hello,

There seems to be a vast difference in the performance of pure vs mixtures.

using Clapeyron, BenchmarkTools

fluid = cPR(["Propane"],idealmodel = ReidIdeal); z = [1.0]
p  =101325.0
T(x) = Clapeyron.PH.temperature(fluid,p,x,z)

h = collect(range(-1,1,100));
@btime T.($h) # 1.243 ms (209 allocations: 4.30 KiB)

fluid2 = cPR(["Propane","R134a"],idealmodel = ReidIdeal); z2 = [1.0,1.0]
T2(x) =  Clapeyron.PH.temperature(fluid2,p,x,z2) # 11.435 ms (120729 allocations: 11.51 MiB)
@btime T2.($h);
nothing

There is almost no difference now between Tproperty and Clapeyron.PH.temperature since 0.6.15.

I am not sure where the bottlenecks might be now. But the difference seems significant.

[abcdvvvv]

I tested the speed of v0.6.14 and 0.6.15 on my computer, Julia 1.10.10. The performance difference between pure components and mixtures has existed since 0.6.14. Furthermore, 0.6.15 has improved speed on the T(x) and T2(x) functions you provided.

For @btime T.($h)
0.6.14 1.936 ms (6108 allocations: 632.39 KiB)
0.6.15 1.685 ms (208 allocations: 4.27 KiB)

for @btime T2.($h);
0.6.14 23.201 ms (111724 allocations: 17.22 MiB)
0.6.15 20.636 ms (73824 allocations: 14.50 MiB)

The PH_property function was changed a lot.

[Sush1090]

Yes indeed. That difference is due to the fixes for #407.
But for this I am unable to find the most significant bottlenecks as of now. I am using ProfileViewer.jl and Cthulhu.jl but generic bottlenecks are not straightforward for me to spot.

[abcdvvvv]

I'm as concerned as you are about the performance of this package. I have already identified the bottleneck. I used PProf, and the code is as follows:

fluid2 = cPR(["Propane", "R134a"], idealmodel=ReidIdeal);
z2 = [1.0, 1.0];
T2(x) = Clapeyron.PH.temperature(fluid2, p, x, z2) # 11.435 ms (120729 allocations: 11.51 MiB)
using PProf, Profile
@pprof for _ = 1:1000
    T2.(h)
end

The result is that the x0_Property function takes up 50% of the time, calling dew_temperature and bubble_temperature, with ForwardDiff taking up the most time internally. I think the author is also aware of this.

Btw, I don't know what the other 50% on the left side is.

[longemen3000]

One path to impeove the performance is improving the bubbledew solvers. At the moment, we use nested ForwardDiff.gradient calls, that, while composable and really generic, leave some performance on the table:

Clapeyron.jl/src/methods/property_solvers/multicomponent/multicomponent.jl

Line 136 in 9ccb7be

function μp_equality2(models::NTuple{2,M}, F, PT::TPspec, v, w, short_view) where M <: EoSModel

Another path to make the operations faster is to merge the initial point generation for the bubbledew solvers in one routine. To generate the initial points, we:

• split the model in pure models
• calculate the saturation conditions for each pure (or critical extrapolation if needed)
• solve the bubble or dew Raoult problem

The first two steps could be done just once.

[Sush1090]

There might also be a third option.

We make Tproperty and Pproperty free of bubble and dew computations and rely purely on FindEdge for phase separation.

I am not sure on the robustness of this. But possibly worth a try.

[Sush1090]

Hello,

Here is the proof of concept for my previous comment. It is very premature hence difficult to make any conclusion.

using Clapeyron, Roots

function Tproperty2_enthalpy(model::EoSModel,p,h,z)

f(x) = enthalpy(model,p,x,z) 
f0(x) = enthalpy(model,p,x,z) - h
edge_val,prop_edge_left, prop_edge_right = Clapeyron.FindEdge(f,1,500)

if h ∈ range(prop_edge_left,prop_edge_right)
    @warn "In two phase"
    return edge_val
end

if h < prop_edge_left
    prob = ZeroProblem(f0,edge_val - 10)
    return solve(prob)
end
if h > prop_edge_right
    prob = ZeroProblem(f0,edge_val + 10)
    return solve(prob)
end

return NaN
end

model = cPR(["R134a","propane"],idealmodel=ReidIdeal);
p = 101325*10; z = [1.0,1.0];
h = -50000;
@time Tproperty2_enthalpy(model,p,h,z) #  0.000104 seconds (275 allocations: 14.438 KiB)
@time Tproperty(model,p,h,z,enthalpy) # 0.000745 seconds (1.26 k allocations: 120.688 KiB)

With this we evade computation of bubble and dew points for initialization. As of now this is not robust but just an idea. Let me know what you think.

[Sush1090]

I have an additional comment on performance of Tproperty on mixtures.

The internal computation of critical point takes about 55% of memory.

julia> model = cPR(["R134a","propane"],idealmodel=ReidIdeal);

julia> p = 101325*10; z = [1.0,1.0];

julia> @time crit_mix(model,z);
  0.000338 seconds (588 allocations: 65.188 KiB)

julia> @time Tproperty(model,p,1000,z,enthalpy)
  0.000398 seconds (1.21 k allocations: 118.062 KiB)
319.6393678738907

I am thinking of a way to by-pass this. This was implemented because of #309.

[longemen3000]

Hmmm, i think i got an idea on how to implement this.
Also, the find_edge function is equivalent to finding a liquid and gas phase of the same composition, with equal gibbs energy. This can be calculated reliably (if the spinodals fail, there is a high chance we are over the mechanical critical point, and the procedure is reduced to 1-phase root finding) using mechanical spinodals for the initial points.

If the property is outside the edge, you still need at least one bubble/dew calculation to confirm that your property point is outside the saturation dome, but in any case, this approach reduces the amount of bubble/dew calculations from two to one/zero

[Sush1090]

• solve the bubble or dew Raoult problem

If we do this will there be inconsistency with the current flashes? Because the current one are based over the default dew and bubble solvers ChemPot

I am not sure but this might break the code?

[longemen3000]

i think i got something going on, on the fast_tp_property branch.

[Sush1090]

I was going through the code of the fast_tp_property branch. Is the idea now to use Raolut's law along with FindEdge?
The idea mentioned in your comment before.

[longemen3000]

Yes, i was looking at a good initial point, and it seems like the bubble/dew raoult pressures are good enough (specially for MultiFluid where there could be multiple spinodals)

[Sush1090]

On the latest master now for mixtures if we run the code on the issue we get:

julia> using Clapeyron, BenchmarkTools

julia> fluid2 = cPR(["Propane","R134a"],idealmodel = ReidIdeal); z2 = [1.0,1.0];

julia> h = collect(range(-1,1,100));

julia> T2(x) =  Clapeyron.PH.temperature(fluid2,p,x,z2)

julia> @btime T2.($h);
  6.976 ms (97089 allocations: 7.52 MiB) 

The computation speed has increased by almost two times and the allocations have dropped 20%.

So should we close this for now?

[longemen3000]

I'm gonna move this into a discussion, because there is probably more things to do

[Sush1090]

Hello,

Is it possible that there is a bottleneck in flash computation when in two-phase?

julia> using Clapeyron

julia> p = 101325.0; z = [1.0,1.0];

julia> fluid = cPR(["Propane","R134a"],idealmodel = ReidIdeal);

julia> Tb = bubble_temperature(fluid,p,z)[1];

julia> Td = dew_temperature(fluid,p,z)[1];

julia> hb = enthalpy(fluid,p,Tb,z);

julia> hd = enthalpy(fluid,p,Td,z);

julia> h_ = collect(range(hb,hd,100));

julia> T2(x) =  Clapeyron.PH.temperature(fluid,p,x,z);

julia> @btime T2.($h_);
  31.784 ms (319636 allocations: 33.31 MiB)

I am not able to completely understand the algorithm GeneralizedXYFlash when in two phase

[abcdvvvv]

I will begin implementing some extreme performance optimizations for this discussion. We need a continuous performance log; please consider the following options:

• PkgBenchmark.jl
• Directly record local bench results within this discussion
• AirspeedVelocity.jl

[Sush1090]

I will check this thank you

[longemen3000]

I've been slowly working towards an unified cache system, around dlnphi_cache. The next step is to allow GeneralizedXYFlash to use those caches, the same goes to the dew and bubble functions.
in the meantime, I can add that algorithm in a proper way in Clapeyron.

[abcdvvvv]

Great!

[longemen3000]

i just added a RRQXFlash method, with a qp_flash implementation. with some minor touches, it could also work for activity models and non-condensable/non-volatile specifications. the qt_flash implementation is just copying the qp_flash and changing temperatures to pressures

[abcdvvvv]

Thank you, it works fine on my computer. A simplest example is as follows:

mixture = PR(["IsoButane", "n-Butane", "n-Pentane", "n-Hexane"]; idealmodel=ReidIdeal)
z = [0.25, 0.25, 0.25, 0.25]
Clapeyron.qp_flash_impl(mixture, 1, 1e5, z, RRQXFlash(equilibrium=:vle))
#=
Flash result at T = 312.935, p = 100000.0 with 2 phases:
 (x = [0.0515604, 0.0715463, 0.219938, 0.656956], β = 0.0, v = 0.000124542)
 (x = [0.25, 0.25, 0.25, 0.25], β = 1.0, v = 0.0251559) 

92.200 μs  Memory estimate: 61.35 KiB, allocs estimate: 692.
=#

I plan to use it as a replacement for the dew point calculation function, as it is a full 50% faster.

Clapeyron.dew_temperature(mixture, 1e5, z)
#=
(312.93460566879094, 0.0001245423558751178, 0.0251559096118607, 
[0.05156042304823789, 0.07154633770179127, 0.2199376042154166, 0.6569556350345542])

173.200 μs Memory estimate: 101.74 KiB, allocs estimate: 981.
=#

[longemen3000]

One big performance improvement could be achieved calculating ∂lnϕ∂P and ∂lnϕ∂T in one gradient call instead of calculating the hessian of eos_res and deriving ∂lnϕ∂P/∂lnϕ∂T from there. I am testing a second_order::Bool on Fug-based bubble/dew functions, and most of the slowdown is explained by those methods.

[abcdvvvv]

Thanks for the second_order support.

I’ve confirmed that I can use it like this:

Clapeyron.tp_flash2(mix, p, T, n, MichelsenTPFlash(equilibrium = :vle, second_order = true))

However, in my tests most flash problems already converge in Stage 1 (successive substitution), so the speed-up from enabling second_order is not always visible in practice.

Are there any other higher-level APIs in Clapeyron that already expose this second_order option, or are intended to benefit from it?

Any pointers or examples would be very helpful.😇

[longemen3000]

to be fair, maybe the second_order name in fugacity solvers is not the most appropiate name haha, the thing is that most of our solvers are non-linear solvers, so there is a jacobian instead.

[longemen3000]

in the bubbledew2 branch, there is a set of solvers for bubble/dew conditions, that work in terms of lnK-volumes. some limited benchmarking shows around a 50% speed improvements in comparison with the current solvers, but it loses on reliability near the mixture critical point.

What is important is that is those implementations rely on some new functions (Clapeyron.lnf and friends) that are already implemented, so at this stage one would need to just rewrite the current solver (nc-1 fractions - volumes) using those functions instead.

[Sush1090]

Hello,

With the implementation of the new methods for solving ( I think since 0.6.20) lot has changed.

The core MWE of this thread had gained significant performance somewhere around v0.6.17

See : comment in this thread

Now it seems to me that with the new implementations the performance has dropped back almost to 0.6.15 the same as beginning of this thread.

using Clapeyron, BenchmarkTools
p = 101325.0
h = collect(range(-1,1,100));
fluid2 = cPR(["Propane","R134a"],idealmodel = ReidIdeal); z2 = [1.0,1.0]
T2(x) =  Clapeyron.PH.temperature(fluid2,p,x,z2) # 11.435 ms (120729 allocations: 11.51 MiB) 
@btime T2.($h); # 9.125 ms (103969 allocations: 10.82 MiB)
nothing

I had not checked the performance on this after the unified cache system was implemented as well.