Try to implement the method of Aasen

[ianhbell]

Basically you find the r/sigma where the value is 1, and chop the integrand at that value, so the quadrature only covers the domain where it is needed.

Described in Aasen and in this post from @longemen3000 : https://github.com/ClapeyronThermo/Clapeyron.jl/issues/152#issuecomment-1480324192

[ianhbell]

The tricky thing as pointed out in the Clapeyron.jl code, is that you need to implement the rootfinding with autodiff types. That in principle is not too tricky, but in practice isn't very straightforward. You could use autodiff for the derivative too, but that is probably overkill, and a bit of derivative-taking is likely worth it.

[longemen3000]

i don't know on the VRQ Mie case, but maybe a good function can be found analytically? i mean

C = (λr/(λr-λa))*(λr/λa)^(λa/(λr-λa)) 
u(j) = C*ϵ*(j^-λr - j^-λa), j = r/σ
EPS = exp(-u(j)/T)
-T*log(EPS)/(C*ϵ) = (j^-λr - j^-λa) # * j^λr
-T*log(EPS)/(C*ϵ)*j^λr = 1 - j^(λr - λa)

we now consider logK = -T*log(EPS)/(C*ϵ) and K = log(-T*log(EPS)/(C*ϵ))

logK*j^λr = 1 - j^(λr - λa)
K + λr*log(j) = log(1 - j^(λr - λa))
f(j) = K+ λr*log(j) - log(1 - j^(λr - λa)) = 0

f(j) seems better for rootfinding that naively doing exp(-u(j)/T)

[longemen3000]

a good initial point is j0 = exp(-K/λr), a better point is doing a quadratic expansion for f around j0:

T: 50
solution (rootfinding): 0.8832075676655532
solution (initial point): 0.9355323202328244
solution (initial point + quadratic taylor interp): 0.8832075676655532
exp(-u/T) (rootfinding): 2.9954715767218476e-18
exp(-u/T) (initial point): 3.67276723019937e-6
exp(-u/T) (initial point + quadratic taylor interp): 2.9954715767218476e-18

T: 100
solution (rootfinding): 0.852864230458969
solution (initial point): 0.884671798931776
solution (initial point + quadratic taylor interp): 0.852864230458969
exp(-u/T) (rootfinding): 1.785485725746699e-16
exp(-u/T) (initial point): 3.102541108581542e-9
exp(-u/T) (initial point + quadratic taylor interp): 1.785485725746699e-16

T: 300
solution (rootfinding): 0.7930333214007922
solution (initial point): 0.8096636887283977
solution (initial point + quadratic taylor interp): 0.7930333214007922
exp(-u/T) (rootfinding): 2.2006266890486548e-16
exp(-u/T) (initial point): 2.507915367849278e-12
exp(-u/T) (initial point + quadratic taylor interp): 2.2006266890486548e-16

T: 1000
solution (rootfinding): 0.7265743254617227
solution (initial point): 0.7347456473933294
solution (initial point + quadratic taylor interp): 0.7265743254617227
exp(-u/T) (rootfinding): 2.2191962696951637e-16
exp(-u/T) (initial point): 3.33886433197667e-14
exp(-u/T) (initial point + quadratic taylor interp): 2.2191962696951637e-16

function:

using Clapeyron, Clapeyron.Roots
model = SAFTVRMie(["ethane"])
function d_test(model, T)
    ϵ = model.params.epsilon[1]
    σ = model.params.sigma[1]
    λa = model.params.lambda_a[1]
    λr = model.params.lambda_r[1]
    C  = (λr/(λr-λa))*(λr/λa)^(λa/(λr-λa)) 
    EPS = eps(typeof(C*T))
    K = log(-T*log(EPS)/(C*ϵ))
    #initial point
    j0 = exp(-K/λr)
   #taylor expansion 
   Δλ = λr - λa
    b = 1 - λr
    j0Δλ = j0^Δλ
    c = -log(1 - j0Δλ) + λr*log(j0) + K
    b = (λr - (Δλ*j0Δλ)/(j0Δλ - 1))/j0 
    a = ((Δλ*j0Δλ*(j0Δλ + Δλ - 1))/(j0Δλ - 1)^2 - λr)/(2*j0*j0) 
    y = (-b + sqrt(b*b - 4*a*c))/(2*a)
    j = y + j0
    _u(x) = C*ϵ*(x^-λr - x^-λa)#, j = r/σ
    _f(x) = exp(-_u(x)/T)
    xinf0 = Roots.find_zero(_f,j)
    println("T: $T")
    println("solution (rootfinding): $xinf0")
    println("solution (initial point): $j0")
    println("solution (initial point + quadratic taylor interp): $j")
    println("exp(-u/T) (rootfinding): $(_f(xinf0))")
    println("exp(-u/T) (initial point): $(_f(j0))")
    println("exp(-u/T) (initial point + quadratic taylor interp): $(_f(j))")
