From 2561a7488ae1fa1f4ee4bed298dce88975c2a6e7 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Wed, 13 Nov 2024 20:57:57 +0000 Subject: [PATCH] benchmarked Decker calibration --- burnman/calibrants/Decker_1971.py | 28 +++++--------------- misc/benchmarks/calibrant_benchmarks.py | 35 +++++++++++++++++++++++++ misc/ref/calibrant_benchmarks.py.out | 28 ++++++++++++++++++++ misc/ref/example_calibrants.py.out | 16 +++++------ 4 files changed, 77 insertions(+), 30 deletions(-) diff --git a/burnman/calibrants/Decker_1971.py b/burnman/calibrants/Decker_1971.py index 60a73018..5a2dddef 100644 --- a/burnman/calibrants/Decker_1971.py +++ b/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, } diff --git a/misc/benchmarks/calibrant_benchmarks.py b/misc/benchmarks/calibrant_benchmarks.py index 7ecb40d6..b22bfaf5 100644 --- a/misc/benchmarks/calibrant_benchmarks.py +++ b/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() diff --git a/misc/ref/calibrant_benchmarks.py.out b/misc/ref/calibrant_benchmarks.py.out index 7cffaede..c6797f7f 100644 --- a/misc/ref/calibrant_benchmarks.py.out +++ b/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... diff --git a/misc/ref/example_calibrants.py.out b/misc/ref/example_calibrants.py.out index d0f96700..e0a7f54c 100644 --- a/misc/ref/example_calibrants.py.out +++ b/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)