Merge pull request #26 from pathsim/fix/audit-findings

Fix four audit findings: flash holdup, UNIQUAC NaN, Redlich-Kister, TCAP stub

Author: milanofthe

src/pathsim_chem/process/flash_drum.py

@@ -43,11 +43,20 @@ class FlashDrum(DynamicalSystem):
 
     Dynamic States
     ---------------
-    The holdup moles of each component in the liquid phase:
+    The drum holds well-mixed liquid; vapor exits immediately at feed-flash
+    composition. Liquid component holdup follows a CSTR balance with the
+    feed-flash liquid as inlet:
 
     .. math::
 
-        \\frac{dN_i}{dt} = F z_i - V y_i - L x_i
+        \\frac{dN_i}{dt} = L_{in} \\, (x_{eq,i} - x_{drum,i})
+
+    where :math:`L_{in} = (1-\\beta) F` is the liquid feed to the drum,
+    :math:`x_{eq}` the Rachford-Rice liquid composition for the feed, and
+    :math:`x_{drum} = N / \\sum_j N_j` the drum liquid composition. Total
+    holdup :math:`M = \\sum_i N_i` is preserved exactly. The output ``x_1``
+    is :math:`x_{drum,1}` (dynamic), while ``y_1``, ``V_rate``, ``L_rate``
+    follow the instantaneous feed flash.
 
     Parameters
     ----------
@@ -147,20 +156,24 @@ def _pad_u(u):
                 return padded
             return u
 
-        # rhs of flash drum ode
+        # rhs of flash drum ode: liquid CSTR fed by feed-flash liquid
         def _fn_d(x, u, t):
             u = _pad_u(u)
             F_in, z_1, T, P = u
             z = np.array([z_1, 1.0 - z_1])
 
-            beta, y, x_eq = _solve_vle(z, T, P)
+            if F_in <= 0.0:
+                return np.zeros(2)
 
-            V_rate = beta * F_in
-            L_rate = (1.0 - beta) * F_in
+            beta, _, x_eq = _solve_vle(z, T, P)
+            L_in = (1.0 - beta) * F_in
+
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
 
-            return F_in * z - V_rate * y - L_rate * x_eq
+            return L_in * (x_eq - x_drum)
 
-        # algebraic output
+        # algebraic output: V/L/y from feed flash (instantaneous), x from drum
         def _fn_a(x, u, t):
             u = _pad_u(u)
             F_in, z_1, T, P = u
@@ -171,7 +184,10 @@ def _fn_a(x, u, t):
             V_rate = beta * F_in
             L_rate = (1.0 - beta) * F_in
 
