Merge 92eaf71 into 6744a6f

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
 
 
@@ -629,6 +630,7 @@ def Rachford_Rice_err_fprime2(V_over_F, zs_k_minus_1, zs_k_minus_1_2, zs_k_minus
         err1 += num1*t2
         err2 += num2*t2*VF_kim1_1_inv
 #            print(err0, err1, err2)
+    # print(err0, V_over_F)
     return err0, err1, err2
 
 def Rachford_Rice_err(V_over_F, zs_k_minus_1, K_minus_1):
@@ -880,6 +882,150 @@ 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, d=0):
+    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
+    i_Kmin = Ks.index(Kmin)
+    i_Kmax = Ks.index(Kmax)# numba: delete
+    z_of_Kmax = zs[i_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 d == 0:
+        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, d=1)
+
+    xs = zs_k_minus_1
+    ys = K_minus_1
+    for i in range(N):
+        if d == 0:
+            switch = -1.01 < V_over_F*K_minus_1[i] < -0.99
+            if switch:
+                xs[i] = zs[i]/(LF + (1.0 - LF)*Ks[i])
+            else:
+                # Attempt to avoid truncation error
+                xs[i] = zs[i]/(1. + V_over_F*K_minus_1[i])
+        else:
+            xs[i] = zs[i]/(1. + V_over_F*K_minus_1[i])
+
+
+
+    i_max_x = xs.index(max(xs))
+    x_sum = sum(xs)
+    xs[i_max_x] = xs[i_max_x]  + (1.0-x_sum)
+
+    for i in range(N):
+        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
@@ -1157,7 +1303,7 @@ def Rachford_Rice_err_LN2(y, zs, cis_ys, x0, V_over_F_min, N):
         x7 = zix5*x5x1x6
         dF0 += x7
         ddF0 += x7*(t51 + x5x1x6 + x5x1x6)
-
+    # print(F0, y)
     return F0, -dF0, ddF0
 
 @mark_numba_uncacheable
