Merge 60af9c2 into 68c694d

chemicals/rachford_rice.py

@@ -43,6 +43,7 @@
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_LN2
 .. autofunction:: chemicals.rachford_rice.Li_Johns_Ahmadi_solution
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_polynomial
+.. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_mpmath
 
 
 Three Phase
@@ -69,13 +70,14 @@
            'Rachford_Rice_solution2', 'Rachford_Rice_solutionN',
            'Rachford_Rice_flashN_f_jac', 'Rachford_Rice_flash2_f_jac',
            'Li_Johns_Ahmadi_solution', 'flash_inner_loop',
-           'flash_inner_loop_all_methods', 'flash_inner_loop_methods']
+           'flash_inner_loop_all_methods', 'flash_inner_loop_methods',
+           'Rachford_Rice_solution_mpmath']
 
 from fluids.constants import R
 from fluids.numerics import IS_PYPY, one_epsilon_larger, one_epsilon_smaller, NotBoundedError, numpy as np
 from fluids.numerics import newton_system, roots_cubic, roots_quartic, secant, horner, brenth, newton, linspace, horner_and_der, halley, solve_2_direct, py_solve, solve_3_direct, solve_4_direct
 from chemicals.utils import exp, log
-from chemicals.utils import normalize, mark_numba_uncacheable
+from chemicals.utils import normalize, mark_numba_uncacheable, mark_numba_incompatible
 from chemicals.exceptions import PhaseCountReducedError
 try:
     from itertools import combinations
@@ -877,6 +879,157 @@ def Rachford_Rice_solution_numpy(zs, Ks, guess=None):
 #    return V_over_F, xs, ys # numba: uncomment
     return float(V_over_F), xs.tolist(), ys.tolist() # numba: delete
 
