Merge branch 'RRLN' into RRLF

chemicals/numba.py

@@ -43,6 +43,7 @@
      'rachford_rice.Rachford_Rice_polynomial_5',
      'rachford_rice.Rachford_Rice_solution_polynomial',
      'rachford_rice.Rachford_Rice_numpy_err_fprime2',
+     'rachford_rice.Rachford_Rice_solution_Leibovici_Neoschil',
      'rachford_rice.Li_Johns_Ahmadi_solution',
      'rachford_rice._Rachford_Rice_analytical_3',
      'rachford_rice.flash_inner_loop',

chemicals/rachford_rice.py

@@ -42,6 +42,7 @@
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_LN2
 .. autofunction:: chemicals.rachford_rice.Li_Johns_Ahmadi_solution
+.. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_Leibovici_Neoschil
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_polynomial
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_mpmath
 .. autofunction:: chemicals.rachford_rice.Rachford_Rice_solution_binary_dd
@@ -71,10 +72,11 @@
            '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',
-           'Rachford_Rice_solution_mpmath', 'Rachford_Rice_solution_binary_dd']
+           'Rachford_Rice_solution_mpmath', 'Rachford_Rice_solution_binary_dd',
+           'Rachford_Rice_solution_Leibovici_Neoschil']
 
 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 IS_PYPY, one_epsilon_larger, one_epsilon_smaller, one_10_epsilon_larger, one_10_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 fluids.numerics import add_dd, div_dd, mul_dd, mul_noerrors_dd
 from chemicals.utils import exp, log
@@ -85,7 +87,6 @@
 except:
     pass
 
-
 R_inv = 1.0/R
 
 
@@ -881,6 +882,136 @@ 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
 
