added Fe benchmarks for SLB2024

misc/benchmarks/slb_2024_benchmarks.py

@@ -0,0 +1,81 @@
+from __future__ import absolute_import
+
+from burnman import Composite, equilibrate
+from burnman.tools.eos import check_eos_consistency
+from burnman.minerals import SLB_2024
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+
+fe_bcc = SLB_2024.alpha_bcc_iron()
+fe_fcc = SLB_2024.gamma_fcc_iron()
+fe_hcp = SLB_2024.epsilon_hcp_iron()
+
+check_eos_consistency(fe_bcc, verbose=True)
+
+
+temperatures = np.linspace(1.0, 2000.0, 101)
+pressures = temperatures * 0.0 + 1.0e5
+
+fig = plt.figure(figsize=(10, 5))
+ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
+
+fig1 = mpimg.imread("figures/SLB_2024_Fe_Cp_V.png")
+ax[0].imshow(fig1, extent=[0.0, 2000.0, -15.0, 70.0], aspect="auto")
+ax[1].imshow(fig1, extent=[0.0, 2000.0, 6.9, 8.4], aspect="auto")
+
+for m in [fe_bcc, fe_fcc]:
+    Cp, V = m.evaluate(["C_p", "V"], pressures, temperatures)
+    ax[0].plot(temperatures, Cp, linestyle=":", linewidth=3.0)
+    ax[1].plot(temperatures, V * 1.0e6, linestyle=":", linewidth=3.0)
+
+ax[0].set_ylim(0.0, 70.0)
+ax[1].set_ylim(6.9, 7.8)
+plt.show()
+
+
+fig = plt.figure(figsize=(10, 5))
+ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)]
+fig1 = mpimg.imread("figures/SLB_2024_Fe_phase_diagram.png")
+ax[0].imshow(fig1, extent=[0.0, 250.0, 300.0, 5800.0], aspect="auto")
+
+# Calculate the invariant point
+irons = Composite([fe_bcc, fe_fcc, fe_hcp], [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0])
+irons.set_state(5.0e9, 600.0)
+sol = equilibrate(
+    {"Fe": 1.0},
+    irons,
+    equality_constraints=[
+        ["phase_fraction", [fe_bcc, 0.0]],
+        ["phase_fraction", [fe_hcp, 0.0]],
+    ],
+)
+P_invariant = irons.pressure
+T_invariant = irons.temperature
+
+print(f"The bcc-fcc-hcp invariant is at {P_invariant/1.e9} GPa, {T_invariant} K")
+
+for Tmin, Tmax, m1, m2 in [
+    [300.0, T_invariant, fe_bcc, fe_hcp],
+    [T_invariant, 2000.0, fe_bcc, fe_fcc],
+    [T_invariant, 4800.0, fe_fcc, fe_hcp],
+]:
+
+    temperatures = np.linspace(Tmin, Tmax, 101)
+    irons = Composite([m1, m2], [1.0 / 2.0, 1.0 / 2.0])
+    sols = equilibrate(
+        {"Fe": 1.0},
+        irons,
+        equality_constraints=[["phase_fraction", [m1, 0.0]], ["T", temperatures]],
+    )[0]
+    Ps, Ts = np.array(
+        [
+            [sol.assemblage.pressure, sol.assemblage.temperature]
+            for sol in sols
+            if sol.success
+        ]
+    ).T
+    ax[0].plot(Ps / 1.0e9, Ts, linestyle=":", linewidth=3.0)
+
+plt.show()