@@ -1235,7 +1381,7 @@ def Rachford_Rice_solution_LN2(zs, Ks, guess=None):
     Examples
     --------
     >>> Rachford_Rice_solution_LN2(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
-    (0.6907302627738541, [0.3394086969663436, 0.3650560590371706, 0.29553524399648573], [0.571903654388289, 0.27087159580558057, 0.1572247498061304])
+    (0.6907302627738, [0.3394086969663, 0.3650560590371, 0.29553524399648], [0.571903654388, 0.27087159580558, 0.1572247498061])
 
     References
     ----------
@@ -1279,19 +1425,40 @@ def Rachford_Rice_solution_LN2(zs, Ks, guess=None):
 
     # Suggests guess V_over_F_min, not using
     log_transform = (V_over_F_max-guess)/(guess-V_over_F_min)
+    if abs(log_transform-1) < 1e-10:
+        # The derivative is huge around a log transform of 1, better to start at a hard-coded location
+        log_transform = 0.8
     if log_transform > 0.0:
         guess = -log(log_transform)
     else:
         # Case where guess was less than V_over_F_min - nasty
         guess = 0.5*(V_over_F_min + V_over_F_max)
         guess = -log((V_over_F_max-guess)/(guess-V_over_F_min))
     # Should always converge - no poles
+
+    # The function cannot be evaluated in some regions due to a zero division however
+    # ci should be the minimum
+    if V_over_F_max > 0.0:
+        near_high = V_over_F_max*one_epsilon_smaller
+    else:
+        near_high = V_over_F_max*one_epsilon_larger
+    solver_high = -log((V_over_F_max-near_high)/(near_high-V_over_F_min))
+
+    if V_over_F_min < 0.0:
+        near_low = V_over_F_min*one_epsilon_smaller
+    elif V_over_F_min > 0.0:
+        near_low = V_over_F_min*one_epsilon_larger
+    else:
+        # V_over_F_min equals zero case, cannot evaluate there
+        near_low = min(1e-20, V_over_F_max*1e-15)
+    solver_low = -log((V_over_F_max-near_low)/(near_low-V_over_F_min))
+
     try:
 #        V_over_F = halley(Rachford_Rice_err_LN2, guess, xtol=1e-10, args=(zs, cis_ys, x0, V_over_F_min, N)) # numba: uncomment
         if one_m_Kmin == 1.0: # numba: delete
-            V_over_F = newton(Rachford_Rice_err_LN2, guess, fprime=True, fprime2=True, xtol=1e-10, args=(zs, cis_ys, x0, V_over_F_min, N)) # numba: delete
+            V_over_F = newton(Rachford_Rice_err_LN2, guess, fprime=True, fprime2=True, xtol=1e-10, low=solver_low, high=solver_high, args=(zs, cis_ys, x0, V_over_F_min, N)) # numba: delete
         else: # numba: delete
-            V_over_F = newton(Rachford_Rice_err_LN2, guess, fprime=True, fprime2=True, ytol=1e-8, args=(zs, cis_ys, x0, V_over_F_min, N)) # numba: delete
+            V_over_F = newton(Rachford_Rice_err_LN2, guess, fprime=True, fprime2=True, xtol=1.48e-12, low=solver_low, high=solver_high, bisection=True, args=(zs, cis_ys, x0, V_over_F_min, N)) # numba: delete
     except:
 #        return Rachford_Rice_solution(zs=zs, Ks=Ks, fprime=True) # numba: delete
 #        raise ValueError("Could not solve") # numba: uncomment
@@ -1317,12 +1484,23 @@ def Rachford_Rice_solution_LN2(zs, Ks, guess=None):
             #     V_over_F = -log((V_over_F_max - V_over_F_max) / (V_over_F_max - V_over_F_min))
 
     V_over_F = (V_over_F_min + (V_over_F_max - V_over_F_min)/(1.0 + exp(-V_over_F)))
-
-    xs = [0.0]*N
-    for i in range(N):
-        xi = zs[i]/(1.0 + V_over_F*(Ks[i] - 1.0))
-        xs[i] = xi
-        cis_ys[i] = Ks[i]*xi
+
+    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 = cis_ys
+        for i in range(N):
+            Ks_inv[i] = 1.0/Ks[i]
+        LF, xs, ys = Rachford_Rice_solution_LN2(zs, Ks_inv, guess=1.0-V_over_F)
+        return V_over_F, ys, xs
+
+    else:
+        xs = [0.0]*N
+        for i in range(N):
+            xi = zs[i]/(1.0 + V_over_F*(Ks[i] - 1.0))
+            xs[i] = xi
+            cis_ys[i] = Ks[i]*xi
     return V_over_F, xs, cis_ys
 
 def LJA_err(x1, t1, terms_2, terms_3, N2):
@@ -1692,7 +1870,7 @@ def flash_inner_loop(zs, Ks, method=None, guess=None, check=False):
     Examples
     --------
     >>> flash_inner_loop(zs=[0.5, 0.3, 0.2], Ks=[1.685, 0.742, 0.532])
-    (0.6907302627738541, [0.3394086969663436, 0.3650560590371706, 0.29553524399648573], [0.571903654388289, 0.27087159580558057, 0.1572247498061304])
+    (0.6907302627738, [0.3394086969663, 0.3650560590371, 0.29553524399648], [0.571903654388, 0.27087159580558, 0.1572247498061])
     '''
     l = len(zs)
     if method is None:
@@ -1753,6 +1931,8 @@ def flash_inner_loop(zs, Ks, method=None, guess=None, check=False):
         return Rachford_Rice_solution(zs, Ks, fprime=False, fprime2=False, guess=guess)
     elif method2 == FLASH_INNER_ANALYTICAL:
         if l == 2:
+            # _, VF, xs, ys = Rachford_Rice_solution_binary_dd(zs, Ks)
+            # return VF, xs, ys
             z1, z2 = zs
             K1, K2 = Ks
             if (K1 < 1.0 and K2 < 1.0) or (K1 > 1.0 and K2 > 1.0):
@@ -1794,7 +1974,7 @@ def flash_inner_loop(zs, Ks, method=None, guess=None, check=False):
         ys = [0.0]*l
         for i in range(l):
             try:
-                xs[i] = zs[i]/(1.+V_over_F*(Ks[i] - 1.))
+                xs[i] = zs[i]/(1.0 + V_over_F*Ks[i]-V_over_F)
             except:
                 pass
             ys[i] = xs[i]*Ks[i]

requirements.txt

@@ -1,4 +1,4 @@
 numpy
 scipy
 pandas
-fluids>=1.0.9
+fluids>=1.0.11

requirements_docs.txt

@@ -3,7 +3,7 @@ scipy
 numpydoc
 pint
 nbsphinx
-fluids>=1.0.9
+fluids>=1.0.11
 IPython
 ipython
 numba

requirements_test.txt

@@ -2,7 +2,7 @@ numpy
 scipy
 pandas
 sympy
-fluids>=1.0.9
+fluids>=1.0.11
 pytest
 pytest-cov
 coveralls

setup.py

@@ -45,6 +45,7 @@
     'Programming Language :: Python :: 3.8',
     'Programming Language :: Python :: 3.9',
     'Programming Language :: Python :: 3.10',
+    'Programming Language :: Python :: 3.11',
     'Programming Language :: Python :: Implementation :: CPython',
     'Programming Language :: Python :: Implementation :: PyPy',
     'Topic :: Education',
@@ -61,7 +62,7 @@
   version = '1.0.12',
   description = 'Chemical properties component of Chemical Engineering Design Library (ChEDL)',
   author = 'Caleb Bell',
-  install_requires=['fluids>=1.0.9', 'scipy', 'numpy', 'pandas'],
+  install_requires=['fluids>=1.0.11', 'scipy', 'numpy', 'pandas'],
   extras_require = {
       'Coverage documentation':  ['wsgiref>=0.1.2', 'coverage>=4.0.3']
   },