-            return np.array([V_rate, L_rate, y[0], x_eq[0]])
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
+
+            return np.array([V_rate, L_rate, y[0], x_drum[0]])
 
         super().__init__(
             func_dyn=_fn_d,

src/pathsim_chem/process/multicomponent_flash.py

@@ -35,11 +35,20 @@ class MultiComponentFlash(DynamicalSystem):
 
     Dynamic States
     ---------------
-    The holdup moles of each component in the liquid phase:
+    The drum holds well-mixed liquid; vapor exits immediately at feed-flash
+    composition. Liquid component holdup follows a CSTR balance with the
+    feed-flash liquid as inlet:
 
     .. math::
 
-        \\frac{dN_i}{dt} = F z_i - V y_i - L x_i
+        \\frac{dN_i}{dt} = L_{in} \\, (x_{eq,i} - x_{drum,i})
+
+    where :math:`L_{in} = (1-\\beta) F` is the liquid feed to the drum,
+    :math:`x_{eq}` the Rachford-Rice liquid composition for the feed, and
+    :math:`x_{drum} = N / \\sum_j N_j` the drum liquid composition. Total
+    holdup :math:`M = \\sum_i N_i` is preserved exactly. Outputs ``x_i``
+    are :math:`x_{drum,i}` (dynamic); ``V_rate``, ``L_rate``, ``y_i`` come
+    from the instantaneous feed flash.
 
     Parameters
     ----------
@@ -207,22 +216,26 @@ def _extract_z(u):
             z_last = 1.0 - np.sum(z_partial)
             return np.append(z_partial, z_last)
 
-        # rhs of flash drum ode
+        # rhs of flash drum ode: liquid CSTR fed by feed-flash liquid
         def _fn_d(x, u, t):
             u = _pad_u(u)
             F_in = u[0]
             z = _extract_z(u)
             T = u[nc]
             P = u[nc + 1]
 
-            beta, y, x_eq = _solve_vle(z, T, P)
+            if F_in <= 0.0:
+                return np.zeros(nc)
 
-            V_rate = beta * F_in
-            L_rate = (1.0 - beta) * F_in
+            beta, _, x_eq = _solve_vle(z, T, P)
+            L_in = (1.0 - beta) * F_in
+
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
 
-            return F_in * z - V_rate * y - L_rate * x_eq
+            return L_in * (x_eq - x_drum)
 
-        # algebraic output
+        # algebraic output: V/L/y from feed flash (instant), x from drum state
         def _fn_a(x, u, t):
             u = _pad_u(u)
             F_in = u[0]
@@ -235,12 +248,15 @@ def _fn_a(x, u, t):
             V_rate = beta * F_in
             L_rate = (1.0 - beta) * F_in
 
+            M = x.sum()
+            x_drum = x / M if M > 1e-30 else x_eq
+
             # output: V_rate, L_rate, y_1..y_{nc-1}, x_1..x_{nc-1}
             result = np.empty(2 + 2 * (nc - 1))
             result[0] = V_rate
             result[1] = L_rate
             result[2:2 + nc - 1] = y[:nc - 1]
-            result[2 + nc - 1:] = x_eq[:nc - 1]
+            result[2 + nc - 1:] = x_drum[:nc - 1]
 
             return result
 

src/pathsim_chem/thermodynamics/activity_coefficients.py

@@ -306,17 +306,22 @@ def _eval(self, T):
         tau = np.exp(self.a + self.b_param / T + self.c_param * np.log(T) + self.d_param * T)
 
         # volume and surface fractions
-        V = self.r * x / np.dot(self.r, x)
-        F = self.q * x / np.dot(self.q, x)
+        sum_rx = np.dot(self.r, x)
+        sum_qx = np.dot(self.q, x)
         Fp = self.q_prime * x / np.dot(self.q_prime, x)
 
         # combinatorial part
+        # V_i / x_i = r_i / sum_j(r_j x_j) is well-defined even at x_i=0
+        # F_i / V_i = (q_i / r_i) * sum_rx / sum_qx (algebraic identity)
+        V_over_x = self.r / sum_rx
+        F_over_V = (self.q / self.r) * (sum_rx / sum_qx)
+
         ln_gamma_C = np.zeros(n)
         sum_xl = np.dot(x, self.l)
         for i in range(n):
-            ln_gamma_C[i] = (np.log(V[i] / x[i])
-                             + self.z / 2 * self.q[i] * np.log(F[i] / V[i])
-                             + self.l[i] - V[i] / x[i] * sum_xl)
+            ln_gamma_C[i] = (np.log(V_over_x[i])
+                             + self.z / 2 * self.q[i] * np.log(F_over_V[i])
+                             + self.l[i] - V_over_x[i] * sum_xl)
 
         # residual part
         ln_gamma_R = np.zeros(n)

src/pathsim_chem/thermodynamics/enthalpy.py

@@ -301,22 +301,24 @@ class ExcessEnthalpyRedlichKister(Function):
     Computes the molar excess enthalpy using a Redlich-Kister polynomial
     expansion. This is a flexible, empirical model for representing binary
     excess properties. For a multi-component mixture, the excess enthalpy
-    is computed as a sum of binary pair contributions:
+    is summed over the unique binary pairs supplied:
 
     .. math::
 
-        h^E = \frac{1}{2} \sum_i \sum_j h^E_{ij}
+        h^E = \sum_{(i,j)} h^E_{ij}
 
-    Each binary pair contribution:
+    Each binary pair contribution follows the standard Redlich-Kister form
 
     .. math::
 
         h^E_{ij} = \frac{x_i x_j}{x_i + x_j}
-        \left(A(T) x_d + B(T) x_d^2 + C(T) x_d^3 + \ldots\right)
+        \left(A_0(T) + A_1(T) x_d + A_2(T) x_d^2 + \ldots\right)
 
-    where :math:`x_d = x_i - x_j` and the coefficients :math:`A(T)`,
-    :math:`B(T)`, ... are temperature-dependent polynomials:
-    :math:`A(T) = a_0 + a_1 T + a_2 T^2 + \ldots`
+    where :math:`x_d = x_i - x_j` and each :math:`A_k(T)` is a temperature
+    polynomial :math:`A_k(T) = a_{k,0} + a_{k,1} T + a_{k,2} T^2 + \ldots`.
+    For a true binary system :math:`x_i + x_j = 1` and the prefactor
+    reduces to :math:`x_i x_j`; in multicomponent mixtures the
+    :math:`/(x_i+x_j)` term provides the symmetric extension.
 
     **Input port:** ``T`` -- temperature [K].
 
@@ -328,12 +330,15 @@ class ExcessEnthalpyRedlichKister(Function):
         Mole fractions [N].
     coeffs : dict
         Redlich-Kister coefficients keyed by binary pair ``(i, j)`` as tuples.
-        Each value is a list of polynomial coefficient arrays, one per
-        Redlich-Kister term. Example for a single pair (0,1) with 3 terms::
+        Provide each unordered pair only once (e.g. ``(0, 1)``, not also
+        ``(1, 0)``). Each value is a list of polynomial coefficient arrays,
+        one per Redlich-Kister term. Example for a single pair (0,1) with
+        3 terms::
 
             {(0, 1): [[a0, a1], [b0, b1], [c0]]}
 
-        This gives A(T)=a0+a1*T, B(T)=b0+b1*T, C(T)=c0.
+        This gives A_0(T)=a0+a1*T, A_1(T)=b0+b1*T, A_2(T)=c0, where
+        :math:`A_k` multiplies :math:`x_d^k`.
     """
 
     input_port_labels = {"T": 0}
@@ -358,9 +363,9 @@ def _eval(self, T):
 
             xd = xi - xj
 
-            # evaluate each Redlich-Kister term
+            # evaluate Redlich-Kister polynomial: A_0 + A_1*xd + A_2*xd^2 + ...
             pair_hE = 0.0
-            xd_power = xd
+            xd_power = 1.0
             for poly_coeffs in terms:
                 # temperature polynomial: c0 + c1*T + c2*T^2 + ...
                 coeff_val = 0.0
@@ -373,7 +378,7 @@ def _eval(self, T):
 
             hE += xi * xj / x_sum * pair_hE
 
-        return 0.5 * hE
+        return hE
 
 
 # ENTHALPY DEPARTURE FROM EOS (9.7) =====================================================

src/pathsim_chem/tritium/tcap.py

@@ -12,15 +12,16 @@
 # BLOCKS ================================================================================
 
 class TCAP1D(ODE):
-    """This block models the Thermal Cycle Absorption Process (TCAP) in 1d. 
+    """This block models the Thermal Cycle Absorption Process (TCAP) in 1d.
 
-    The model uses a 1d finite difference spatial discretization to construct 
+    The model uses a 1d finite difference spatial discretization to construct
     a nonlinear ODE internally as proposed in
 
         https://doi.org/10.1016/j.ijhydene.2023.03.101
 
-
     """
-    raise NotImplementedError("TCAP1D block is currently not impolemented!")
+
+    def __init__(self, *args, **kwargs):
+        raise NotImplementedError("TCAP1D block is not yet implemented")
 
 

tests/process/test_flash_drum.py

@@ -136,6 +136,62 @@ def test_no_feed(self):
         self.assertAlmostEqual(F.outputs[0], 0.0)  # V_rate
         self.assertAlmostEqual(F.outputs[1], 0.0)  # L_rate
 
+    def test_holdup_dynamics_drives_to_equilibrium(self):
+        """At steady state the drum liquid composition must equal the RR
+        equilibrium liquid composition for the feed."""
+        F = FlashDrum(holdup=100.0, N0=[80.0, 20.0])  # off-equilibrium init
+        F.set_solver(EUF, parent=None)
+
+        # T=370 K gives two-phase region for benzene/toluene defaults at 1 atm
+        T, P = 370.0, 101325.0
+        u = np.array([10.0, 0.5, T, P])
+
+        # at the RR equilibrium x_eq, dN/dt must vanish
+        # compute x_eq via direct VLE (binary RR with same Antoine defaults)
+        Psat = np.exp(F.antoine_A - F.antoine_B / (T + F.antoine_C))
+        K = Psat / P
+        z = np.array([0.5, 0.5])
+        d1, d2 = K[0] - 1, K[1] - 1
+        beta = -(z[0]*d1 + z[1]*d2) / (d1*d2)
+        x_eq = z / (1.0 + beta * (K - 1.0))
+        x_eq = x_eq / x_eq.sum()
+
+        # state at equilibrium with same total holdup
+        N_eq = 100.0 * x_eq
+        dN = F.op_dyn(N_eq, u, 0.0)
+        self.assertTrue(np.allclose(dN, 0.0, atol=1e-10))
+
+        # state away from equilibrium: dN must be non-zero
+        dN_off = F.op_dyn(np.array([80.0, 20.0]), u, 0.0)
+        self.assertGreater(np.linalg.norm(dN_off), 1e-3)
+
+    def test_holdup_total_moles_conserved(self):
+        """dM/dt = sum(dN/dt) must be exactly zero (perfect level control)."""
+        F = FlashDrum(holdup=100.0, N0=[70.0, 30.0])
+        F.set_solver(EUF, parent=None)
+
+        u = np.array([5.0, 0.4, 355.0, 101325.0])
+        for state in (np.array([70.0, 30.0]),
+                      np.array([20.0, 80.0]),
+                      np.array([1.0, 99.0])):
+            dN = F.op_dyn(state, u, 0.0)
+            self.assertAlmostEqual(dN.sum(), 0.0, places=10,
+                                    msg=f"dM/dt != 0 for state {state}")
+
+    def test_x_output_uses_drum_state(self):
+        """x_1 output must reflect drum state, not feed composition."""
+        F = FlashDrum(holdup=100.0, N0=[90.0, 10.0])
+        F.set_solver(EUF, parent=None)
+
+        F.inputs[0] = 10.0
+        F.inputs[1] = 0.3       # different from drum
+        F.inputs[2] = 360.0
+        F.inputs[3] = 101325.0
+
+        F.update(None)
+        # drum state x_1 = 90/100 = 0.9
+        self.assertAlmostEqual(F.outputs[3], 0.9, places=8)
+
 
 # RUN TESTS LOCALLY ====================================================================
 

tests/thermodynamics/test_activity_coefficients.py

@@ -174,6 +174,9 @@ def test_pure_component_gamma_is_one(self):
         )
         gamma1, gamma2 = eval_block(U, 350)
         self.assertAlmostEqual(gamma1, 1.0, places=10)
