Merge pull request geodynamics#667 from bobmyhill/E_to_F

Enhancement: Isothermal equations of state now require F_0, not E_0

Author: bobmyhill

burnman/calibrants/Dorogokupets_2007.py

@@ -15,7 +15,7 @@
 
 _materials_data = {
     "Ag": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -45,7 +45,7 @@
         "S": 0.732,
     },
     "Al": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -75,7 +75,7 @@
         "S": 0.998,
     },
     "Au": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -105,7 +105,7 @@
         "S": 1.067,
     },
     "Cu": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -135,7 +135,7 @@
         "S": 1.407,
     },
     "Pt": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -165,7 +165,7 @@
         "S": 0.631,
     },
     "Ta": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -195,7 +195,7 @@
         "S": 4.910,
     },
     "W": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,
@@ -225,7 +225,7 @@
         "S": 0.672,
     },
     "MgO": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 2.0,
@@ -255,7 +255,7 @@
         "S": 0,
     },
     "Diamond": {
-        "E_0": 0.0,
+        "F_0": 0.0,
         "P_0": 0.0,
         "T_0": 298.15,
         "n": 1.0,

burnman/eos/birch_murnaghan.py

@@ -276,9 +276,9 @@ def shear_modulus(self, pressure, temperature, volume, params):
         elif self.order == 3:
             return shear_modulus_third_order(volume, params)
 
-    def _molar_internal_energy(self, pressure, temperature, volume, params):
+    def _molar_helmholtz_energy(self, pressure, temperature, volume, params):
         """
-        Returns the internal energy :math:`\\mathcal{E}`
+        Returns the Helmholtz energy :math:`\\mathcal{F}`
         of the mineral. :math:`[J/mol]`
         """
         x = np.power(volume / params["V_0"], -1.0 / 3.0)
@@ -299,7 +299,7 @@ def _molar_internal_energy(self, pressure, temperature, volume, params):
             )
         )
 
-        return -intPdV + params["E_0"]
+        return -intPdV + params["F_0"]
 
     def gibbs_energy(self, pressure, temperature, volume, params):
         """
@@ -309,7 +309,7 @@ def gibbs_energy(self, pressure, temperature, volume, params):
         # G = int VdP = [PV] - int PdV = E + PV
 
         return (
-            self._molar_internal_energy(pressure, temperature, volume, params)
+            self._molar_helmholtz_energy(pressure, temperature, volume, params)
             + volume * pressure
         )
 
@@ -318,11 +318,18 @@ def validate_parameters(self, params):
         Check for existence and validity of the parameters
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 0.0
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # If G and Gprime are not included this is presumably deliberate,
         # as we can model density and bulk modulus just fine without them,
         # so just add them to the dictionary as nans
@@ -400,9 +407,9 @@ def shear_modulus(self, pressure, temperature, volume, params):
         """
         return 0.0
 
-    def _molar_internal_energy(self, pressure, temperature, volume, params):
+    def _molar_helmholtz_energy(self, pressure, temperature, volume, params):
         """
-        Returns the internal energy :math:`\\mathcal{E}` of the mineral. :math:`[J/mol]`
+        Returns the Helmholtz free energy :math:`\\mathcal{F}` of the mineral. :math:`[J/mol]`
         """
         x = np.power(volume / params["V_0"], -1.0 / 3.0)
         x2 = x * x
@@ -433,7 +440,7 @@ def _molar_internal_energy(self, pressure, temperature, volume, params):
             )
         )
 
-        return -intPdV + params["E_0"]
+        return -intPdV + params["F_0"]
 
     def validate_parameters(self, params):
         """

burnman/eos/dks_solid.py

@@ -294,12 +294,17 @@ def validate_parameters(self, params):
         if "E_0" not in params:
             params["E_0"] = float("nan")
 
-        # First, let's check the EoS parameters for Tref
-        bm.BirchMurnaghanBase.validate_parameters(bm.BirchMurnaghanBase(), params)
-
         # Now check all the required keys for the
         # thermal part of the EoS are in the dictionary
-        expected_keys = ["Cv", "grueneisen_0", "q_0", "eta_s_0"]
+        expected_keys = [
+            "V_0",
+            "K_0",
+            "Kprime_0",
+            "Cv",
+            "grueneisen_0",
+            "q_0",
+            "eta_s_0",
+        ]
         for k in expected_keys:
             if k not in params:
                 raise KeyError("params object missing parameter : " + k)

burnman/eos/macaw.py

@@ -92,9 +92,9 @@ def pressure(self, temperature, volume, params):
         term3 = C * np.power(Vrel, 1.5) + B
         return term1 * term2 * term3 - A * (B + C) + params["P_0"]
 
-    def _molar_internal_energy(self, pressure, temperature, volume, params):
+    def _molar_helmholtz_energy(self, pressure, temperature, volume, params):
         """
