improvements to reactive equilibrium

tests/test_reaction.py

@@ -278,17 +278,17 @@ def rate(self, stream):
             LaEt, La, H2O, EtOH = stream.mol / stream.F_mol
             return 3600 * (kf * La * EtOH - kr * LaEt * H2O) # kmol / kg-catalyst / hr
 
-    rxn = Esterification('LacticAcid + Ethanol -> H2O + EthylLactate')
+    rxn = Esterification('LacticAcid + Ethanol -> H2O + EthylLactate', reactant='LacticAcid')
     stream = tmo.Stream(
         H2O=2, Ethanol=5, LacticAcid=1, T=355,
     )
     rxn(stream)
     assert_allclose(
         stream.mol,
         [0.0015876828181456534,
-         1.0015876828181456,
+         0.9984123171818543,
          2.001587682818146,
-         5.001587682818146],
+         4.998412317181854],
         atol=1e-3,
         rtol=1e-3,
     )
@@ -300,19 +300,16 @@ def rate(self, stream):
     stream.vle(T=T, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
-        [0.026722998037919946,
-         1.02672299803792,
-         0.9054900513291582,
-         1.9033162674693367],
+        [0.026512250430257022,
+         0.9451332614822996,
+         0.8872670426652305,
+         1.7832800372000892],
         rtol=1e-3,
         atol=1e-3,
     )
     assert_allclose(
         stream.imol['g'],
-        [0.0,
-         0.0, 
-         1.1212329467087616, 
-         3.1234067305685835],
+        [0.0, 0.028354488087443397, 1.1392452077650264, 3.1902077123696535],
         rtol=1e-3,
         atol=1e-3,
     )
@@ -324,19 +321,16 @@ def rate(self, stream):
     stream.vle(V=V, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
-        [0.026426291229780553,
-         1.0264262912297806,
-         0.9176810961603543,
-         1.932662367552195],
+        [0.0265122504353963,
+         0.9451332614841681,
+         0.8872670426420364,
+         1.7832800356800504],
         rtol=1e-3,
         atol=1e-3,
     )
     assert_allclose(
         stream.imol['g'],
-        [0.0, 
-         0.0,
-         1.1087451950694263, 
-         3.0937639236775856],
+        [0.0, 0.028354488080435617, 1.13924520779336, 3.1902077138845533],
         rtol=1e-3,
         atol=1e-3,
     )
@@ -346,16 +340,16 @@ def rate(self, stream):
     stream.vle(V=V, T=T, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
-        [0.026556271540719167,
-         1.0265562715407193,
-         0.9176843861342889,
-         1.9328505665839004],
+        [0.026512250408482565,
+         0.9451332615716925,
+         0.8872670442912585,
+         1.7832800377855502],
         rtol=1e-3,
         atol=1e-3,
     )
     assert_allclose(
         stream.imol['g'],
-        [0.0, 0.0, 1.1088718854064301, 3.0937057049568186],
+        [0.0, 0.02835448801982482, 1.1392452061172242, 3.190207711805967],
         rtol=1e-3,
         atol=1e-3,
     )
@@ -365,19 +359,19 @@ def rate(self, stream):
     stream.vle(H=H, P=P, liquid_conversion=rxn)
     assert_allclose(
         stream.imol['l'],
-        [0.026722998039327987,
-         1.026722997950048,
-         0.90549005089524, 
-         1.9033162679063462],
+        [0.026512250412710874,
+         0.945133261339238,
+         0.8872670433265364,
+         1.7832800367594983],
         rtol=1e-3,
         atol=1e-3,
     )
     assert_allclose(
         stream.imol['g'],
-        [2.3858290549839817e-12,
-         9.166582118983998e-11,
-         1.1212329471464737,
-         3.1234067301353674],
+        [4.9752893699702054e-12,
+         0.0283544882430759,
+         1.1392452070911496,
+         3.190207712822816],
         rtol=1e-3,
         atol=1e-3,
     )

thermosteam/_chemical.py