+def Rachford_Rice_err_fprime_Leibovici_Neoschil(V_over_F, zs_k_minus_1, zs_k_minus_1_2, K_minus_1, V_over_F_min, V_over_F_max):
+    plain_err, plan_diff = 0.0, 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.0/(1. + V_over_F*Kim1)
+        plain_err += num0*VF_kim1_1_inv
+        plan_diff += num1*VF_kim1_1_inv*VF_kim1_1_inv
+    err = (V_over_F - V_over_F_min)*(V_over_F_max - V_over_F)*plain_err
+    fprime = (plan_diff*(-V_over_F + V_over_F_max)*(V_over_F - V_over_F_min) 
+              + plain_err*(-V_over_F + V_over_F_max)
+              + plain_err*(-V_over_F + V_over_F_min))
+    # print(err, V_over_F)
+    return err, fprime
+
+
+@mark_numba_uncacheable
+def Rachford_Rice_solution_Leibovici_Neoschil(zs, Ks, guess=None):
+    r'''Solves the objective function of the Rachford-Rice flash equation as
+    modified by Leibovici and Neoschil. This modification helps
+    convergence near the vapor fraction boundaries only; it slows
+    convergence in other regions.
+
+    .. math::
+        \left(\frac{V}{F} - \alpha_L\right)\left(\alpha_R - \frac{V}{F}\right)
+        \sum_i \frac{z_i(K_i-1)}{1 + \frac{V}{F}(K_i-1)} = 0
+        
+    .. math::
+        \alpha_L = - \frac{1}{K_{max} - 1}
+    
+    .. math::
+        \alpha_R = \frac{1}{1 - K_{min}}
+
+    Parameters
+    ----------
+    zs : list[float]
+        Overall mole fractions of all species, [-]
+    Ks : list[float]
+        Equilibrium K-values, [-]
+    guess : float, optional
+        Optional initial guess for vapor fraction, [-]
+
+    Returns
+    -------
+    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
+    -----
+    The initial guess is the average of the following, as described in [2]_.
+
+    .. math::
+        \left(\frac{V}{F}\right)_{min} = \frac{(K_{max}-K_{min})z_{of\;K_{max}}
+        - (1-K_{min})}{(1-K_{min})(K_{max}-1)}
+
+    .. math::
+        \left(\frac{V}{F}\right)_{max} = \frac{1}{1-K_{min}}
+
+    Examples
+    --------
+    >>> Rachford_Rice_solution_Leibovici_Neoschil(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
+    (0.69073026277385, [0.339408696966343, 0.36505605903717, 0.29553524399648], [0.57190365438828, 0.270871595805580, 0.157224749806130])
+
+    References
+    ----------
+    .. [1] Leibovici, ClaudeF., and Jean Neoschil. "A New Look at the 
+       Rachford-Rice Equation." Fluid Phase Equilibria 74 (July 15, 1992): 
+       303-8. https://doi.org/10.1016/0378-3812(92)85069-K.
+    '''
+    N = len(Ks)
+    Kmin = min(Ks) # numba: delete
+    Kmax = max(Ks)# numba: delete
+    z_of_Kmax = zs[Ks.index(Kmax)]# numba: delete
+
+#    Kmin, Kmax, z_of_Kmax = Ks[0], Ks[0], zs[0] # numba: uncomment
+#    for i in range(N): # numba: uncomment
+#        if Ks[i] > Kmax: # numba: uncomment
+#            z_of_Kmax = zs[i] # numba: uncomment
+#            Kmax = Ks[i] # numba: uncomment
+#        if Ks[i] < Kmin: # numba: uncomment
+#            Kmin = Ks[i] # numba: uncomment
+    if Kmin > 1.0 or Kmax < 1.0:
+        raise PhaseCountReducedError("For provided K values, there is no positive-composition solution; Ks=%s" % (Ks))  # numba: delete
+#        raise PhaseCountReducedError("For provided K values, there is no positive-composition solution") # numba: uncomment
+    V_over_F_min = ((Kmax-Kmin)*z_of_Kmax - (1.- Kmin))/((1.- Kmin)*(Kmax- 1.))
+    V_over_F_min_LN = -1.0/(Kmax-1)
+    V_over_F_max = 1./(1.-Kmin)
+
+    V_over_F_min2 = V_over_F_min
+    V_over_F_max2 = V_over_F_max
+    if guess is not None and guess > V_over_F_min and guess < V_over_F_max:
+        x0 = guess
+    else:
+        x0 = (V_over_F_min2 + V_over_F_max2)*0.5
+
+    K_minus_1 = [0.0]*N
+    zs_k_minus_1 = [0.0]*N
+    for i in range(N):
+        Kim1 = Ks[i] - 1.0
+        K_minus_1[i] = Kim1
+        zs_k_minus_1[i] = zs[i]*Kim1
+
+    zs_k_minus_1_2 = [0.0]*N
+    for i in range(N):
+        zs_k_minus_1_2[i] = -zs_k_minus_1[i]*K_minus_1[i]
+
+    # Right the boundaries, the derivative goes very large and microscopic steps are made and the newton solver switches
+    # The boundaries need to be handled with bisection-style solvers
+    low, high = V_over_F_min*one_10_epsilon_larger, V_over_F_max*one_10_epsilon_smaller
+    V_over_F = newton(Rachford_Rice_err_fprime_Leibovici_Neoschil, x0, xtol=1e-15, fprime=True, high=high,
+                        low=low, bisection=True, args=(zs_k_minus_1, zs_k_minus_1_2, K_minus_1, V_over_F_min_LN, V_over_F_max))
+    if 1-1e-4 < V_over_F < 1+1e-4:
+        # Re-calculate while inverting phase compositions and inverting the K values
+        # This is necessary for the correct calculation for x compositions when the
+        # liquid fraction is very near zero.
+        Ks_inv = [0.0]*N
+        for i in range(N):
+            Ks_inv[i] = 1.0/Ks[i]
+        LF, xs, ys = Rachford_Rice_solution_Leibovici_Neoschil(zs, Ks_inv, guess=1.0-V_over_F)
+        return V_over_F, ys, xs
+
+    xs = zs_k_minus_1
+    ys = K_minus_1
+    for i in range(N):
+        xs[i] = zs[i]/(1. + V_over_F*K_minus_1[i])
+        ys[i] = xs[i]*Ks[i]
+    return V_over_F, xs, ys
+
 def Rachford_Rice_solution_binary_dd(zs, Ks):
     r'''Solves the the Rachford-Rice flash equation for a binary system using
     double-double math. This increases the range in which the

tests/test_rachford_rice.py

@@ -28,6 +28,7 @@
 import sys
 from chemicals.exceptions import PhaseCountReducedError
 from fluids.constants import calorie, R
+from fluids.numerics import derivative
 from chemicals.rachford_rice import *
 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
@@ -112,6 +113,47 @@ def assert_same_RR_results(zs, Ks, f0, f1, rtol=1e-9):
     assert_close1d(zs, zs_recalc2, rtol=rtol)
 
 
+
+def test_Rachford_Rice_solution_Leibovici_Neoschil():
+    # Check the derivative of the objective function
+    from chemicals.rachford_rice import Rachford_Rice_err_fprime_Leibovici_Neoschil
+    args = ([0.3425, -0.0774, -0.0936], [-0.23461250000000003, -0.0199692, -0.0438048], [0.685, -0.258, -0.46799999999999997], 0.3384490610767984, 2.1367521367521367)
+    point = 1 # vapor fraction
+    num_der = derivative(lambda VF, *args: Rachford_Rice_err_fprime_Leibovici_Neoschil(VF, *args)[0], point, order=3, dx=point*8e-6, args=args)
+    implemented_der = Rachford_Rice_err_fprime_Leibovici_Neoschil(point, *args)[1]
+    assert_close(num_der, implemented_der, rtol=1e-10)
+
+    # Single test point
+    VF, xs, ys = Rachford_Rice_solution_Leibovici_Neoschil(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
+    assert_close(VF, 0.6907302627738544, rtol=1e-15)
+    assert_close1d(xs, [0.33940869696634357, 0.3650560590371706, 0.29553524399648584], rtol=1e-15)
+    assert_close1d(ys, [0.5719036543882889, 0.27087159580558057, 0.15722474980613046], rtol=1e-15)
+
+
+    # This test corrected the limit one one side
+    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]    
+
+    VF, xs, ys = Rachford_Rice_solution_Leibovici_Neoschil(zs=zs, Ks=Ks)
+    assert_close(VF, -0.19984637297628843, rtol=1e-15)
+    xs_expect = [0.0029141328046673174, 0.06800825195999576, 0.07134616627394522, 0.07361365636097347, 0.028539335400666326, 0.04294614179555682, 0.03639035607104222, 0.029942400390529558, 0.057109473261243694, 0.05806676028080419, 0.012302868485612144, 0.06309490688185262, 0.0436078461819446, 0.04480997284023247, 0.018671571073170388, 0.04099575210629456, 0.008088060495929661, 0.055913202622192716, 0.06034108330284007, 0.029444673394774787, 0.1538533880157314]
+    ys_expect = [4.68221649024986e-07, 0.001248416966212587, 0.0067383016390009645, 0.008631243534556089, 0.001531434267154102, 0.02310086508481239, 0.0009501207862215337, 0.00319150938900673, 0.0007994620291340757, 0.009431051221051083, 9.299619787569067e-05, 0.004071770597033691, 0.03187762234314891, 0.00024383091621489136, 0.0002912837043954206, 0.0004927981084267735, 0.0019316497925475297, 0.0012707923472371917, 0.00010728738256074635, 0.00022369237878699152, 0.9037734030929736]
+    assert_close1d(xs, xs_expect, rtol=1e-15)
+    assert_close1d(ys, ys_expect, rtol=1e-15)
+
+    # This test avoids the derivative getting very high at the boundary
+    # Probably introducing an issue with very-close-to-solution points now
+    zs = [0.0031837223319269936, 0.0010664802407109979, 0.0477453142702529, 0.0044117514817209914, 0.0018066278358659963, 0.04562721609715891, 0.09478037830858681, 0.011239515827540977, 0.04637631678894391, 0.07213704782208787, 0.01912611083738096, 0.05746885708647488, 0.09174883552566682, 0.03285749312824994, 0.04287268287821692, 0.029304384114075942, 0.02629105012965095, 0.08706123918262182, 0.062288212638427876, 0.07050856380648586, 0.0043157065296399915, 0.02150562264794196, 0.07459766410965485, 0.051679206380714895]
+    Ks = [0.230879675310901, 34.6613844104893, 2.36200912834507, 3.56139443114014, 253.767128654543, 27.4982543972241, 6.49752569557429, 4.86413631545473, 11.4295807067618, 70.922105955787, 39.3998758557032, 15.9210754842061, 6.30574129622118, 1.78196413354547, 4.52681059087712, 292.041527320902, 0.968296775957929, 39.850923110968, 2.1745799861653, 4.54565379548083, 7.65425289667786, 0.834631636599059, 24.7143417566771, 581.818086414345]
+
+    VF, xs, ys = Rachford_Rice_solution_Leibovici_Neoschil(zs=zs, Ks=Ks)
+    assert_close(VF, 1.2951827117920327, rtol=1e-14)
+    xs_expect = [0.8272304659719063, 2.3913376706490657e-05, 0.017273675435485428, 0.0010218363101008369, 5.501644293848875e-06, 0.0012918208231946849, 0.011672028807351659, 0.0018717669067864103, 0.003196556184776442, 0.0007878500294484951, 0.00037698167675803083, 0.0028274238440099183, 0.011655227358251828, 0.016324381309082817, 0.007700023217015088, 7.753468068480375e-05, 0.02741682521219134, 0.0016964704032019513, 0.02470484237272209, 0.012608219942347239, 0.0004486893491095034, 0.027367188538055545, 0.002352169716456059, 6.860689006322322e-05]
+    ys_expect = [0.19099070139087906, 0.0008288707425765132, 0.04080057905868659, 0.0036391621443299096, 0.0013961364753286797, 0.0355228176318389, 0.07583930709525073, 0.009104529385366147, 0.03653529689760093, 0.055875983265816075, 0.014853031264141248, 0.04501562844632608, 0.07349484846977544, 0.029089461995105627, 0.034856546648783605, 0.022643346567528526, 0.026547623459966935, 0.06760591159803384, 0.053722655785089925, 0.057312602835187826, 0.00343438175012992, 0.022841521358632308, 0.05813232624220132, 0.03991672949142387]
+    assert_close1d(xs, xs_expect, rtol=1e-11)
+    assert_close1d(ys, ys_expect, rtol=1e-11)
+
+
 def test_RR_numpy():
     def Wilson_K_value(T, P, Tc, Pc, omega):
         return Pc/P*exp((5.37*(1.0 + omega)*(1.0 - Tc/T)))
@@ -346,7 +388,7 @@ def test_flash_inner_loop_methods():
               flash_inner_loop_NR, flash_inner_loop_halley,
               flash_inner_loop_numpy, flash_inner_loop_LJA,
               flash_inner_loop_poly, flash_inner_loop_LN2,
-              RR_solution_mpmath]
+              RR_solution_mpmath, Rachford_Rice_solution_Leibovici_Neoschil]
 
 @pytest.mark.parametrize("algorithm", algorithms)
 @pytest.mark.parametrize("array", [False])