end

[ianhbell]

i don't know on the VRQ Mie case, but maybe a good function can be found analytically? i mean

C = (λr/(λr-λa))*(λr/λa)^(λa/(λr-λa)) 
u(j) = C*ϵ*(j^-λr - j^-λa), j = r/σ
EPS = exp(-u(j)/T)
-T*log(EPS)/(C*ϵ) = (j^-λr - j^-λa) # * j^λr
-T*log(EPS)/(C*ϵ)*j^λr = 1 - j^(λr - λa)

we now consider logK = -T*log(EPS)/(C*ϵ) and K = log(-T*log(EPS)/(C*ϵ))

logK*j^λr = 1 - j^(λr - λa)
K + λr*log(j) = log(1 - j^(λr - λa))
f(j) = K+ λr*log(j) - log(1 - j^(λr - λa)) = 0

f(j) seems better for rootfinding that naively doing exp(-u(j)/T)

I wonder if the code in Clapeyron is doing something analogous, and approximating the cutoff by something along these lines.

[ianhbell]

a good initial point is j0 = exp(-K/λr), a better point is doing a quadratic expansion for f around j0:

T: 50
solution (rootfinding): 0.8832075676655532
solution (initial point): 0.9355323202328244
solution (initial point + quadratic taylor interp): 0.8832075676655532
exp(-u/T) (rootfinding): 2.9954715767218476e-18
exp(-u/T) (initial point): 3.67276723019937e-6
exp(-u/T) (initial point + quadratic taylor interp): 2.9954715767218476e-18

T: 100
solution (rootfinding): 0.852864230458969
solution (initial point): 0.884671798931776
solution (initial point + quadratic taylor interp): 0.852864230458969
exp(-u/T) (rootfinding): 1.785485725746699e-16
exp(-u/T) (initial point): 3.102541108581542e-9
exp(-u/T) (initial point + quadratic taylor interp): 1.785485725746699e-16

T: 300
solution (rootfinding): 0.7930333214007922
solution (initial point): 0.8096636887283977
solution (initial point + quadratic taylor interp): 0.7930333214007922
exp(-u/T) (rootfinding): 2.2006266890486548e-16
exp(-u/T) (initial point): 2.507915367849278e-12
exp(-u/T) (initial point + quadratic taylor interp): 2.2006266890486548e-16

T: 1000
solution (rootfinding): 0.7265743254617227
solution (initial point): 0.7347456473933294
solution (initial point + quadratic taylor interp): 0.7265743254617227
exp(-u/T) (rootfinding): 2.2191962696951637e-16
exp(-u/T) (initial point): 3.33886433197667e-14
exp(-u/T) (initial point + quadratic taylor interp): 2.2191962696951637e-16

function:

using Clapeyron, Clapeyron.Roots
model = SAFTVRMie(["ethane"])
function d_test(model, T)
    ϵ = model.params.epsilon[1]
    σ = model.params.sigma[1]
    λa = model.params.lambda_a[1]
    λr = model.params.lambda_r[1]
    C  = (λr/(λr-λa))*(λr/λa)^(λa/(λr-λa)) 
    EPS = eps(typeof(C*T))
    K = log(-T*log(EPS)/(C*ϵ))
    #initial point
    j0 = exp(-K/λr)
   #taylor expansion 
   Δλ = λr - λa
    b = 1 - λr
    j0Δλ = j0^Δλ
    c = -log(1 - j0Δλ) + λr*log(j0) + K
    b = (λr - (Δλ*j0Δλ)/(j0Δλ - 1))/j0 
    a = ((Δλ*j0Δλ*(j0Δλ + Δλ - 1))/(j0Δλ - 1)^2 - λr)/(2*j0*j0) 
    y = (-b + sqrt(b*b - 4*a*c))/(2*a)
    j = y + j0
    _u(x) = C*ϵ*(x^-λr - x^-λa)#, j = r/σ
    _f(x) = exp(-_u(x)/T)
    xinf0 = Roots.find_zero(_f,j)
    println("T: $T")
    println("solution (rootfinding): $xinf0")
    println("solution (initial point): $j0")
    println("solution (initial point + quadratic taylor interp): $j")
    println("exp(-u/T) (rootfinding): $(_f(xinf0))")
    println("exp(-u/T) (initial point): $(_f(j0))")
    println("exp(-u/T) (initial point + quadratic taylor interp): $(_f(j))")