+@mark_numba_incompatible
+def Rachford_Rice_solution_mpmath(zs, Ks, dps=200, tol=1e-100):
+    r'''Solves the the Rachford-Rice flash equation using numerical 
+    root-finding to a high precision using the `mpmath` library.
+    
+    .. math::
+        \sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0
+
+    Parameters
+    ----------
+    zs : list[float]
+        Overall mole fractions of all species, [-]
+    Ks : list[float]
+        Equilibrium K-values, [-]
+    dps : int, optional
+        Number of decimal places to use in the intermediate values of the 
+        calculation, [-]
+    tol : float, optional
+        The tolerance of the solver used in `mpmath`, [-]
+
+    Returns
+    -------
+    L_over_F : float
+        Liquid fraction solution [-]
+    V_over_F : float
+        Vapor fraction solution [-]
+    xs : list[float]
+        Mole fractions of each species in the liquid phase, [-]
+    ys : list[float]
+        Mole fractions of each species in the vapor phase, [-]
+
+    Notes
+    -----
+    This function is written solely for development purposes with the aim
+    of returning bit-accurate solutions.
+    
+    Note that the liquid fraction is also returned; it is insufficient to 
+    compute it as :math:`\frac{L}{F} = 1 - \frac{V}{F}`.
+
+    Examples
+    --------
+    >>> Rachford_Rice_solution_mpmath(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
+    (0.3092697372261456, 0.6907302627738544, [0.33940869696634357, 0.3650560590371706, 0.29553524399648584], [0.5719036543882889, 0.27087159580558057, 0.15722474980613046])
+    >>> Rachford_Rice_solution_mpmath(zs=[0.999999999999, 1e-12], Ks=[2.0, 1e-12])
+    (1e-12, 0.999999999999, [0.49999999999975003, 0.50000000000025], [0.9999999999995001, 5.0000000000025e-13])
+    '''
+    # extremely important to validate high decimal precision with mpmath
+    # numerical issues make this an open research problem with respect to maintaining speed
+    from mpmath import mp, mpf, findroot
+
+    mp.dps = dps
+    N = len(zs)
+
+    def make_objf(zs_k_minus_1, K_minus_1):
+        def Rachford_Rice_err(V_over_F):
+            # Compute the objective function
+            err = 0
+            for i in range(len(zs_k_minus_1)):
+                err += zs_k_minus_1[i]/(1 + V_over_F*K_minus_1[i])
+            return err
+        return Rachford_Rice_err
+
+    def make_objf_der(zs_k_minus_1, zs_k_minus_1_2, K_minus_1):
+        def Rachford_Rice_err_fprime(V_over_F):
+            # Compute the objective function and its derivative w.r.t. V_over_F
+            err0, err1 = 0, 0
+            for num0, num1, Kim1 in zip(zs_k_minus_1, zs_k_minus_1_2, K_minus_1):
+                VF_kim1_1_inv = 1/(1 + V_over_F*Kim1)
+                err0 += num0*VF_kim1_1_inv
+                err1 += num1*VF_kim1_1_inv*VF_kim1_1_inv
+            return err0, err1
+        return Rachford_Rice_err_fprime
+
+    zs_orig, Ks_orig = zs, Ks
+
+    solve_order = 1 # 1 for newton, 0 for secant - development only, should always return the correct answer
+
+    Kmin = mpf(min(Ks))
+    Kmax = mpf(max(Ks))
+    if Kmin > 1.0 or Kmax < 1.0:
+        raise PhaseCountReducedError("For provided K values, there is no positive-composition solution; Ks=%s" % (Ks))
+
+    z_of_Kmax = mpf(zs[Ks.index(Kmax)])
+    V_over_F_min = ((Kmax-Kmin)*z_of_Kmax - (1- Kmin))/((1- Kmin)*(Kmax- 1))
+    V_over_F_max = 1/(1 - Kmin)
+
+    # Attempt to avoid zero division errors at either boundary
+    V_over_F_min += V_over_F_min*tol
+    V_over_F_max -= V_over_F_max*tol
+
+    zs_mp = [mpf(i) for i in zs]
+    Ks_mp = [mpf(i) for i in Ks]
+    Ks_minus_1 = [Ki - 1 for Ki in Ks_mp]
+
+    zs_k_minus_1 = [zs_mp[i]*Ks_minus_1[i] for i in range(N)]
+
+    zs_k_minus_1_2 = [0.0]*N
+    for i in range(N):
+        zs_k_minus_1_2[i] = -zs_k_minus_1[i]*Ks_minus_1[i]
+
+    if solve_order == 0:
+        objf = make_objf(zs_k_minus_1, Ks_minus_1)
+    elif solve_order == 1:
+        objf = make_objf_der(zs_k_minus_1, zs_k_minus_1_2, Ks_minus_1)
+
+    guess = None
+    try:
+        # Try to obtain an initial guess for speed and convergence from the
+        # other solvers
+        guess, _, _ = flash_inner_loop(zs_orig, Ks_orig)
+    except:
+        pass
+    if guess is None or (guess < V_over_F_min or guess > V_over_F_max):
+        # Check the guess to ensure it is within bounds
+        guess = 0.5*(V_over_F_min + V_over_F_max)
+
+    # mpmath does not have an option to respect boundaries, so its solvers
+    # often go outside the correct range
+    # try:
+    #     V_over_F = findroot(objf, guess, tol=tol, solver='newton')
+    # except:
+    if N == 2:
+        # A peculiar amount of numerical issues come up in binaries with the solvers
+        # Fortunately, the analytical solution is easy to compute and we have lots of
+        # decimals!
+        z1, z2 = zs_mp
+        K1, K2 = Ks_mp
+        z1z2 = z1 + z2
+        K1z1 = K1*z1
+        K2z2 = K2*z2
+        t1 = z1z2 - K1z1 - K2z2
+        den = t1 + K2*K1z1 + K1*K2z2 - K1*z2 - K2*z1
+        V_over_F = t1/den
+    else:
+        # Requiring a low ytol seems to guarantee a correct solution
+        if solve_order == 0:
+            p1 = 0.4*V_over_F_min + 0.6*V_over_F_max
+            V_over_F = secant(objf, guess, low=V_over_F_min, high=V_over_F_max, xtol=tol, maxiter=1000, ytol=1e-50, x1=p1)
+        elif solve_order == 1:
+            V_over_F = newton(objf, guess, low=V_over_F_min, high=V_over_F_max, xtol=tol, maxiter=1000, ytol=1e-50, fprime=True)
+
+    xs = [0]*N
+    ys = [0]*N
+
+    for i in range(N):
+        x_calc = zs_mp[i]/(1 + V_over_F*Ks_minus_1[i])
+        xs[i] = float(x_calc)
+        ys[i] = float(x_calc*Ks_mp[i])
+    return float(1-V_over_F), float(V_over_F), xs, ys
+
+
 
 def Rachford_Rice_err_LN2(y, zs, cis_ys, x0, V_over_F_min, N):
     x1 = exp(-y)

