Even the unbounded solver Rachford_Rice_solution_LN2 apparently still…

… neeeds limits

chemicals/rachford_rice.py

@@ -1157,7 +1157,6 @@ 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
@@ -1289,12 +1288,18 @@ def Rachford_Rice_solution_LN2(zs, Ks, guess=None):
         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
+    near_high = V_over_F_max*one_epsilon_smaller
+    solver_high = -log((V_over_F_max-near_high)/(near_high-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, 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, xtol=1.48e-12, bisection=True, 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, 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

tests/test_rachford_rice.py

@@ -186,6 +186,19 @@ def test_Rachford_Rice_solution_LN2_points():
     assert_close1d(ys[:-1], ys_expect[:-1], rtol=1e-15)
     assert_close(V_over_F, -87.47697657610895, rtol=1e-15)
 
+    # The following case with the first component being lightest has a small composition error
+    # The issue with the solver was that the newton solver was jumping so high the math thought it was at VF max, and zero divided
+    zs = [0.012692054090137975, 0.010666797278186979, 0.0009194857929749981, 0.017470790744137966, 0.016843066321609965, 0.0002561964223459995, 0.05279574102694189, 0.027585130658382945, 0.0030041555848069937, 0.01006221306055498, 0.008441411298762982, 0.02350560912243095, 0.04162860619876592, 0.05564310727820389, 0.0497144837663429, 0.05409552150533089, 0.010985515989031977, 0.040266824975770915, 0.0470163214248899, 0.03554837372222793, 0.045399591363595906, 0.023349127390622955, 0.038940262794703924, 0.012157361905424977, 0.00507766294761099, 0.05467583955106789, 0.014954684178677971, 0.0534356116296629, 0.051497730943570894, 0.00519826703202799, 0.0535400616269459, 0.014088047967291972, 0.02114528671286696, 0.022257393945042954, 0.02665248264602195, 0.007507531190584985, 0.030981649912439936] 
+    Ks = [0.723182996154432, 1.23519391187981, 3.64790350249987, 18.5637643987892, 3.82730273431659, 9.55998243548818, 7.28123344668702, 8.3768614499351, 5.43212719859162, 3.27302690798611, 2.49190317901262, 9.76956309032054, 5.78718965354511, 7.36863057377774, 2.0532435682409, 0.976329915082406, 24.3687952929195, 1.81112520953641, 1.52874821640877, 8.58329677042206, 0.826908749936913, 4.87507235159708, 5.40262156026164, 5.97055327322089, 2.94555106080988, 15.3811087714507, 15.1963390052198, 12.605673455515, 5.252476057943, 1.37013848462909, 7.84987161014184, 15.7522352359933, 3.82689447066015, 3.51649883363176, 6.78102108048812, 3.40186445663176, 40.2867569231891]
+    V_over_F, xs, ys = Rachford_Rice_solution_LN2(zs=zs, Ks=Ks)
+    xs_expect = [0.7434579686842309, 0.005812550126126229, 8.839307186387247e-05, 0.0002757131246052513, 0.001525743263942414, 8.160424371206078e-06, 0.0022655662585562113, 0.0010143854482568462, 0.0001794843527914173, 0.0011092583644856495, 0.0013404416645140496, 0.0007313694772945749, 0.002312897454584283, 0.002356372542758535, 0.010488548558731955, 0.05905936009245475, 0.00013081348124969965, 0.01037761294651069, 0.016339342461795813, 0.001272905735045562, 0.11780367052125561, 0.0015819515293255151, 0.002341154375732744, 0.0006518847251555405, 0.000642066867929061, 0.0010501501995641331, 0.0002908979504556291, 0.001265955755191919, 0.003198658084369941, 0.0022461543311921585, 0.0021143121274689357, 0.0002639071891884155, 0.0019157150712173644, 0.0022401562677665074, 0.0012380733331850689, 0.0007878947124000653, 0.00022050942443154514]
+    ys_expect = [0.5376561613079499, 0.007179626528287341, 0.0003224493964489431, 0.005118273486825894, 0.0058394813659519195, 7.801351365485979e-05, 0.016496116817485057, 0.008497366356877912, 0.0009749818345198717, 0.0036306324748701948, 0.0033402508450835282, 0.007145160250764105, 0.013385176218880983, 0.017363238761780933, 0.02153554486839875, 0.057661420023887584, 0.003187766946128094, 0.018795156422236937, 0.02497874064576243, 0.01092572768466829, 0.09741288592871145, 0.007712128162181535, 0.012648371106234603, 0.0038921124795401126, 0.0018912407439393229, 0.01615247444585659, 0.004420583871047374, 0.01595822485907922, 0.016800875005698938, 0.003077542491582691, 0.016597078744396997, 0.004157128124565709, 0.007331239413402048, 0.0078775069027538, 0.008395401371518145, 0.0026803110176818854, 0.008883609581345995]
+    assert_close1d(xs, xs_expect, rtol=1e-13)
+    assert_close1d(ys, ys_expect, rtol=1e-13)
+    assert_close1d(xs[1:], xs_expect[1:], rtol=1e-15)
+    assert_close1d(ys[1:], ys_expect[1:], rtol=1e-15)
+    assert_close(V_over_F, 3.5508236001931017, rtol=1e-15)
+
 
 def test_Rachford_Rice_solution_LN2_near_VF1():