Merge pull request #724 from bobmyhill/fix_H_elasticsolution

fix enthalpy calculation in ElasticSolution class, test consistency

Author: bobmyhill

burnman/classes/elasticsolution.py

@@ -554,25 +554,15 @@ def excess_enthalpy(self):
         Returns excess molar enthalpy [J/mol].
         Property specific to solutions.
         """
-        return self.solution_model.excess_enthalpy(
-            self.molar_volume, self.temperature, self.molar_fractions
-        )
+        return NotImplementedError
 
     @material_property
     def molar_enthalpy(self):
         """
         Returns molar enthalpy of the solution [J/mol].
         Aliased with self.H.
         """
-        return (
-            sum(
-                [
-                    self.solution_model.endmembers[i][0].H * self.molar_fractions[i]
-                    for i in range(self.n_endmembers)
-                ]
-            )
-            + self.excess_enthalpy
-        )
+        return self.molar_gibbs + self.temperature * self.molar_entropy
 
     @material_property
     def isothermal_bulk_modulus_reuss(self):

burnman/classes/elasticsolutionmodel.py

@@ -121,31 +121,6 @@ def excess_entropy(self, volume, temperature, molar_fractions):
             self.excess_partial_entropies(volume, temperature, molar_fractions),
         )
 
-    def excess_enthalpy(self, volume, temperature, molar_fractions):
-        """
-        Given a list of molar fractions of different phases,
-        compute the excess enthalpy of the solution.
-        The base class implementation assumes that the excess enthalpy is zero.
-
-        :param volume: Volume at which to evaluate the solution model. [m^3/mol]
-        :type volume: float
-
-        :param temperature: Temperature at which to evaluate the solution. [K]
-        :type temperature: float
-
-        :param molar_fractions: List of molar fractions of the different
-            endmembers in solution.
-        :type molar_fractions: list of floats
-
-        :returns: The excess enthalpy of the solution.
-        :rtype: float
-        """
-        return (
-            self.excess_helmholtz_energy(volume, temperature, molar_fractions)
-            + temperature * self.excess_entropy(volume, temperature, molar_fractions)
-            - volume * self.excess_pressure(volume, temperature, molar_fractions)
-        )
-
     def excess_partial_helmholtz_energies(self, volume, temperature, molar_fractions):
         """
         Given a list of molar fractions of different phases,
@@ -595,7 +570,7 @@ def __init__(
         self.Wijks = np.zeros_like(self.Wijke)
         self.Wijkp = np.zeros_like(self.Wijke)
 
-        # setup excess enthalpy interaction matrix
+        # setup excess interaction matrices
         for i in range(self.n_endmembers):
             for j in range(i + 1, self.n_endmembers):
                 w0 = energy_interaction[i][j - i - 1][0] / 2.0
@@ -647,7 +622,7 @@ def __init__(
 
         if pressure_ternary_terms is not None:
             for i, j, k, v in pressure_ternary_terms:
-                self.Wijkv[i, j, k] += v
+                self.Wijkp[i, j, k] += v
 
         # initialize ideal solution model
         ElasticIdealSolution.__init__(self, endmembers)

tests/test_elasticsolution.py

@@ -345,14 +345,6 @@ def test_set_state_with_volume(self):
         m.set_state(P1, T1)  # return to new state
         self.assertFloatEqual(V, m.V)
 
-    def test_gt_Wh(self):
-        gt_ss = garnet_ss()
-        H_excess = gt_ss.solution_model.excess_enthalpy(
-            1.0e5, 1000.0, [0.5, 0.5]
-        )  # Hxs = Exs if Vxs=0
-        We = gt_ss.solution_model.We[0][1]
-        self.assertArraysAlmostEqual([We / 4.0], [H_excess])
-
     def test_order_disorder(self):
         opx = orthopyroxene()
         opx.set_composition(np.array([0.0, 1.0]))
@@ -469,6 +461,14 @@ def test_subregular_model_ternary_partial_volume_multicomponent_change(self):
 
         self.assertArraysAlmostEqual(dVdx, dVdx2)
 
+    def test_subregular_model_ternary_eos_consistency(self):
+        ss = two_site_ss_subregular_ternary()
+        f0 = np.array([0.25, 0.35, 0.4])
+        ss.set_composition(f0)
+        ss.set_state(1.0e9, 1000.0)
+
+        self.assertTrue(burnman.tools.eos.check_eos_consistency(ss))
+
     def test_subregular(self):
         ss0 = two_site_ss()
         ss1 = two_site_ss_subregular()