From e5b2226f3837ce08620f15d1f824c727b28a3b78 Mon Sep 17 00:00:00 2001 From: cortespea Date: Sun, 26 Nov 2023 12:21:14 -0600 Subject: [PATCH] better integration and code reduction --- .../equilibrium/binary_phase_fraction.py | 2 +- .../equilibrium/fugacity_coefficients.py | 2 +- thermosteam/equilibrium/vle.py | 79 ++++++++----------- 3 files changed, 35 insertions(+), 48 deletions(-) diff --git a/thermosteam/equilibrium/binary_phase_fraction.py b/thermosteam/equilibrium/binary_phase_fraction.py index 4d971dbc..1d485c9a 100644 --- a/thermosteam/equilibrium/binary_phase_fraction.py +++ b/thermosteam/equilibrium/binary_phase_fraction.py @@ -54,7 +54,7 @@ def phase_composition(zs, Ks, phi): return zs * Ks / (phi * Ks + (1. - phi)) # @njit(cache=True) -def solve_phase_fraction_Rashford_Rice(zs, Ks, guess, za, zb): +def solve_phase_fraction_Rashford_Rice(zs, Ks, guess, za=0, zb=0): """ Return phase fraction for N-component binary equilibrium by numerically solving the Rashford Rice equation. diff --git a/thermosteam/equilibrium/fugacity_coefficients.py b/thermosteam/equilibrium/fugacity_coefficients.py index 7723f044..a914b451 100644 --- a/thermosteam/equilibrium/fugacity_coefficients.py +++ b/thermosteam/equilibrium/fugacity_coefficients.py @@ -91,7 +91,7 @@ def __new__(cls, chemicals): @classmethod def subclass(cls, EOS, name=None): - if name is None: name = EOS.__name__[:-3] + 'FugacityCoefficients' + if name is None: name = EOS.__name__.replace('MIX', '') + 'FugacityCoefficients' return type(name, (cls,), dict(EOS=EOS, cache={})) @property diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py index c9574866..00df23ef 100644 --- a/thermosteam/equilibrium/vle.py +++ b/thermosteam/equilibrium/vle.py @@ -25,49 +25,44 @@ __all__ = ('VLE', 'VLECache') -# @njit(cache=True) -def xV_iter(xV, pcf_Psat_over_P_phi, T, P, z, z_light, z_heavy, f_gamma, gamma_args): - xV = xV.copy() +@njit(cache=True) +def xy(xV, Ks): x = xV[:-1] - V = xV[-1] x[x < 0.] = 0. - x = fn.normalize(x) - Ks = pcf_Psat_over_P_phi * f_gamma(x, T, *gamma_args) - V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, V, z_light, z_heavy) + x /= x.sum() + return x, Ks * x + +@njit(cache=True) +def update_xV(xV, V, Ks, z): if V < 0.: V = 0. elif V > 1.: V = 1. xV[-1] = V xV[:-1] = z/(1. + V * (Ks - 1.)) + +def xV_iter(xV, pcf_Psat_over_P, T, P, + z, z_light, z_heavy, f_gamma, gamma_args, + Ks, f_phi): + xV = xV.copy() + x, y = xy(xV, Ks) + Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) + V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, xV[-1], z_light, z_heavy) + update_xV(xV, V, Ks, z) return xV -@njit(cache=True) -def xV_iter_2n(xV, pcf_Psat_over_P_phi, T, P, z, f_gamma, gamma_args): +def xV_iter_2n(xV, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, Ks, f_phi): xV = xV.copy() - x = xV[:-1] - V = xV[-1] - x[x < 0.] = 0. - x = fn.normalize(x) - Ks = pcf_Psat_over_P_phi * f_gamma(x, T, *gamma_args) + x, y = xy(xV, Ks) + Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) V = binary.compute_phase_fraction_2N(z, Ks) - if V < 0.: V = 0. - elif V > 1.: V = 1. - xV[-1] = V - xV[:-1] = z/(1. + V * (Ks - 1.)) + update_xV(xV, V, Ks, z) return xV -@njit(cache=True) -def xV_iter_3n(xV, pcf_Psat_over_P_phi, T, P, z, f_gamma, gamma_args): +def xV_iter_3n(xV, pcf_Psat_over_P, T, P, z, f_gamma, gamma_args, Ks, f_phi): xV = xV.copy() - x = xV[:-1] - V = xV[-1] - x[x < 0.] = 0. - x = fn.normalize(x) - Ks = pcf_Psat_over_P_phi * f_gamma(x, T, *gamma_args) + x, y = xy(xV, Ks) + Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P) V = binary.compute_phase_fraction_3N(z, Ks) - if V < 0.: V = 0. - elif V > 1.: V = 1. - xV[-1] = V - xV[:-1] = z/(1. + V * (Ks - 1.)) + update_xV(xV, V, Ks, z) return xV def set_flows(vapor_mol, liquid_mol, index, vapor_data, total_data): @@ -1196,23 +1191,23 @@ def _V_err_at_P(self, P, V): def _V_err_at_T(self, T, V): return self._solve_v(T, self._P).sum()/self._F_mol_vle - V - def _y_iter(self, y, pcf_Psat_over_P, T, P): - phi = self._phi(y, T, P) + def _solve_y(self, y, pcf_Psat_over_P, T, P): gamma = self._gamma x = self._x - pcf_Psat_over_P_phi = pcf_Psat_over_P / phi + pcf_Psat_over_P = pcf_Psat_over_P N = self._N z = self._z + Ks = pcf_Psat_over_P.copy() if N > 3 or self._z_light or self._z_heavy: f = xV_iter - args = (pcf_Psat_over_P_phi, T, P, z, self._z_light, - self._z_heavy, gamma.f, gamma.args) + args = (pcf_Psat_over_P, T, P, z, self._z_light, + self._z_heavy, gamma.f, gamma.args, Ks, self._phi) elif N == 2: f = xV_iter_2n - args = (pcf_Psat_over_P_phi, T, P, z, gamma.f, gamma.args) + args = (pcf_Psat_over_P, T, P, z, gamma.f, gamma.args, Ks, self._phi) elif N == 3: f = xV_iter_3n - args = (pcf_Psat_over_P_phi, T, P, z, gamma.f, gamma.args) + args = (pcf_Psat_over_P, T, P, z, gamma.f, gamma.args, Ks, self._phi) xV = np.zeros(x.size + 1) xV[:-1] = x xV[-1] = self._V @@ -1228,7 +1223,7 @@ def _y_iter(self, y, pcf_Psat_over_P, T, P): else: Ks = (z / x - 1) / V + 1. self._z_last = z - v = self._F_mol * V * x * Ks + v = self._F_mol * V * x * Ks return fn.normalize(v, v.sum() + self._F_mol_light) def _solve_v(self, T, P): @@ -1252,15 +1247,7 @@ def _solve_v(self, T, P): self._bubble_point.Psats]) pcf_Psats_over_P = self._pcf(T, P, Psats) * Psats / P self._T = T - if isinstance(self._phi, IdealFugacityCoefficients): - y = self._y_iter(self._y, pcf_Psats_over_P, T, P) - else: - y = flx.aitken(self._y_iter, self._y, self.y_tol, - args=(pcf_Psats_over_P, T, P), - checkiter=False, - checkconvergence=False, - convergenceiter=5, - maxiter=self.maxiter) + y = self._solve_y(self._y, pcf_Psats_over_P, T, P) self._v = v = self._F_mol * self._V * y mask = v > self._mol_vle v[mask] = self._mol_vle[mask]