end

This is pretty much what happens after one step of Halley's method (the rootfinding can be expressed as a rootfinding problem for a Taylor expansion). So I think your original approach is probably best; use the guess value, do a step or two of Halley's method, you're done.

[ianhbell]

This is what I settled on below, to be implemented in teqp. The key point is to use your starting value, and do root finding in $log(f(j))-log(EPS)$, which turns into $-u/T-log(EPS)$. I used my multicomplex derivative package here for derivatives, but I plan to switch to analytical ones in teqp.

The only mystery I have now is how you got the initial value for j0. It does indeed seem to work well.

Code:

# Here is a small example of using adaptive quadrature 
# to obtain the quasi-exact value of d for ethane
# according to the pure-fluid parameters given in 
# Lafitte et al.

epskB = 206.12 # [K]
sigma_m = 3.7257e-10 # [m]
lambda_r = 12.4 
lambda_a = 6.0
C = lambda_r/(lambda_r-lambda_a)*(lambda_r/lambda_a)**(lambda_a/(lambda_r-lambda_a))
T = 50.0 # [K]

# The classical method based on adaptive quadrature
def integrand(r_m):
    u = C*epskB*((sigma_m/r_m)**(lambda_r) - (sigma_m/r_m)**(lambda_a))
    return 1-np.exp(-u/T)

print('quasi-exact; (value, error estimate):')
exact, exact_error = scipy.integrate.quad(integrand, 0.0, sigma_m)
print(exact, exact_error)

EPS = np.finfo(float).eps

def f(j):
    u = C*epskB*((j)**(-lambda_r) - (j)**(-lambda_a))
    return np.exp(-u/T)

def res(j):
    return np.log(f(j))-np.log(EPS)

import multicomplex
def ders(j):
    return multicomplex.diff_mcx1(res, j, 2, True)

x = np.linspace(0.0000001*sigma_m, sigma_m, 10000)
y = integrand(x)
plt.plot(x/sigma_m, 1-y)
plt.yscale('log')

K = np.log(-T/epskB*np.log(EPS)/C)
j0 = np.exp(-K/lambda_r)

res = scipy.optimize.root_scalar(ders, x0=j0,fprime=True,fprime2=True, method='Halley')
print('Rootfinding')
print(res)
print(f(j0), j0)
print(f(res.root), res.root)

plt.axvline(j0, dashes=[2,2])
plt.axvline(res.root)
plt.close()

r0 = res.root
print('and now with fixed quadrature')
for n in [5,6,7,8,9,10,20,40]:
    val = r0*sigma_m + scipy.integrate.fixed_quad(integrand, r0*sigma_m, sigma_m, n=n)[0]
    print(n, val, val/exact-1)

yielding

quasi-exact; (value, error estimate):
3.6940193410704346e-10 5.8624831636284655e-12
Rootfinding
      converged: True
           flag: 'converged'
 function_calls: 4
     iterations: 4
           root: 0.88880166306088
3.67276723019937e-06 0.9355323202328244
2.22044604925035e-16 0.88880166306088
and now with fixed quadrature
5 3.6939698784422905e-10 -1.3389921269313376e-05
6 3.6940240893171906e-10 1.2853876272256315e-06
7 3.6940193069488347e-10 -9.236984621630029e-09
8 3.694019315295026e-10 -6.977605160329858e-09
9 3.6940193571500047e-10 4.3528656235025665e-09
10 3.694019351113992e-10 2.7188697959701358e-09
20 3.694019351651498e-10 2.8643767358005334e-09
40 3.694019351651499e-10 2.8643771798897433e-09

[ianhbell]

@longemen3000 now I see the logic; you consider just the repulsive part of the potential to get the j0

[ianhbell]

All set now. I did as we discussed, thanks for your help @longemen3000 : https://github.com/usnistgov/teqp/blob/72ed1188c30d47d5d06a62fe672d935020dcc375/include/teqp/models/saftvrmie.hpp#L255-L336

That works well, and converges nicely, but the implementation in Clapeyron.jl is still very competitive on accuracy without as many points in the quadrature and also ignoring the problem of breaking the integral into two pieces.

[longemen3000]

reading the paper that contains the Clapeyron.jl quadrature, it seems that the best approach would be: if T/ϵ < threshold then: use gauss-laguerre quadrature (Clapeyron.jl current) else use gauss-legendre + assen (teqp current)

