diff --git a/thermo/eos.py b/thermo/eos.py index 91a71b3e..29bc83fe 100644 --- a/thermo/eos.py +++ b/thermo/eos.py @@ -867,7 +867,7 @@ class GCEOS: N = 1 """The number of components in the EOS""" - scalar = True + vectorized = False multicomponent = False """Whether or not the EOS is multicomponent or not""" @@ -1028,7 +1028,7 @@ def as_json(self): ''' # vaguely jsonpickle compatible d = object_data(self) - if not self.scalar: + if self.vectorized: d = serialize.arrays_to_lists(d) # TODO: delete kwargs and reconstruct it # Need to add all kwargs attributes diff --git a/thermo/eos_alpha_functions.py b/thermo/eos_alpha_functions.py index 9c6d04f2..f7077c07 100644 --- a/thermo/eos_alpha_functions.py +++ b/thermo/eos_alpha_functions.py @@ -2007,9 +2007,9 @@ def a_alphas_vectorized(self, T): Tr = T/Tcs[i] a_alpha = ais[i]*(Tr**(coeffs[2]*(coeffs[1] - 1.0))*exp(coeffs[0]*(1.0 - (Tr)**(coeffs[1]*coeffs[2])))) a_alphas.append(a_alpha) - if self.scalar: - return a_alphas - return array(a_alphas) + if self.vectorized: + return array(a_alphas) + return a_alphas def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -2064,9 +2064,9 @@ def a_alpha_and_derivatives_vectorized(self, T): da_alpha_dTs[i] = x8*(x1 - x7)*T_inv d2a_alpha_dT2s[i] = d2a_alpha_dT2 - if self.scalar: - return a_alphas, da_alpha_dTs, d2a_alpha_dT2s - return array(a_alphas), array(da_alpha_dTs), array(d2a_alpha_dT2s) + if self.vectorized: + return array(a_alphas), array(da_alpha_dTs), array(d2a_alpha_dT2s) + return a_alphas, da_alpha_dTs, d2a_alpha_dT2s class Soave_1993_a_alpha(a_alpha_base): @@ -2350,9 +2350,9 @@ def a_alphas_vectorized(self, T): Tcs, omegas, ais = self.Tcs, self.omegas, self.ais a_alphas = [TWU_a_alpha_common(T, Tcs[i], omegas[i], ais[i], full=False, method='SRK') for i in range(self.N)] - if self.scalar: - return a_alphas - return array(a_alphas) + if self.vectorized: + return array(a_alphas) + return a_alphas def a_alpha_and_derivatives_vectorized(self, T): Tcs, omegas, ais = self.Tcs, self.omegas, self.ais @@ -2362,9 +2362,9 @@ def a_alpha_and_derivatives_vectorized(self, T): r0.append(v0) r1.append(v1) r2.append(v2) - if self.scalar: - return r0, r1, r2 - return array(r0), array(r1), array(r2) + if self.vectorized: + return array(r0), array(r1), array(r2) + return r0, r1, r2 @@ -2473,9 +2473,9 @@ def a_alphas_vectorized(self, T): Tcs, omegas, ais = self.Tcs, self.omegas, self.ais a_alphas = [TWU_a_alpha_common(T, Tcs[i], omegas[i], ais[i], full=False, method='PR') for i in range(self.N)] - if self.scalar: - return a_alphas - return array(a_alphas) + if self.vectorized: + return array(a_alphas) + return a_alphas def a_alpha_and_derivatives_vectorized(self, T): Tcs, omegas, ais = self.Tcs, self.omegas, self.ais @@ -2485,9 +2485,9 @@ def a_alpha_and_derivatives_vectorized(self, T): r0.append(v0) r1.append(v1) r2.append(v2) - if self.scalar: - return r0, r1, r2 - return array(r0), array(r1), array(r2) + if self.vectorized: + return array(r0), array(r1), array(r2) + return r0, r1, r2 class Soave_1979_a_alpha(a_alpha_base): diff --git a/thermo/eos_mix.py b/thermo/eos_mix.py index 9361abfa..a6903218 100644 --- a/thermo/eos_mix.py +++ b/thermo/eos_mix.py @@ -328,9 +328,9 @@ class GCEOSMIX(GCEOS): multicomponent = True """All inherited classes of GCEOSMIX are multicomponent. """ - scalar = True - """Whether the model is implemented using pure-Python lists of floats, - or numpy arrays of float64. + vectorized = False + """Whether the model is implemented using numpy arrays of float64, + or pure-Python lists of floats. """ translated = False """Whether or not the model implements volume translation. @@ -457,7 +457,7 @@ def from_json(cls, json_repr): eos_name = d['py/object'] del d['py/object'] del d['json_version'] - if not d['scalar']: + if d['vectorized']: d = serialize.naive_lists_to_arrays(d) try: @@ -532,9 +532,9 @@ def to_TP_zs_fast(self, T, P, zs, only_l=False, only_g=False, full_alphas=True): RKMIX(Tcs=[126.1, 190.6], Pcs=[3394000.0, 4604000.0], omegas=[0.04, 0.011], kijs=[[0.0, 0.0], [0.0, 0.0]], zs=[0.6, 0.4], T=300, P=100000.0) ''' new = self.__class__.__new__(self.__class__) # potentially also object.__new__(self.__class__) - new.N, new.Tcs, new.Pcs, new.omegas, new.kijs, new.one_minus_kijs, new.kwargs, new.ais, new.bs, new.scalar = ( + new.N, new.Tcs, new.Pcs, new.omegas, new.kijs, new.one_minus_kijs, new.kwargs, new.ais, new.bs, new.vectorized = ( self.N, self.Tcs, self.Pcs, self.omegas, self.kijs, self.one_minus_kijs, self.kwargs, self.ais, self.bs, - self.scalar + self.vectorized ) # new.N = self.N # new.Tcs = self.Tcs @@ -545,7 +545,7 @@ def to_TP_zs_fast(self, T, P, zs, only_l=False, only_g=False, full_alphas=True): # new.kwargs = self.kwargs # new.ais = self.ais # new.bs = self.bs - # new.scalar = self.scalar + # new.vectorized = self.vectorized if T == self.T: new.a_alphas = self.a_alphas @@ -1061,10 +1061,10 @@ def a_alpha_and_derivatives(self, T, full=True, quick=True, else: self.a_alphas = a_alphas = self.a_alphas_vectorized(T) da_alpha_dTs = d2a_alpha_dT2s = None - if self.scalar: - self.a_alpha_roots = [sqrt(i) for i in a_alphas] - else: + if self.vectorized: self.a_alpha_roots = npsqrt(a_alphas) + else: + self.a_alpha_roots = [sqrt(i) for i in a_alphas] else: try: a_alphas, da_alpha_dTs, d2a_alpha_dT2s, = self.a_alphas, self.da_alpha_dTs, self.d2a_alpha_dT2s @@ -1076,12 +1076,12 @@ def a_alpha_and_derivatives(self, T, full=True, quick=True, self.a_alphas = a_alphas = self.a_alphas_vectorized(T) da_alpha_dTs = d2a_alpha_dT2s = None - zs, one_minus_kijs, scalar, N, a_alpha_roots = self.zs, self.one_minus_kijs, self.scalar, self.N, self.a_alpha_roots + zs, one_minus_kijs, vectorized, N, a_alpha_roots = self.zs, self.one_minus_kijs, self.vectorized, self.N, self.a_alpha_roots if full: - if scalar: - a_alpha_j_rows, da_alpha_dT_j_rows = [0.0]*N, [0.0]*N - else: + if vectorized: a_alpha_j_rows, da_alpha_dT_j_rows = zeros(N), zeros(N) + else: + a_alpha_j_rows, da_alpha_dT_j_rows = [0.0]*N, [0.0]*N a_alpha, da_alpha_dT, d2a_alpha_dT2, self.a_alpha_j_rows, self.da_alpha_dT_j_rows = ( a_alpha_and_derivatives_quadratic_terms(a_alphas, a_alpha_roots, da_alpha_dTs, d2a_alpha_dT2s, T, zs, one_minus_kijs, @@ -1089,7 +1089,7 @@ def a_alpha_and_derivatives(self, T, full=True, quick=True, return a_alpha, da_alpha_dT, d2a_alpha_dT2 else: - a_alpha_j_rows = [0.0]*N if scalar else zeros(N) + a_alpha_j_rows = zeros(N) if vectorized else [0.0]*N a_alpha, self.a_alpha_j_rows = a_alpha_quadratic_terms(a_alphas, a_alpha_roots, T, zs, one_minus_kijs, a_alpha_j_rows) return a_alpha @@ -1111,23 +1111,23 @@ def a_alpha_and_derivatives(self, T, full=True, quick=True, # assert same_T # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = self.a_alpha_ijs, self.a_alpha_roots, self.a_alpha_ij_roots_inv # except (AttributeError, AssertionError): - # if self.scalar: - # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = [[0.0]*N for _ in range(N)], [0.0]*N, [[0.0]*N for _ in range(N)] - # else: + # if self.vectorized: # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = zeros((N, N)), zeros(N), zeros((N, N)) + # else: + # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = [[0.0]*N for _ in range(N)], [0.0]*N, [[0.0]*N for _ in range(N)] # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = a_alpha_aijs_composition_independent(a_alphas, one_minus_kijs, a_alpha_ijs=a_alpha_ijs, a_alpha_roots=a_alpha_roots, a_alpha_ij_roots_inv=a_alpha_ij_roots_inv) # self.a_alpha_ijs, self.a_alpha_roots, self.a_alpha_ij_roots_inv = a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv # if full: # try: # a_alpha, da_alpha_dT, d2a_alpha_dT2, a_alpha_ijs, da_alpha_dT_ijs, d2a_alpha_dT2_ijs = a_alpha_and_derivatives_full(a_alphas, da_alpha_dTs, d2a_alpha_dT2s, T, zs, one_minus_kijs, # a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv) - # if not self.scalar: + # if self.vectorized: # d2a_alpha_dT2_ijs = array(d2a_alpha_dT2_ijs) # except: # if self.N == 1: # a_alpha, da_alpha_dT, d2a_alpha_dT2 = a_alphas[0], da_alpha_dTs[0], d2a_alpha_dT2s[0] # d2a_alpha_dT2_ijs, da_alpha_dT_ijs, a_alpha_ijs = [[d2a_alpha_dT2s[0]]], [[da_alpha_dTs[0]]], [[a_alphas[0]]] - # if not self.scalar: + # if self.vectorized: # d2a_alpha_dT2_ijs, da_alpha_dT_ijs, a_alpha_ijs = array(d2a_alpha_dT2_ijs), array(da_alpha_dT_ijs), array(a_alpha_ijs) # self.d2a_alpha_dT2_ijs, self.da_alpha_dT_ijs, self.a_alpha_ijs = d2a_alpha_dT2_ijs, da_alpha_dT_ijs, a_alpha_ijs @@ -1143,12 +1143,12 @@ def a_alpha_and_derivatives_py(self, a_alphas, da_alpha_dTs, d2a_alpha_dT2s, T, # 4 ms pypy for 44*4, 1.3 ms for pythran, 10 ms python with numpy # 2 components 1.89 pypy, pythran 1.75 us, regular python 12.7 us. # 10 components - regular python 148 us, 9.81 us PyPy, 8.37 pythran in PyPy (flags have no effect; 14.3 us in regular python) - zs, one_minus_kijs, scalar, N, a_alpha_roots = self.zs, self.one_minus_kijs, self.scalar, self.N, self.a_alpha_roots + zs, one_minus_kijs, vectorized, N, a_alpha_roots = self.zs, self.one_minus_kijs, self.vectorized, self.N, self.a_alpha_roots if full: - if scalar: - a_alpha_j_rows, da_alpha_dT_j_rows = [0.0]*N, [0.0]*N - else: + if vectorized: a_alpha_j_rows, da_alpha_dT_j_rows = zeros(N), zeros(N) + else: + a_alpha_j_rows, da_alpha_dT_j_rows = [0.0]*N, [0.0]*N a_alpha, da_alpha_dT, d2a_alpha_dT2, self.a_alpha_j_rows, self.da_alpha_dT_j_rows = ( a_alpha_and_derivatives_quadratic_terms(a_alphas, a_alpha_roots, da_alpha_dTs, d2a_alpha_dT2s, T, zs, one_minus_kijs, @@ -1156,7 +1156,7 @@ def a_alpha_and_derivatives_py(self, a_alphas, da_alpha_dTs, d2a_alpha_dT2s, T, return a_alpha, da_alpha_dT, d2a_alpha_dT2 else: - a_alpha_j_rows = [0.0]*N if scalar else zeros(N) + a_alpha_j_rows = zeros(N) if vectorized else [0.0]*N a_alpha, self.a_alpha_j_rows = a_alpha_quadratic_terms(a_alphas, a_alpha_roots, T, zs, one_minus_kijs, a_alpha_j_rows=a_alpha_j_rows) return a_alpha @@ -1174,10 +1174,10 @@ def a_alpha_and_derivatives_numpy(self, a_alphas, da_alpha_dTs, d2a_alpha_dT2s, z_products = np.einsum('i,j', zs, zs) a_alpha = np.einsum('ij,ji', a_alpha_ijs, z_products) - if self.scalar: - self.a_alpha_ijs = a_alpha_ijs.tolist() - else: + if self.vectorized: self.a_alpha_ijs = a_alpha_ijs + else: + self.a_alpha_ijs = a_alpha_ijs.tolist() if full: term0 = np.einsum('j,i', a_alphas, da_alpha_dTs) @@ -1486,30 +1486,30 @@ def fugacities(self, only_l=False, only_g=False): .. [2] Walas, Stanley M. Phase Equilibria in Chemical Engineering. Butterworth-Heinemann, 1985. ''' - P, zs, scalar = self.P, self.zs, self.scalar + P, zs, vectorized = self.P, self.zs, self.vectorized if not only_g and hasattr(self, 'V_l'): self.lnphis_l = lnphis_l = self.fugacity_coefficients(self.Z_l) - if scalar: + if vectorized: + self.phis_l = phis_l = npexp(lnphis_l) + self.fugacities_l = zs*P*phis_l + else: try: self.phis_l = [exp(i) for i in lnphis_l] except: self.phis_l = [trunc_exp(i, trunc=1e308) for i in lnphis_l] self.fugacities_l = [phi*x*P for phi, x in zip(self.phis_l, zs)] - else: - self.phis_l = phis_l = npexp(lnphis_l) - self.fugacities_l = zs*P*phis_l if not only_l and hasattr(self, 'V_g'): self.lnphis_g = lnphis_g = self.fugacity_coefficients(self.Z_g) - if scalar: + if vectorized: + self.phis_g = phis_g = npexp(lnphis_g) + self.fugacities_g = zs*P*phis_g + else: try: self.phis_g = phis_g = [exp(i) for i in lnphis_g] except: self.phis_g = phis_g = [trunc_exp(i, trunc=1e308) for i in lnphis_g] self.fugacities_g = [phi*y*P for phi, y in zip(phis_g, zs)] - else: - self.phis_g = phis_g = npexp(lnphis_g) - self.fugacities_g = zs*P*phis_g @@ -2470,10 +2470,10 @@ def _a_alpha_j_rows(self): pass zs, N = self.zs, self.N a_alpha_ijs = self.a_alpha_ijs - if self.scalar: - a_alpha_j_rows = [0.0]*N - else: + if self.vectorized: a_alpha_j_rows = zeros(N) + else: + a_alpha_j_rows = [0.0]*N for i in range(N): l = a_alpha_ijs[i] @@ -2488,10 +2488,10 @@ def _set_alpha_matrices(self): N = self.N if not hasattr(self, 'd2a_alpha_dT2s'): self.a_alpha_and_derivatives(self.T, full=True, pure_a_alphas=True) - if self.scalar: - a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = [[0.0]*N for _ in range(N)], [0.0]*N, [[0.0]*N for _ in range(N)] - else: + if self.vectorized: a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = zeros((N, N)), zeros(N), zeros((N, N)) + else: + a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = [[0.0]*N for _ in range(N)], [0.0]*N, [[0.0]*N for _ in range(N)] a_alpha_ijs, a_alpha_roots, a_alpha_ij_roots_inv = a_alpha_aijs_composition_independent(self.a_alphas, self.one_minus_kijs, a_alpha_ijs=a_alpha_ijs,a_alpha_roots=a_alpha_roots, a_alpha_ij_roots_inv=a_alpha_ij_roots_inv) @@ -2499,7 +2499,7 @@ def _set_alpha_matrices(self): _, _, _, a_alpha_ijs, da_alpha_dT_ijs, d2a_alpha_dT2_ijs = a_alpha_and_derivatives_full( self.a_alphas, self.da_alpha_dTs, self.d2a_alpha_dT2s, self.T, self.zs, self.one_minus_kijs, a_alpha_ijs, self.a_alpha_roots, a_alpha_ij_roots_inv) - if not self.scalar: + if self.vectorized: d2a_alpha_dT2_ijs = array(d2a_alpha_dT2_ijs) da_alpha_dT_ijs = array(da_alpha_dT_ijs) self._d2a_alpha_dT2_ijs = d2a_alpha_dT2_ijs @@ -2616,7 +2616,7 @@ def _da_alpha_dT_j_rows(self): return self.da_alpha_dT_j_rows except: pass - zs, N, scalar = self.zs, self.N, self.scalar + zs, N, vectorized = self.zs, self.N, self.vectorized da_alpha_dT_ijs = self.da_alpha_dT_ijs # Handle the case of attempting to avoid a full alpha derivative matrix evaluation @@ -2624,10 +2624,10 @@ def _da_alpha_dT_j_rows(self): self.resolve_full_alphas() da_alpha_dT_ijs = self.da_alpha_dT_ijs - if scalar: - da_alpha_dT_j_rows = [0.0]*N - else: + if vectorized: da_alpha_dT_j_rows = zeros(N) + else: + da_alpha_dT_j_rows = [0.0]*N for i in range(N): l = da_alpha_dT_ijs[i] @@ -2645,7 +2645,7 @@ def _d2a_alpha_dT2_j_rows(self): return self.d2a_alpha_dT2_j_rows except AttributeError: pass - d2a_alpha_dT2_ijs, N, scalar = self.d2a_alpha_dT2_ijs, self.N, self.scalar + d2a_alpha_dT2_ijs, N, vectorized = self.d2a_alpha_dT2_ijs, self.N, self.vectorized # Handle the case of attempting to avoid a full alpha derivative matrix evaluation if d2a_alpha_dT2_ijs is None: @@ -2653,10 +2653,10 @@ def _d2a_alpha_dT2_j_rows(self): d2a_alpha_dT2_ijs = self.d2a_alpha_dT2_ijs zs = self.zs - if scalar: - d2a_alpha_dT2_j_rows = [0.0]*N - else: + if vectorized: d2a_alpha_dT2_j_rows = zeros(N) + else: + d2a_alpha_dT2_j_rows = [0.0]*N for i in range(N): l = d2a_alpha_dT2_ijs[i] for j in range(i): @@ -2708,10 +2708,10 @@ def db_dns(self): This derivative is checked numerically. ''' b = self.b - if self.scalar: - return [bi - b for bi in self.bs] - else: + if self.vectorized: return self.bs - b + else: + return [bi - b for bi in self.bs] @property def dnb_dns(self): @@ -2754,9 +2754,7 @@ def d2b_dzizjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[0.0]*N for i in range(N)] - return zeros((N, N)) + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d2b_dninjs(self): @@ -2779,16 +2777,16 @@ def d2b_dninjs(self): ''' bb = 2.0*self.b bs = self.bs - if self.scalar: - d2b_dninjs = [] - for bi in bs: - d2b_dninjs.append([bb - bi - bj for bj in bs]) - else: + if self.vectorized: N = self.N d2b_dninjs = np.full((N, N), bb, float) d2b_dninjs -= bs d2b_dninjs = d2b_dninjs.transpose() d2b_dninjs -= bs + else: + d2b_dninjs = [] + for bi in bs: + d2b_dninjs.append([bb - bi - bj for bj in bs]) return d2b_dninjs @property @@ -2812,11 +2810,8 @@ def d3b_dzizjzks(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: - return zeros((N, N, N)) - + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] + @property def d3b_dninjnks(self): r'''Helper method for calculating the third partial mole number @@ -2839,16 +2834,7 @@ def d3b_dninjnks(self): ''' bs = self.bs n6b = -6.0*self.b - if self.scalar: - bs2 = [bi + bi for bi in bs] - d3b_dninjnks = [] - for bi2 in bs2: - d3b_dnjnks = [] - for bj2 in bs2: - base = n6b + bi2 + bj2 - d3b_dnjnks.append([base + bk2 for bk2 in bs2]) - d3b_dninjnks.append(d3b_dnjnks) - else: + if self.vectorized: bs2 = 2.0*self.bs N = self.N d3b_dninjnks = np.full((N, N, N), n6b) @@ -2857,6 +2843,15 @@ def d3b_dninjnks(self): d3b_dninjnks += bs2 d3b_dninjnks = d3b_dninjnks.transpose((0, 2, 1)) d3b_dninjnks += bs2 + else: + bs2 = [bi + bi for bi in bs] + d3b_dninjnks = [] + for bi2 in bs2: + d3b_dnjnks = [] + for bj2 in bs2: + base = n6b + bi2 + bj2 + d3b_dnjnks.append([base + bk2 for bk2 in bs2]) + d3b_dninjnks.append(d3b_dnjnks) return d3b_dninjnks @property @@ -2878,11 +2873,8 @@ def d3epsilon_dzizjzks(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: - return zeros((N, N, N)) - + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] + @property def d3delta_dzizjzks(self): r'''Helper method for calculating the third composition derivatives @@ -2903,11 +2895,8 @@ def d3delta_dzizjzks(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: - return zeros((N, N, N)) - + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] + @property def da_alpha_dzs(self): r'''Helper method for calculating the composition derivatives of @@ -2931,9 +2920,7 @@ def da_alpha_dzs(self): a_alpha_j_rows = self.a_alpha_j_rows except: a_alpha_j_rows = self._a_alpha_j_rows - if self.scalar: - return [i + i for i in a_alpha_j_rows] - return 2.0*a_alpha_j_rows + return 2.0*a_alpha_j_rows if self.vectorized else [i + i for i in a_alpha_j_rows] @property def da_alpha_dns(self): @@ -2959,9 +2946,7 @@ def da_alpha_dns(self): except: a_alpha_j_rows = self._a_alpha_j_rows a_alpha_n_2 = -2.0*self.a_alpha - if self.scalar: - return [2.0*t + a_alpha_n_2 for t in a_alpha_j_rows] - return 2.0*a_alpha_j_rows + a_alpha_n_2 + return 2.0*a_alpha_j_rows + a_alpha_n_2 if self.vectorized else [2.0*t + a_alpha_n_2 for t in a_alpha_j_rows] @property def dna_alpha_dns(self): @@ -2987,9 +2972,7 @@ def dna_alpha_dns(self): except: a_alpha_j_rows = self._a_alpha_j_rows a_alpha = self.a_alpha - if self.scalar: - return [t + t - a_alpha for t in a_alpha_j_rows] - return 2.0*a_alpha_j_rows - a_alpha + return 2.0*a_alpha_j_rows - a_alpha if self.vectorized else [t + t - a_alpha for t in a_alpha_j_rows] @property def d2a_alpha_dzizjs(self): @@ -3012,10 +2995,7 @@ def d2a_alpha_dzizjs(self): This derivative is checked numerically. ''' a_alpha_ijs = self.a_alpha_ijs - if self.scalar: - return [[i+i for i in row] for row in a_alpha_ijs] - else: - return 2.0*a_alpha_ijs + return 2.0*a_alpha_ijs if self.vectorized else [[i+i for i in row] for row in a_alpha_ijs] @property def d2a_alpha_dninjs(self): @@ -3055,10 +3035,10 @@ def d2a_alpha_dninjs(self): zs = self.zs a_alpha3 = 3.0*a_alpha - if self.scalar: - hessian = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: hessian = zeros((N, N)) + else: + hessian = [[0.0]*N for _ in range(N)] for i in range(N): for j in range(i+1): if i == j: @@ -3094,7 +3074,7 @@ def d3a_alpha_dzizjzks(self): This derivative is checked numerically. ''' N = self.N - return [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] @property def d3a_alpha_dninjnks(self): @@ -3126,7 +3106,7 @@ def d3a_alpha_dninjnks(self): N = self.N zs = self.zs a_alpha6 = -6.0*a_alpha - matrix = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + matrix = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] for i in range(N): l = [] for j in range(N): @@ -3165,10 +3145,8 @@ def da_alpha_dT_dzs(self): da_alpha_dT_j_rows = self.da_alpha_dT_j_rows except: da_alpha_dT_j_rows = self._da_alpha_dT_j_rows - if self.scalar: - return [i + i for i in da_alpha_dT_j_rows] - return 2.0*da_alpha_dT_j_rows - + return 2.0*da_alpha_dT_j_rows if self.vectorized else [i + i for i in da_alpha_dT_j_rows] + @property def da_alpha_dT_dns(self): r'''Helper method for calculating the mole number derivatives of @@ -3197,9 +3175,7 @@ def da_alpha_dT_dns(self): except: da_alpha_dT_j_rows = self._da_alpha_dT_j_rows da_alpha_dT = self.da_alpha_dT - if self.scalar: - return [2.0*(t - da_alpha_dT) for t in da_alpha_dT_j_rows] - return 2.0*(da_alpha_dT_j_rows - da_alpha_dT) + return 2.0*(da_alpha_dT_j_rows - da_alpha_dT) if self.vectorized else [2.0*(t - da_alpha_dT) for t in da_alpha_dT_j_rows] @property def dna_alpha_dT_dns(self): @@ -3229,10 +3205,7 @@ def dna_alpha_dT_dns(self): except: da_alpha_dT_j_rows = self._da_alpha_dT_j_rows da_alpha_dT = self.da_alpha_dT - if self.scalar: - return [t + t - da_alpha_dT for t in da_alpha_dT_j_rows] - return 2.0*da_alpha_dT_j_rows - da_alpha_dT - + return 2.0*da_alpha_dT_j_rows - da_alpha_dT if self.vectorized else [t + t - da_alpha_dT for t in da_alpha_dT_j_rows] @property def d2a_alpha_dT2_dzs(self): @@ -3258,9 +3231,7 @@ def d2a_alpha_dT2_dzs(self): d2a_alpha_dT2_j_rows = self.d2a_alpha_dT2_j_rows except: d2a_alpha_dT2_j_rows = self._d2a_alpha_dT2_j_rows - if self.scalar: - return [i + i for i in d2a_alpha_dT2_j_rows] - return 2.0*d2a_alpha_dT2_j_rows + return 2.0*d2a_alpha_dT2_j_rows if self.vectorized else [i + i for i in d2a_alpha_dT2_j_rows] @property def d2a_alpha_dT2_dns(self): @@ -3288,9 +3259,7 @@ def d2a_alpha_dT2_dns(self): except: d2a_alpha_dT2_j_rows = self._d2a_alpha_dT2_j_rows d2a_alpha_dT2 = self.d2a_alpha_dT2 - if self.scalar: - return [2.0*(t - d2a_alpha_dT2) for t in d2a_alpha_dT2_j_rows] - return 2.0*(d2a_alpha_dT2_j_rows - d2a_alpha_dT2) + return 2.0*(d2a_alpha_dT2_j_rows - d2a_alpha_dT2) if self.vectorized else [2.0*(t - d2a_alpha_dT2) for t in d2a_alpha_dT2_j_rows] def dV_dzs(self, Z): r'''Calculates the molar volume composition derivative @@ -3350,7 +3319,7 @@ def dV_dzs(self, Z): [(-R*T*(V(x)**2 + V(x)*delta(x) + epsilon(x))**3*Derivative(b(x), x) + (V(x) - b(x))**2*(V(x)**2 + V(x)*delta(x) + epsilon(x))**2*Derivative(a \alpha(x), x) - (V(x) - b(x))**2*V(x)**3*a \alpha(x)*Derivative(delta(x), x) - (V(x) - b(x))**2*V(x)**2*a \alpha(x)*delta(x)*Derivative(delta(x), x) - (V(x) - b(x))**2*V(x)**2*a \alpha(x)*Derivative(epsilon(x), x) - (V(x) - b(x))**2*V(x)*a \alpha(x)*delta(x)*Derivative(epsilon(x), x) - (V(x) - b(x))**2*V(x)*a \alpha(x)*epsilon(x)*Derivative(delta(x), x) - (V(x) - b(x))**2*a \alpha(x)*epsilon(x)*Derivative(epsilon(x), x))/(-R*T*(V(x)**2 + V(x)*delta(x) + epsilon(x))**3 + 2*(V(x) - b(x))**2*V(x)**3*a \alpha(x) + 3*(V(x) - b(x))**2*V(x)**2*a \alpha(x)*delta(x) + (V(x) - b(x))**2*V(x)*a \alpha(x)*delta(x)**2 + 2*(V(x) - b(x))**2*V(x)*a \alpha(x)*epsilon(x) + (V(x) - b(x))**2*a \alpha(x)*delta(x)*epsilon(x))] ''' N = self.N - out = [0.0]*N if self.scalar else zeros(N) + out = zeros(N) if self.vectorized else [0.0]*N return eos_mix_dV_dzs(self.T, self.P, Z, self.b, self.delta, self.epsilon, self.a_alpha, self.db_dzs, self.ddelta_dzs, self.depsilon_dzs, self.da_alpha_dzs, N, out) @@ -3376,7 +3345,7 @@ def dV_dns(self, Z): dV_dns : float Molar volume mole number derivatives, [m^3/mol^2] ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N dV_dns = dxs_to_dns(self.dV_dzs(Z), self.zs, out) return dV_dns @@ -3404,7 +3373,7 @@ def dnV_dns(self, Z): [m^3/mol] ''' V = Z*R*self.T/self.P - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dn_partials(self.dV_dzs(Z), self.zs, V, out) def _d2V_dij_wrapper(self, V, d_Vs, dbs, d2bs, d_epsilons, d2_epsilons, @@ -3436,7 +3405,7 @@ def _d2V_dij_wrapper(self, V, d_Vs, dbs, d2bs, d_epsilons, d2_epsilons, x39 = x10*x38 N = self.N - hessian = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + hessian = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): for j in range(N): # TODO optimize this - symmetric, others @@ -3598,7 +3567,7 @@ def dZ_dzs(self, Z): ''' factor = self.P/(self.T*R) dV_dzs =self.dV_dzs(Z) - return [dV*factor for dV in dV_dzs] if self.scalar else dV_dzs*factor + return dV_dzs*factor if self.vectorized else [dV*factor for dV in dV_dzs] def dZ_dns(self, Z): r'''Calculates the compressibility mole number derivatives @@ -3621,7 +3590,7 @@ def dZ_dns(self, Z): dZ_dns : float Compressibility number derivatives, [1/mol] ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dZ_dzs(Z), self.zs, out) def dnZ_dns(self, Z): @@ -3647,7 +3616,7 @@ def dnZ_dns(self, Z): Partial compressibility of the mixture of the specified phase, [-] ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dn_partials(self.dZ_dzs(Z), self.zs, Z, out) def dH_dep_dzs(self, Z): @@ -3737,7 +3706,7 @@ def dH_dep_dzs(self, Z): t1 = x10*t0*x13 t2 = 2.0*x10*x13/(x13*x3*x3 - 1.0) x3_x13 = x3*x13 - dH_dzs = [0.0]*N if self.scalar else zeros(N) + dH_dzs = zeros(N) if self.vectorized else [0.0]*N for i in range(self.N): x1 = dV_dzs[i] x11 = ddelta_dzs[i] @@ -3780,9 +3749,9 @@ def dS_dep_dzs(self, Z): dH_dep_dzs = self.dH_dep_dzs(Z) dG_dep_dzs = self.dG_dep_dzs(Z) T_inv = 1.0/self.T - if self.scalar: - return [T_inv*(dH_dep_dzs[i] - dG_dep_dzs[i]) for i in range(self.N)] - return T_inv*(dH_dep_dzs - dG_dep_dzs) + if self.vectorized: + return T_inv*(dH_dep_dzs - dG_dep_dzs) + return [T_inv*(dH_dep_dzs[i] - dG_dep_dzs[i]) for i in range(self.N)] def dS_dep_dns(self, Z): r'''Calculates the molar departure entropy mole number derivatives @@ -3805,7 +3774,7 @@ def dS_dep_dns(self, Z): dS_dep_dns : float Departure entropy mole number derivatives, [J/mol^2/K] ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dS_dep_dzs(Z), self.zs, out) def dP_dns_Vt(self, phase): @@ -3845,7 +3814,7 @@ def dP_dns_Vt(self, phase): t3 = a_alpha*t2*t2 t4 = t1*Vt -t3*(Vt*delta + Vt2 + Vt2) - dP_dns_Vt = [0.0]*N if self.scalar else zeros(N) + dP_dns_Vt = zeros(N) if self.vectorized else [0.0]*N for i in range(N): v = (t4 + t1*db_dns[i] + t3*(Vt*ddelta_dns[i] + depsilon_dns[i]) - t2*da_alpha_dns[i]) dP_dns_Vt[i] = v @@ -3899,7 +3868,7 @@ def d2P_dninjs_Vt(self, phase): t4 = 2.0*x12*x9_inv3 t5 = 2.0*x0*x7_inv*x7_inv*x7_inv - hess = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): x15 = ddelta_dns[i] x17 = -x15*Vt + x16 - depsilon_dns[i] @@ -3949,7 +3918,7 @@ def d3P_dninjnks_Vt(self, phase): d3a_alpha_dninjnks = self.d3a_alpha_dninjnks d3b_dninjnks = self.d3b_dninjnks - mat = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + mat = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] for i in range(N): for j in range(N): @@ -4041,7 +4010,7 @@ def dH_dep_dns(self, Z): dH_dep_dns : float Departure enthalpy mole number derivatives, [J/mol^2] ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dH_dep_dzs(Z), self.zs, out) def dnH_dep_dns(self, Z): @@ -4072,12 +4041,12 @@ def dnH_dep_dns(self, Z): F = self.H_dep_g except: F = self.H_dep_g - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dn_partials(self.dH_dep_dzs(Z), self.zs, F, out) def _G_dep_lnphi_d_helper(self, Z, dbs, depsilons, ddelta, dVs, da_alphas, G=True): - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return G_dep_lnphi_d_helper(self.T, self.P, self.b, self.delta, self.epsilon, self.a_alpha, self.N, Z, dbs, depsilons, ddelta, dVs, da_alphas, @@ -4269,7 +4238,7 @@ def dnG_dep_dns(self, Z): except: F = self.G_dep_g dG_dns = self.dG_dep_dns(Z) - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dns_to_dn_partials(dG_dns, F, out) def fugacity_coefficients(self, Z): @@ -4315,7 +4284,7 @@ def fugacity_coefficients(self, Z): logF = log(F) except: logF = -690.7755278982137 - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N log_phis = dns_to_dn_partials(self.dlnphi_dns(Z), logF, out) return log_phis @@ -4326,7 +4295,7 @@ def _d2_G_dep_lnphi_d2_helper(self, V, d_Vs, d2Vs, dbs, d2bs, d_epsilons, d2_eps N = self.N RT = T*R RT_inv = 1.0/RT - hess = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): for j in range(N): # x1: i @@ -4559,7 +4528,7 @@ def dlnphis_dns(self, Z): dns = self.dlnphi_dns(Z) d2ns = self.d2lnphi_dninjs(Z) ans = d2ns_to_dn2_partials(d2ns, dns) - if not self.scalar: ans = array(ans) + if self.vectorized: ans = array(ans) return ans def dlnfugacities_dns(self, phase): @@ -4664,7 +4633,7 @@ def dfugacities_dns(self, phase): P = self.P N = self.N - matrix = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + matrix = zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] for i in range(N): phi_P = P*phis[i] ziPphi = phi_P*zs[i] @@ -4741,7 +4710,7 @@ def _d2_A_dep_d2_helper(self, V, d_Vs, d2Vs, dbs, d2bs, d_epsilons, b = self.b N = self.N RT = T*R - hess = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): for j in range(N): @@ -4878,7 +4847,7 @@ def dA_dep_dns_Vt(self, phase): x18 = x0*x10 x19 = x14*catanh(x11*x8**-0.5).real - jac = [0.0]*N if self.scalar else zeros(N) + jac = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x20 = ddelta_dns[i] x21 = x20*x4 - 2*depsilon_dns[i] @@ -4916,7 +4885,7 @@ def d2A_dep_dninjs_Vt(self, phase): dP_dns_Vt = self.dP_dns_Vt(phase) d2P_dninjs_Vt = self.d2P_dninjs_Vt(phase) - hess = [[0.0]*N for i in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] for i in range(N): for j in range(i+1): @@ -5031,7 +5000,7 @@ def dScomp_dns(self, phase): tot += zs[i]*logzs[i] const = R*self.T/self.P - out = [0.0]*N if self.scalar else zeros(N) + out = zeros(N) if self.vectorized else [0.0]*N for i in range(N): out[i] = mRT*(tot - logzs[i]) + const*dP_dns_Vt[i] return out @@ -5050,7 +5019,7 @@ def d2Scomp_dninjs(self, phase): logzs = [log(zi) for zi in zs] - hess = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): row = [] for j in range(N): @@ -5087,7 +5056,7 @@ def d2A_dninjs_Vt(self, phase): d2A_dep_dninjs_Vt = self.d2A_dep_dninjs_Vt(phase) d2Scomp_dninjs = self.d2Scomp_dninjs(phase) - hess = [[0.0]*N for i in range(N)] if self.scalar else zeros((N, N)) + hess = zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] for i in range(N): for j in range(N): hess[i][j] = d2Scomp_dninjs[i][j] + d2A_dep_dninjs_Vt[i][j] @@ -5098,12 +5067,12 @@ def d2nA_dninjs_Vt(self, phase): d2ns = [[i+j for i, j in zip(r1, r2)] for r1, r2 in zip(self.d2A_dep_dninjs_Vt(phase), self.d2Scomp_dninjs(phase))] dns = [i+j for i, j in zip(self.dA_dep_dns_Vt(phase), self.dScomp_dns(phase))] ans = d2ns_to_dn2_partials(d2ns, dns) - if not self.scalar: ans = array(ans) + if self.vectorized: ans = array(ans) return ans def d2A_dninjs_Vt_another(self, phase): d2ns = [[i+j for i, j in zip(r1, r2)] for r1, r2 in zip(self.d2A_dep_dninjs_Vt(phase), self.d2Scomp_dninjs(phase))] - if not self.scalar: d2ns = array(d2ns) + if self.vectorized: d2ns = array(d2ns) return d2ns # dns = [i+j for i, j in zip(self.dA_dep_dns_Vt(phase), self.dScomp_dns(phase))] # return d2ns_to_dn2_partials(d2ns, dns) @@ -5566,7 +5535,7 @@ def dlnphis_dP(self, phase): t50 = 1.0/(x0*x0) N = self.N - dlnphis_dPs = [0.0]*N if self.scalar else zeros(N) + dlnphis_dPs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): # number dependent calculations x1 = dV_dns[i] # Derivative(x0, n) @@ -5688,7 +5657,7 @@ def dlnphis_dT(self, phase): x34 = x7*self.da_alpha_dT x35 = 8*x13*x29*x5/x17**2 - dlnphis_dTs = [0.0]*N if self.scalar else zeros(N) + dlnphis_dTs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x2 = d2V_dTdns[i] x8 = x2*x7 @@ -5737,9 +5706,7 @@ def dlnphis_dzs(self, Z): ''' d2dxs = self.d2lnphi_dzizjs(Z) d2ns = d2xs_to_dxdn_partials(d2dxs, self.zs) - if self.scalar: - return d2ns - return array(d2ns) + return array(d2ns) if self.vectorized else d2ns class EpsilonZeroMixingRules: @property @@ -5760,9 +5727,7 @@ def depsilon_dzs(self): ----- This derivative is checked numerically. ''' - if self.scalar: - return [0.0]*self.N - return zeros(self.N) + return zeros(self.N) if self.vectorized else [0.0]*self.N @property def depsilon_dns(self): @@ -5782,9 +5747,7 @@ def depsilon_dns(self): ----- This derivative is checked numerically. ''' - if self.scalar: - return [0.0]*self.N - return zeros(self.N) + return zeros(self.N) if self.vectorized else [0.0]*self.N @property def d2epsilon_dzizjs(self): @@ -5805,9 +5768,7 @@ def d2epsilon_dzizjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[0.0]*N for i in range(N)] - return zeros((N, N)) + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d2epsilon_dninjs(self): @@ -5828,9 +5789,7 @@ def d2epsilon_dninjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[0.0]*N for i in range(N)] - return zeros((N, N)) + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d3epsilon_dninjnks(self): @@ -5853,9 +5812,7 @@ def d3epsilon_dninjnks(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[[0.0]*N for _ in range(N)] for _ in range(N)] - return zeros((N, N, N)) + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] # # Python 2/3 compatibility # try: @@ -6142,10 +6099,7 @@ def _zeros2d(self): def _zeros3d(self): N = self.N - if self.scalar: - return [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: - return zeros((N, N, N)) + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] @property def a_alpha_roots(self): @@ -6329,11 +6283,11 @@ def __init__(self, zs, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list - if scalar: - self.zeros2d = zeros2d = [[0.0]*N for _ in range(N)] - else: + self.vectorized = vectorized = type(zs) is ndarray + if vectorized: self.zeros2d = zeros2d = zeros((N, N)) + else: + self.zeros2d = zeros2d = [[0.0]*N for _ in range(N)] if kijs is None: kijs = zeros2d self.kijs = kijs @@ -6584,12 +6538,12 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -6598,16 +6552,16 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, self.V = V c1R2_c2R, c2R = self.c1R2_c2R, self.c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - b = float((bs*zs).sum()) self.b = self.delta = b @@ -6618,11 +6572,11 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, def _fast_init_specific(self, other): b = 0.0 - if self.scalar: + if self.vectorized: + b = float((self.bs*self.zs).sum()) + else: for bi, zi in zip(self.bs, self.zs): b += bi*zi - else: - b = float((self.bs*self.zs).sum()) self.b = self.delta = b def a_alphas_vectorized(self, T): @@ -6649,7 +6603,7 @@ def a_alphas_vectorized(self, T): [0.1449810919468, 0.30019773677] ''' return RK_a_alphas_vectorized(T, self.Tcs, self.ais, - a_alphas=[0.0]*self.N if self.scalar else zeros(self.N)) + a_alphas=zeros(self.N) if self.vectorized else [0.0]*self.N) def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -6688,10 +6642,10 @@ def a_alpha_and_derivatives_vectorized(self, T): ([0.1449810919468, 0.30019773677], [-0.000630352573681, -0.00130520755121], [8.2219900915e-06, 1.7024446320e-05]) ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return RK_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) @@ -6749,9 +6703,7 @@ def ddelta_dns(self): This derivative is checked numerically. ''' b = self.b - if self.scalar: - return [(bi - b) for bi in self.bs] - return self.bs - b + return self.bs - b if self.vectorized else [(bi - b) for bi in self.bs] @property def d2delta_dzizjs(self): @@ -6772,10 +6724,7 @@ def d2delta_dzizjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[0.0]*N for i in range(N)] - else: - return zeros((N, N)) + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d2delta_dninjs(self): @@ -6818,7 +6767,7 @@ def d3delta_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N) ] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N) ] for _ in range(N)] return RK_d3delta_dninjnks(self.b, self.bs, N, out) @@ -6921,7 +6870,7 @@ class PRMIX(GCEOSMIX, PR): 'V_g', 'Z_g', 'PIP_g', 'dP_dT_g', 'dP_dV_g', 'dV_dT_g', 'dV_dP_g', 'dT_dV_g', 'dT_dP_g', 'd2P_dT2_g', 'd2P_dV2_g', 'd2P_dTdV_g', 'H_dep_g', 'S_dep_g', 'G_dep_g', 'Cp_dep_g', 'Cv_dep_g', 'phase', 'kwargs',) - eos_mix_slots = ('lnphis_l', 'phis_l', 'fugacities_l', 'N', 'Tcs', 'Pcs', 'omegas', 'zs', 'scalar', 'kijs', 'one_minus_kijs', 'bs', + eos_mix_slots = ('lnphis_l', 'phis_l', 'fugacities_l', 'N', 'Tcs', 'Pcs', 'omegas', 'zs', 'vectorized', 'kijs', 'one_minus_kijs', 'bs', 'bs', 'ais', 'a_alphas', 'da_alpha_dTs', 'd2a_alpha_dT2s', 'a_alpha_roots', 'a_alpha_j_rows', 'da_alpha_dT_j_rows', 'lnphis_g', 'phis_g', 'fugacities_g') @@ -6935,12 +6884,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -6951,19 +6900,19 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, # optimization, unfortunately c1R2_c2R, c2R = self.c1R2_c2R, self.c2R # Also tried to store the inverse of Pcs, without success - slows it down - self.scalar = scalar = type(Tcs) is list - if scalar: + self.vectorized = vectorized = type(Tcs) is ndarray + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + self.kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] self.kappas = [omega*(-0.26992*omega + 1.54226) + 0.37464 for omega in omegas] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - self.kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 - b = float((bs*zs).sum()) self.b = b @@ -6975,12 +6924,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.fugacities() def _fast_init_specific(self, other): - if self.scalar: + if self.vectorized: + b = float(prodsum(self.bs, self.zs)) + else: b = 0.0 for bi, zi in zip(self.bs, self.zs): b += bi*zi - else: - b = float(prodsum(self.bs, self.zs)) self.kappas, self.b, self.delta, self.epsilon = other.kappas, b, 2.0*b, -b*b def a_alphas_vectorized(self, T): @@ -7002,7 +6951,7 @@ def a_alphas_vectorized(self, T): Coefficient calculated by EOS-specific method, [J^2/mol^2/Pa] ''' return PR_a_alphas_vectorized(T, self.Tcs, self.ais, self.kappas, - a_alphas=[0.0]*self.N if self.scalar else zeros(self.N)) + a_alphas=zeros(self.N) if self.vectorized else [0.0]*self.N) def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -7039,10 +6988,10 @@ def a_alpha_and_derivatives_vectorized(self, T): EOS-specific method, [J^2/mol^2/Pa/K**2] ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return PR_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, self.kappas, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) @@ -7069,11 +7018,11 @@ def d3a_alpha_dT3(self): tot = 0.0 zs = self.zs vs = self.d3a_alpha_dT3_vectorized(self.T) - if self.scalar: + if self.vectorized: + tot += float(dot(zs, vs)) + else: for i in range(self.N): tot += zs[i]*vs[i] - else: - tot += float(dot(zs, vs)) self._d3a_alpha_dT3 = tot return tot @@ -7097,10 +7046,9 @@ def d3a_alpha_dT3_vectorized(self, T): T_inv = 1.0/T N = self.N - d3a_alpha_dT3s = [0.0]*N if self.scalar else zeros(N) + d3a_alpha_dT3s = zeros(N) if self.vectorized else [0.0]*N for i in range(N): kappa = kappas[i] - x0 = 1.0/Tcs[i] x1 = sqrt(T*x0) v = (-ais[i]*0.75*kappa*(kappa*x0 - x1*(kappa*(x1 - 1.0) - 1.0)*T_inv)*T_inv*T_inv) @@ -7137,7 +7085,7 @@ def fugacity_coefficients(self, Z): ''' N = self.N return PR_lnphis(self.T, self.P, Z, self.b, self.a_alpha, self.bs, self.a_alpha_j_rows, N, - lnphis=[0.0]*N if self.scalar else zeros(N)) + lnphis=zeros(N) if self.vectorized else [0.0]*N) a_alpha = self.a_alpha # a_alpha_ijs = self.a_alpha_ijs T_inv = 1.0/self.T @@ -7173,11 +7121,11 @@ def fugacity_coefficients(self, Z): return [0.0]*self.N t51 = (x4 + (Z - 1.0)*two_root_two_B)/(b*two_root_two_B) - if self.scalar: + if self.vectorized: + return bs*t51 - x0 - t50*a_alpha_j_rows + else: return [bs[i]*t51 - x0 - t50*a_alpha_j_rows[i] for i in range(self.N)] - else: - return bs*t51 - x0 - t50*a_alpha_j_rows def dlnphis_dT(self, phase): r'''Formula for calculating the temperature derivaitve of @@ -7264,7 +7212,7 @@ def dlnphis_dT(self, phase): a_alpha_j_rows = self._a_alpha_j_rows da_alpha_dT_j_rows = self._da_alpha_dT_j_rows - d_lnphis_dTs = [0.0]*N if self.scalar else zeros(N) + d_lnphis_dTs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): d_lnphis_dTs[i] = x52 + bs[i]*x58 + x50*(x59*a_alpha_j_rows[i] + da_alpha_dT_j_rows[i]) return d_lnphis_dTs @@ -7316,7 +7264,7 @@ def dlnphis_dP(self, phase): x50 = -2.0/a_alpha N = self.N - d_lnphi_dPs = [0.0]*N if self.scalar else zeros(N) + d_lnphi_dPs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x3 = bs[i]*x2 x10 = x50*a_alpha_j_rows[i] @@ -7608,7 +7556,7 @@ def dlnphis_dzs(self, Z): G_inv = 1.0/G logG = log(G) C_inv = 1.0/C - dlnphis_dxs = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + dlnphis_dxs = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] t61s = [C_inv*dC_dxi for dC_dxi in dC_dxs] for i in range(N): dD_dxs_i = dD_dxs[i] @@ -7643,7 +7591,7 @@ def ddelta_dzs(self): This derivative is checked numerically. ''' N = self.N - return PR_ddelta_dzs(self.bs, N, out=[0.0]*N if self.scalar else zeros(N)) + return PR_ddelta_dzs(self.bs, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def ddelta_dns(self): @@ -7664,7 +7612,7 @@ def ddelta_dns(self): This derivative is checked numerically. ''' N = self.N - return PR_ddelta_dns(self.bs, self.b, N, out=[0.0]*N if self.scalar else zeros(N)) + return PR_ddelta_dns(self.bs, self.b, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def d2delta_dzizjs(self): @@ -7685,9 +7633,7 @@ def d2delta_dzizjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - return [[0.0]*N for i in range(N)] - return zeros((N, N)) + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d2delta_dninjs(self): @@ -7708,7 +7654,7 @@ def d2delta_dninjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_d2delta_dninjs(self.b, self.bs, N, out) @property @@ -7732,7 +7678,7 @@ def d3delta_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N) ] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N) ] for _ in range(N)] return PR_d3delta_dninjnks(self.b, self.bs, N, out) @property @@ -7754,7 +7700,7 @@ def depsilon_dzs(self): This derivative is checked numerically. ''' N = self.N - return PR_depsilon_dzs(self.b, self.bs, N, out=[0.0]*N if self.scalar else zeros(N)) + return PR_depsilon_dzs(self.b, self.bs, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def depsilon_dns(self): @@ -7775,7 +7721,7 @@ def depsilon_dns(self): This derivative is checked numerically. ''' N = self.N - return PR_depsilon_dns(self.b, self.bs, N, out=[0.0]*N if self.scalar else zeros(N)) + return PR_depsilon_dns(self.b, self.bs, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def d2epsilon_dzizjs(self): @@ -7796,7 +7742,7 @@ def d2epsilon_dzizjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_d2epsilon_dzizjs(self.b, self.bs, N, out) @property @@ -7819,7 +7765,7 @@ def d2epsilon_dninjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_d2epsilon_dninjs(self.b, self.bs, N, out) @@ -7846,7 +7792,7 @@ def d3epsilon_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N) ] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N) ] for _ in range(N)] return PR_d3epsilon_dninjnks(self.b, self.bs, N, out) def solve_T(self, P, V, quick=True, solution=None): @@ -7979,12 +7925,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.T = T @@ -7992,30 +7938,30 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.V = V c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: - b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] - self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] - else: + if vectorized: b0s = c2R*Tcs/Pcs self.ais = c1R2_c2R*Tcs*b0s + else: + b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] + self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] if cs is None: - if scalar: - cs = [0.0]*N - else: + if vectorized: cs = zeros(N) - if scalar: + else: + cs = [0.0]*N + if vectorized: + self.kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + bs = b0s - cs + else: self.kappas = [omega*(-0.26992*omega + 1.54226) + 0.37464 for omega in omegas] b0, c = 0.0, 0.0 for i in range(N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] bs = [b0s[i] - cs[i] for i in range(N)] - else: - self.kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) - bs = b0s - cs self.kwargs = {'kijs': kijs, 'cs': cs} self.cs = cs @@ -8035,14 +7981,14 @@ def _fast_init_specific(self, other): self.kappas = other.kappas zs = self.zs self.b0s = b0s = other.b0s - if self.scalar: + if self.vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + else: b0, c = 0.0, 0.0 for i in range(self.N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) self.c = c self.b = b0 - c self.delta = 2.0*(c + b0) @@ -8069,7 +8015,7 @@ def ddelta_dzs(self): ''' N = self.N return PR_translated_ddelta_dzs(self.b0s, self.cs, N, - [0.0]*N if self.scalar else zeros(N)) + zeros(N) if self.vectorized else [0.0]*N) # Zero in both cases d2delta_dzizjs = PRMIX.d2delta_dzizjs @@ -8095,7 +8041,7 @@ def ddelta_dns(self): This derivative is checked numerically. ''' N = self.N - return PR_translated_ddelta_dns(self.b0s, self.cs, self.delta, N, [0.0]*N if self.scalar else zeros(N)) + return PR_translated_ddelta_dns(self.b0s, self.cs, self.delta, N, zeros(N) if self.vectorized else [0.0]*N) @property @@ -8118,7 +8064,7 @@ def d2delta_dninjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_translated_d2delta_dninjs(self.b0s, self.cs, self.b, self.c, self.delta, N, out) @property @@ -8144,7 +8090,7 @@ def d3delta_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] return PR_translated_d3delta_dninjnks(self.b0s, self.cs, self.delta, N, out) @property @@ -8168,7 +8114,7 @@ def depsilon_dzs(self): ''' N = self.N return PR_translated_depsilon_dzs(self.c, self.b, self.b0s, self.cs, N, - [0.0]*N if self.scalar else zeros(N)) + zeros(N) if self.vectorized else [0.0]*N) @property def depsilon_dns(self): @@ -8192,7 +8138,7 @@ def depsilon_dns(self): ''' c, b = self.c, self.b N, b0s, cs = self.N, self.b0s, self.cs - return PR_translated_depsilon_dns(b, c, b0s, cs, N, out=([0.0]*N if self.scalar else zeros(N))) + return PR_translated_depsilon_dns(b, c, b0s, cs, N, out=(zeros(N) if self.vectorized else [0.0]*N)) @property def d2epsilon_dzizjs(self): @@ -8213,7 +8159,7 @@ def d2epsilon_dzizjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_translated_d2epsilon_dzizjs(self.b0s, self.cs, N=N, out=out) d3epsilon_dzizjzks = GCEOSMIX.d3epsilon_dzizjzks # Zeros @@ -8243,7 +8189,7 @@ def d2epsilon_dninjs(self): ''' # Not trusted yet - numerical check does not have enough digits N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return PR_translated_d2epsilon_dninjs(self.b0s, self.cs, self.b, self.c, N, out=out) @property @@ -8277,7 +8223,7 @@ def d3epsilon_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] return PR_translated_d3epsilon_dninjnks(self.b0s, self.cs, self.b, self.c, self.epsilon, N, out) @@ -8382,12 +8328,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.T = T @@ -8395,31 +8341,31 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.V = V c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: - b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] - self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] - else: + if vectorized: b0s = c2R*Tcs/Pcs self.ais = c1R2_c2R*Tcs*b0s + else: + b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] + self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] if cs is None: - if scalar: - cs = [0.0]*N - else: + if vectorized: cs = zeros(N) + else: + cs = [0.0]*N - if scalar: + if vectorized: + self.kappas = omegas*(omegas*(0.1063*omegas - 0.2721) + 1.4996) + 0.3919 + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + bs = b0s - cs + else: self.kappas = [omega*(omega*(0.1063*omega - 0.2721) + 1.4996) + 0.3919 for omega in omegas] b0, c = 0.0, 0.0 for i in range(N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] bs = [b0s[i] - cs[i] for i in range(N)] - else: - self.kappas = omegas*(omegas*(0.1063*omegas - 0.2721) + 1.4996) + 0.3919 - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) - bs = b0s - cs self.kwargs = {'kijs': kijs, 'cs': cs} self.cs = cs @@ -8553,12 +8499,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.T = T @@ -8567,19 +8513,19 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, c1R2_c2R, c2R = self.c1R2_c2R, self.c2R - if scalar: - b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] - self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] - else: + if vectorized: b0s = c2R*Tcs/Pcs self.ais = c1R2_c2R*Tcs*b0s + else: + b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] + self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] if cs is None: - if scalar: + if vectorized: + cs = R*Tcs/Pcs*(0.0198*npmin(npmax(omegas, -0.01), 1.48) - 0.0065) + else: cs = [R*Tcs[i]/Pcs[i]*(0.0198*min(max(omegas[i], -0.01), 1.48) - 0.0065) for i in range(N)] - else: - cs = R*Tcs/Pcs*(0.0198*npmin(npmax(omegas, -0.01), 1.48) - 0.0065) if alpha_coeffs is None: alpha_coeffs = [] for i in range(N): @@ -8592,16 +8538,16 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.alpha_coeffs = alpha_coeffs self.cs = cs - if scalar: + if vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + bs = b0s - cs + else: b0, c = 0.0, 0.0 for i in range(N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] bs = [b0s[i] - cs[i] for i in range(N)] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) - bs = b0s - cs self.b0s = b0s self.bs = bs @@ -8619,14 +8565,14 @@ def _fast_init_specific(self, other): zs = self.zs self.b0s = b0s = other.b0s - if self.scalar: + if self.vectorized: + b0 = float((zs*b0s).sum()) + c = float((zs*cs).sum()) + else: b0, c = 0.0, 0.0 for i in range(self.N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] - else: - b0 = float((zs*b0s).sum()) - c = float((zs*cs).sum()) self.c = c self.b = b0 - c @@ -8740,29 +8686,29 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} self.T = T self.P = P self.V = V - if self.scalar: - self.ais = [self.c1*R2*Tc*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] - self.bs = [self.c2*R*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] - ms = [omega*(1.574 - 0.176*omega) + 0.480 for omega in omegas] - b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) - else: + if self.vectorized: Tc_Pc_ratio = Tcs/Pcs self.ais = self.c1R2*Tcs*Tc_Pc_ratio self.bs = bs = self.c2R*Tc_Pc_ratio ms = omegas*(1.574 - 0.176*omegas) + 0.480 b = float((bs*zs).sum()) + else: + self.ais = [self.c1*R2*Tc*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] + self.bs = [self.c2*R*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] + ms = [omega*(1.574 - 0.176*omega) + 0.480 for omega in omegas] + b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) self.b = b self.ms = ms @@ -8774,10 +8720,10 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, def _fast_init_specific(self, other): self.ms = other.ms - if self.scalar: - self.b = b = sum([bi*zi for bi, zi in zip(self.bs, self.zs)]) - else: + if self.vectorized: self.b = b = float((self.bs*self.zs).sum()) + else: + self.b = b = sum([bi*zi for bi, zi in zip(self.bs, self.zs)]) self.delta = b def a_alphas_vectorized(self, T): @@ -8799,7 +8745,7 @@ def a_alphas_vectorized(self, T): Coefficient calculated by EOS-specific method, [J^2/mol^2/Pa] ''' return SRK_a_alphas_vectorized(T, self.Tcs, self.ais, self.ms, - a_alphas=[0.0]*self.N if self.scalar else zeros(self.N)) + a_alphas=zeros(self.N) if self.vectorized else [0.0]*self.N) def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -8835,10 +8781,10 @@ def a_alpha_and_derivatives_vectorized(self, T): EOS-specific method, [J^2/mol^2/Pa/K**2] ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return SRK_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, self.ms, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) @@ -8872,7 +8818,7 @@ def fugacity_coefficients(self, Z): ''' N = self.N return SRK_lnphis(self.T, self.P, Z, self.b, self.a_alpha, self.bs, self.a_alpha_j_rows, N, - lnphis=[0.0]*N if self.scalar else zeros(N)) + lnphis=zeros(N) if self.vectorized else [0.0]*N) def dlnphis_dT(self, phase): @@ -8936,7 +8882,7 @@ def dlnphis_dT(self, phase): x52 = (dZ_dT + x2*x4)/(x6 - Z) # Composition stuff - d_lnphis_dTs = [0.0]*N if self.scalar else zeros(N) + d_lnphis_dTs = zeros(N) if self.vectorized else [0.0]*N a_alpha_j_rows = self.a_alpha_j_rows @@ -8999,7 +8945,7 @@ def dlnphis_dP(self, phase): x10 = a_alpha*x9*(self.P*dZ_dP*x9 - 1.0)*RT_inv*RT_inv/(x5*x9 + 1.0) x50 = 2.0/a_alpha - d_lnphi_dPs = [0.0]*N if self.scalar else zeros(N) + d_lnphi_dPs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x8 = x50*a_alpha_j_rows[i] x3 = bs[i]*x2 @@ -9109,17 +9055,17 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, T=None, P=None, V=N self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if cs is None: - if scalar: - cs = [0.0]*N - else: + if vectorized: cs = zeros(N) - if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] else: + cs = [0.0]*N + if kijs is None: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs, 'cs': cs} @@ -9130,26 +9076,26 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, T=None, P=None, V=N c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: - b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] - self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] - self.ms = [0.480 + omega*(1.574 - 0.176*omega) for omega in omegas] - else: + if vectorized: b0s = c2R*Tcs/Pcs self.ais = c1R2_c2R*Tcs*b0s self.ms = 0.480 + omegas*(1.574 - 0.176*omegas) + else: + b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] + self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] + self.ms = [0.480 + omega*(1.574 - 0.176*omega) for omega in omegas] self.cs = cs - if scalar: + if vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + bs = b0s - cs + else: b0, c = 0.0, 0.0 for i in range(N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] bs = [b0s[i] - cs[i] for i in range(N)] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) - bs = b0s - cs self.b0s = b0s self.bs = bs @@ -9167,14 +9113,14 @@ def _fast_init_specific(self, other): zs = self.zs self.b0s = b0s = other.b0s - if self.scalar: + if self.vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + else: b0, c = 0.0, 0.0 for i in range(self.N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) self.c = c self.b = b0 - c @@ -9201,9 +9147,7 @@ def ddelta_dzs(self): This derivative is checked numerically. ''' b0s, cs = self.b0s, self.cs - if self.scalar: - return [(2.0*cs[i] + b0s[i]) for i in range(self.N)] - return 2.0*cs + b0s + return 2.0*cs + b0s if self.vectorized else [(2.0*cs[i] + b0s[i]) for i in range(self.N)] # Zero in both cases d2delta_dzizjs = PRMIX.d2delta_dzizjs @@ -9229,7 +9173,7 @@ def ddelta_dns(self): This derivative is checked numerically. ''' N = self.N - return SRK_translated_ddelta_dns(self.b0s, self.cs, self.delta, N, out=[0.0]*N if self.scalar else zeros(N)) + return SRK_translated_ddelta_dns(self.b0s, self.cs, self.delta, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def d2delta_dninjs(self): @@ -9251,7 +9195,7 @@ def d2delta_dninjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return SRK_translated_d2delta_dninjs(self.b0s, self.cs, self.b, self.c, self.delta, N, out) @property @@ -9277,7 +9221,7 @@ def d3delta_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] return SRK_translated_d3delta_dninjnks(self.b0s, self.cs, self.b, self.c, self.delta, N, out) @property @@ -9300,7 +9244,7 @@ def depsilon_dzs(self): This derivative is checked numerically. ''' N = self.N - out = [0.0]*N if self.scalar else zeros(N) + out = zeros(N) if self.vectorized else [0.0]*N return SRK_translated_depsilon_dzs(self.b0s, self.cs, self.b, self.c, N, out) @property @@ -9324,7 +9268,7 @@ def depsilon_dns(self): This derivative is checked numerically. ''' N = self.N - return SRK_translated_depsilon_dns(self.b, self.c, self.b0s, self.cs, N, out=[0.0]*N if self.scalar else zeros(N)) + return SRK_translated_depsilon_dns(self.b, self.c, self.b0s, self.cs, N, out=zeros(N) if self.vectorized else [0.0]*N) @property def d2epsilon_dzizjs(self): @@ -9345,7 +9289,7 @@ def d2epsilon_dzizjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return SRK_translated_d2epsilon_dzizjs(self.b0s, self.cs, self.b, self.c, N, out=out) d3epsilon_dzizjzks = GCEOSMIX.d3epsilon_dzizjzks # Zeros @@ -9370,7 +9314,7 @@ def d2epsilon_dninjs(self): This derivative is checked numerically. ''' N = self.N - out = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + out = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] return SRK_translated_d2epsilon_dninjs(self.b0s, self.cs, self.b, self.c, N, out) @property @@ -9405,7 +9349,7 @@ def d3epsilon_dninjnks(self): This derivative is checked numerically. ''' N = self.N - out = [[[0.0]*N for _ in range(N)] for _ in range(N)] if self.scalar else zeros((N, N, N)) + out = zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N)] for _ in range(N)] return SRK_translated_d3epsilon_dninjnks(self.b0s, self.cs, self.b, self.c, self.epsilon, N, out) @@ -9527,12 +9471,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.T = T @@ -9540,19 +9484,19 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.V = V c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: - b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] - self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] - else: + if vectorized: b0s = c2R*Tcs/Pcs self.ais = c1R2_c2R*Tcs*b0s + else: + b0s = [c2R*Tcs[i]/Pcs[i] for i in cmps] + self.ais = [c1R2_c2R*Tcs[i]*b0s[i] for i in cmps] if cs is None: - if scalar: + if vectorized: + cs = R*Tcs/Pcs*(0.0172*npmin(npmax(omegas, -0.01), 1.46) + 0.0096) + else: cs = [R*Tcs[i]/Pcs[i]*(0.0172*min(max(omegas[i], -0.01), 1.46) + 0.0096) for i in range(N)] - else: - cs = R*Tcs/Pcs*(0.0172*npmin(npmax(omegas, -0.01), 1.46) + 0.0096) if alpha_coeffs is None: alpha_coeffs = [] @@ -9566,16 +9510,16 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.alpha_coeffs = alpha_coeffs self.cs = cs - if scalar: + if vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + bs = b0s - cs + else: b0, c = 0.0, 0.0 for i in range(N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] bs = [b0s[i] - cs[i] for i in range(N)] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) - bs = b0s - cs self.b0s = b0s self.bs = bs @@ -9593,14 +9537,14 @@ def _fast_init_specific(self, other): zs = self.zs self.b0s = b0s = other.b0s - if self.scalar: + if self.vectorized: + b0 = float((b0s*zs).sum()) + c = float((cs*zs).sum()) + else: b0, c = 0.0, 0.0 for i in range(self.N): b0 += b0s[i]*zs[i] c += cs[i]*zs[i] - else: - b0 = float((b0s*zs).sum()) - c = float((cs*zs).sum()) self.c = c self.b = b0 - c @@ -9737,7 +9681,7 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, cs=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: kijs = [[0.0]*N for i in cmps] self.kijs = kijs @@ -9886,7 +9830,7 @@ def __init__(self, Tcs, Pcs, omegas, zs, alpha_coeffs, ge_model, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: kijs = [[0.0]*N for i in cmps] if cs is None: @@ -10048,12 +9992,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -10063,7 +10007,15 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, c1R2_c2R, c2R = self.c1R2_c2R, self.c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + self.kappas = kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 + b = float((bs*zs).sum()) + high_omega_idxs = npwhere(omegas > 0.491) + high_omegas = omegas[high_omega_idxs] + kappas[high_omega_idxs] = high_omegas*(high_omegas*(0.016666*high_omegas - 0.164423) + 1.48503) + 0.379642 + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] self.kappas = kappas = [omega*(-0.26992*omega + 1.54226) + 0.37464 for omega in omegas] @@ -10073,15 +10025,6 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, b = 0.0 for i in cmps: b += bs[i]*zs[i] - - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - self.kappas = kappas = omegas*(-0.26992*omegas + 1.54226) + 0.37464 - b = float((bs*zs).sum()) - high_omega_idxs = npwhere(omegas > 0.491) - high_omegas = omegas[high_omega_idxs] - kappas[high_omega_idxs] = high_omegas*(high_omegas*(0.016666*high_omegas - 0.164423) + 1.48503) + 0.379642 self.b = b self.delta = 2.*b @@ -10180,12 +10123,12 @@ def __init__(self, Tcs, Pcs, zs, kijs=None, T=None, P=None, V=None, self.Tcs = Tcs self.Pcs = Pcs self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*self.N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*self.N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -10194,15 +10137,15 @@ def __init__(self, Tcs, Pcs, zs, kijs=None, T=None, P=None, V=None, self.V = V c1R2, c2R = self.c1R2, self.c2R - if self.scalar: - self.ais = [c1R2*Tc*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] - self.bs = [c2R*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] - self.b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) - else: + if self.vectorized: Tc_Pc_ratio = Tcs/Pcs self.ais = c1R2*Tcs*Tc_Pc_ratio self.bs = bs = c2R*Tc_Pc_ratio self.b = float((bs*zs).sum()) + else: + self.ais = [c1R2*Tc*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] + self.bs = [c2R*Tc/Pc for Tc, Pc in zip(Tcs, Pcs)] + self.b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) self.omegas = omegas self.solve(only_l=only_l, only_g=only_g) @@ -10210,10 +10153,10 @@ def __init__(self, Tcs, Pcs, zs, kijs=None, T=None, P=None, V=None, self.fugacities() def _fast_init_specific(self, other): - if self.scalar: - self.b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) - else: + if self.vectorized: self.b = float((self.bs*self.zs).sum()) + else: + self.b = sum(bi*zi for bi, zi in zip(self.bs, self.zs)) def a_alphas_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` for the VDW EOS. @@ -10264,10 +10207,10 @@ def a_alpha_and_derivatives_vectorized(self, T): Second temperature derivative of coefficient calculated by EOS-specific method, [J^2/mol^2/Pa/K**2] ''' - if self.scalar: - zero_array = [0.0]*self.N - else: + if self.vectorized: zero_array = zeros(self.N) + else: + zero_array = [0.0]*self.N return self.ais, zero_array, zero_array def fugacity_coefficients(self, Z): @@ -10297,7 +10240,7 @@ def fugacity_coefficients(self, Z): ''' N = self.N return VDW_lnphis(self.T, self.P, Z, self.b, self.a_alpha, self.bs, self.a_alpha_roots, N, - lnphis=[0.0]*N if self.scalar else zeros(N)) + lnphis=zeros(N) if self.vectorized else [0.0]*N) def dlnphis_dT(self, phase): r'''Formula for calculating the temperature derivaitve of @@ -10353,7 +10296,7 @@ def dlnphis_dT(self, phase): x15 = x4*(P*x13*(T_inv + x4*dZ_dT) - x14*dZ_dT)/x14 # Composition stuff - d_lnphis_dTs = [0.0]*N if self.scalar else zeros(N) + d_lnphis_dTs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x1 = (ais[i]*x0)**0.5 d_lhphi_dT = -bs[i]*x11 + x1*x5 + x1*x8 - x1*x9 + x15 @@ -10409,7 +10352,7 @@ def dlnphis_dP(self, phase): x14 = x12*x13 - 1.0 x15 = -x5*(-x13*(x12*dZ_dP - 1.0) + x14*dZ_dP)/x14 - d_lnphi_dPs = [0.0]*N if self.scalar else zeros(N) + d_lnphi_dPs = zeros(N) if self.vectorized else [0.0]*N for i in range(N): x1 = (ais[i]*a_alpha)**0.5 d_lnphi_dP = -bs[i]*x11 - x1*x6 + x1*x8 + x15 @@ -10434,11 +10377,7 @@ def ddelta_dzs(self): ----- This derivative is checked numerically. ''' - if self.scalar: - zero_array = [0.0]*self.N - else: - zero_array = zeros(self.N) - return zero_array + return zeros(self.N) if self.vectorized else [0.0]*self.N @property def ddelta_dns(self): @@ -10458,11 +10397,7 @@ def ddelta_dns(self): ----- This derivative is checked numerically. ''' - if self.scalar: - zero_array = [0.0]*self.N - else: - zero_array = zeros(self.N) - return zero_array + return zeros(self.N) if self.vectorized else [0.0]*self.N @property @@ -10484,11 +10419,7 @@ def d2delta_dzizjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - zero_array = [[0.0]*N for i in range(N)] - else: - zero_array = zeros((N, N)) - return zero_array + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d2delta_dninjs(self): @@ -10509,11 +10440,7 @@ def d2delta_dninjs(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - zero_array = [[0.0]*N for i in range(N)] - else: - zero_array = zeros((N, N)) - return zero_array + return zeros((N, N)) if self.vectorized else [[0.0]*N for i in range(N)] @property def d3delta_dninjnks(self): @@ -10536,11 +10463,7 @@ def d3delta_dninjnks(self): This derivative is checked numerically. ''' N = self.N - if self.scalar: - zero_array = [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: - zero_array = zeros((N, N, N)) - return zero_array + return zeros((N, N, N)) if self.vectorized else [[[0.0]*N for _ in range(N) ] for _ in range(N)] class PRSVMIX(PRMIX, PRSV): @@ -10669,21 +10592,21 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0]*self.N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0]*self.N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) if kappa1s is None: - if scalar: - kappa1s = [0.0 for i in range(N)] - else: + if vectorized: kappa1s = zeros(N) + else: + kappa1s = [0.0 for i in range(N)] self.kwargs = {'kijs': kijs, 'kappa1s': kappa1s} self.T = T self.P = P @@ -10692,18 +10615,18 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, c1R2_c2R, c2R = self.c1R2_c2R, self.c2R - if scalar: + if vectorized: + self.kappa0s = omegas*(omegas*(0.0196554*omegas - 0.17131848) + 1.4897153) + 0.378893 + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + b = float((bs*zs).sum()) + else: self.kappa0s = [omega*(omega*(0.0196554*omega - 0.17131848) + 1.4897153) + 0.378893 for omega in omegas] self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.kappa0s = omegas*(omegas*(0.0196554*omegas - 0.17131848) + 1.4897153) + 0.378893 - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - b = float((bs*zs).sum()) self.b = b @@ -10718,11 +10641,10 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.T = self.solve_T(self.P, self.V, solution=solution) else: self.kappa1s = [(0.0 if (T/Tc > 0.7 and self.kappa1_Tr_limit) else kappa1) for kappa1, Tc in zip(kappa1s, Tcs)] - if not scalar: - self.kappa1s = array(self.kappa1s) + if vectorized: self.kappa1s = array(self.kappa1s) self.kappas = [kappa0 + kappa1*(1 + (self.T/Tc)**0.5)*(0.7 - (self.T/Tc)) for kappa0, kappa1, Tc in zip(self.kappa0s, self.kappa1s, self.Tcs)] - if not scalar: + if vectorized: self.kappas = array(self.kappas) self.solve(only_l=only_l, only_g=only_g) @@ -10735,12 +10657,12 @@ def _fast_init_specific(self, other): self.kappa0s = other.kappa0s self.kappa1s = other.kappa1s self.kappas = other.kappas - if self.scalar: + if self.vectorized: + b = float((self.bs*self.zs).sum()) + else: b = 0.0 for bi, zi in zip(self.bs, self.zs): b += bi*zi - else: - b = float((self.bs*self.zs).sum()) self.b = b self.delta = 2.0*b self.epsilon = -b*b @@ -10765,7 +10687,7 @@ def a_alphas_vectorized(self, T): Coefficient calculated by EOS-specific method, [J^2/mol^2/Pa] ''' return PRSV_a_alphas_vectorized(T, self.Tcs, self.ais, self.kappa0s, self.kappa1s, - a_alphas=[0.0]*self.N if self.scalar else zeros(self.N)) + a_alphas=zeros(self.N) if self.vectorized else [0.0]*self.N) def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -10794,10 +10716,10 @@ def a_alpha_and_derivatives_vectorized(self, T): EOS-specific method, [J^2/mol^2/Pa/K**2] ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return PRSV_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, self.kappa0s, self.kappa1s, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) @@ -10917,30 +10839,30 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) - if scalar: + if vectorized: if kappa1s is None: - kappa1s = [0.0]*N + kappa1s = zeros(N) if kappa2s is None: - kappa2s = [0.0]*N + kappa2s = zeros(N) if kappa3s is None: - kappa3s = [0.0]*N + kappa3s = zeros(N) else: if kappa1s is None: - kappa1s = zeros(N) + kappa1s = [0.0]*N if kappa2s is None: - kappa2s = zeros(N) + kappa2s = [0.0]*N if kappa3s is None: - kappa3s = zeros(N) + kappa3s = [0.0]*N self.kwargs = {'kijs': kijs, 'kappa1s': kappa1s, 'kappa2s': kappa2s, 'kappa3s': kappa3s} self.kappa1s = kappa1s @@ -10953,18 +10875,18 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + self.kappa0s = kappa0s = omegas*(omegas*(0.0196554*omegas - 0.17131848) + 1.4897153) + 0.378893 + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] self.kappa0s = kappa0s = [omega*(omega*(0.0196554*omega - 0.17131848) + 1.4897153) + 0.378893 for omega in omegas] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - self.kappa0s = kappa0s = omegas*(omegas*(0.0196554*omegas - 0.17131848) + 1.4897153) + 0.378893 - b = float((bs*zs).sum()) self.b = b self.delta = 2.0*b @@ -10976,16 +10898,16 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.T = T = self.solve_T(self.P, self.V, solution=solution) - if scalar: + if vectorized: + Trs = T/Tcs + sqrtTrs = npsqrt(Trs) + kappas = kappa0s + ((kappa1s + kappa2s*(kappa3s - Trs)*(1. - sqrtTrs))*(1. + sqrtTrs)*(0.7 - Trs)) + else: kappas = [0.0]*N for i in cmps: Tr = T/Tcs[i] sqrtTr = sqrt(Tr) kappas[i] = kappa0s[i] + ((kappa1s[i] + kappa2s[i]*(kappa3s[i] - Tr)*(1. - sqrtTr))*(1. + sqrtTr)*(0.7 - Tr)) - else: - Trs = T/Tcs - sqrtTrs = npsqrt(Trs) - kappas = kappa0s + ((kappa1s + kappa2s*(kappa3s - Trs)*(1. - sqrtTrs))*(1. + sqrtTrs)*(0.7 - Trs)) self.kappas = kappas @@ -10999,12 +10921,12 @@ def _fast_init_specific(self, other): self.kappa2s = other.kappa2s self.kappa3s = other.kappa3s self.kappas = other.kappas - if self.scalar: + if self.vectorized: + b = float((self.bs*self.zs).sum()) + else: b = 0.0 for bi, zi in zip(self.bs, self.zs): b += bi*zi - else: - b = float((self.bs*self.zs).sum()) self.b = b self.delta = b + b self.epsilon = -b*b @@ -11030,7 +10952,7 @@ def a_alphas_vectorized(self, T): [0.0860568595, 0.20174345803] ''' return PRSV2_a_alphas_vectorized(T, self.Tcs, self.ais, self.kappa0s, self.kappa1s, self.kappa2s, self.kappa3s, - a_alphas=([0.0]*self.N if self.scalar else zeros(self.N))) + a_alphas=(zeros(self.N) if self.vectorized else [0.0]*self.N)) def a_alpha_and_derivatives_vectorized(self, T): r'''Method to calculate the pure-component `a_alphas` and their first @@ -11054,10 +10976,10 @@ def a_alpha_and_derivatives_vectorized(self, T): EOS-specific method, [J^2/mol^2/Pa/K**2] ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return PRSV2_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, self.kappa0s, self.kappa1s, self.kappa2s, self.kappa3s, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) @@ -11171,12 +11093,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -11185,16 +11107,16 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.V = V c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - b = float((bs*zs).sum()) self.b = b self.delta = 2.*b @@ -11206,12 +11128,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.fugacities() def _fast_init_specific(self, other): - if self.scalar: + if self.vectorized: + b = float((self.bs*self.zs).sum()) + else: b = 0.0 for bi, zi in zip(self.bs, self.zs): b += bi*zi - else: - b = float((self.bs*self.zs).sum()) self.b = b self.delta = 2.0*b self.epsilon = -b*b @@ -11327,12 +11249,12 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in range(N)] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in range(N)] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -11341,16 +11263,16 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, self.V = V c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - b = float((bs*zs).sum()) self.delta = self.b = b self.check_sufficient_inputs() @@ -11362,11 +11284,11 @@ def __init__(self, Tcs, Pcs, omegas, zs, kijs=None, T=None, P=None, V=None, def _fast_init_specific(self, other): b = 0.0 bs, zs = self.bs, self.zs - if self.scalar: + if self.vectorized: + b = float((bs*zs).sum()) + else: for i in range(self.N): b += bs[i]*zs[i] - else: - b = float((bs*zs).sum()) self.delta = self.b = b class APISRKMIX(SRKMIX, APISRK): @@ -11472,12 +11394,12 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, self.Pcs = Pcs self.omegas = omegas self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if kijs is None: - if scalar: - kijs = [[0.0]*N for i in cmps] - else: + if vectorized: kijs = zeros((N, N)) + else: + kijs = [[0.0]*N for i in cmps] self.kijs = kijs self.one_minus_kijs = one_minus_kijs(kijs) self.kwargs = {'kijs': kijs} @@ -11491,30 +11413,30 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, if S1s is None and omegas is None: raise ValueError('Either acentric factor of S1 is required') if S1s is None: - if scalar: - self.S1s = [omega*(1.55171 - 0.15613*omega) + 0.48508 for omega in omegas] - else: + if vectorized: self.S1s = omegas*(1.55171 - 0.15613*omegas) + 0.48508 + else: + self.S1s = [omega*(1.55171 - 0.15613*omega) + 0.48508 for omega in omegas] else: self.S1s = S1s if S2s is None: - if scalar: - S2s = [0.0]*N - else: + if vectorized: S2s = zeros(N) + else: + S2s = [0.0]*N self.S2s = S2s self.kwargs = {'S1s': self.S1s, 'S2s': self.S2s} c2R, c1R2_c2R = self.c2R, self.c1R2_c2R - if scalar: + if vectorized: + self.bs = bs = c2R*Tcs/Pcs + self.ais = c1R2_c2R*Tcs*bs + b = float((bs*zs).sum()) + else: self.bs = bs = [c2R*Tcs[i]/Pcs[i] for i in cmps] self.ais = [c1R2_c2R*Tcs[i]*bs[i] for i in cmps] b = 0.0 for i in cmps: b += bs[i]*zs[i] - else: - self.bs = bs = c2R*Tcs/Pcs - self.ais = c1R2_c2R*Tcs*bs - b = float((bs*zs).sum()) self.b = self.delta = b self.solve(only_l=only_l, only_g=only_g) if fugacities: @@ -11523,13 +11445,13 @@ def __init__(self, Tcs, Pcs, zs, omegas=None, kijs=None, T=None, P=None, V=None, def _fast_init_specific(self, other): self.S1s = other.S1s self.S2s = other.S2s - if self.scalar: - self.delta = self.b = sum([bi*zi for bi, zi in zip(self.bs, self.zs)]) - else: + if self.vectorized: self.delta = self.b = float((self.bs*self.zs).sum()) + else: + self.delta = self.b = sum([bi*zi for bi, zi in zip(self.bs, self.zs)]) def a_alphas_vectorized(self, T): - a_alphas = [0.0]*self.N if self.scalar else zeros(self.N) + a_alphas = zeros(self.N) if self.vectorized else [0.0]*self.N return APISRK_a_alphas_vectorized(T, self.Tcs, self.ais, self.S1s, self.S2s, a_alphas=a_alphas) def a_alpha_and_derivatives_vectorized(self, T): @@ -11572,10 +11494,10 @@ def a_alpha_and_derivatives_vectorized(self, T): EOS-specific method, [J^2/mol^2/Pa/K**2] ''' N = self.N - if self.scalar: - a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alphas, da_alpha_dTs, d2a_alpha_dT2s = zeros(N), zeros(N), zeros(N) + else: + a_alphas, da_alpha_dTs, d2a_alpha_dT2s = [0.0]*N, [0.0]*N, [0.0]*N return APISRK_a_alpha_and_derivatives_vectorized(T, self.Tcs, self.ais, self.S1s, self.S2s, a_alphas=a_alphas, da_alpha_dTs=da_alpha_dTs, d2a_alpha_dT2s=d2a_alpha_dT2s) diff --git a/thermo/equilibrium.py b/thermo/equilibrium.py index d541015e..58580e77 100644 --- a/thermo/equilibrium.py +++ b/thermo/equilibrium.py @@ -2351,15 +2351,15 @@ def zs_no_water(self, phase=None): water_index = self.water_index if water_index is None: return phase.zs - scalar = self.flasher.scalar - zs = list(phase.zs) if scalar else array(phase.zs) + vectorized = self.flasher.vectorized + zs = array(phase.zs) if vectorized else list(phase.zs) z_water = zs[water_index] m = 1/(1.0 - z_water) - if scalar: + if vectorized: + zs *= m + else: for i in range(self.N): zs[i] *= m - else: - zs *= m zs[water_index] = 0.0 return zs @@ -2383,15 +2383,15 @@ def ws_no_water(self, phase=None): water_index = self.water_index if water_index is None: return phase.ws() - scalar = self.flasher.scalar - ws = list(phase.ws()) if scalar else array(phase.ws()) + vectorized = self.flasher.vectorized + ws = array(phase.ws()) if vectorized else list(phase.ws()) z_water = ws[water_index] m = 1/(1.0 - z_water) - if scalar: + if vectorized: + ws *= m + else: for i in range(self.N): ws[i] *= m - else: - ws *= m ws[water_index] = 0.0 return ws @@ -2439,10 +2439,10 @@ def Ks(self, phase, ref_phase=None): ref_phase = self.gas ref_zs = ref_phase.zs zs = phase.zs - if self.flasher.scalar: - Ks = [g/l for l, g in zip(ref_zs, zs)] - else: + if self.flasher.vectorized: Ks = zs/ref_zs + else: + Ks = [g/l for l, g in zip(ref_zs, zs)] return Ks def Hc(self, phase=None): @@ -2779,7 +2779,7 @@ def Psats(self): pass T = self.T self._Psats = [o.T_dependent_property(T) for o in self.VaporPressures] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Psats = array(self._Psats) return self._Psats @@ -2807,7 +2807,7 @@ def Psubs(self): pass T = self.T self._Psubs = [o.T_dependent_property(T) for o in self.SublimationPressures] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Psubs = array(self._Psubs) return self._Psubs @@ -2835,7 +2835,7 @@ def Hsubs(self): pass T = self.T self._Hsubs = [o.T_dependent_property(T) for o in self.EnthalpySublimations] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Hsubs = array(self._Hsubs) return self._Hsubs @@ -2863,7 +2863,7 @@ def Hvaps(self): pass T = self.T self._Hvaps = [o.T_dependent_property(T) for o in self.EnthalpyVaporizations] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Hvaps = array(self._Hvaps) return self._Hvaps @@ -2887,7 +2887,7 @@ def sigmas(self): pass T = self.T self._sigmas = [o.T_dependent_property(T) for o in self.SurfaceTensions] - if not self.flasher.scalar: + if self.flasher.vectorized: self._sigmas = array(self._sigmas) return self._sigmas @@ -2911,7 +2911,7 @@ def Cpgs(self): pass T = self.T self._Cpgs = [o.T_dependent_property(T) for o in self.HeatCapacityGases] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Cpgs = array(self._Cpgs) return self._Cpgs @@ -2939,7 +2939,7 @@ def Cpls(self): pass T = self.T self._Cpls = [o.T_dependent_property(T) for o in self.HeatCapacityLiquids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Cpls = array(self._Cpls) return self._Cpls @@ -2962,7 +2962,7 @@ def Cpss(self): pass T = self.T self._Cpss = [o.T_dependent_property(T) for o in self.HeatCapacitySolids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Cpss = array(self._Cpss) return self._Cpss @@ -2990,7 +2990,7 @@ def kls(self): pass T = self.T self._kls = [o.T_dependent_property(T) for o in self.ThermalConductivityLiquids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._kls = array(self._kls) return self._kls @@ -3042,7 +3042,7 @@ def kgs(self): pass T = self.T self._kgs = [o.T_dependent_property(T) for o in self.ThermalConductivityGases] - if not self.flasher.scalar: + if self.flasher.vectorized: self._kgs = array(self._kgs) return self._kgs @@ -3068,7 +3068,7 @@ def muls(self): pass T = self.T self._muls = [o.T_dependent_property(T) for o in self.ViscosityLiquids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._muls = array(self._muls) return self._muls @@ -3094,7 +3094,7 @@ def mugs(self): pass T = self.T self._mugs = [o.T_dependent_property(T) for o in self.ViscosityGases] - if not self.flasher.scalar: + if self.flasher.vectorized: self._mugs = array(self._mugs) return self._mugs @@ -3120,7 +3120,7 @@ def Vls(self): pass T = self.T self._Vls = [o.T_dependent_property(T) for o in self.VolumeLiquids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Vls = array(self._Vls) return self._Vls @@ -3144,7 +3144,7 @@ def Vss(self): pass T = self.T self._Vss = [o.T_dependent_property(T) for o in self.VolumeSolids] - if not self.flasher.scalar: + if self.flasher.vectorized: self._Vss = array(self._Vss) return self._Vss diff --git a/thermo/flash/flash_base.py b/thermo/flash/flash_base.py index fad94fb9..ca4c5500 100644 --- a/thermo/flash/flash_base.py +++ b/thermo/flash/flash_base.py @@ -223,7 +223,7 @@ def flash(self, zs=None, T=None, P=None, VF=None, SF=None, V=None, H=None, # print('flashing',T, P, H, S, VF) if zs is None: if self.N == 1: - zs = [1.0] if self.scalar else ones(1) + zs = ones(1) if self.vectorized else [1.0] else: raise ValueError("Composition missing for flash") constants, correlations = self.constants, self.correlations @@ -734,7 +734,7 @@ def grid_flash(self, zs, Ts=None, Ps=None, Vs=None, spec_iters.append(SFs) do_props = props is not None - scalar_props = isinstance(props, str) + vectorized_props = isinstance(props, str) calc_props = [] for n0, spec0 in enumerate(spec_iters[0]): @@ -753,8 +753,8 @@ def grid_flash(self, zs, Ts=None, Ps=None, Vs=None, if store: row_flashes.append(state) if do_props: - if scalar_props: - state_props = state.value(props)if state is not None else None + if vectorized_props: + state_props = state.value(props) if state is not None else None else: state_props = [state.value(s) for s in props] if state is not None else [None for s in props] @@ -817,19 +817,13 @@ def _finish_initialization_base(self): self.unique_phases += solids self.unique_phase_count = (1 if gas is not None else 0) + self.unique_liquid_count + len(solids) - - - - - self.T_MIN_FLASH = max(p.T_MIN_FLASH for p in self.phases) - self.T_MAX_FLASH = min(p.T_MAX_FLASH for p in self.phases) - self.T_MIN_FLASH_ANY = min(p.T_MIN_FLASH for p in self.phases) - self.T_MAX_FLASH_ANY = max(p.T_MAX_FLASH for p in self.phases) - scalar = True - scalar_statuses = {i.scalar for i in self.phases} - if len(scalar_statuses) > 1: + self.T_MIN_FLASH = self.T_MIN_FLASH_ANY = max(p.T_MIN_FLASH for p in self.phases) + self.T_MAX_FLASH = self.T_MAX_FLASH_ANY = min(p.T_MAX_FLASH for p in self.phases) + vectorized = False + vectorized_statuses = {i.vectorized for i in self.phases} + if len(vectorized_statuses) > 1: raise ValueError("Can only perform flashes with all phases in a numpy basis or all phases in a pure Python basis") - self.scalar = scalar_statuses.pop() + self.vectorized = vectorized_statuses.pop() self.supports_lnphis_args = all(p.supports_lnphis_args for p in self.phases) diff --git a/thermo/phases/ceos.py b/thermo/phases/ceos.py index aa6099ab..31071730 100644 --- a/thermo/phases/ceos.py +++ b/thermo/phases/ceos.py @@ -80,9 +80,9 @@ class CEOSPhase(IdealGasDeparturePhase): ''' - __slots__ = ('eos_class', 'eos_kwargs', 'scalar', 'HeatCapacityGases', 'N', + __slots__ = ('eos_class', 'eos_kwargs', 'vectorized', 'HeatCapacityGases', 'N', 'Hfs', 'Gfs', 'Sfs', 'Cpgs_poly_fit', '_Cpgs_data', 'composition_independent', - 'eos_mix', 'T', 'P', 'zs', '_model_hash_ignore_phase', '_model_hash') + 'eos_mix', 'T', 'P' 'zs', '_model_hash_ignore_phase', '_model_hash') ideal_gas_basis = True pure_references = ('HeatCapacityGases',) @@ -127,14 +127,13 @@ def __repr__(self): return base def __init__(self, eos_class, eos_kwargs, HeatCapacityGases=None, Hfs=None, - Gfs=None, Sfs=None, - T=None, P=None, zs=None): + Gfs=None, Sfs=None, T=None, P=None, zs=None): self.eos_class = eos_class self.eos_kwargs = eos_kwargs - - self.scalar = scalar = not (any(type(v) is ndarray for v in eos_kwargs.values()) or any(type(v) is ndarray for v in (zs, Hfs, Gfs, Sfs))) - - + self.vectorized = vectorized = ( + any(type(v) is ndarray for v in eos_kwargs.values()) + or any(type(v) is ndarray for v in (zs, Hfs, Gfs, Sfs)) + ) self.HeatCapacityGases = HeatCapacityGases if HeatCapacityGases is not None: self.N = N = len(HeatCapacityGases) @@ -143,31 +142,23 @@ def __init__(self, eos_class, eos_kwargs, HeatCapacityGases=None, Hfs=None, raise ValueError("A HeatCapacityGas object is required") elif 'Tcs' in eos_kwargs: self.N = N = len(eos_kwargs['Tcs']) - self.Hfs = Hfs self.Gfs = Gfs self.Sfs = Sfs self.Cpgs_poly_fit, self._Cpgs_data = self._setup_Cpigs(HeatCapacityGases) - self.composition_independent = ideal_gas = eos_class is IGMIX - if ideal_gas: - self.force_phase = 'g' - - - if T is not None and P is not None and zs is not None: - self.T = T - self.P = P - self.zs = zs - self.eos_mix = eos_mix = self.eos_class(T=T, P=P, zs=zs, **self.eos_kwargs) - else: - if scalar: - zs = [1.0/N]*N - else: - v = 1.0/N + self.composition_independent = eos_class is IGMIX + if T is None: T = 298.15 + if P is None: P = 101325.0 + if zs is None: + if vectorized: + v = 1.0 / N zs = full(N, v) - self.eos_mix = eos_mix = self.eos_class(T=298.15, P=101325.0, zs=zs, **self.eos_kwargs) - self.T = 298.15 - self.P = 101325.0 - self.zs = zs + else: + zs = [1.0 / N] * N + self.T = T + self.P = P + self.zs = zs + self.eos_mix = self.eos_class(T=T, P=P, zs=zs, **self.eos_kwargs) def to_TP_zs(self, T, P, zs, other_eos=None): r'''Method to create a new Phase object with the same constants as the @@ -215,10 +206,7 @@ def to_TP_zs(self, T, P, zs, other_eos=None): ''' # Why so slow new = self.__class__.__new__(self.__class__) - new.T = T - new.P = P - new.zs = zs - new.scalar = self.scalar + new.vectorized = self.vectorized if other_eos is not None: other_eos.solve_missing_volumes() new.eos_mix = other_eos @@ -230,20 +218,20 @@ def to_TP_zs(self, T, P, zs, other_eos=None): # 1 hour on this because the heat capacity calculation was wrong except AttributeError: new.eos_mix = self.eos_class(T=T, P=P, zs=zs, **self.eos_kwargs) - new.eos_class = self.eos_class new.eos_kwargs = self.eos_kwargs - new.HeatCapacityGases = self.HeatCapacityGases new._Cpgs_data = self._Cpgs_data new.Cpgs_poly_fit = self.Cpgs_poly_fit new.composition_independent = self.composition_independent - if new.composition_independent: - new.force_phase = 'g' - new.Hfs = self.Hfs new.Gfs = self.Gfs new.Sfs = self.Sfs + new.T = T + new.P = P + new.zs = zs + if new.composition_independent: + new.force_phase = 'g' try: new.N = self.N @@ -255,9 +243,8 @@ def to_TP_zs(self, T, P, zs, other_eos=None): def to(self, zs, T=None, P=None, V=None): new = self.__class__.__new__(self.__class__) - new.zs = zs # temporary TODO remove this statement - if not self.scalar and type(zs) is not ndarray: + if self.vectorized and type(zs) is not ndarray: zs = array(zs) if T is not None: @@ -281,17 +268,15 @@ def to(self, zs, T=None, P=None, V=None): T = new.eos_mix.T else: raise ValueError("Two of T, P, or V are needed") - new.P, new.T = P, T - new.eos_class, new.eos_kwargs = self.eos_class, self.eos_kwargs - new.HeatCapacityGases, new._Cpgs_data, new.Cpgs_poly_fit = self.HeatCapacityGases, self._Cpgs_data, self.Cpgs_poly_fit - new.composition_independent, new.scalar = self.composition_independent, self.scalar + new.composition_independent, new.vectorized = self.composition_independent, self.vectorized + new.Hfs, new.Gfs, new.Sfs = self.Hfs, self.Gfs, self.Sfs + new.T = T + new.P = P + new.zs = zs if new.composition_independent: new.force_phase = 'g' - - new.Hfs, new.Gfs, new.Sfs = self.Hfs, self.Gfs, self.Sfs - try: new.N = self.N except: @@ -318,10 +303,10 @@ def lnphis_args(self, most_stable=False): # Could save time by allowing T, P as an argument, and getting a new eos_mix at that N = self.N eos_mix = self.eos_mix - if self.scalar: - a_alpha_j_rows, vec0, lnphis = [0.0]*N, [0.0]*N, [0.0]*N - else: + if self.vectorized: a_alpha_j_rows, vec0, lnphis = zeros(N), zeros(N), zeros(N) + else: + a_alpha_j_rows, vec0, lnphis = [0.0]*N, [0.0]*N, [0.0]*N l, g = (self.is_liquid, self.is_gas) #if not most_stable else (True, True) if eos_mix.translated: return (self.eos_class.model_id, self.T, self.P, self.N, eos_mix.one_minus_kijs, l, g, @@ -339,10 +324,10 @@ def lnphis_at_zs(self, zs, most_stable=False): def fugacities_at_zs(self, zs, most_stable=False): P = self.P lnphis = lnphis_direct(zs, *self.lnphis_args(most_stable)) - if self.scalar: - return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(len(zs))] - else: + if self.vectorized: return trunc_exp_numpy(lnphis)*P*zs + else: + return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(len(zs))] def T_max_at_V(self, V): T_max = self.eos_mix.T_max_at_V(V) @@ -395,7 +380,6 @@ def k(self): return k def _set_mechanical_critical_point(self): - zs = self.zs new = self.eos_mix.to_mechanical_critical_point() self._mechanical_critical_T = new.T self._mechanical_critical_P = new.P @@ -577,10 +561,10 @@ def fugacities_lowest_Gibbs(self): P = self.P zs = self.zs lnphis = self.lnphis_lowest_Gibbs() - if self.scalar: - return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(len(zs))] - else: + if self.vectorized: return trunc_exp_numpy(lnphis)*P*zs + else: + return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(len(zs))] def V_iter(self, force=False): # Can be some severe issues in the very low pressure/temperature range @@ -608,6 +592,7 @@ def V_iter(self, force=False): class CEOSGas(CEOSPhase): is_gas = True is_liquid = False + __slots__ = () @property def phase(self): @@ -659,7 +644,7 @@ def phi_pures(self): return self._phi_pures except: pass - phis_pure = [0.0]*self.N if self.scalar else zeros(self.N) + phis_pure = zeros(self.N) if self.vectorized else [0.0]*self.N for i, o in enumerate(self.eos_mix.pures()): try: phis_pure[i] = o.phi_g @@ -894,6 +879,7 @@ def dS_dep_dzs(self): class CEOSLiquid(CEOSPhase): is_gas = False is_liquid = True + __slots__ = () @property def phase(self): @@ -945,7 +931,7 @@ def phi_pures(self): return self._phi_pures except: pass - phis_pure = [0.0]*self.N if self.scalar else zeros(self.N) + phis_pure = zeros(self.N) if self.vectorized else [0.0]*self.N for i, o in enumerate(self.eos_mix.pures()): try: phis_pure[i] = (o.phi_l) diff --git a/thermo/phases/gibbs_excess.py b/thermo/phases/gibbs_excess.py index 03e42298..194410ee 100644 --- a/thermo/phases/gibbs_excess.py +++ b/thermo/phases/gibbs_excess.py @@ -364,10 +364,10 @@ def __init__(self, VaporPressures, VolumeLiquids=None, self._Psats_data = Psats_data - if self.scalar: - zero_coeffs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: zero_coeffs = zeros((N, N)) + else: + zero_coeffs = [[0.0]*N for _ in range(N)] self.HeatCapacityGases = HeatCapacityGases self.Cpgs_poly_fit, self._Cpgs_data = self._setup_Cpigs(HeatCapacityGases) @@ -503,20 +503,20 @@ def __init__(self, VaporPressures, VolumeLiquids=None, if input_count_henry > 1: raise ValueError("Input only one of henry_abcdef, or (henry_as...henry_fs)") if henry_abcdef is not None: - if self.scalar: - self.henry_as = [[i[0] for i in l] for l in henry_abcdef] - self.henry_bs = [[i[1] for i in l] for l in henry_abcdef] - self.henry_cs = [[i[2] for i in l] for l in henry_abcdef] - self.henry_ds = [[i[3] for i in l] for l in henry_abcdef] - self.henry_es = [[i[4] for i in l] for l in henry_abcdef] - self.henry_fs = [[i[5] for i in l] for l in henry_abcdef] - else: + if self.vectorized: self.henry_as = array(henry_abcdef[:,:,0], order='C', copy=True) self.henry_bs = array(henry_abcdef[:,:,1], order='C', copy=True) self.henry_cs = array(henry_abcdef[:,:,2], order='C', copy=True) self.henry_ds = array(henry_abcdef[:,:,3], order='C', copy=True) self.henry_es = array(henry_abcdef[:,:,4], order='C', copy=True) self.henry_fs = array(henry_abcdef[:,:,5], order='C', copy=True) + else: + self.henry_as = [[i[0] for i in l] for l in henry_abcdef] + self.henry_bs = [[i[1] for i in l] for l in henry_abcdef] + self.henry_cs = [[i[2] for i in l] for l in henry_abcdef] + self.henry_ds = [[i[3] for i in l] for l in henry_abcdef] + self.henry_es = [[i[4] for i in l] for l in henry_abcdef] + self.henry_fs = [[i[5] for i in l] for l in henry_abcdef] else: if henry_abcdef is None: henry_abcdef = multiple_henry_inputs @@ -720,7 +720,7 @@ def lnphis_args(self): Poyntings = self.Poyntings() phis_sat = self.phis_sat() activity_args = self.GibbsExcessModel.gammas_args() - lnphis = [0.0]*self.N if self.scalar else zeros(self.N) + lnphis = zeros(self.N) if self.vectorized else [0.0]*self.N self._lnphis_args = (self.model_id, self.T, self.P, self.N, lnPsats, Poyntings, phis_sat) + activity_args +(lnphis,) return self._lnphis_args @@ -747,10 +747,10 @@ def lnHenry_matrix(self): except: pass N = self.N - if self.scalar: - lnHenry_matrix = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: lnHenry_matrix = zeros((N, N)) + else: + lnHenry_matrix = [[0.0]*N for _ in range(N)] lnHenry_matrix = ln_henries(self.T, N, self.henry_as, self.henry_bs, self.henry_cs, self.henry_ds, self.henry_es, self.henry_fs, lnHenry_matrix) self._lnHenry_matrix = lnHenry_matrix @@ -775,10 +775,10 @@ def dlnHenry_matrix_dT(self): except: pass N = self.N - if self.scalar: - dlnHenry_matrix_dT = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: dlnHenry_matrix_dT = zeros((N, N)) + else: + dlnHenry_matrix_dT = [[0.0]*N for _ in range(N)] dlnHenry_matrix_dT = dln_henries_dT(self.T, N, self.henry_bs, self.henry_cs, self.henry_ds, self.henry_es, self.henry_fs, dlnHenry_matrix_dT) self._dlnHenry_matrix_dT = dlnHenry_matrix_dT return dlnHenry_matrix_dT @@ -802,39 +802,39 @@ def d2lnHenry_matrix_dT2(self): except: pass N = self.N - if self.scalar: - d2lnHenry_matrix_dT2 = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d2lnHenry_matrix_dT2 = zeros((N, N)) + else: + d2lnHenry_matrix_dT2 = [[0.0]*N for _ in range(N)] d2lnHenry_matrix_dT2 = d2ln_henries_dT2(self.T, N, self.henry_bs, self.henry_cs, self.henry_es, self.henry_fs, d2lnHenry_matrix_dT2) self._d2lnHenry_matrix_dT2 = d2lnHenry_matrix_dT2 return d2lnHenry_matrix_dT2 def Henry_constants(self): - zs, scalar, N, henry_components, henry_mode = self.zs, self.scalar, self.N, self.henry_components, self.henry_mode + zs, vectorized, N, henry_components, henry_mode = self.zs, self.vectorized, self.N, self.henry_components, self.henry_mode solvents_with_parameters = henry_mode == 'solvents_with_parameters' lnHenry_matrix = self.lnHenry_matrix() - Hs = [0.0]*N if scalar else zeros(N) + Hs = zeros(N) if vectorized else [0.0]*N Henry_constants(lnHenry_matrix, zs, henry_components, solvents_with_parameters, Hs) return Hs dHenry_constants_dT def dHenry_constants_dT(self): - zs, scalar, N, henry_components, henry_mode = self.zs, self.scalar, self.N, self.henry_components, self.henry_mode + zs, vectorized, N, henry_components, henry_mode = self.zs, self.vectorized, self.N, self.henry_components, self.henry_mode solvents_with_parameters = henry_mode == 'solvents_with_parameters' lnHenry_matrix = self.lnHenry_matrix() dlnHenry_matrix_dT = self.dlnHenry_matrix_dT() - dHs = [0.0]*N if scalar else zeros(N) + dHs = zeros(N) if vectorized else [0.0]*N dHenry_constants_dT(lnHenry_matrix, dlnHenry_matrix_dT, zs, henry_components, solvents_with_parameters, dHs) return dHs def d2Henry_constants_dT2(self): - zs, scalar, N, henry_components, henry_mode = self.zs, self.scalar, self.N, self.henry_components, self.henry_mode + zs, vectorized, N, henry_components, henry_mode = self.zs, self.vectorized, self.N, self.henry_components, self.henry_mode solvents_with_parameters = henry_mode == 'solvents_with_parameters' lnHenry_matrix = self.lnHenry_matrix() dlnHenry_matrix_dT = self.dlnHenry_matrix_dT() d2lnHenry_matrix_dT2 = self.d2lnHenry_matrix_dT2() - d2Hs = [0.0]*N if scalar else zeros(N) + d2Hs = zeros(N) if vectorized else [0.0]*N d2Henry_constants_dT2(lnHenry_matrix, dlnHenry_matrix_dT, d2lnHenry_matrix_dT2, zs, henry_components, solvents_with_parameters, d2Hs) return d2Hs diff --git a/thermo/phases/ideal_gas.py b/thermo/phases/ideal_gas.py index 128c5139..da6167d5 100644 --- a/thermo/phases/ideal_gas.py +++ b/thermo/phases/ideal_gas.py @@ -92,7 +92,7 @@ def __init__(self, HeatCapacityGases=None, Hfs=None, Gfs=None, T=None, P=None, z self.HeatCapacityGases = HeatCapacityGases self.Hfs = Hfs self.Gfs = Gfs - self.scalar = scalar = not any(type(v) is ndarray for v in (zs, Hfs, Gfs)) + self.vectorized = vectorized = any(type(v) is ndarray for v in (zs, Hfs, Gfs)) if Hfs is not None and Gfs is not None and None not in Hfs and None not in Gfs: T_ref_inv = 1.0/298.15 @@ -102,12 +102,12 @@ def __init__(self, HeatCapacityGases=None, Hfs=None, Gfs=None, T=None, P=None, z if zs is not None: self.N = N = len(zs) - self.zeros1d = [0.0]*N if scalar else zeros(N) - self.ones1d = [1.0]*N if scalar else ones(N) + self.zeros1d = zeros(N) if vectorized else [0.0]*N + self.ones1d = ones(N) if vectorized else [1.0]*N elif HeatCapacityGases is not None: self.N = N = len(HeatCapacityGases) - self.zeros1d = [0.0]*N if scalar else zeros(N) - self.ones1d = [1.0]*N if scalar else ones(N) + self.zeros1d = zeros(N) if vectorized else [0.0]*N + self.ones1d = ones(N) if vectorized else [1.0]*N if zs is not None: self.zs = zs if T is not None: @@ -267,7 +267,7 @@ def to_TP_zs(self, T, P, zs): new.N = self.N new.zeros1d = self.zeros1d new.ones1d = self.ones1d - new.scalar = self.scalar + new.vectorized = self.vectorized new.HeatCapacityGases = self.HeatCapacityGases new.Hfs = self.Hfs @@ -298,7 +298,7 @@ def to(self, zs, T=None, P=None, V=None): new.Hfs = self.Hfs new.Gfs = self.Gfs new.Sfs = self.Sfs - new.scalar = self.scalar + new.vectorized = self.vectorized return new diff --git a/thermo/phases/phase.py b/thermo/phases/phase.py index b03b748b..af506385 100644 --- a/thermo/phases/phase.py +++ b/thermo/phases/phase.py @@ -142,7 +142,7 @@ class Phase: Psats_poly_fit = False Cpgs_poly_fit = False composition_independent = False - scalar = True + vectorized = False supports_lnphis_args = False @@ -210,7 +210,7 @@ def as_json(self): >>> assert phase == new_phase ''' d = object_data(self) - if not self.scalar: + if self.vectorized: d = arrays_to_lists(d) for obj_name in self.obj_references: o = d[obj_name] @@ -230,8 +230,6 @@ def as_json(self): d['json_version'] = 1 return d - - @classmethod def from_json(cls, json_repr): r'''Method to create a phase from a JSON @@ -279,9 +277,8 @@ def from_json(cls, json_repr): for ref_name, ref_lookup in zip(new.pointer_references, new.pointer_reference_dicts): d[ref_name] = ref_lookup[d[ref_name]] - for k, v in d.items(): + for k, v in d.items(): setattr(new, k, v) - # new.__dict__ = d return new def __hash__(self): @@ -681,11 +678,11 @@ def S_phi_consistency(self): lnphis = self.lnphis() dlnphis_dT = self.dlnphis_dT() T, zs = self.T, self.zs - if self.scalar: + if self.vectorized: + S0 -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) + else: for i in range(self.N): S0 -= zs[i]*(R*lnphis[i] + R*T*dlnphis_dT[i]) - else: - S0 -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) self._S_phi_consistency = abs(1.0 - S0/self.S()) return self._S_phi_consistency @@ -733,10 +730,10 @@ def G_dep_phi_consistency(self): zs, T = self.zs, self.T G_dep_RT = 0.0 lnphis = self.lnphis() - if self.scalar: - G_dep_RT = sum(zs[i]*lnphis[i] for i in range(self.N)) - else: + if self.vectorized: G_dep_RT = float(dot(zs, lnphis)) + else: + G_dep_RT = sum(zs[i]*lnphis[i] for i in range(self.N)) G_dep = G_dep_RT*R*T self._G_dep_phi_consistency = abs(1.0 - G_dep/self.G_dep()) return self._G_dep_phi_consistency @@ -763,10 +760,10 @@ def H_dep_phi_consistency(self): H_dep_RT2 = 0.0 dlnphis_dTs = self.dlnphis_dT() zs, T = self.zs, self.T - if self.scalar: - H_dep_RT2 = sum([zs[i]*dlnphis_dTs[i] for i in range(self.N)]) - else: + if self.vectorized: H_dep_RT2 = float(dot(zs, dlnphis_dTs)) + else: + H_dep_RT2 = sum([zs[i]*dlnphis_dTs[i] for i in range(self.N)]) H_dep_recalc = -H_dep_RT2*R*T*T H_dep = self.H_dep() self._H_dep_phi_consistency = abs(1.0 - H_dep/H_dep_recalc) @@ -796,11 +793,11 @@ def S_dep_phi_consistency(self): dlnphis_dT = self.dlnphis_dT() T, zs = self.T, self.zs S_dep = 0.0 - if self.scalar: + if self.vectorized: + S_dep -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) + else: for i in range(self.N): S_dep -= zs[i]*(R*lnphis[i] + R*T*dlnphis_dT[i]) - else: - S_dep -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) self._S_dep_phi_consistency = abs(1.0 - S_dep/self.S_dep()) return self._S_dep_phi_consistency @@ -826,10 +823,10 @@ def V_phi_consistency(self): pass zs, P = self.zs, self.P dlnphis_dP = self.dlnphis_dP() - if self.scalar: - lhs = sum(zs[i]*dlnphis_dP[i] for i in range(self.N)) - else: + if self.vectorized: lhs = float(dot(zs, dlnphis_dP)) + else: + lhs = sum(zs[i]*dlnphis_dP[i] for i in range(self.N)) Z_calc = lhs*P + 1.0 V_calc = Z_calc*self.R*self.T/P V = self.V() @@ -854,11 +851,11 @@ def H_from_phi(self): H0 = self.H_ideal_gas() dlnphis_dT = self.dlnphis_dT() T, zs = self.T, self.zs - if self.scalar: + if self.vectorized: + H0 -= R*T*T*float(dot(zs, dlnphis_dT)) + else: for i in range(self.N): H0 -= R*T*T*zs[i]*dlnphis_dT[i] - else: - H0 -= R*T*T*float(dot(zs, dlnphis_dT)) return H0 def S_from_phi(self): @@ -880,11 +877,11 @@ def S_from_phi(self): lnphis = self.lnphis() dlnphis_dT = self.dlnphis_dT() T, zs = self.T, self.zs - if self.scalar: + if self.vectorized: + S0 -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) + else: for i in range(self.N): S0 -= zs[i]*(R*lnphis[i] + R*T*dlnphis_dT[i]) - else: - S0 -= float((zs*(R*lnphis + R*T*dlnphis_dT)).sum()) return S0 def V_from_phi(self): @@ -903,10 +900,10 @@ def V_from_phi(self): ''' zs, P = self.zs, self.P dlnphis_dP = self.dlnphis_dP() - if self.scalar: - obj = sum(zs[i]*dlnphis_dP[i] for i in range(self.N)) - else: + if self.vectorized: obj = float(dot(zs, dlnphis_dP)) + else: + obj = sum(zs[i]*dlnphis_dP[i] for i in range(self.N)) Z = P*obj + 1.0 return Z*self.R*self.T/P @@ -930,11 +927,11 @@ def G_min_criteria(self): zs = self.zs log_zs = self.log_zs() G_crit = 0.0 - if self.scalar: + if self.vectorized: + G_crit += float(dot(zs, log_zs)) + else: for i in range(self.N): G_crit += zs[i]*log_zs[i] - else: - G_crit += float(dot(zs, log_zs)) G_crit = G_crit*R*self.T + self.G_dep() return G_crit @@ -978,10 +975,10 @@ def fugacities_at_zs(self, zs, most_stable=False): # P = self.P # fugacities_lowest_Gibbs # lnphis = self.lnphis_at_zs(zs, most_stable) - # if self.scalar: - # return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(self.N)] - # else: + # if self.vectorized: # return P*zs*trunc_exp_numpy(lnphis) + # else: + # return [P*zs[i]*trunc_exp(lnphis[i]) for i in range(self.N)] def lnphi(self): r'''Method to calculate and return the log of fugacity coefficient of @@ -1067,10 +1064,10 @@ def fugacities(self): P = self.P zs = self.zs lnphis = self.lnphis() - if self.scalar: - self._fugacities = [P*zs[i]*trunc_exp(lnphis[i]) for i in range(self.N)] - else: + if self.vectorized: self._fugacities = P*zs*trunc_exp_numpy(lnphis) + else: + self._fugacities = [P*zs[i]*trunc_exp(lnphis[i]) for i in range(self.N)] return self._fugacities def lnfugacities(self): @@ -1093,10 +1090,10 @@ def lnfugacities(self): lnphis = self.lnphis() logP = log(P) log_zs = self.log_zs() - if self.scalar: - lnfugacities = [logP + log_zs[i] + lnphis[i] for i in range(self.N)] - else: + if self.vectorized: lnfugacities = logP + log_zs + lnphis + else: + lnfugacities = [logP + log_zs[i] + lnphis[i] for i in range(self.N)] self._lnfugacities = lnfugacities return lnfugacities @@ -1110,10 +1107,10 @@ def lnphis_lowest_Gibbs(self): P = self.P zs = self.zs fugacities_lowest_Gibbs = self.fugacities_lowest_Gibbs() - if self.scalar: - lnphis_lowest_Gibbs = [trunc_log(fi/(zi*P)) for fi, zi in zip(fugacities_lowest_Gibbs, zs)] - else: + if self.vectorized: lnphis_lowest_Gibbs = trunc_log_numpy(fugacities_lowest_Gibbs/(zs*P)) + else: + lnphis_lowest_Gibbs = [trunc_log(fi/(zi*P)) for fi, zi in zip(fugacities_lowest_Gibbs, zs)] self._lnphis_lowest_Gibbs = lnphis_lowest_Gibbs return lnphis_lowest_Gibbs @@ -1140,10 +1137,10 @@ def dfugacities_dT(self): pass dphis_dT = self.dphis_dT() P, zs = self.P, self.zs - if self.scalar: - self._dfugacities_dT = [P*zs[i]*dphis_dT[i] for i in range(self.N)] - else: + if self.vectorized: self._dfugacities_dT = P*zs*dphis_dT + else: + self._dfugacities_dT = [P*zs[i]*dphis_dT[i] for i in range(self.N)] return self._dfugacities_dT def lnphis_G_min(self): @@ -1175,10 +1172,10 @@ def phis(self): return self._phis except: pass - if self.scalar: - self._phis = [trunc_exp(i) for i in self.lnphis()] - else: + if self.vectorized: self._phis = trunc_exp_numpy(self.lnphis()) + else: + self._phis = [trunc_exp(i) for i in self.lnphis()] return self._phis def dphis_dT(self): @@ -1211,10 +1208,10 @@ def dphis_dT(self): phis = self._phis except AttributeError: phis = self.phis() - if self.scalar: - self._dphis_dT = [dlnphis_dT[i]*phis[i] for i in range(self.N)] - else: + if self.vectorized: self._dphis_dT = dlnphis_dT*phis + else: + self._dphis_dT = [dlnphis_dT[i]*phis[i] for i in range(self.N)] return self._dphis_dT def dphis_dP(self): @@ -1247,10 +1244,10 @@ def dphis_dP(self): phis = self._phis except AttributeError: phis = self.phis() - if self.scalar: - self._dphis_dP = [dlnphis_dP[i]*phis[i] for i in range(self.N)] - else: + if self.vectorized: self._dphis_dP = dlnphis_dP*phis + else: + self._dphis_dP = [dlnphis_dP[i]*phis[i] for i in range(self.N)] return self._dphis_dP @@ -1288,7 +1285,7 @@ def dphis_dzs(self): N = self.N self._dphis_dzs = [[dlnphis_dzs[i][j]*phis[i] for j in range(N)] for i in range(N)] - if not self.scalar: + if self.vectorized: self._dphis_dzs = array(self._dphis_dzs) return self._dphis_dzs @@ -1324,10 +1321,10 @@ def dfugacities_dP(self): phis = self.phis() P, zs = self.P, self.zs - if self.scalar: - return [zs[i]*(P*dphis_dP[i] + phis[i]) for i in range(self.N)] - else: + if self.vectorized: return zs*(P*dphis_dP + phis) + else: + return [zs[i]*(P*dphis_dP[i] + phis[i]) for i in range(self.N)] def dfugacities_dns(self): r'''Method to calculate and return the mole number derivative of the @@ -1359,7 +1356,7 @@ def dfugacities_dns(self): phis = self.phis() dlnphis_dns = self.dlnphis_dns() P, zs, N = self.P, self.zs, self.N - matrix = [[0.0]*N for _ in range(N)] if self.scalar else zeros((N, N)) + matrix = zeros((N, N)) if self.vectorized else [[0.0]*N for _ in range(N)] for i in range(N): phi_P = P*phis[i] ziPphi = phi_P*zs[i] @@ -1388,7 +1385,7 @@ def dlnfugacities_dns(self): ''' fugacities = self.fugacities() dlnfugacities_dns = [list(i) for i in self.dfugacities_dns()] - if not self.scalar: + if self.vectorized: dlnfugacities_dns = array(dlnfugacities_dns) fugacities_inv = [1.0/fi for fi in fugacities] cmps = range(self.N) @@ -1444,10 +1441,10 @@ def log_zs(self): return self._log_zs except AttributeError: pass - if self.scalar: - self._log_zs = [trunc_log(zi) for zi in self.zs] - else: + if self.vectorized: self._log_zs = trunc_log_numpy(self.zs) + else: + self._log_zs = [trunc_log(zi) for zi in self.zs] # except ValueError: # self._log_zs = _log_zs = [] # for zi in self.zs: @@ -1547,7 +1544,7 @@ def dH_dns(self): Notes ----- ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dH_dzs(), self.zs, out) def dS_dns(self): @@ -1566,7 +1563,7 @@ def dS_dns(self): Notes ----- ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dS_dzs(), self.zs, out) def dG_dT(self): @@ -2108,11 +2105,11 @@ def H_reactive(self): except AttributeError: pass H = self.H() - if self.scalar: + if self.vectorized: + H += float(dot(self.zs, self.Hfs)) + else: for zi, Hf in zip(self.zs, self.Hfs): H += zi*Hf - else: - H += float(dot(self.zs, self.Hfs)) self._H_reactive = H return H @@ -2136,11 +2133,11 @@ def S_reactive(self): except: pass S = self.S() - if self.scalar: + if self.vectorized: + S += float(dot(self.zs, self.Sfs)) + else: for zi, Sf in zip(self.zs, self.Sfs): S += zi*Sf - else: - S += float(dot(self.zs, self.Sfs)) self._S_reactive = S return S @@ -2219,11 +2216,11 @@ def H_formation_ideal_gas(self): except AttributeError: pass Hf_ideal_gas = 0.0 - if self.scalar: + if self.vectorized: + Hf_ideal_gas = float(self.zs, self.Hfs) + else: for zi, Hf in zip(self.zs, self.Hfs): Hf_ideal_gas += zi*Hf - else: - Hf_ideal_gas = float(self.zs, self.Hfs) self._H_formation_ideal_gas = Hf_ideal_gas return Hf_ideal_gas @@ -2248,11 +2245,11 @@ def S_formation_ideal_gas(self): except: pass Sf_ideal_gas = 0.0 - if self.scalar: + if self.vectorized: + Sf_ideal_gas = float(dot(self.zs, self.Sfs)) + else: for zi, Sf in zip(self.zs, self.Sfs): Sf_ideal_gas += zi*Sf - else: - Sf_ideal_gas = float(dot(self.zs, self.Sfs)) self._S_formation_ideal_gas = Sf_ideal_gas return Sf_ideal_gas @@ -2497,11 +2494,11 @@ def chemical_potential(self): dS_dzs = self.dS_dzs() dH_dzs = self.dH_dzs() T, Hfs, Sfs = self.T, self.Hfs, self.Sfs - if self.scalar: - dG_reactive_dzs = [Hfs[i] - T*(Sfs[i] + dS_dzs[i]) + dH_dzs[i] for i in range(self.N)] - else: + if self.vectorized: dG_reactive_dzs = Hfs - T*(Sfs + dS_dzs) + dH_dzs - dG_reactive_dns = [0.0]*self.N if self.scalar else zeros(self.N) + else: + dG_reactive_dzs = [Hfs[i] - T*(Sfs[i] + dS_dzs[i]) + dH_dzs[i] for i in range(self.N)] + dG_reactive_dns = zeros(self.N) if self.vectorized else [0.0]*self.N dG_reactive_dns = dxs_to_dns(dG_reactive_dzs, self.zs, dG_reactive_dns) chemical_potentials = dns_to_dn_partials(dG_reactive_dns, self.G_reactive()) self._chemical_potentials = chemical_potentials @@ -2532,9 +2529,9 @@ def activities(self): # CORRECT DO NOT CHANGE fugacities = self.fugacities() fugacities_std = self.fugacities_std() # TODO implement fugacities_std - if self.scalar: - return [fugacities[i]/fugacities_std[i] for i in range(self.N)] - return fugacities/fugacities_std + if self.vectorized: + return fugacities/fugacities_std + return [fugacities[i]/fugacities_std[i] for i in range(self.N)] def gammas(self): r'''Method to calculate and return the activity coefficients of the @@ -2566,10 +2563,10 @@ def gammas(self): # the most generally used one for EOSs; and activity methods # override this phis = self.phis() - T, P, N, scalar = self.T, self.P, self.N, self.scalar - self._gammas = gammas = [0.0]*N if scalar else zeros(N) + T, P, N, vectorized = self.T, self.P, self.N, self.vectorized + self._gammas = gammas = zeros(N) if vectorized else [0.0]*N for i in range(N): - comp = [0.0]*N if scalar else zeros(N) + comp = zeros(N) if vectorized else [0.0]*N comp[i] = 1.0 phi = self.to_TP_zs(T=T, P=P, zs=comp).phis()[i] gammas[i] = phis[i]/phi @@ -2599,12 +2596,12 @@ def gammas_infinite_dilution(self): zs_base = self.zs x_infinite_dilution = self._x_infinite_dilution # x_infinite_dilution = 1e-7 - if self.scalar: - gammas_inf = [0.0]*N - copy_fun = list - else: + if self.vectorized: gammas_inf = zeros(N) copy_fun = array + else: + gammas_inf = [0.0]*N + copy_fun = list phis = self.phis() for i in range(N): zs = copy_fun(zs_base) @@ -3039,9 +3036,9 @@ def dZ_dzs(self): ----- ''' factor = self.P/(self.T*self.R) - if self.scalar: - return [dV*factor for dV in self.dV_dzs()] - return factor*self.dV_dzs() + if self.vectorized: + return factor*self.dV_dzs() + return [dV*factor for dV in self.dV_dzs()] def dZ_dns(self): r'''Method to calculate and return the mole number derivatives of the @@ -3059,7 +3056,7 @@ def dZ_dns(self): Notes ----- ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dZ_dzs(), self.zs, out) def dV_dzs(self): @@ -3094,7 +3091,7 @@ def dV_dns(self): Notes ----- ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dxs_to_dns(self.dV_dzs(), self.zs, out) def dnV_dns(self): @@ -3113,7 +3110,7 @@ def dnV_dns(self): Notes ----- ''' - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N return dns_to_dn_partials(self.dV_dns(), self.V(), out) # Derived properties @@ -3661,7 +3658,7 @@ def _setup_Cpigs(self, HeatCapacityGases): def _Cp_pure_fast(self, Cps_data): T, N = self.T, self.N - Cps = [0.0]*N if self.scalar else zeros(N) + Cps = zeros(N) if self.vectorized else [0.0]*N Tmins, Tmaxs, coeffs = Cps_data[0], Cps_data[3], Cps_data[12] Tmin_slopes = Cps_data[1] Tmin_values = Cps_data[2] @@ -3682,7 +3679,7 @@ def _Cp_pure_fast(self, Cps_data): def _dCp_dT_pure_fast(self, Cps_data): T, N = self.T, self.N - dCps = [0.0]*N if self.scalar else zeros(N) + dCps = zeros(N) if self.vectorized else [0.0]*N Tmins, Tmaxs, coeffs = Cps_data[0], Cps_data[3], Cps_data[12] Tmin_slopes = Cps_data[1] Tmax_slopes = Cps_data[4] @@ -3701,7 +3698,7 @@ def _dCp_dT_pure_fast(self, Cps_data): def _Cp_integrals_pure_fast(self, Cps_data): T, N = self.T, self.N - Cp_integrals_pure = [0.0]*N if self.scalar else zeros(N) + Cp_integrals_pure = zeros(N) if self.vectorized else [0.0]*N Tmins, Tmaxes, int_coeffs = Cps_data[0], Cps_data[3], Cps_data[13] for i in range(N): # If indeed everything is working here, need to optimize to decide what to store @@ -3757,7 +3754,7 @@ def _Cp_integrals_pure_fast(self, Cps_data): def _Cp_integrals_over_T_pure_fast(self, Cps_data): T, N = self.T, self.N Tmins, Tmaxes, T_int_T_coeffs = Cps_data[0], Cps_data[3], Cps_data[14] - Cp_integrals_over_T_pure = [0.0]*N if self.scalar else zeros(N) + Cp_integrals_over_T_pure = zeros(N) if self.vectorized else [0.0]*N logT = log(T) for i in range(N): Tmin = Tmins[i] @@ -3807,7 +3804,7 @@ def Cpigs_pure(self): T = self.T Cpigs = [i.T_dependent_property(T) for i in self.HeatCapacityGases] - if not self.scalar: + if self.vectorized: Cpigs = array(Cpigs) self._Cpigs = Cpigs return Cpigs @@ -3840,7 +3837,7 @@ def Cpig_integrals_pure(self): T, T_REF_IG, HeatCapacityGases = self.T, self.T_REF_IG, self.HeatCapacityGases Cpig_integrals_pure = [obj.T_dependent_property_integral(T_REF_IG, T) for obj in HeatCapacityGases] - if not self.scalar: + if self.vectorized: Cpig_integrals_pure = array(Cpig_integrals_pure) self._Cpig_integrals_pure = Cpig_integrals_pure return Cpig_integrals_pure @@ -3876,7 +3873,7 @@ def Cpig_integrals_over_T_pure(self): T, T_REF_IG, HeatCapacityGases = self.T, self.T_REF_IG, self.HeatCapacityGases Cpig_integrals_over_T_pure = [obj.T_dependent_property_integral_over_T(T_REF_IG, T) for obj in HeatCapacityGases] - if not self.scalar: + if self.vectorized: Cpig_integrals_over_T_pure = array(Cpig_integrals_over_T_pure) self._Cpig_integrals_over_T_pure = Cpig_integrals_over_T_pure return Cpig_integrals_over_T_pure @@ -3907,7 +3904,7 @@ def dCpigs_dT_pure(self): T = self.T dCpigs_dT = [i.T_dependent_property_derivative(T) for i in self.HeatCapacityGases] - if not self.scalar: + if self.vectorized: dCpigs_dT = array(dCpigs_dT) self._dCpigs_dT = dCpigs_dT return dCpigs_dT @@ -3924,7 +3921,7 @@ def _Cpls_pure(self): T = self.T Cpls = [i.T_dependent_property(T) for i in self.HeatCapacityLiquids] - if not self.scalar: + if self.vectorized: Cpls = array(Cpls) self._Cpls = Cpls return Cpls @@ -3949,7 +3946,7 @@ def _Cpl_integrals_pure(self): T, T_REF_IG, HeatCapacityLiquids = self.T, self.T_REF_IG, self.HeatCapacityLiquids Cpl_integrals_pure = [obj.T_dependent_property_integral(T_REF_IG, T) for obj in HeatCapacityLiquids] - if not self.scalar: + if self.vectorized: Cpl_integrals_pure = array(Cpl_integrals_pure) self._Cpl_integrals_pure = Cpl_integrals_pure return Cpl_integrals_pure @@ -3975,7 +3972,7 @@ def _Cpl_integrals_over_T_pure(self): T, T_REF_IG, HeatCapacityLiquids = self.T, self.T_REF_IG, self.HeatCapacityLiquids Cpl_integrals_over_T_pure = [obj.T_dependent_property_integral_over_T(T_REF_IG, T) for obj in HeatCapacityLiquids] - if not self.scalar: + if self.vectorized: Cpl_integrals_over_T_pure = array(Cpl_integrals_over_T_pure) self._Cpl_integrals_over_T_pure = Cpl_integrals_over_T_pure return Cpl_integrals_over_T_pure @@ -4009,12 +4006,12 @@ def H_ideal_gas(self): return self._H_ideal_gas except AttributeError: pass - if self.scalar: + if self.vectorized: + H = float(dot(self.zs, self.Cpig_integrals_pure())) + else: H = 0.0 for zi, Cp_int in zip(self.zs, self.Cpig_integrals_pure()): H += zi*Cp_int - else: - H = float(dot(self.zs, self.Cpig_integrals_pure())) self._H_ideal_gas = H return H @@ -4041,13 +4038,13 @@ def S_ideal_gas(self): S = 0.0 S -= R*log(P*P_REF_IG_INV) - if self.scalar: + if self.vectorized: + S -= R*float(dot(zs, log_zs)) + S += float(dot(zs, Cpig_integrals_over_T_pure)) + else: S -= R*sum([zs[i]*log_zs[i] for i in cmps]) # ideal composition entropy composition for i in cmps: S += zs[i]*Cpig_integrals_over_T_pure[i] - else: - S -= R*float(dot(zs, log_zs)) - S += float(dot(zs, Cpig_integrals_over_T_pure)) self._S_ideal_gas = S return S @@ -4069,11 +4066,11 @@ def Cp_ideal_gas(self): pass Cpigs_pure = self.Cpigs_pure() Cp, zs = 0.0, self.zs - if self.scalar: + if self.vectorized: + Cp = float(dot(zs, Cpigs_pure)) + else: for i in range(self.N): Cp += zs[i]*Cpigs_pure[i] - else: - Cp = float(dot(zs, Cpigs_pure)) self._Cp_ideal_gas = Cp return Cp @@ -4715,12 +4712,12 @@ def MW(self): except AttributeError: pass zs, MWs = self.zs, self.constants.MWs - if self.scalar: + if self.vectorized: + MW = float(dot(zs, MWs)) + else: MW = 0.0 for i in range(self.N): MW += zs[i]*MWs[i] - else: - MW = float(dot(zs, MWs)) self._MW = MW return MW @@ -5362,15 +5359,15 @@ def ws(self): pass MWs = self.constants.MWs zs, cmps = self.zs, range(self.N) - if self.scalar: + if self.vectorized: + ws = zs*MWs + Mavg = 1.0/ws.sum() + ws *= Mavg + else: ws = [zs[i]*MWs[i] for i in cmps] Mavg = 1.0/sum(ws) for i in cmps: ws[i] *= Mavg - else: - ws = zs*MWs - Mavg = 1.0/ws.sum() - ws *= Mavg self._ws = ws return ws @@ -5657,10 +5654,10 @@ def ns(self): return self._ns except: n = self.result.n*self.beta - if self.scalar: - self._ns = [n*zi for zi in self.zs] - else: + if self.vectorized: self._ns = self.zs*n + else: + self._ns = [n*zi for zi in self.zs] return self._ns except: return None @@ -5687,10 +5684,10 @@ def ms(self): except: pass m = self.result.m*self.beta_mass - if self.scalar: - self._ms = [m*wi for wi in self.ws()] - else: + if self.vectorized: self._ms = m*self.ws() + else: + self._ms = [m*wi for wi in self.ws()] return self._ms except: return None @@ -5721,10 +5718,10 @@ def Qgs(self): V = R*settings.T_gas_ref/settings.P_gas_ref n = self.n Vn = V*n - if self.scalar: - self._Qgs = [zi*Vn for zi in self.zs] - else: + if self.vectorized: self._Qgs = self.zs*Vn + else: + self._Qgs = [zi*Vn for zi in self.zs] return self._Qgs Qgs_calc = Qgs @@ -5779,10 +5776,10 @@ def Qls(self): pass ns = self.ns Vmls = self.result.V_liquids_ref() - if self.scalar: - self._Qls = [ns[i]*Vmls[i] for i in range(self.N)] - else: + if self.vectorized: self._Qls = ns*Vmls + else: + self._Qls = [ns[i]*Vmls[i] for i in range(self.N)] return self._Qls Qls_calc = Qls @@ -5922,10 +5919,10 @@ def concentrations(self): pass rho = self.rho() zs = self.zs - if self.scalar: - self._concentrations = concentrations = [rho*zi for zi in zs] - else: + if self.vectorized: self._concentrations = concentrations = rho*zs + else: + self._concentrations = concentrations = [rho*zi for zi in zs] return concentrations def concentrations_gas(self): @@ -5943,10 +5940,10 @@ def concentrations_gas(self): ''' rho = self.rho_gas() zs = self.zs - if self.scalar: - concentrations = [rho*zi for zi in zs] - else: + if self.vectorized: concentrations = rho*zs + else: + concentrations = [rho*zi for zi in zs] return concentrations def concentrations_gas_normal(self): @@ -5964,10 +5961,10 @@ def concentrations_gas_normal(self): ''' rho = self.rho_gas_normal() zs = self.zs - if self.scalar: - concentrations = [rho*zi for zi in zs] - else: + if self.vectorized: concentrations = rho*zs + else: + concentrations = [rho*zi for zi in zs] return concentrations def concentrations_gas_standard(self): @@ -5985,10 +5982,10 @@ def concentrations_gas_standard(self): ''' rho = self.rho_gas_standard() zs = self.zs - if self.scalar: - concentrations = [rho*zi for zi in zs] - else: + if self.vectorized: concentrations = rho*zs + else: + concentrations = [rho*zi for zi in zs] return concentrations def concentrations_mass(self): @@ -6009,10 +6006,10 @@ def concentrations_mass(self): pass rho_mass = self.rho_mass() ws = self.ws() - if self.scalar: - self._concentrations_mass = [rho_mass*wi for wi in ws] - else: + if self.vectorized: self._concentrations_mass = rho_mass*ws + else: + self._concentrations_mass = [rho_mass*wi for wi in ws] return self._concentrations_mass def concentrations_mass_gas(self): @@ -6030,10 +6027,10 @@ def concentrations_mass_gas(self): ''' rho_mass = self.rho_mass_gas() ws = self.ws() - if self.scalar: - concentrations_mass = [rho_mass*wi for wi in ws] - else: + if self.vectorized: concentrations_mass = rho_mass*ws + else: + concentrations_mass = [rho_mass*wi for wi in ws] return concentrations_mass def concentrations_mass_gas_normal(self): @@ -6051,10 +6048,10 @@ def concentrations_mass_gas_normal(self): ''' rho_mass = self.rho_mass_gas_normal() ws = self.ws() - if self.scalar: - concentrations_mass = [rho_mass*wi for wi in ws] - else: + if self.vectorized: concentrations_mass = rho_mass*ws + else: + concentrations_mass = [rho_mass*wi for wi in ws] return concentrations_mass def concentrations_mass_gas_standard(self): @@ -6072,10 +6069,10 @@ def concentrations_mass_gas_standard(self): ''' rho_mass = self.rho_mass_gas_standard() ws = self.ws() - if self.scalar: - concentrations_mass = [rho_mass*wi for wi in ws] - else: + if self.vectorized: concentrations_mass = rho_mass*ws + else: + concentrations_mass = [rho_mass*wi for wi in ws] return concentrations_mass def partial_pressures(self): @@ -6099,10 +6096,10 @@ def partial_pressures(self): except: pass P = self.P - if self.scalar: - self._partial_pressures = [zi*P for zi in self.zs] - else: + if self.vectorized: self._partial_pressures = self.zs*P + else: + self._partial_pressures = [zi*P for zi in self.zs] return self._partial_pressures @@ -6115,11 +6112,11 @@ def H(self): except AttributeError: pass H = self.H_dep() - if self.scalar: + if self.vectorized: + H += float(dot(self.zs, self.Cpig_integrals_pure())) + else: for zi, Cp_int in zip(self.zs, self.Cpig_integrals_pure()): H += zi*Cp_int - else: - H += float(dot(self.zs, self.Cpig_integrals_pure())) self._H = H return H @@ -6136,20 +6133,18 @@ def S(self): R = self.R S = 0.0 - if self.scalar: + if self.vectorized: + S -= R*float(dot(zs, log_zs)) + S += float(dot(zs, Cpig_integrals_over_T_pure)) + else: S -= R*sum([zs[i]*log_zs[i] for i in cmps]) # ideal composition entropy composition for i in cmps: S += zs[i]*Cpig_integrals_over_T_pure[i] - else: - S -= R*float(dot(zs, log_zs)) - S += float(dot(zs, Cpig_integrals_over_T_pure)) S -= R*log(P*P_REF_IG_INV) S += self.S_dep() self._S = S return S - - def Cp(self): try: return self._Cp @@ -6157,11 +6152,11 @@ def Cp(self): pass Cpigs_pure = self.Cpigs_pure() Cp, zs = 0.0, self.zs - if self.scalar: + if self.vectorized: + Cp = float(dot(zs, Cpigs_pure)) + else: for i in range(self.N): Cp += zs[i]*Cpigs_pure[i] - else: - Cp = float(dot(zs, Cpigs_pure)) Cp += self.Cp_dep() self._Cp = Cp return Cp @@ -6178,7 +6173,6 @@ def dH_dP(self): dH_dP_T = dH_dP - def dH_dT_V(self): dH_dT_V = self.Cp_ideal_gas() dH_dT_V += self.dH_dep_dT_V() @@ -6204,11 +6198,11 @@ def d2H_dT2(self): pass dCpigs_pure = self.dCpigs_dT_pure() dCp, zs = 0.0, self.zs - if self.scalar: + if self.vectorized: + dCp = float(dot(zs, dCpigs_pure)) + else: for i in range(self.N): dCp += zs[i]*dCpigs_pure[i] - else: - dCp = float(dot(zs, dCpigs_pure)) dCp += self.d2H_dep_dT2() self._d2H_dT2 = dCp return dCp @@ -6216,14 +6210,13 @@ def d2H_dT2(self): def d2H_dT2_V(self): dCpigs_pure = self.dCpigs_dT_pure() dCp, zs = 0.0, self.zs - if self.scalar: + if self.vectorized: + dCp = float(dot(zs, dCpigs_pure)) + else: for i in range(self.N): dCp += zs[i]*dCpigs_pure[i] - else: - dCp = float(dot(zs, dCpigs_pure)) return dCp + self.d2H_dep_dT2_V() - def dH_dzs(self): try: return self._dH_dzs @@ -6231,20 +6224,14 @@ def dH_dzs(self): pass dH_dep_dzs = self.dH_dep_dzs() Cpig_integrals_pure = self.Cpig_integrals_pure() - if self.scalar: - self._dH_dzs = [dH_dep_dzs[i] + Cpig_integrals_pure[i] for i in range(self.N)] - else: + if self.vectorized: self._dH_dzs = dH_dep_dzs + Cpig_integrals_pure + else: + self._dH_dzs = [dH_dep_dzs[i] + Cpig_integrals_pure[i] for i in range(self.N)] return self._dH_dzs def dS_dT(self): - HeatCapacityGases = self.HeatCapacityGases - cmps = range(self.N) - T, zs = self.T, self.zs - T_REF_IG = self.T_REF_IG - P_REF_IG_INV = self.P_REF_IG_INV - - dS_dT = self.Cp_ideal_gas()/T + dS_dT = self.Cp_ideal_gas() / self.T dS_dT += self.dS_dep_dT() return dS_dT @@ -6308,10 +6295,10 @@ def dS_dzs(self): integrals = self.Cpig_integrals_over_T_pure() dS_dep_dzs = self.dS_dep_dzs() R = self.R - if self.scalar: - self._dS_dzs = [integrals[i] - R*(log_zs[i] + 1.0) + dS_dep_dzs[i] for i in cmps] - else: + if self.vectorized: self._dS_dzs = integrals - R*(log_zs + 1.0) + dS_dep_dzs + else: + self._dS_dzs = [integrals[i] - R*(log_zs[i] + 1.0) + dS_dep_dzs[i] for i in cmps] return self._dS_dzs def gammas(self): @@ -6321,10 +6308,10 @@ def gammas(self): pass phis = self.phis() phi_pures = self.phi_pures() - if self.scalar: - self._gammas = [phis[i]/phi_pures[i] for i in range(self.N)] - else: + if self.vectorized: self._gammas = phis/phi_pures + else: + self._gammas = [phis[i]/phi_pures[i] for i in range(self.N)] return self._gammas diff --git a/thermo/phases/virial_phase.py b/thermo/phases/virial_phase.py index 8e993ae9..ce870d9d 100644 --- a/thermo/phases/virial_phase.py +++ b/thermo/phases/virial_phase.py @@ -68,12 +68,11 @@ from fluids.constants import R from fluids.numerics import log, newton from fluids.numerics import numpy as np - from thermo.heat_capacity import HeatCapacityGas from thermo.phases.phase import IdealGasDeparturePhase, Phase try: - array, zeros, ones, delete, npsum, nplog = np.array, np.zeros, np.ones, np.delete, np.sum, np.log + ndarray, array, zeros, ones, delete, npsum, nplog = np.ndarray, np.array, np.zeros, np.ones, np.delete, np.sum, np.log except (ImportError, AttributeError): pass @@ -361,40 +360,40 @@ def __init__(self, Tcs, Pcs, Vcs, omegas, self.Vcs = Vcs self.omegas = omegas self.N = N = len(Tcs) - self.scalar = scalar = type(Tcs) is list + self.vectorized = vectorized = type(Tcs) is ndarray self.T = T self.B_model = B_model self.cross_B_model = cross_B_model if cross_B_model_kijs is None: - if scalar: - cross_B_model_kijs = [[0.0]*N for i in range(N)] - else: + if vectorized: cross_B_model_kijs = zeros((N, N)) + else: + cross_B_model_kijs = [[0.0]*N for i in range(N)] self.cross_B_model_kijs = cross_B_model_kijs # Parameters specific to `B` model if B_model_Meng_as is None: - if scalar: - B_model_Meng_as = [[0.0]*N for i in range(N)] - else: + if vectorized: B_model_Meng_as = zeros((N, N)) + else: + B_model_Meng_as = [[0.0]*N for i in range(N)] B_model_Meng_as_pure = [B_model_Meng_as[i][i] for i in range(N)] if B_model_Tsonopoulos_extended_as is None: - if scalar: - B_model_Tsonopoulos_extended_as = [[0.0]*N for i in range(N)] - else: + if vectorized: B_model_Tsonopoulos_extended_as = zeros((N, N)) + else: + B_model_Tsonopoulos_extended_as = [[0.0]*N for i in range(N)] B_model_Tsonopoulos_extended_as_pure = [B_model_Tsonopoulos_extended_as[i][i] for i in range(N)] if B_model_Tsonopoulos_extended_bs is None: - if scalar: - B_model_Tsonopoulos_extended_bs = [[0.0]*N for i in range(N)] - else: + if vectorized: B_model_Tsonopoulos_extended_bs = zeros((N, N)) + else: + B_model_Tsonopoulos_extended_bs = [[0.0]*N for i in range(N)] B_model_Tsonopoulos_extended_bs_pure = [B_model_Tsonopoulos_extended_bs[i][i] for i in range(N)] self.B_model_Meng_as = B_model_Meng_as @@ -446,7 +445,7 @@ def to(self, T=None): new.Vcs = self.Vcs new.omegas = self.omegas new.N = self.N - new.scalar = self.scalar + new.vectorized = self.vectorized new.B_model = self.B_model # Parameters specific to `B` model new.B_model_Meng_as = self.B_model_Meng_as @@ -477,17 +476,16 @@ def to(self, T=None): def B_interactions_at_T(self, T): N = self.N Tcijs, Pcijs, Vcijs, omegaijs = self.cross_B_model_Tcijs, self.cross_B_model_Pcijs, self.cross_B_model_Vcijs, self.cross_B_model_omegaijs - if self.scalar: - Bs = [[0.0]*N for _ in range(N)] - dB_dTs = [[0.0]*N for _ in range(N)] - d2B_dT2s = [[0.0]*N for _ in range(N)] - d3B_dT3s = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: Bs = zeros((N, N)) dB_dTs = zeros((N, N)) d2B_dT2s = zeros((N, N)) d3B_dT3s = zeros((N, N)) - + else: + Bs = [[0.0]*N for _ in range(N)] + dB_dTs = [[0.0]*N for _ in range(N)] + d2B_dT2s = [[0.0]*N for _ in range(N)] + d3B_dT3s = [[0.0]*N for _ in range(N)] if self.B_model == VIRIAL_B_ZERO: return Bs, dB_dTs, d2B_dT2s, d3B_dT3s @@ -532,16 +530,16 @@ def _set_B_and_der_interactions(self): def B_pures_at_T(self, T): N = self.N Tcs, Pcs, Vcs, omegas = self.Tcs, self.Pcs, self.Vcs, self.omegas - if self.scalar: - Bs = [0.0]*N - dB_dTs = [0.0]*N - d2B_dT2s = [0.0]*N - d3B_dT3s = [0.0]*N - else: + if self.vectorized: Bs = zeros(N) dB_dTs = zeros(N) d2B_dT2s = zeros(N) d3B_dT3s = zeros(N) + else: + Bs = [0.0]*N + dB_dTs = [0.0]*N + d2B_dT2s = [0.0]*N + d3B_dT3s = [0.0]*N if self.B_model == VIRIAL_B_ZERO: return Bs, dB_dTs, d2B_dT2s, d3B_dT3s @@ -694,16 +692,16 @@ def d3B_dT3_interactions(self): def C_interactions_at_T(self, T): N = self.N Tcijs, Pcijs, Vcijs, omegaijs = self.cross_C_model_Tcijs, self.cross_C_model_Pcijs, self.cross_C_model_Vcijs, self.cross_C_model_omegaijs - if self.scalar: - Cs = [[0.0]*N for _ in range(N)] - dC_dTs = [[0.0]*N for _ in range(N)] - d2C_dT2s = [[0.0]*N for _ in range(N)] - d3C_dT3s = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: Cs = zeros((N, N)) dC_dTs = zeros((N, N)) d2C_dT2s = zeros((N, N)) d3C_dT3s = zeros((N, N)) + else: + Cs = [[0.0]*N for _ in range(N)] + dC_dTs = [[0.0]*N for _ in range(N)] + d2C_dT2s = [[0.0]*N for _ in range(N)] + d3C_dT3s = [[0.0]*N for _ in range(N)] if self.C_model == VIRIAL_C_ZERO: return Cs, dC_dTs, d2C_dT2s, d3C_dT3s @@ -720,16 +718,16 @@ def C_interactions_at_T(self, T): def C_pures_at_T(self, T): N = self.N Tcs, Pcs, Vcs, omegas = self.Tcs, self.Pcs, self.Vcs, self.omegas - if self.scalar: - Cs = [0.0]*N - dC_dTs = [0.0]*N - d2C_dT2s = [0.0]*N - d3C_dT3s = [0.0]*N - else: + if self.vectorized: Cs = zeros(N) dC_dTs = zeros(N) d2C_dT2s = zeros(N) d3C_dT3s = zeros(N) + else: + Cs = [0.0]*N + dC_dTs = [0.0]*N + d2C_dT2s = [0.0]*N + d3C_dT3s = [0.0]*N if self.C_model == VIRIAL_C_ZERO: return Cs, dC_dTs, d2C_dT2s, d3C_dT3s @@ -984,7 +982,7 @@ def __init__(self, model, HeatCapacityGases=None, Hfs=None, Gfs=None, break if zs is not None: self.zs = zs - self.scalar = scalar = type(zs) is list + self.vectorized = vectorized = type(zs) is ndarray if T is not None: self.T = T self.model.T = T @@ -1056,10 +1054,10 @@ def dV_dzs(self): C = self.C() V = self._V N = self.N - if self.scalar: - dV_dzs = [0.0]*N - else: + if self.vectorized: dV_dzs = zeros(N) + else: + dV_dzs = [0.0]*N self._dV_dzs = dV_dzs_virial(B=B, C=C, V=V, dB_dzs=dB_dzs, dC_dzs=dC_dzs, dV_dzs=dV_dzs) return dV_dzs @@ -1087,10 +1085,10 @@ def d2V_dzizjs(self): V = self._V dV_dzs = self.dV_dzs() N = self.N - if self.scalar: - d2V_dzizjs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d2V_dzizjs = zeros((N,N)) + else: + d2V_dzizjs = [[0.0]*N for _ in range(N)] self._d2V_dzizjs = d2V_dzizjs_virial(B=B, C=C, V=V, dB_dzs=dB_dzs, dC_dzs=dC_dzs, dV_dzs=dV_dzs, d2B_dzizjs=d2B_dzizjs, d2C_dzizjs=d2C_dzizjs, d2V_dzizjs=d2V_dzizjs) return d2V_dzizjs @@ -1121,10 +1119,10 @@ def dG_dep_dzs(self): C = self.C() V = self._V N = self.N - if self.scalar: - dG_dep_dzs = [0.0]*N - else: + if self.vectorized: dG_dep_dzs = zeros(N) + else: + dG_dep_dzs = [0.0]*N for i in range(N): x0 = V x1 = x0*x0 @@ -1150,10 +1148,10 @@ def dG_dep_dns(self): except: pass N = self.N - if self.scalar: - dG_dep_dns = [0.0]*N - else: + if self.vectorized: dG_dep_dns = zeros(N) + else: + dG_dep_dns = [0.0]*N self._dG_dep_dns = dG_dep_dns = dxs_to_dns(self.dG_dep_dzs(), self.zs, dG_dep_dns) return dG_dep_dns @@ -1165,10 +1163,10 @@ def dnG_dep_dns(self): except: pass N = self.N - if self.scalar: - dnG_dep_dns = [0.0]*N - else: + if self.vectorized: dnG_dep_dns = zeros(N) + else: + dnG_dep_dns = [0.0]*N self._dnG_dep_dns = dnG_dep_dns = dxs_to_dn_partials(self.dG_dep_dzs(), self.zs, self.G_dep(), dnG_dep_dns) return dnG_dep_dns @@ -1191,12 +1189,12 @@ def lnphis(self): RT_inv = 1.0/(R*T) lnphi = self.G_dep()*RT_inv dG_dep_dns = self.dG_dep_dns() - if self.scalar: - dG_dep_dns_RT = [v*RT_inv for v in dG_dep_dns] - else: + if self.vectorized: dG_dep_dns_RT = RT_inv*dG_dep_dns + else: + dG_dep_dns_RT = [v*RT_inv for v in dG_dep_dns] - out = [0.0]*self.N if self.scalar else zeros(self.N) + out = zeros(self.N) if self.vectorized else [0.0]*self.N log_phis = dns_to_dn_partials(dG_dep_dns_RT, lnphi, out) return log_phis @@ -1809,7 +1807,7 @@ def to_TP_zs(self, T, P, zs): new.P = P new.zs = zs new.N = self.N - new.scalar = self.scalar + new.vectorized = self.vectorized new.cross_B_coefficients = self.cross_B_coefficients new.cross_C_coefficients = self.cross_C_coefficients new.cross_B_model = self.cross_B_model @@ -1827,7 +1825,7 @@ def to_TP_zs(self, T, P, zs): def to(self, zs, T=None, P=None, V=None): new = self.__class__.__new__(self.__class__) new.zs = zs - new.scalar = self.scalar + new.vectorized = self.vectorized new.N = self.N new.cross_B_coefficients = self.cross_B_coefficients new.cross_C_coefficients = self.cross_C_coefficients @@ -2150,10 +2148,10 @@ def dB_dzs(self): self._dB_dzs = dB_dzs = Bs return dB_dzs N = self.N - if self.scalar: - dB_dzs = [0.0]*N - else: + if self.vectorized: dB_dzs = zeros(N) + else: + dB_dzs = [0.0]*N B_interactions = self.model.B_interactions() self._dB_dzs = dB_dzs = dBVirial_mixture_dzs(zs, B_interactions, dB_dzs) @@ -2182,10 +2180,10 @@ def d2B_dTdzs(self): return d2B_dTdzs N = self.N - if self.scalar: - d2B_dTdzs = [0.0]*N - else: + if self.vectorized: d2B_dTdzs = zeros(N) + else: + d2B_dTdzs = [0.0]*N B_interactions = self.model.dB_dT_interactions() self._d2B_dTdzs = d2B_dTdzs = dBVirial_mixture_dzs(zs, B_interactions, d2B_dTdzs) return d2B_dTdzs @@ -2213,10 +2211,10 @@ def d3B_dT2dzs(self): return d3B_dT2dzs N = self.N - if self.scalar: - d3B_dT2dzs = [0.0]*N - else: + if self.vectorized: d3B_dT2dzs = zeros(N) + else: + d3B_dT2dzs = [0.0]*N B_interactions = self.model.d2B_dT2_interactions() self._d3B_dT2dzs = d3B_dT2dzs = dBVirial_mixture_dzs(zs, B_interactions, d3B_dT2dzs) return d3B_dT2dzs @@ -2243,10 +2241,10 @@ def d4B_dT3dzs(self): self._d4B_dT3dzs = d4B_dT3dzs = Bs return d4B_dT3dzs N = self.N - if self.scalar: - d4B_dT3dzs = [0.0]*N - else: + if self.vectorized: d4B_dT3dzs = zeros(N) + else: + d4B_dT3dzs = [0.0]*N B_interactions = self.model.d3B_dT3_interactions() self._d4B_dT3dzs = d4B_dT3dzs = dBVirial_mixture_dzs(zs, B_interactions, d4B_dT3dzs) @@ -2269,10 +2267,10 @@ def d2B_dzizjs(self): N = self.N zs = self.zs - if self.scalar: - d2B_dzizjs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d2B_dzizjs = zeros((N, N)) + else: + d2B_dzizjs = [[0.0]*N for _ in range(N)] if not self.cross_B_coefficients: self._d2B_dzizjs = d2B_dzizjs return d2B_dzizjs @@ -2300,10 +2298,10 @@ def d3B_dTdzizjs(self): N = self.N zs = self.zs - if self.scalar: - d3B_dTdzizjs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d3B_dTdzizjs = zeros((N, N)) + else: + d3B_dTdzizjs = [[0.0]*N for _ in range(N)] if not self.cross_B_coefficients: self._d3B_dTdzizjs = d3B_dTdzizjs return d3B_dTdzizjs @@ -2331,10 +2329,10 @@ def d4B_dT2dzizjs(self): N = self.N zs = self.zs - if self.scalar: - d4B_dT2dzizjs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d4B_dT2dzizjs = zeros((N, N)) + else: + d4B_dT2dzizjs = [[0.0]*N for _ in range(N)] if not self.cross_B_coefficients: self._d4B_dT2dzizjs = d4B_dT2dzizjs return d4B_dT2dzizjs @@ -2360,10 +2358,10 @@ def d3B_dzizjzks(self): N = self.N zs = self.zs - if self.scalar: - d3B_dzizjzks = [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: + if self.vectorized: d3B_dzizjzks = zeros((N, N, N)) + else: + d3B_dzizjzks = [[[0.0]*N for _ in range(N)] for _ in range(N)] if not self.cross_B_coefficients: self._d3B_dzizjzks = d3B_dzizjzks return d3B_dzizjzks @@ -2391,10 +2389,10 @@ def dB_dns(self): except: pass N = self.N - if self.scalar: - dB_dns = [0.0]*N - else: + if self.vectorized: dB_dns = zeros(N) + else: + dB_dns = [0.0]*N self._dB_dns = dB_dns = dxs_to_dns(self.dB_dzs(), self.zs, dB_dns) return dB_dns @@ -2413,10 +2411,10 @@ def dnB_dns(self): except: pass N = self.N - if self.scalar: - dnB_dns = [0.0]*N - else: + if self.vectorized: dnB_dns = zeros(N) + else: + dnB_dns = [0.0]*N self._dnB_dns = dnB_dns = dxs_to_dn_partials(self.dB_dzs(), self.zs, self.B(), dnB_dns) return dnB_dns @@ -2442,10 +2440,10 @@ def dC_dzs(self): self._dC_dzs = dC_dzs = Cs return dC_dzs N = self.N - if self.scalar: - dC_dzs = [0.0]*N - else: + if self.vectorized: dC_dzs = zeros(N) + else: + dC_dzs = [0.0]*N self._dC_dzs = dC_dzs if not self.model.C_zero: @@ -2473,10 +2471,10 @@ def d2C_dTdzs(self): self._d2C_dTdzs = d2C_dTdzs = self.model.dC_dT_pures() return d2C_dTdzs N = self.N - if self.scalar: - d2C_dTdzs = [0.0]*N - else: + if self.vectorized: d2C_dTdzs = zeros(N) + else: + d2C_dTdzs = [0.0]*N self._d2C_dTdzs = d2C_dTdzs if not self.model.C_zero: @@ -2504,10 +2502,10 @@ def d2C_dzizjs(self): pass N = self.N - if self.scalar: - d2C_dzizjs = [[0.0]*N for _ in range(N)] - else: + if self.vectorized: d2C_dzizjs = zeros((N, N)) + else: + d2C_dzizjs = [[0.0]*N for _ in range(N)] zs = self.zs if not self.cross_C_coefficients: @@ -2540,10 +2538,10 @@ def d3C_dzizjzks(self): pass N = self.N - if self.scalar: - d3C_dzizjzks = [[[0.0]*N for _ in range(N)] for _ in range(N)] - else: + if self.vectorized: d3C_dzizjzks = zeros((N, N, N)) + else: + d3C_dzizjzks = [[[0.0]*N for _ in range(N)] for _ in range(N)] zs = self.zs if not self.cross_C_coefficients: @@ -2572,10 +2570,10 @@ def dC_dns(self): except: pass N = self.N - if self.scalar: - dC_dns = [0.0]*N - else: + if self.vectorized: dC_dns = zeros(N) + else: + dC_dns = [0.0]*N self._dC_dns = dC_dns = dxs_to_dns(self.dC_dzs(), self.zs, dC_dns) return dC_dns @@ -2594,10 +2592,10 @@ def dnC_dns(self): except: pass N = self.N - if self.scalar: - dnC_dns = [0.0]*N - else: + if self.vectorized: dnC_dns = zeros(N) + else: + dnC_dns = [0.0]*N self._dnC_dns = dnC_dns = dxs_to_dn_partials(self.dC_dzs(), self.zs, self.C(), dnC_dns) return dnC_dns