-        Returns the internal energy :math:`\\mathcal{E}` of the mineral. :math:`[J/mol]`
+        Returns the Helmholtz energy :math:`\\mathcal{F}` of the mineral. :math:`[J/mol]`
         """
         A, B, C = make_params(params["K_0"], params["Kprime_0"], params["Kprime_inf"])
         Vrel = volume / params["V_0"]
@@ -103,14 +103,14 @@ def _molar_internal_energy(self, pressure, temperature, volume, params):
             - 1.0
         )
         I0 = (-A * (B + C) + params["P_0"]) * params["V_0"] * (Vrel - 1.0)
-        return -A * I1 - I0
+        return params["F_0"] - A * I1 - I0
 
     def gibbs_energy(self, pressure, temperature, volume, params):
         """
         Returns the Gibbs free energy :math:`\\mathcal{G}` of the mineral. :math:`[J/mol]`
         """
         return (
-            self._molar_internal_energy(pressure, temperature, volume, params)
+            self._molar_helmholtz_energy(pressure, temperature, volume, params)
             + pressure * volume
         )
 
@@ -127,11 +127,18 @@ def validate_parameters(self, params):
         between 5/3 and :math:`K'_0` :cite:`StaceyDavis2004`.
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 1.0e5
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # Check that all the required keys are in the dictionary
         expected_keys = ["V_0", "K_0", "Kprime_0", "Kprime_inf"]
         for k in expected_keys:

burnman/eos/modified_tait.py

@@ -148,18 +148,25 @@ def gibbs_energy(self, pressure, temperature, volume, params):
             + (1.0 - a) * (pressure - params["P_0"])
         )
 
-        return intVdP + params["E_0"] + params["V_0"] * params["P_0"]
+        return intVdP + params["F_0"] + params["V_0"] * params["P_0"]
 
     def validate_parameters(self, params):
         """
         Check for existence and validity of the parameters
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 1.0e5
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # G and Gprime are not defined in this equation of state,
         # We can model density and bulk modulus just fine without them,
         # so just add them to the dictionary as nans

burnman/eos/morse_potential.py

@@ -105,9 +105,9 @@ def shear_modulus(self, pressure, temperature, volume, params):
         """
         return shear_modulus(volume, params)
 
-    def _molar_internal_energy(self, pressure, temperature, volume, params):
+    def _molar_helmholtz_energy(self, pressure, temperature, volume, params):
         """
-        Returns the internal energy :math:`\\mathcal{E}` of the mineral. :math:`[J/mol]`
+        Returns the Helmholtz energy :math:`\\mathcal{F}` of the mineral. :math:`[J/mol]`
         """
 
         x = (params["Kprime_0"] - 1) * (1 - np.power(volume / params["V_0"], 1.0 / 3.0))
@@ -120,14 +120,14 @@ def _molar_internal_energy(self, pressure, temperature, volume, params):
             * (2.0 * np.exp(x) - np.exp(2.0 * x) - 1.0)
         )
 
-        return -intPdV + params["E_0"]
+        return -intPdV + params["F_0"]
 
     def gibbs_energy(self, pressure, temperature, volume, params):
         """
         Returns the Gibbs free energy :math:`\\mathcal{G}` of the mineral. :math:`[J/mol]`
         """
         return (
-            self._molar_internal_energy(pressure, temperature, volume, params)
+            self._molar_helmholtz_energy(pressure, temperature, volume, params)
             + volume * pressure
         )
 
@@ -136,11 +136,18 @@ def validate_parameters(self, params):
         Check for existence and validity of the parameters
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 0.0
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # If G and Gprime are not included this is presumably deliberate,
         # as we can model density and bulk modulus just fine without them,
         # so just add them to the dictionary as nans

burnman/eos/murnaghan.py

@@ -68,13 +68,13 @@ def shear_modulus(self, pressure, temperature, volume, params):
         """
         return 0.0
 
-    def _molar_internal_energy(self, pressure, temperature, volume, params):
+    def _molar_helmholtz_energy(self, pressure, temperature, volume, params):
         """
-        Returns the internal energy :math:`\\mathcal{E}` of the mineral.
+        Returns the Helmholtz energy :math:`\\mathcal{F}` of the mineral.
         :math:`[J/mol]`
         """
         return energy(
-            volume, params["E_0"], params["V_0"], params["K_0"], params["Kprime_0"]
+            volume, params["F_0"], params["V_0"], params["K_0"], params["Kprime_0"]
         )
 
     def gibbs_energy(self, pressure, temperature, volume, params):
@@ -84,7 +84,7 @@ def gibbs_energy(self, pressure, temperature, volume, params):
         """
         # G = E + PV
         return (
-            self._molar_internal_energy(pressure, temperature, volume, params)
+            self._molar_helmholtz_energy(pressure, temperature, volume, params)
             + volume * pressure
         )
 
@@ -93,11 +93,18 @@ def validate_parameters(self, params):
         Check for existence and validity of the parameters
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 0.0
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # G is not included in the Murnaghan EOS so we shall set them to NaN's
         if "G_0" not in params:
             params["G_0"] = float("nan")

burnman/eos/reciprocal_kprime.py

@@ -221,7 +221,7 @@ def gibbs_energy(self, pressure, temperature, volume, params):
         # G = E0 + int VdP (when S = 0)
         K = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
         return (
-            params["E_0"]
+            params["F_0"]
             + params["P_0"] * params["V_0"]
             + self._intVdP((pressure - params["P_0"]) / K, params)
             - self._intVdP(0.0, params)
@@ -234,11 +234,18 @@ def validate_parameters(self, params):
         between 5/3 and :math:`K'_0` :cite:`StaceyDavis2004`.
         """
 
-        if "E_0" not in params:
-            params["E_0"] = 0.0
+        if "F_0" not in params:
+            params["F_0"] = 0.0
         if "P_0" not in params:
             params["P_0"] = 0.0
 
+        if "E_0" in params:
+            raise KeyError(
+                "Isothermal equations of state should be "
+                "defined in terms of Helmholtz free energy "
+                "F_0, not internal energy E_0."
+            )
+
         # If G and Gprime_inf are not included this is presumably deliberate,
         # as we can model density and bulk modulus just fine without them,
         # so just add them to the dictionary as nans