added anharmonic component to ModularMGD

Author: bobmyhill

burnman/eos/modular_mie_grueneisen_debye.py

@@ -24,6 +24,7 @@ def decorator(fn):
 from . import equation_of_state as eos
 from ..utils.math import bracket
 from . import bukowinski_electronic as el
+from .anharmonic_debye_pade import AnharmonicDebyePade as Anharmonic
 
 
 class ModularMGD(eos.EquationOfState):
@@ -91,6 +92,11 @@ def pressure(self, temperature, volume, params):
                 * (temperature * temperature - T_0 * T_0)
                 / volume
             )
+
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            P += Anharmonic.pressure(temperature, volume, params)
+
         return P
 
     def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
@@ -119,6 +125,7 @@ def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
         )
 
         KT = KT_ref + V * d2FthdV2
+
         # If the material is conductive, add the electronic contribution
         if params["bel_0"] is not None:
             KT += volume * el.KToverV(
@@ -129,6 +136,11 @@ def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
                 params["bel_0"],
                 params["gel"],
             )
+
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            KT += Anharmonic.isothermal_bulk_modulus(temperature, volume, params)
+
         return KT
 
     def _molar_heat_capacity_v(self, pressure, temperature, volume, params):
@@ -143,6 +155,11 @@ def _molar_heat_capacity_v(self, pressure, temperature, volume, params):
             bel_0 = params["bel_0"]
             gel = params["gel"]
             C_v += temperature * el.CVoverT(volume, params["V_0"], bel_0, gel)
+
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            C_v += Anharmonic.heat_capacity_v(temperature, volume, params)
+
         return C_v
 
     def thermal_expansivity(self, pressure, temperature, volume, params):
@@ -163,6 +180,10 @@ def thermal_expansivity(self, pressure, temperature, volume, params):
             gel = params["gel"]
             aKT += el.aKT(temperature, volume, params["V_0"], bel_0, gel)
 
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            aKT += Anharmonic.dSdV(temperature, volume, params)
+
         KT = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params)
 
         return aKT / KT
@@ -180,6 +201,11 @@ def entropy(self, pressure, temperature, volume, params):
             S += el.entropy(
                 temperature, volume, params["V_0"], params["bel_0"], params["gel"]
             )
+
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            S += Anharmonic.entropy(temperature, volume, params)
+
         return S
 
     def _helmholtz_energy(self, pressure, temperature, volume, params):
@@ -213,6 +239,10 @@ def _helmholtz_energy(self, pressure, temperature, volume, params):
                 params["gel"],
             )
 
+        # If the material has an anharmonic component, add it
+        if params["a_anh"] is not None:
+            F += Anharmonic.helmholtz_energy(temperature, volume, params)
+
         return F
 
     # Derived properties from here
@@ -296,6 +326,11 @@ def validate_parameters(self, params):
                 raise KeyError("params object missing parameter : " + k)
 
         # Check if material is conductive
-        if not hasattr(params, "bel_0"):
+        if "bel_0" not in params:
             params["bel_0"] = None
             params["gel"] = None
+
+        # Check if material has an anharmonic component
+        if "a_anh" not in params:
+            params["a_anh"] = None
+            params["m_anh"] = None

tests/test_modular_mgd.py

@@ -204,15 +204,46 @@ def test_SLB_electronic_contribution(self):
         params["grueneisen_0"] = 1.2
         params["q_0"] = 1.1
         params["P_0"] = 1.0e5
-        params["bel_0"] = 0.005
-        params["gel"] = 1.5
+        m_without_component = Mineral(params)
 
-        m = Mineral(params)
+        params2 = params.copy()
+        params2["bel_0"] = 0.005
+        params2["gel"] = 1.5
+
+        m = Mineral(params2)
+        consistent = check_eos_consistency(
+            m, 2.0e9, 2000.0, including_shear_properties=False, tol=1.0e-4
+        )
+        self.assertTrue(consistent)
+
+        m.set_state(2.0e9, 2000.0)
+        m_without_component.set_state(2.0e9, 2000.0)
+        self.assertFalse(m_without_component.S == m.S)
+
+    def test_SLB_anharmonic_contribution(self):
+        params = mineral_params.copy()
+        params["reference_eos"] = create("bm3")
+        params["debye_temperature_model"] = theta_SLB()
+        params["Debye_0"] = 1000.0
+        params["grueneisen_0"] = 1.2
+        params["q_0"] = 1.1
+        params["P_0"] = 1.0e5
+        m_without_component = Mineral(params)
+
+        params2 = params.copy()
+        params2["a_anh"] = 1000.0
+        params2["m_anh"] = 3.0
+
+        m = Mineral(params2)
         consistent = check_eos_consistency(
             m, 2.0e9, 2000.0, including_shear_properties=False, tol=1.0e-4
         )
         self.assertTrue(consistent)
 
+        m.set_state(2.0e9, 2000.0)
+        m_without_component.set_state(2.0e9, 2000.0)
+        self.assertFalse(m_without_component.S == m.S)
+
     def test_SPOCK_isothermal_contribution(self):
         params = mineral_params.copy()
         params["reference_eos"] = create("spock")