i need to explore this numerically, to find a proper threshold (maybe T/ϵ = 10?)

[ianhbell]

Check out the checks in : https://github.com/usnistgov/teqp/blob/master/notebooks/Check%20SAFT-VR-Mie%20Implementations.ipynb

@longemen3000 which paper are you referring to?

[longemen3000]

i think i reverse-engineered the integral. basically is a translated gauss-laguerre: if we take:

θ = C*ϵ/T
f(j) = exp(-θ*(j^-λr - j^-λa)) 

#= 
change of variables: 
x = j^(-λr), j = x^(-1/λr)
j -> 0, x -> Inf
j -> 1, x -> 1, we need to to invert the limits, and multiply by -1
dx = -λr*j^(-λr - 1)
dj = -x^(-1/λr - 1)/λr
=# 

f2(x) = x^(-1/λr - 1)/λr * exp(-θ*(x - x^λa/λr)) = x^(-1/λr - 1)/λr * exp(θ*x^λa/λr))*exp(-θx)

that function form conforms to a translated gauss-laguerre with r = 1 and a = θ. the integral in Clapeyron can then be rewrited as:

∑fi = 0.0
θ0 = exp(-θ)/θ #constant factor for quadrature
for j ∈ 1:5
    xj = u[j]/θ + 1 #translated expression
    fuj = xj^(-1/λr - 1)*exp(θ*xj^(λa/λr))/λr #f(x) evaluated at xj
    ∑fi += θ0 * w[j] * fuj #quadrature
end

(the function above is functionally equivalent as the one currently present in Clapeyron.jl) applying the assen strategy here seems a little harder, maybe expressing the integral as this?:

integral(f2,1,inf) - integral(f2,xcut,inf), xcut = rcut^-λr

[ianhbell]

Could we just crank up the number of Gauss-Laguerre nodes a bit? Would that be a good enough solution?

[longemen3000]

I will try. the ideal way would be to use generalized gauss laguerre, that integrates functions of the form x^a exp(-x) f(x), but that would require calculating weights.

On the assen's cut, another way would be to use false position on this function:

rcut(r) = exp((-K + log1p(-r^(λr - λa)))/λr)

i also got some problems of finding the cut when T/ϵ is really low.

another misc data, is the limit of the diameter when T -> 0. assen has an analytical expression for d, maybe it could be useful:

d_T0 = σ*cbrt(1 + 3*T/((λr - λa)*C*ϵ))

Finally, i doubt that this will affect the specific tests for now, but it looks like the notebook on where the tests are running is using an old version of Clapeyron (0.3.11), you could update the local package enviroment with Pkg.update("Clapeyron")

[ianhbell]

Yes I also have problems with getting the cut when T/(epsilon/k) is very small. I would prefer to avoid the cut algorithm for this reason. Also, rootfinding is a generally unreliable method.

[longemen3000]

did a comparison. i used integration over 256bit BigFloats and adaptative quadrature for reference results, also i cranked up the results up until T = 20000 (or T/epsilon = 100, cases that could be present with electrolytes).

While using 10 points improves the results (impressively, getting to Float64 accuracy on T/eps < 0.5), it just only gets marginally better at high temps. so the assen cut seems better

