Merge pull request #185 from Liozou/cubicrealroots

Fix detection of real cubic roots

Project.toml

@@ -1,7 +1,7 @@
 name = "Clapeyron"
 uuid = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
 authors = ["Hon Wa Yew <yewhonwa@gmail.com>", "Pierre Walker <pwalker@mit.edu>", "Andrés Riedemann <andres.riedemann@gmail.com>"]
-version = "0.4.11"
+version = "0.4.12"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"

src/models/cubic/equations.jl

@@ -215,24 +215,14 @@ function crit_pure_tp(model::ABCCubicModel)
 end
 
 function volume_impl(model::ABCubicModel,p,T,z=SA[1.0],phase=:unknown,threaded=false,vol0=nothing)
-    lb_v   =lb_volume(model,z)
+    lb_v = lb_volume(model,z)
     if iszero(p)
         vl,_ = zero_pressure_impl(model,T,z)
         return vl
     end
     nRTp = sum(z)*R̄*T/p
     _poly,c̄ = cubic_poly(model,p,T,z)
-    sols = Solvers.roots3(_poly)
-    function imagfilter(x)
-        absx = abs(imag(x))
-        return absx < 8 * eps(typeof(absx))
-    end
-    x1, x2, x3 = sols
-    sols = (x1, x2, x3)
-    xx = (x1, x2, x3)
-    isreal = imagfilter.(xx)
-    vvv = extrema(real.(xx))
-    zl,zg = vvv
+    num_isreal, zl, zg = Solvers.real_roots3(_poly)
     vvl,vvg = nRTp*zl,nRTp*zg
     #err() = @error("model $model Failed to converge to a volume root at pressure p = $p [Pa], T = $T [K] and compositions = $z")
     if !isfinite(vvl) && !isfinite(vvg) && phase != :unknown