tests/test_rachford_rice.py

@@ -29,7 +29,7 @@
 from chemicals.exceptions import PhaseCountReducedError
 from fluids.constants import calorie, R
 from chemicals.rachford_rice import *
-from chemicals.rachford_rice import Rachford_Rice_solution_numpy
+from chemicals.rachford_rice import Rachford_Rice_solution_numpy, Rachford_Rice_solution_mpmath
 from chemicals.rachford_rice import Rachford_Rice_valid_solution_naive, Rachford_Rice_solution2
 from chemicals.rachford_rice import Rachford_Rice_flash2_f_jac, Rachford_Rice_flashN_f_jac
 from fluids.numerics import isclose, assert_close, assert_close1d, assert_close2d, normalize
@@ -38,38 +38,27 @@
 
 is_pypy = 'PyPy' in sys.version
 
-def RR_solution_mpmath(zs, Ks, dps=200):
-    # extremely important to validate high decimal precision with mpmath
-    # numerical issues make this an open research problem with respect to maintaining speed
-    from mpmath import mp, mpf, findroot
-
-    def make_objf(zs_k_minus_1, K_minus_1):
-        def Rachford_Rice_err(V_over_F):
-            err = 0
-            for i in range(len(zs_k_minus_1)):
-                err += zs_k_minus_1[i]/(1 + V_over_F*K_minus_1[i])
-            return err
-        return Rachford_Rice_err
-
-    mp.dps = dps
-    N = len(zs)
-    zs_mp = [mpf(i) for i in zs]
-    Ks_mp = [mpf(i) for i in Ks]
-    Ks_minus_1 = [Ki - 1 for Ki in Ks_mp]
-
-    zs_k_minus_1 = [zs_mp[i]*Ks_minus_1[i] for i in range(N)]
-    objf = make_objf(zs_k_minus_1, Ks_minus_1)
-
-    V_over_F = findroot(objf, .5, tol=1e-100)
-    xs = [0]*N
-    ys = [0]*N
-
-    for i in range(N):
-        xs[i] = float(zs[i]/(1 + V_over_F*Ks_minus_1[i]))
-        ys[i] = float(xs[i]*Ks[i])
-    return float(V_over_F), xs, ys
+def RR_solution_mpmath(zs, Ks, guess=None, dps=200):
+    LV, VF, xs, ys = Rachford_Rice_solution_mpmath(zs, Ks, dps=dps) 
+    return (VF, xs, ys)
 
 
+def test_RR_mpmath_points():
+    # points from RR contest update
+    zs = [0.003496418103652993, 0.08134996284399883, 0.08425781368698183, 0.08660015587158083, 0.03393676648390093, 0.046912148366909906, 0.04347295855013991, 0.03528846893106793, 0.06836282476823487, 0.06778643033352585, 0.014742967176819971, 0.07489046659005885, 0.04595208887032791, 0.053716354661293896, 0.022344804838081954, 0.0490901205575939, 0.009318396845552981, 0.06683329012632486, 0.07237858894810185, 0.03528438046562893, 0.003984592980220992] 
+    Ks = [0.000160672721667, 0.018356845386159, 0.094445181723264, 0.117250574977987, 0.053660474066903, 0.53790315308842, 0.026109136837427, 0.106588294438016, 0.013998763838654, 0.162417382603121, 0.007558903680426, 0.064534061436341, 0.731006576434576, 0.005441443070815, 0.015600385380209, 0.012020711491011, 0.238827317565149, 0.022727947741144, 0.001778015519249, 0.007597040584824, 5.87425090047782]    
+    LF, VF, xs, ys = Rachford_Rice_solution_mpmath(zs, Ks)
+    assert_close(VF, -0.19984637297628843, rtol=1e-16)
+
+
+    zs = [0.004157284954252984, 0.0014213285047209943, 0.014848717759705941, 0.03044981105958088, 0.11359772999756554, 0.10729135314739957, 0.029101839157508882, 0.06599513817538073, 0.010324751675473958, 0.04473237952961382, 0.06963277351376072, 0.08645767360282165, 0.08930957642684664, 0.0251765398366809, 0.01594252170408694, 0.09571836326377162, 0.04078154182757284, 0.08243486397206067, 0.013175199166126948, 0.029268130397145882, 0.030182482327921877]
+    Ks = [0.723423674540192, 1.16174499507754, 1.17547301996336, 1.03037454150302, 1.13879814918106, 1.70014675310204, 1.05877209838981, 1.28326927930614, 1.16351944743047, 1.41069373097001, 1.00536461101331, 1.77523784049187, 1.51071399916121, 1.14835523974627, 1.60369203193725, 1.29606646034495, 1.17053941125983, 1.03935905847522, 1.15389315235544, 1.09443270044454, 1.93013701123472]
+    xs_expect = [0.44904982468349636, 0.0008999183019312609, 0.009117623147622485, 0.02746178684472478, 0.07587355807545791, 0.030584404973243343, 0.02404054965205764, 0.032756525185209176, 0.006510943046530854, 0.018101680695012657, 0.0683198757920764, 0.02289038644790009, 0.03156415667512449, 0.01643985772271979, 0.005041074553077012, 0.046452621333235224, 0.02531599345252268, 0.07224850132659927, 0.00849316591461233, 0.021870066188561105, 0.006967485988285152]
+    ys_expect = [0.32485327422416393, 0.001045475583247321, 0.01071752001622364, 0.028295926028986965, 0.08640466750811313, 0.05199797681081755, 0.02545346320155348, 0.0420354424669968, 0.00757560885575084, 0.02553592747647521, 0.06868638535017856, 0.040635880205794526, 0.04768441336082832, 0.01887879675656845, 0.008084331093171238, 0.0602056845051105, 0.029633368071373612, 0.07509213431505989, 0.009800205990689797, 0.02393531559764776, 0.013448202581248493]
+    LF, VF, xs, ys = Rachford_Rice_solution_mpmath(zs, Ks)
+    assert_close(VF, 3.582165028608595, rtol=1e-16)
+    assert_close1d(xs, xs_expect, rtol=1e-16)
+    assert_close1d(ys, ys_expect, rtol=1e-16)
 
 def assert_same_RR_results(zs, Ks, f0, f1, rtol=1e-9):
     N = len(zs)
@@ -242,11 +231,13 @@ def test_flash_inner_loop_methods():
 flash_inner_loop_poly = lambda zs, Ks, guess=None: flash_inner_loop(zs=zs, Ks=Ks, guess=guess, method='Rachford-Rice (polynomial)')
 flash_inner_loop_LN2 = lambda zs, Ks, guess=None: flash_inner_loop(zs=zs, Ks=Ks, guess=guess, method='Leibovici and Nichita 2')
 
+
 algorithms = [Rachford_Rice_solution, Li_Johns_Ahmadi_solution,
               flash_inner_loop, flash_inner_loop_secant,
               flash_inner_loop_NR, flash_inner_loop_halley,
               flash_inner_loop_numpy, flash_inner_loop_LJA,
-              flash_inner_loop_poly, flash_inner_loop_LN2]
+              flash_inner_loop_poly, flash_inner_loop_LN2,
+              RR_solution_mpmath]
 
 @pytest.mark.parametrize("algorithm", algorithms)
 @pytest.mark.parametrize("array", [False])