+        # gamma_2 at infinite dilution must be finite (not NaN)
+        self.assertTrue(np.isfinite(gamma2))
+        self.assertGreater(gamma2, 0.0)
 
     def test_symmetric_case(self):
         # Identical r, q, symmetric a => gamma_1 = gamma_2

tests/thermodynamics/test_enthalpy.py

@@ -201,28 +201,40 @@ def test_init(self):
         )
         self.assertEqual(hE.n, 2)
 
-    def test_symmetric_binary(self):
-        # Simplest case: one-term Redlich-Kister with constant A
-        # h^E_ij = x_i*x_j/(x_i+x_j) * A * (x_i-x_j)
-        # For x=[0.5, 0.5]: h^E = 0.5 * 0.5 / 1.0 * A * 0 = 0
+    def test_constant_coeff_only(self):
+        # One-term Redlich-Kister with constant A_0:
+        # h^E = x_i*x_j/(x_i+x_j) * A_0
+        # For x=[0.5, 0.5], A_0=1000: h^E = 0.25/1.0 * 1000 = 250
         hE = ExcessEnthalpyRedlichKister(
             x=[0.5, 0.5],
             coeffs={(0, 1): [[1000.0]]},
         )
         result = eval_block_T(hE, 300)
-        self.assertAlmostEqual(result, 0.0, places=10)
+        self.assertAlmostEqual(result, 250.0, places=10)
 
