added Fig 1c and 1d benchmarks

geodynamics · May 17, 2024 · 52caa1b · 52caa1b
1 parent 7cff9b9
commit 52caa1b
Show file tree

Hide file tree

Showing 3 changed files with 77 additions and 6 deletions.
diff --git a/misc/benchmarks/figures/SLB_2024_fper_spin_20_50_80.png b/misc/benchmarks/figures/SLB_2024_fper_spin_20_50_80.png
diff --git a/misc/benchmarks/figures/SLB_2024_fper_spin_percents.png b/misc/benchmarks/figures/SLB_2024_fper_spin_percents.png
diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py
@@ -25,14 +25,19 @@ def check_fig_1_fper_relaxed():
     pressures = np.linspace(1.0e5, 140.0e9, 141)
     temperatures = pressures * 0.0 + 300.0
 
-    fig = plt.figure(figsize=(10, 5))
+    fig = plt.figure(figsize=(12, 8))
     fig.suptitle("Figure 1 (Stixrude and Lithgow-Bertelloni, 2024)")
-    ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)]
+    ax = [fig.add_subplot(2, 2, i) for i in range(1, 5)]
 
     fig1 = mpimg.imread("figures/SLB_2024_fper_V.png")
     fig2 = mpimg.imread("figures/SLB_2024_fper_K_T.png")
+    fig3 = mpimg.imread("figures/SLB_2024_fper_spin_20_50_80.png")
+    fig4 = mpimg.imread("figures/SLB_2024_fper_spin_percents.png")
+
     ax[0].imshow(fig1, extent=[0.0, 140.0, 50.0, 80.0], aspect="auto")
     ax[1].imshow(fig2, extent=[0.0, 140.0, 100.0, 700.0], aspect="auto")
+    ax[2].imshow(fig3, extent=[0.0, 1.0, 0.0, 100.0], aspect="auto")
+    ax[3].imshow(fig4, extent=[0.0, 200.0, 0.0, 4000.0], aspect="auto")
 
     for x_fe in [0.1, 0.3, 0.5, 1.0]:
         fper.set_composition([1.0 - x_fe, x_fe, 0.0, 0.0])
@@ -43,6 +48,69 @@ def check_fig_1_fper_relaxed():
 
     ax[0].set_ylim(50.0, 80.0)
     ax[1].set_ylim(100.0, 700.0)
+
+    fper = SLB_2024.ferropericlase()
+    assemblage = simplify_composite_with_composition(
+        Composite([fper]), {"Mg": 0.5, "Fe": 0.5, "O": 1.0}
+    )
+    assemblage.set_state(50.0e9, 300.0)
+    fper = assemblage.phases[0]
+    x_MgOs = np.linspace(1.0e-5, 0.9999, 101)
+    pressures = np.empty((3, 101))
+    for i, fLS in enumerate([0.2, 0.5, 0.8]):
+        for j, x_MgO in enumerate(x_MgOs):
+            composition = {"Mg": x_MgO, "Fe": 1.0 - x_MgO, "O": 1.0}
+            fper.set_composition([x_MgO, (1.0 - x_MgO) / 2.0, (1.0 - x_MgO) / 2.0])
+            equality_constraints = [
+                ["T", 300.0],
+                [
+                    "phase_composition",
+                    (
+                        fper,
+                        (
+                            ["Fe_A", "Fels_A", "Fe_B", "Fels_B"],
+                            [0.0, 1.0, 0.0, 1.0],
+                            [1.0, 1.0, 1.0, 1.0],
+                            fLS,
+                        ),
+                    ),
+                ],
+            ]
+            equilibrate(composition, assemblage, equality_constraints)
+            pressures[i, j] = assemblage.pressure
+    ax[2].plot(
+        1.0 - x_MgOs,
+        (pressures[2] - pressures[0]) / 1.0e9,
+        linestyle=":",
+        linewidth=3.0,
+    )
+    ax[2].plot(1.0 - x_MgOs, (pressures[1]) / 1.0e9, linestyle=":", linewidth=3.0)
+
+    composition = {"Mg": 0.75, "Fe": 0.25, "O": 1.0}
+    fractions = np.linspace(0.1, 0.9, 9)
+    temperatures = np.linspace(1.0, 4000.0, 21)
+    equality_constraints = [
+        ["T", temperatures],
+        [
+            "phase_composition",
+            (
+                fper,
+                (
+                    ["Fe_A", "Fels_A", "Fe_B", "Fels_B"],
+                    [0.0, 1.0, 0.0, 1.0],
+                    [1.0, 1.0, 1.0, 1.0],
+                    fractions,
+                ),
+            ),
+        ],
+    ]
+    solss = equilibrate(composition, assemblage, equality_constraints)[0]
+    pressures = np.array([[sol.assemblage.pressure for sol in sols] for sols in solss])
+    for i, f in enumerate(fractions):
+        ax[3].plot(pressures[:, i] / 1.0e9, temperatures, linestyle=":", linewidth=3.0)
+    ax[3].set_xlim(0.0, 200.0)
+    ax[3].set_ylim(0.0, 4000.0)
+
     plt.show()
 
 
@@ -64,13 +132,16 @@ def check_fig_3_fcc_ferric_fper():
         for i, x_MgO in enumerate(x_MgOs):
             composition = {"Mg": x_MgO, "Fe": 1.0 - x_MgO, "O": 1.0}
             fper.set_composition([x_MgO, 1.0 - x_MgO, 0.0, 0.0, 0.0])
-            assemblage = Composite([fper, fcc], [1.0, 0.0])
+            assemblage = Composite([fper, fcc], [0.9, 0.1])
             assemblage = simplify_composite_with_composition(assemblage, composition)
-            equilibrate(composition, assemblage, equality_constraints)
 
-            f = assemblage.phases[0].formula
-            x_Fe3oversumFes[i] = 2.0 * ((f["O"] - f["Mg"]) / f["Fe"] - 1.0)
+            try:
+                equilibrate(composition, assemblage, equality_constraints)
 
+                f = assemblage.phases[0].formula
+                x_Fe3oversumFes[i] = 2.0 * ((f["O"] - f["Mg"]) / f["Fe"] - 1.0)
+            except AssertionError:
+                pass
         ax[0].plot(x_MgOs, x_Fe3oversumFes, linestyle=":", linewidth=3.0)
     plt.show()