Merge pull request #606 from bobmyhill/decker_original

benchmarked Decker calibration

burnman/calibrants/Decker_1971.py

@@ -10,6 +10,7 @@
 
 from burnman.eos.mie_grueneisen_debye import MGDBase
 from burnman.classes.calibrant import Calibrant
+from burnman.eos.birch_murnaghan_4th import birch_murnaghan_fourth
 
 
 class NaCl_B1(Calibrant):
@@ -25,43 +26,26 @@ class NaCl_B1(Calibrant):
 
     def __init__(self):
         def _pressure_Decker_NaCl(volume, temperature, params):
-
             # Isothermal pressure (GPa)
-            a = (3 / 2) * (params["Kprime_0"] - 4)
-            b = (
-                9 * params["K_0"] * params["Kprime_prime_0"]
-                + 9 * params["Kprime_0"] ** 2
-                - 63 * params["Kprime_0"]
-                + 143
-            ) / 6.0
-            f = 0.5 * ((volume / params["V_0"]) ** (-2 / 3) - 1)
-            K_T = (
-                params["K_0"]
-                * ((1 + 2 * f) ** (5 / 2))
-                * (1 + (7 + 2 * a) * f + (9 * a + 3 * b) * f**2 + 11 * b * f**3)
-            )
-            P0 = 3 * f * params["K_0"] * (1 + 2 * f) ** (5 / 2) * (1 + a * f + b * f**2)
+            P0 = birch_murnaghan_fourth(params["V_0"] / volume, params)
 
             # Thermal pressure
             thermal_model = MGDBase()
             Pth0 = thermal_model._thermal_pressure(params["T_0"], volume, params)
             Pth = thermal_model._thermal_pressure(temperature, volume, params)
 
-            # Total pressure
-            P = (P0 * 1e9) + Pth - Pth0
-
-            return P
+            return P0 + Pth - Pth0
 
         _params_Decker_NaCl = {
             "V_0": 2.7013e-05,
-            "K_0": 23.7,
+            "K_0": 23.7e9,
             "Kprime_0": 4.91,
-            "Kprime_prime_0": -0.267,
+            "Kprime_prime_0": -0.267e-9,
             "Debye_0": 279.0,
             "grueneisen_0": 1.59,
             "q_0": 0.93,
             "n": 2.0,
-            "T_0": 300.0,
+            "T_0": 298.15,
             "P_0": 0.0,
             "Z": 4.0,
         }

misc/benchmarks/calibrant_benchmarks.py

@@ -63,5 +63,40 @@ def check_figures():
     )
 
 
+def check_Decker_1971():
+    NaCl = calibrants.Decker_1971.NaCl_B1()
+    V_0 = NaCl.params["V_0"]
+
+    # First column in table
+    V = V_0
+
+    print(
+        "Pressures from Decker 1971 calibrant "
+        "vs. tabulated data in original paper (given in GPa)"
+    )
+    print(f"\nV={V:.4e} m^3/mol (standard state volume):")
+    T_C = [25.0, 100.0, 200.0, 300.0, 500.0, 800.0]
+    P_kbar = [0.00, 2.13, 5.00, 7.89, 13.72, 22.48]
+    for i, T in enumerate(T_C):
+        print(f"{T} C: {NaCl.pressure(V, T+273.15)/1.e9:.2f}, {P_kbar[i]/10.:.2f}")
+
+    # Middle column in table
+    V = (1.0 - 0.1904) * V_0
+    print(f"\nV={V:.4e} m^3/mol (middle row):")
+    T_C = [0.0, 25.0, 100.0, 200.0, 300.0, 500.0, 800.0]
+    P_kbar = [83.24, 83.93, 86.02, 88.89, 91.79, 97.65, 106.52]
+    for i, T in enumerate(T_C):
+        print(f"{T} C: {NaCl.pressure(V, T+273.15)/1.e9:.2f}, {P_kbar[i]/10.:.2f}")
+
+    # Last column in table
+    V = (1.0 - 0.2950) * V_0
+    print(f"\nV={V:.4e} m^3/mol (last row):")
+    P_kbar = [193.45, 194.12, 196.18, 199.03, 201.93, 207.81, 216.73]
+    for i, T in enumerate(T_C):
+        print(f"{T} C: {NaCl.pressure(V, T+273.15)/1.e9:.2f}, {P_kbar[i]/10.:.2f}")
+    print("")
+
+
 if __name__ == "__main__":
+    check_Decker_1971()
     check_figures()

misc/ref/calibrant_benchmarks.py.out

@@ -1,3 +1,31 @@
+Pressures from Decker 1971 calibrant vs. tabulated data in original paper (given in GPa)
+
+V=2.7013e-05 m^3/mol (standard state volume):
+25.0 C: 0.00, 0.00
+100.0 C: 0.21, 0.21
+200.0 C: 0.50, 0.50
+300.0 C: 0.79, 0.79
+500.0 C: 1.37, 1.37
+800.0 C: 2.25, 2.25
+
+V=2.1870e-05 m^3/mol (middle row):
+0.0 C: 8.32, 8.32
+25.0 C: 8.39, 8.39
+100.0 C: 8.60, 8.60
+200.0 C: 8.88, 8.89
+300.0 C: 9.17, 9.18
+500.0 C: 9.76, 9.77
+800.0 C: 10.65, 10.65
+
+V=1.9044e-05 m^3/mol (last row):
+0.0 C: 19.33, 19.34
+25.0 C: 19.40, 19.41
+100.0 C: 19.60, 19.62
+200.0 C: 19.89, 19.90
+300.0 C: 20.18, 20.19
+500.0 C: 20.77, 20.78
+800.0 C: 21.66, 21.67
+
 Checking burnman.calibrants.Fei_2007.Au...
 Checking burnman.calibrants.Fei_2007.Pt...
 Checking burnman.calibrants.Holmes_1989.Pt...

misc/ref/example_calibrants.py.out

@@ -1,26 +1,26 @@
 Experimental observations:
-Volume: 19.0287 +/- 0.1903 m^3/mol
+Volume: 19.0328 +/- 0.1903 m^3/mol
 Temperature: 1073.15 +/- 10.0 K
 
 Pressure estimated using the Decker equation of state:
-21.7355 GPa
+21.7180 GPa
 
 Pressure with propagated uncertainty
 estimated using the Decker equation of state:
-21.7355 +/- 1.0254 GPa
+21.7180 +/- 1.0245 GPa
 PVT correlation matrix:
-[[ 1.         -0.99957671  0.02909288]
- [-0.99957671  1.          0.        ]
- [ 0.02909288  0.          1.        ]]
+[[ 1.         -0.99957603  0.02911621]
+ [-0.99957603  1.          0.        ]
+ [ 0.02911621  0.          1.        ]]
 The uncertainties and correlation matrix only take into account 
 the uncertainties in the measured volume and temperature, 
 not the uncertainties in the calibrated equation of state.
 
 Consistency checks:
-Volume: 19.0287 +/- 0.1903 m^3/mol
+Volume: 19.0328 +/- 0.1903 m^3/mol
 Temperature: 1073.15 +/- 10.0000 K
 V-T correlation: 0.000000
 
 Result of conversion from Decker calibration to Decker calibration:
-21.7355 +/- 1.0254 GPa
+21.7180 +/- 1.0245 GPa
 (This should be the same as the pressures above)