-    def test_asymmetric_composition(self):
-        # x=[0.3, 0.7], one-term A=2000
-        # h^E = 0.5 * 0.3*0.7/1.0 * 2000 * (0.3-0.7) = 0.5 * 0.21 * 2000 * (-0.4) = -84
+    def test_two_term_polynomial(self):
+        # x=[0.3, 0.7], coeffs A_0=2000, A_1=500
+        # h^E = 0.21 * (2000 + 500*(-0.4)) = 0.21 * 1800 = 378
         hE = ExcessEnthalpyRedlichKister(
             x=[0.3, 0.7],
-            coeffs={(0, 1): [[2000.0]]},
+            coeffs={(0, 1): [[2000.0], [500.0]]},
         )
         result = eval_block_T(hE, 300)
-        expected = 0.5 * 0.3 * 0.7 / 1.0 * 2000.0 * (0.3 - 0.7)
+        expected = 0.3 * 0.7 / 1.0 * (2000.0 + 500.0 * (0.3 - 0.7))
         self.assertAlmostEqual(result, expected, places=5)
 
+    def test_antisymmetric_under_pair_swap(self):
+        # Odd-order RK terms should flip sign when (i,j) <-> (j,i),
+        # since x_d = x_i - x_j flips. Verify by swapping x values.
+        hE_a = ExcessEnthalpyRedlichKister(
+            x=[0.3, 0.7], coeffs={(0, 1): [[0.0], [1000.0]]},  # A_1 only
+        )
+        hE_b = ExcessEnthalpyRedlichKister(
+            x=[0.7, 0.3], coeffs={(0, 1): [[0.0], [1000.0]]},
+        )
+        self.assertAlmostEqual(eval_block_T(hE_a, 300),
+                               -eval_block_T(hE_b, 300), places=8)
+
     def test_temperature_dependent_coeffs(self):
         # A(T) = 1000 + 2*T
         hE = ExcessEnthalpyRedlichKister(