T     d(exact)            d(laguerre 5)       d(laguerre 10)      err%(laguerre 5) err%(laguerre 10)
   50 3.69401935165149808 3.69401935371413526 3.69401935165149808 5.584e-08         0.000e+00
   60 3.68859851870539179 3.68859851960730989 3.68859851870539179 2.445e-08         0.000e+00
   70 3.68340042216842489 3.68340042257674849 3.68340042216842489 1.109e-08         0.000e+00
   80 3.67840285128657563 3.67840285147834445 3.67840285128657563 5.213e-09         0.000e+00
   90 3.67358718756225944 3.67358718765633752 3.67358718756225944 2.561e-09         0.000e+00
  100 3.66893761045343059 3.66893761050484679 3.66893761045343059 1.401e-09         0.000e+00
  120 3.66008411989695759 3.66008411997046190 3.66008411989696603 2.008e-09         2.220e-13
  140 3.65175326575490100 3.65175326619121554 3.65175326575500980 1.195e-08         2.975e-12
  160 3.64387588507065185 3.64387588711970967 3.64387588507150184 5.623e-08         2.334e-11
  180 3.63639690468588173 3.63639691196848602 3.63639690469049004 2.003e-07         1.267e-10
  195 3.63102187995729553 3.63102189634933437 3.63102187997095838 4.514e-07         3.763e-10
  200 3.62927155930897483 3.62927158032592967 3.62927155932805068 5.791e-07         5.256e-10
  205 3.62754099900295346 3.62754102568787395 3.62754099902924976 7.356e-07         7.249e-10
  250 3.61277687034389849 3.61277703120265281 3.61277687064171049 4.452e-06         8.243e-09
  300 3.59783859272094952 3.59783929656464307 3.59783859494307112 1.956e-05         6.176e-08
  350 3.58415640219504894 3.58415858784754393 3.58415641268620888 6.098e-05         2.927e-07
  400 3.57151392390723776 3.57151932743030498 3.57151396035398694 1.513e-04         1.020e-06
  500 3.54873703186315659 3.54875831750520465 3.54873727434672226 5.998e-04         6.833e-06
  600 3.52859476394662597 3.52865278580458153 3.52859573989969055 1.644e-03         2.766e-05
  700 3.51049915509707899 3.51062501937698324 3.51050202546950807 3.585e-03         8.177e-05
  800 3.49404617578916987 3.49428043063848737 3.49405301000310420 6.704e-03         1.956e-04
  900 3.47894468746539287 3.47933579620028066 3.47895868858363366 1.124e-02         4.025e-04
 1000 3.46497713222883164 3.46557971204631876 3.46500279148203649 1.739e-02         7.405e-04
 2000 3.36300309864212377 3.36911230405509299 3.36367374976556155 1.817e-01         1.994e-02
 5000 3.20596949733352599 3.25130034457504546 3.21763689018254784 1.414e+00         3.639e-01
 7000 3.14365501060091512 3.21989229481234052 3.16828308216813204 2.425e+00         7.834e-01
10000 3.07571140772415985 3.19608033897958554 3.12321929457141367 3.914e+00         1.545e+00
15000 2.99678065129743443 3.18063616540862215 3.08399046947168021 6.135e+00         2.910e+00
20000 2.94006254901040842 3.17627436025302590 3.06465771021096467 8.034e+00         4.238e+00

Code (Clapeyron, master branch):

using Printf #pretty printing, available in base
using QuadGK #for integration with arbitrary precision
function d_vrmie_experiment(λa,λr,σ,ϵ)
    C = Clapeyron.Cλ_mie(λa, λr)
    λrinv = 1/λr
    λaλr = λa/λr
    f_laguerre(x,θ) = x^(-λrinv - 1)*exp(θ*x^(λaλr))*λrinv
    f_normal(r,θ) = exp(-θ*(r^-λr - r^-λa))

    println("T     d(exact)            d(laguerre 5)       d(laguerre 10)      err%(laguerre 5) err%(laguerre 10)")
    Ts = [50, 60, 70, 80, 90, 100, 120, 140, 160, 180, 195, 200, 205, 250, 300, 350, 400, 500, 600, 700, 800, 900, 1000,2000,5000,7000,10000,15000,20000]
    for T in Ts
        Tx = T/ϵ
        θ = C/Tx
        fi_normal = r -> f_normal(r,θ)
        ∑fi_exact = Float64(QuadGK.quadgk(fi_normal, BigFloat(0.), BigFloat(1.))[1])
        fi = x -> f_laguerre(x,θ)
        ∑fi5 = Clapeyron.Solvers.laguerre5(fi,θ,1.)
        ∑fi10 = Clapeyron.Solvers.laguerre10(fi,θ,1.)
        di_5 = σ*(1-∑fi5)*1e10
        di_10 = σ*(1-∑fi10)*1e10
        di_exact = σ*(1-∑fi_exact)*1e10
        err5 = abs(di_5/di_exact-1.0)*100
        err10 = abs(di_10/di_exact-1.0)*100
        Printf.@printf("%5i %.17f %.17f %.17f %.3e         %.3e", T,di_exact,di_5,di_10,err5,err10)
        println()
    end
end

function d_vrmie_experiment(model::EoSModel)
    ϵ = model.params.epsilon[1]
    σ = model.params.sigma[1]
    λa = model.params.lambda_a[1]
    λr = model.params.lambda_r[1]
    return d_vrmie_experiment(λa,λr,σ,ϵ)
end

model = SAFTVRMie(["ethane"])
d_vrmie_experiment(model)

[ianhbell]

Thanks for the test! Seems like 10 point Gauss-Laguerre below T^ of 100, and split + Gauss-Legendre with 10 points above T^=100 would be a good choice.