@@ -241,19 +231,12 @@ function volume_impl(model::ABCubicModel,p,T,z=SA[1.0],phase=:unknown,threaded=f
         #isnan(v) && err()
         return v
     end
-    if sum(isreal) == 3 #3 roots
+    if num_isreal == 3 # 3 real roots
         vg = vvg
         _vl = vvl
         vl = ifelse(_vl > lb_v, _vl, vg) #catch case where solution is unphysical
-    elseif sum(isreal) == 1
-        i = findfirst(imagfilter, sols)::Int
-        vl = real(sols[i]) * nRTp
-        vg = real(sols[i]) * nRTp
-    elseif sum(isreal) == 0
-        V0 = x0_volume(model, p, T, z; phase)
-        v = _volume_compress(model, p, T, z, V0)
-        #isnan(v) && err()
-        return v
+    else # 1 real root (or 2 with the second one degenerate)
+        vg = vl = zl * nRTp
     end
 
     function gibbs(v)

src/modules/solvers/poly.jl

@@ -33,7 +33,7 @@ function solve_cubic_eq(poly::NTuple{4,T}) where {T<:Real}
 # copied from PolynomialRoots.jl, adapted to be AD friendly
 # Cubic equation solver for complex polynomial (degree=3)
 # http://en.wikipedia.org/wiki/Cubic_function   Lagrange's method
-# poly = (a,b,c,d) that represents ax3 + bx2 + cx2 + d 
+# poly = (a,b,c,d) that represents a + bx + cx3 + dx4
 
 _1 = one(T)
 third = _1/3
@@ -77,4 +77,44 @@ end
 function roots3(a,b,c,d) 
 x = (a,b,c,d)
 return roots3(x)
-end
+end
+
+"""
+    real_roots3(pol::NTuple{4,T}) where {T<:Real}
+
+Given a cubic real polynom of the form `pol[1] + pol[2]*x + pol[3]*x^2 + pol[4]*x^3`,
+return `(n, zl, zg)` where `n` is the number of real roots and:
+- if `n == 1`, `zl` and `zg` are equal to the only real root, the other two are complex.
+- if `n == 2`, `zl` is the single real root and `zg` is the double (degenerate) real root.
+- if `n == 3`, `zl` is the lowest real root and `zg` the greatest real root.
+
+!!! info
+    If there is a single root triply degenerate, e.g. with `pol == (1,3,3,1)` corresponding
+    to `(x+1)^3`, this will return `n == 2` and `zl == zg` equal to the root.
+"""
+function real_roots3(pol::NTuple{4,T}) where {T<:Real}
+    z1, z2, z3 = sort(roots3(pol); by=real)
+    x1, x2, x3 = real(z1), real(z2), real(z3)
+    mid12 = (x1 + x2)/2
+    mid23 = (x2 + x3)/2
+    between1 = evalpoly(mid12, pol)
+    between2 = evalpoly(mid23, pol)
+    if abs(between1) <  eps(typeof(between1))
+        (2, x3, mid12) # first the single root, then the double root
+    elseif abs(between2) < eps(typeof(between2))
+        (2, x1, mid23) # first the single root, then the double root
+    else
+        sign1 = signbit(between1)
+        sign2 = signbit(between2)
+        if sign1 == sign2 # only one root
+            if sign1 ⊻ (pol[4] > 0)
+                (1, x1, x1)
+            else
+                (1, x3, x3)
+            end
+        else # three distinct roots
+            (3, x1, x3)
+        end
+    end
+end
+real_roots3(a,b,c,d) = real_roots3((a,b,c,d))

test/test_methods_api.jl

@@ -435,7 +435,8 @@ end
         T = 202.694
         v0 = [-4.136285855713797, -4.131888756537859, 0.9673991775701574, 0.014192499147585259, 0.014746430039492817, 0.003661893242764558]
         model = PCSAFT(["methane","butane","isobutane","pentane"])
-        @test_broken bubble_pressure(model,T,x;v0 = v0)[1] ≈ 5.913118531569793e6 rtol = 1e-4
+        # @test_broken bubble_pressure(model,T,x;v0 = v0)[1] ≈ 5.913118531569793e6 rtol = 1e-4
+        # FIXME: The test does not yield the same value depending on the OS and the julia version
     end
 end
 

test/test_methods_eos.jl

@@ -298,6 +298,13 @@ end
     @test vc/system.params.Vc[1] ≈ 1.168 rtol = 1e-4 #if Zc_exp < 0.29, this should hold, by definition
 end
 
+@testset "EPPR78, single component" begin
+    system = EPPR78(["carbon dioxide"])
+    T = 400u"K"
+    @test Clapeyron.volume(system, 3311.0u"bar", T) ≈ 3.363139140634349e-5u"m^3"
+    @test Clapeyron.molar_density(system, 3363.1u"bar", T) ≈ 29810.118127751106u"mol*m^-3"
+end
+
 @testset "Cubic methods, multi-components" begin
     system = RK(["ethane","undecane"])
     p = 1e7

test/test_solvers.jl

@@ -43,7 +43,13 @@ end
         #@test [@inferred(SOL.roots3(SA[-6im, -(3 + 4im), 2im-2, 1.0]))...] ≈ [3, -2im, -1]
         @test @inferred(SOL.roots3(1.0, -3.0, 3.0, -1.0)) ≈ [1, 1, 1]
         @test @inferred(SOL.roots3([1.0, -3.0, 3.0, -1.0])) ≈ [1, 1, 1]
+        @test @inferred(SOL.real_roots3((0.4179831841886642, -16.86256511617483, 1.6508521434627874, 1.0))) == (3, -5.0238891763914015, 3.3481880327202824)
+        @test @inferred(SOL.real_roots3(0.7336510424849684, -17.464843881306653, 1.6925644348171853, 1.0)) == (3, -5.126952070514365, 3.392203575902995)
+        @test @inferred(SOL.real_roots3(-1, 3, -3, 1)) == (2, 1.0, 1.0) # triply degenerate root
+        @test collect(SOL.real_roots3(1, -1, -1, 1)) ≈ [2, -1.0, 1.0] # one doubly degenerate
+        @test SOL.real_roots3(1, 0, 1, 2) == (1, -1.0, -1.0) # only one real root
     end
+
     # A difficult MCP.
     #
     # Presented in Miranda and Fackler (2002): "Applied Computational Economics and