From c25413200276982af2ff8e9fcd21a0f5e3cc6fec Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Sat, 11 Oct 2025 00:53:05 +0100 Subject: [PATCH] improve SLB24 benchmarking --- burnman/tools/equilibration.py | 12 ++++--- misc/benchmarks/slb_2024_benchmarks.py | 43 +++++++++++++++++++------- 2 files changed, 39 insertions(+), 16 deletions(-) diff --git a/burnman/tools/equilibration.py b/burnman/tools/equilibration.py index 8149ae85..60d51025 100644 --- a/burnman/tools/equilibration.py +++ b/burnman/tools/equilibration.py @@ -631,22 +631,22 @@ def get_equilibration_parameters(assemblage, composition, free_compositional_vec # Process parameter names prm.parameter_names = ["Pressure (Pa)", "Temperature (K)"] - prm.default_tolerances = [1.0, 1.0e-3] + prm.default_tolerances = [1.0, 1.0e-6] for i, n in enumerate(assemblage.endmembers_per_phase): prm.parameter_names.append("x({0})".format(assemblage.phases[i].name)) - prm.default_tolerances.append(1.0e-6) + prm.default_tolerances.append(1.0e-9) if n > 1: p_names = [ "p({0} in {1})".format(n, assemblage.phases[i].name) for n in assemblage.phases[i].endmember_names[1:] ] prm.parameter_names.extend(p_names) - prm.default_tolerances.extend([1.0e-6] * len(p_names)) + prm.default_tolerances.extend([1.0e-9] * len(p_names)) n_free_compositional_vectors = len(free_compositional_vectors) for i in range(n_free_compositional_vectors): prm.parameter_names.append(f"v_{i}") - prm.default_tolerances.append(1.0e-6) + prm.default_tolerances.append(1.0e-9) prm.default_tolerances = np.array(prm.default_tolerances) @@ -1036,6 +1036,10 @@ def equilibrate( if sol.success and len(assemblage.reaction_affinities) > 0.0: maxres = np.max(np.abs(assemblage.reaction_affinities)) + 1.0e-5 assemblage.equilibrium_tolerance = maxres + assert maxres < 10.0, ( + "Equilibrium not satisfied to high enough precision. " + f"Max reaction affinity is {maxres} J/mol. Try reducing tol." + ) if store_assemblage: sol.assemblage = assemblage.copy() diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py index 30ed5d46..e0bdc267 100644 --- a/misc/benchmarks/slb_2024_benchmarks.py +++ b/misc/benchmarks/slb_2024_benchmarks.py @@ -60,6 +60,7 @@ def check_fper_entropy(): def check_fig_1_fper_relaxed(): print("\nChecking Figure 1...") fper = SLB_2024.ferropericlase_relaxed() + fper_unrelaxed = SLB_2024.ferropericlase() c = molar_volume_from_unit_cell_volume(1.0, SLB_2024.periclase().params["Z"]) pressures = np.linspace(1.0e5, 140.0e9, 141) @@ -80,7 +81,26 @@ def check_fig_1_fper_relaxed(): 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]) + # First, we need to equilibrate the system at 1 bar and 1473 K + # according to the caption of SLB2024 Figure 1 + assemblage = Composite([fper_unrelaxed, SLB_2024.gamma_fcc_iron()], [0.5, 0.5]) + composition = { + "Mg": (1.0 - x_fe), + "Fe": 1.0, + "O": 1.0, + "Na": 1.0e-9, + "Al": 1.0e-9, + } + fper_unrelaxed.set_composition( + [1.0 - x_fe, x_fe * 0.98, x_fe * 0.01, 0.0, x_fe * 0.01] + ) + equality_constraints = [["P", 1.0e5], ["T", 1473.0]] + _, _ = equilibrate(composition, assemblage, equality_constraints) + f = fper_unrelaxed.molar_fractions + molar_fractions = [f[0], f[1] + f[2], f[3], f[4]] # combine HS and LS Fe + + # Now use that composition of ferropericlase in the relaxed solution + fper.set_composition(molar_fractions) Vs, K_Ss = fper.evaluate(["V", "K_S"], pressures, temperatures) ax[0].plot(pressures / 1.0e9, Vs / c, linestyle=":", linewidth=3.0) @@ -353,6 +373,7 @@ def check_fig_7_fO2(): wu = simplify_composite_with_composition( Composite([wu]), {"Fe": 1, "O": 1.01} ).phases[0] + assert len(wu.endmembers) == 3 # check that the simplification worked hepv = SLB_2024.hepv() hppv = SLB_2024.hppv() @@ -435,17 +456,15 @@ def check_fig_7_fO2(): T = 1500.0 for i, P in enumerate(pressures): assemblage.set_state(P, T) - try: - sol = equilibrate(assemblage.formula, assemblage, [["P", P], ["T", T]])[ - 0 - ] - if sol.success: - mu_O2 = assemblage.chemical_potential([{"O": 2.0}])[0] - O2_gas.set_state(1.0e5, T) - O2_gibbs = O2_gas.gibbs + f - logfO2[i] = (mu_O2 - O2_gibbs) / (gas_constant * T) / np.log(10.0) - except AssertionError: - pass + sol, _ = equilibrate( + assemblage.formula, assemblage, [["P", P], ["T", T]], tol=1.0e-10 + ) + if sol.success: + mu_O2 = sol.assemblage.chemical_potential([{"O": 2.0}])[0] + O2_gas.set_state(1.0e5, T) + O2_gibbs = O2_gas.gibbs + f + logfO2[i] = (mu_O2 - O2_gibbs) / (gas_constant * T) / np.log(10.0) + ax[0].plot(pressures / 1.0e9, logfO2, linestyle=":", linewidth=3.0) ax[0].set_ylim(-14, 22)