@@ -619,7 +619,8 @@ def new(cls, ID, CAS, eos=None, phase_ref=None, phase=None, **data):
 
     @classmethod
     def blank(cls, ID, CAS=None, phase_ref=None, phase=None,
-              formula=None, aliases=None, synonyms=None, free_energies=True, **data):
+              formula=None, aliases=None, synonyms=None, free_energies=True, 
+              atoms=None, **data):
         """
         Return a new Chemical object without any thermodynamic models or data 
         (unless provided).
@@ -698,6 +699,7 @@ def blank(cls, ID, CAS=None, phase_ref=None, phase=None,
         self._CAS = CAS or ID
         for i,j in data.items(): setfield(self, '_' + i , j)
         if formula: self.formula = formula
+        if atoms: self.atoms = atoms
         self._eos = create_eos(self.EOS_default, self._Tc, self._Pc, self._omega)
         TDependentProperty.RAISE_PROPERTY_CALCULATION_ERROR = False
         self._estimate_missing_properties(None)
@@ -959,6 +961,8 @@ def formula(self):
         return self._formula
     @formula.setter
     def formula(self, formula):
+        if isinstance(formula, dict):
+            formula = ''.join([f'{j}{i}' for i, j in formula.items()])
         self._formula = str(formula)
         self._MW = compute_molecular_weight(self.atoms)
         if self._Hf: self.reset_combustion_data()

thermosteam/_chemicals.py

@@ -231,6 +231,7 @@ def compile(self, skip_checks=False):
         setattr(self, '__class__', CompiledChemicals)
         try: self._compile(chemicals, skip_checks)
         except Exception as error:
+            raise error
             setattr(self, '__class__', Chemicals)
             setattr(self, '__dict__', {i.ID: i for i in chemicals})
             raise error

thermosteam/equilibrium/lle.py

@@ -165,11 +165,11 @@ class LLE(Equilibrium, phases='lL'):
                                       'popsize': 12,
                                       'tol': 1e-6}
     pseudo_equilibrium_outer_loop_options = dict(
-        xtol=1e-9, maxiter=100, checkiter=False, 
+        xtol=1e-12, maxiter=100, checkiter=False, 
         checkconvergence=False, convergenceiter=10,
     )
     pseudo_equilibrium_inner_loop_options = dict(
-        xtol=1e-6, maxiter=20, checkiter=False,
+        xtol=1e-9, maxiter=20, checkiter=False,
         checkconvergence=False, convergenceiter=5,
     )
     default_composition_cache_tolerance = 1e-5

thermosteam/equilibrium/vle.py

@@ -67,12 +67,16 @@ def xV_iter(
     V = xV[n]
     Ks = np.exp(xVlogK[-n:])
     x, y = xy(xV, Ks)
+    Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P)
     if gas_conversion: 
         z = z + gas_conversion(material=y * V, T=T, P=P, phase='g')
+        z[z <= 0] = 1e-16
+        z /= z.sum()
     if liquid_conversion:
         z = z + liquid_conversion(material=x * (1 - V), T=T, P=P, phase='l')
-    Ks[:] = pcf_Psat_over_P * f_gamma(x, T, *gamma_args) / f_phi(y, T, P)
-    V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, xV[-1], z_light, z_heavy)
+        z[z <= 0] = 1e-16
+        z /= z.sum()
+    V = binary.solve_phase_fraction_Rashford_Rice(z, Ks, V, z_light, z_heavy)
     update_xV(xV, V, Ks, z)
     Ks[Ks < 1e-16] = 1e-16
     xVlogK[-n:] = np.log(Ks)
@@ -1359,6 +1363,7 @@ def _solve_y(self, y, pcf_Psat_over_P, T, P, gas_conversion, liquid_conversion):
         xVlogK = np.zeros(2*x.size + 1)
         xVlogK[:n] = x
         xVlogK[n] = self._V
+        Ks[Ks < 1e-16] = 1e-16
         xVlogK[-n:] = np.log(Ks)
         xVlogK = flx.aitken(
             f, xVlogK, self.x_tol, args, checkiter=False, 

thermosteam/reaction/_reaction.py

@@ -1879,8 +1879,7 @@ def conversion_bounds(self, material):
         )
         return bounds
 
-    def _equilibrium_objective(self, conversion, data):
-        stream = self._stream
+    def _equilibrium_objective(self, conversion, stream, data):
         stream.mol += conversion * self._stoichiometry
         rate = self.rate(stream)
         stream.set_data(data)
@@ -1889,22 +1888,29 @@ def _equilibrium_objective(self, conversion, data):
     def equilibrium(self, stream):
         initial_condition = stream.get_data()
         f = self._equilibrium_objective
-        args = (initial_condition,)
-        conversion = flx.IQ_interpolation(f, *self.conversion_bounds(stream.imol.data), args=args)
-        return conversion * self._stoichiometry
+        args = (stream, initial_condition)
+        conversion = flx.IQ_interpolation(f, *self.conversion_bounds(stream.imol.data), xtol=1e-16, args=args)
+        return conversion
 
     def conversion_handle(self, stream):
         return KineticConversion(self, stream)
 
     def _rate(self, stream):
         conversion = self.volume(stream) * self.rate(stream)
         mol = stream.mol
-        excess = mol + conversion
-        mask = excess < 0        
-        if excess[mask].any() and (mol[mask] > 0).all():
-            x = (-mol[mask] / conversion[mask]).min()
-            conversion *= x
-            assert 0 <= x <= 1
+        full_conversion = conversion * self._stoichiometry
+        new_mol = mol + full_conversion
+        mask = new_mol < 0
+        if mask.any(): conversion = (1 - 1e-12) * (-mol[mask] / full_conversion[mask]).min() * conversion
+        # conversion_eq = self.equilibrium(stream)
+        # if conversion * conversion_eq < 0:
+        #     mol = stream.mol
+        #     full_conversion = conversion * self._stoichiometry
+        #     new_mol = mol + full_conversion
+        #     mask = new_mol < 0
+        #     if mask.any(): conversion = (1 - 1e-12) * (-mol[mask] / full_conversion[mask]).min() * conversion
+        # elif conversion_eq == 0 or conversion / conversion_eq > 1:
+        #     conversion = conversion_eq
         return conversion
 
     def conversion(self, stream, time=None):
@@ -1915,15 +1921,15 @@ def conversion(self, stream, time=None):
             conversion = self.equilibrium(stream)
         else: # Plug flow reaction (PFR)
             raise NotImplementedError('kinetic integration not yet implemented')
-        return conversion
+        return conversion * self._stoichiometry
 
     def __call__(self, stream, time=None):
         if stream.chemicals is not self.chemicals: self.reset_chemicals(stream.chemicals)
         values = stream.imol.data
         if time is None and self.volume: # Ideal continuous stirred tank CSTR
-            values += self._rate(stream)
+            values += self._rate(stream) * self._stoichiometry
         elif time is None: # Instantaneous reaction
-            values += self.equilibrium(stream)
+            values += self.equilibrium(stream) * self._stoichiometry
         else: # Plug flow reaction (PFR)
             raise NotImplementedError('kinetic integration not yet implemented')
         if tmo.reaction.CHECK_FEASIBILITY: