add anharmonic parameter validation

Author: bobmyhill

burnman/eos/anharmonic_debye.py

@@ -1,4 +1,5 @@
-# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
+# the Earth and Planetary Sciences
 # Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
@@ -31,7 +32,7 @@ def helmholtz_energy(temperature, volume, params):
         A = params["anharmonic_prefactor_model"].value(x, params)
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
         F_a = anharmonic_model.nondimensional_helmholtz_energy(
             temperature, debye_T, params
         )
@@ -46,7 +47,7 @@ def entropy(temperature, volume, params):
         A = params["anharmonic_prefactor_model"].value(x, params)
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
         S_a = anharmonic_model.nondimensional_entropy(temperature, debye_T, params)
         return A * S_a
 
@@ -56,7 +57,7 @@ def heat_capacity_v(temperature, volume, params):
         A = params["anharmonic_prefactor_model"].value(x, params)
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
         Cv_a = anharmonic_model.nondimensional_heat_capacity(
             temperature, debye_T, params
         )
@@ -69,7 +70,7 @@ def pressure(temperature, volume, params):
         dAdV = params["anharmonic_prefactor_model"].dVrel(x, params) / params["V_0"]
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
         F_a = anharmonic_model.nondimensional_helmholtz_energy(
             temperature, debye_T, params
         )
@@ -93,12 +94,11 @@ def isothermal_bulk_modulus(temperature, volume, params):
         A = params["anharmonic_prefactor_model"].value(x, params)
         dAdV = params["anharmonic_prefactor_model"].dVrel(x, params) / params["V_0"]
         d2AdV2 = (
-            params["anharmonic_prefactor_model"].d2dVrel2(x, params)
-            / params["V_0"] ** 2
+            params["anharmonic_prefactor_model"].dVrel2(x, params) / params["V_0"] ** 2
         )
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
 
         F_a = anharmonic_model.nondimensional_helmholtz_energy(
             temperature, debye_T, params
@@ -123,7 +123,7 @@ def isothermal_bulk_modulus(temperature, volume, params):
             d2AdV2 * (F_a - F_a0)
             + 2 * dAdV * (F_ad - F_ad0) * theta_model.dVrel(x, params) / params["V_0"]
             + A * (F_add - F_add0) * (theta_model.dVrel(x, params) / params["V_0"]) ** 2
-            + A * (F_ad - F_ad0) * theta_model.d2dVrel2(x, params) / params["V_0"] ** 2
+            + A * (F_ad - F_ad0) * theta_model.dVrel2(x, params) / params["V_0"] ** 2
         )
 
     @staticmethod
@@ -133,7 +133,7 @@ def dSdV(temperature, volume, params):
         dAdV = params["anharmonic_prefactor_model"].dVrel(x, params) / params["V_0"]
         theta_model = params["debye_temperature_model"]
         anharmonic_model = params["anharmonic_thermal_model"]
-        debye_T = theta_model(x, params)
+        debye_T = theta_model.value(x, params)
 
         S_a = anharmonic_model.nondimensional_entropy(temperature, debye_T, params)
         S_ad = anharmonic_model.nondimensional_dentropy_dTheta(
@@ -143,3 +143,66 @@ def dSdV(temperature, volume, params):
         aK_T = dAdV * S_a + A * (theta_model.dVrel(x, params) / params["V_0"]) * S_ad
 
         return aK_T
+
+    @staticmethod
+    def validate_parameters(params):
+        # Check for all required keys
+        expected_keys = [
+            "debye_temperature_model",
+            "anharmonic_prefactor_model",
+            "anharmonic_thermal_model",
+            "V_0",
+            "T_0",
+        ]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain an '{key}' key")
+
+        # Validate the three models:
+        models = [
+            params["debye_temperature_model"],
+            params["anharmonic_prefactor_model"],
+            params["anharmonic_thermal_model"],
+        ]
+        for model in models:
+            model.validate_parameters(params)
+
+        # Check that the required methods are present in the
+        # debye_temperature_model
+        expected_methods = ["value", "dVrel", "dVrel2"]
+        for method in expected_methods:
+            if not hasattr(params["debye_temperature_model"], method) or not callable(
+                getattr(params["debye_temperature_model"], method)
+            ):
+                raise AttributeError(
+                    f"params['debye_temperature_model'] must have a {method} method that takes arguments Vrel and params"
+                )
+
+        # Check that the required methods are present in the
+        # anharmonic_prefactor_model
+        expected_methods = ["value", "dVrel", "dVrel2"]
+        for method in expected_methods:
+            if not hasattr(
+                params["anharmonic_prefactor_model"], method
+            ) or not callable(getattr(params["anharmonic_prefactor_model"], method)):
+                raise AttributeError(
+                    f"params['anharmonic_prefactor_model'] must have a {method} method that takes arguments Vrel and params"
+                )
+
+        # Check that the required methods are present in the
+        # anharmonic_thermal_model
+        expected_methods = [
+            "nondimensional_helmholtz_energy",
+            "nondimensional_dhelmholtz_dTheta",
+            "nondimensional_d2helmholtz_dTheta2",
+            "nondimensional_entropy",
+            "nondimensional_dentropy_dTheta",
+            "nondimensional_heat_capacity",
+        ]
+        for method in expected_methods:
+            if not hasattr(params["anharmonic_thermal_model"], method) or not callable(
+                getattr(params["anharmonic_thermal_model"], method)
+            ):
+                raise AttributeError(
+                    f"params['anharmonic_thermal_model'] must have a {method} method that takes arguments temperature, debye_T, and params"
+                )

burnman/eos/anharmonic_prefactor_models.py

@@ -16,8 +16,13 @@ def value(self, Vrel, params):
     def dVrel(self, Vrel, params):
         raise NotImplementedError("need to implement dVrel() in derived class!")
 
-    def d2dVrel2(self, Vrel, params):
-        raise NotImplementedError("need to implement d2dVrel2() in derived class!")
+    def dVrel2(self, Vrel, params):
+        raise NotImplementedError("need to implement dVrel2() in derived class!")
+
+    def validate_parameters(self, params):
+        raise NotImplementedError(
+            "need to implement validate_parameters() in derived class!"
+        )
 
 
 class PowerLaw(AnharmonicPrefactorModel):
@@ -38,14 +43,26 @@ def value(self, Vrel, params):
     def dVrel(self, Vrel, params):
         return params["a_anh"] * params["m_anh"] * np.power(Vrel, params["m_anh"] - 1)
 
-    def d2dVrel2(self, Vrel, params):
+    def dVrel2(self, Vrel, params):
         return (
             params["a_anh"]
             * params["m_anh"]
             * (params["m_anh"] - 1)
             * np.power(Vrel, params["m_anh"] - 2)
         )
 
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["a_anh", "m_anh"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")
+
+        # Check that the required parameters are valid numbers
+        for key in expected_keys:
+            if not isinstance(params[key], (int, float)):
+                raise TypeError(f"params['{key}'] must be a number")
+
 
 class Sigmoid(AnharmonicPrefactorModel):
     """
@@ -74,10 +91,22 @@ def dVrel(self, Vrel, params):
         denom = (1 + (Vrel + b) ** c) ** 2
         return a * c * (Vrel + b) ** (c - 1) / denom
 
-    def d2dVrel2(self, Vrel, params):
+    def dVrel2(self, Vrel, params):
         a = params["a_anh"]
         b = params["b_anh"]
         c = params["c_anh"]
         num = (c - 1) - (c + 1) * (Vrel + b) ** c
         denom = (1 + (Vrel + b) ** c) ** 3
         return a * c * (Vrel + b) ** (c - 2) * num / denom
+
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["a_anh", "b_anh", "c_anh"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")
+
+        # Check that the required parameters are valid numbers
+        for key in expected_keys:
+            if not isinstance(params[key], (int, float)):
+                raise TypeError(f"params['{key}'] must be a number")

burnman/eos/anharmonic_thermal_models.py

@@ -305,6 +305,9 @@ def nondimensional_dentropy_dTheta(self, T, debye_T, params=None):
         t = T / debye_T
         return _d2helmholtzdt2_pade(t) * t / debye_T
 
+    def validate_parameters(self, params):
+        pass
+
 
 @singleton
 class LogNormal(object):
@@ -481,3 +484,15 @@ def nondimensional_dentropy_dTheta(self, T, debye_T, params):
         return (
             self._d2helmholtzdt2(t, params["mu_anh"], params["sigma_anh"]) * t / debye_T
         )
+
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["mu_anh", "sigma_anh"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")
+
+        # Check that the required parameters are valid numbers
+        for key in expected_keys:
+            if not isinstance(params[key], (int, float)):
+                raise TypeError(f"params['{key}'] must be a number")

burnman/eos/debye_temperature_models.py

@@ -22,7 +22,7 @@ class DebyeTemperatureModelBase(ABC):
     """Abstract base class for Debye temperature models."""
 
     @abstractmethod
-    def __call__(self, Vrel: float, params: dict) -> float:
+    def value(self, Vrel: float, params: dict) -> float:
         """Return the Debye temperature for the given relative volume."""
 
     @abstractmethod
@@ -35,13 +35,13 @@ def _gamma_prime(self, Vrel: float, params: dict) -> float:
 
     def dVrel(self, Vrel: float, params: dict) -> float:
         """First derivative of Debye temperature wrt relative volume."""
-        theta = self.__call__(Vrel, params)
+        theta = self.value(Vrel, params)
         gamma = self._gamma(Vrel, params)
         return -theta * gamma / Vrel
 
-    def d2dVrel2(self, Vrel: float, params: dict) -> float:
+    def dVrel2(self, Vrel: float, params: dict) -> float:
         """Second derivative of Debye temperature wrt relative volume."""
-        theta = self.__call__(Vrel, params)
+        theta = self.value(Vrel, params)
         gamma = self._gamma(Vrel, params)
         gamma_prime = self._gamma_prime(Vrel, params)
         return theta / Vrel**2 * (gamma * (gamma + 1.0) - Vrel * gamma_prime)
@@ -68,7 +68,7 @@ class SLB(DebyeTemperatureModelBase):
     :type params: dict
     """
 
-    def __call__(self, Vrel, params):
+    def value(self, Vrel, params):
         # This method computes the Debye temperature.
         a1_ii, a2_iikk, f = self.a1_a2_and_f(Vrel, params)
         nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f
@@ -108,6 +108,13 @@ def a1_a2_and_f(self, Vrel, params):
         f = 0.5 * (pow(Vrel, -2.0 / 3.0) - 1.0)
         return a1_ii, a2_iikk, f
 
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["grueneisen_0", "q_0", "Debye_0"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")
+
 
 @singleton
 class PowerLawGammaSimple(DebyeTemperatureModelBase):
@@ -131,7 +138,7 @@ class PowerLawGammaSimple(DebyeTemperatureModelBase):
     :type params: dict
     """
 
-    def __call__(self, Vrel, params):
+    def value(self, Vrel, params):
         # This method computes the Debye temperature.
         return params["Debye_0"] * np.exp(
             -params["grueneisen_0"]
@@ -147,6 +154,13 @@ def _gamma_prime(self, Vrel, params):
             params["grueneisen_0"] * params["q_0"] * np.power(Vrel, params["q_0"] - 1.0)
         )
 
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["grueneisen_0", "q_0", "Debye_0"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")
+
 
 @singleton
 class PowerLawGamma(DebyeTemperatureModelBase):
@@ -170,7 +184,7 @@ class PowerLawGamma(DebyeTemperatureModelBase):
     :type params: dict
     """
 
-    def __call__(self, Vrel, params):
+    def value(self, Vrel, params):
         # This method computes the Debye temperature.
         return (
             params["Debye_0"]
@@ -192,3 +206,10 @@ def _gamma_prime(self, Vrel, params):
         return params["c_1"] * params["q_1"] * np.power(
             Vrel, params["q_1"] - 1.0
         ) + params["c_2"] * params["q_2"] * np.power(Vrel, params["q_2"] - 1.0)
+
+    def validate_parameters(self, params):
+        # Check for all required keys
+        expected_keys = ["Debye_0", "grueneisen_0", "c_1", "c_2", "q_1", "q_2"]
+        for key in expected_keys:
+            if key not in params:
+                raise AttributeError(f"params dictionary must contain a '{key}' key")