diff --git a/burnman/classes/solutionmodel.py b/burnman/classes/solutionmodel.py index fe6ef84f..d7736696 100644 --- a/burnman/classes/solutionmodel.py +++ b/burnman/classes/solutionmodel.py @@ -108,7 +108,7 @@ def _non_ideal_interactions_subreg(p, n_endmembers, Wijk): return Wint -def logish(x, eps=1.0e-5): +def logish(x, eps=1.0e-7): """ 2nd order series expansion of log(x) about eps: log(eps) - sum_k=1^infty (f_eps)^k / k @@ -434,7 +434,13 @@ def __init__(self, endmembers): def _calculate_endmember_configurational_entropies(self): S_conf = -( constants.gas_constant - * (self.endmember_noccupancies * logish(self.endmember_occupancies)).sum(-1) + * ( + self.endmember_noccupancies + * ( + logish(self.endmember_noccupancies) + - logish(self.site_multiplicities) + ) + ).sum(-1) ) self.endmember_configurational_entropies = S_conf diff --git a/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py b/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py index 6be7776e..bec3e914 100644 --- a/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py +++ b/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py @@ -61,7 +61,7 @@ def rfloat(x, m=1.0): hardcoded_site_occupancies = { "wu": "[Fe]2[Fe]2", "wuls": "[Fels]2[Fels]2", - "mag": "[Vac1/2Fe1/2]2[Fef]2", + "mag": "[Fef]0[Fef]2", # first site not counted "smag": "[Fe][Fef][Fef]", "hepv": "[Fef][Fef]", "hlpv": "[Fef][Fefls]", @@ -79,7 +79,6 @@ def rfloat(x, m=1.0): "mag": {"p": 0.4}, "smag": {"p": 0.4}, "hmag": {"p": 0.4}, - "wu": {"p": 0.6}, } hefesto_path = "HeFESTo_parameters_010123" @@ -424,6 +423,8 @@ def rfloat(x, m=1.0): "from __future__ import absolute_import\n" "\n\n" "from ..eos.slb import SLB3\n" + "from ..eos import debye\n" + "from ..eos import birch_murnaghan as bm\n" "from ..classes.mineral import Mineral\n" "from ..classes.solution import Solution, RelaxedSolution\n" "from ..classes.solutionmodel import (\n" @@ -437,22 +438,40 @@ def rfloat(x, m=1.0): st_class = ( "class SLB3Stishovite(SLB3):\n" " def shear_modulus(self, pressure, temperature, volume, params):\n" - " G_bare = SLB3.shear_modulus(self, pressure, temperature, volume, params)\n" + ' """\n' + " Returns the shear modulus. :math:`[Pa]`\n" + " SLB2024 first calculate the shear modulus at T_0 using\n" + " the third order Birch-Murnaghan equation of state,\n" + " then apply the softening relevant to the temperature\n" + " and pressure of interest, and then apply the usual thermal term.\n" + ' """\n' "\n" - " Pc = 51.6e9 + 11.1*(temperature - 300.)*1.e6 # Eq B4\n" - " PcQ = Pc + 50.7e9 # From the HeFESTo code in stishtran.f. PcQ is T-independent in paper.\n" - " Cs0 = 128.e9*PcQ/Pc\n" - " Delta = 0. # 100.e9 in the HeFESTo code in stishtran.f, 0. in paper.\n" + " G_bare = bm.shear_modulus_third_order(volume, params)\n" + " Pc = 51.6e9 + 11.1 * (temperature - 300.0) * 1.0e6 # Eq B4\n" + " PcQ = (\n" + " Pc + 50.7e9\n" + " ) # From the HeFESTo code in stishtran.f. PcQ is T-independent in paper.\n" + " Cs0 = 128.0e9 * PcQ / Pc\n" + " Delta = 0.0 # 100.e9 in the HeFESTo code in stishtran.f, 0. in paper.\n" " PminusPc = pressure - Pc\n" - " if (PminusPc < 0.):\n" - " Cs = Cs0*PminusPc/(pressure - PcQ) # Eq B2\n" + " if PminusPc < 0.0:\n" + " Cs = Cs0 * PminusPc / (pressure - PcQ) # Eq B2\n" " else:\n" - " Ast = (Cs0 - Delta)*(PcQ - Pc)\n" - " Cs = Cs0 - Ast/(PcQ - pressure + 3.*PminusPc) # Eq B3\n" + " Ast = (Cs0 - Delta) * (PcQ - Pc)\n" + " Cs = Cs0 - Ast / (PcQ - pressure + 3.0 * PminusPc) # Eq B3\n" "\n" " # VRH-like average of the bare shear modulus with softened (C11-C12)/2.\n" - " G = 0.5*(13./15.*G_bare + 2./15.*Cs) + 0.5/(13./15./G_bare + 2./15./Cs) # Eq B5\n" - " return G\n" + " G = 0.5 * (13.0 / 15.0 * G_bare + 2.0 / 15.0 * Cs) + 0.5 / (\n" + " 13.0 / 15.0 / G_bare + 2.0 / 15.0 / Cs\n" + " ) # Eq B5\n" + "\n" + " debye_T = self._debye_temperature(params['V_0'] / volume, params)\n" + " eta_s = self._isotropic_eta_s(params['V_0'] / volume, params)\n" + "\n" + " E_th = debye.thermal_energy(temperature, debye_T, params['n'])\n" + " E_th_ref = debye.thermal_energy(params['T_0'], debye_T, params['n'])\n" + "\n" + " return G - eta_s * (E_th - E_th_ref) / volume\n" ) fper_relaxed_class = ( @@ -463,7 +482,8 @@ def rfloat(x, m=1.0): " - pe ([Mg]2[Mg]2)\n" " - wu and wuls ([Fe,Fels]2[Fe,Fels]2)\n" " - anao ([Na]2[Al]2)\n" - " - mag ([Vac1/2Fe1/2]2[Fef]2)\n" + " - mag ([Fe]0[Fef]2)\n" + " The entropy from the first site in magnetite is not counted.\n" " The equilibrium spin state is calculated automatically.\n" ' """\n' " solution = ferropericlase()\n" @@ -486,7 +506,7 @@ def rfloat(x, m=1.0): " - mgpv ([Mg][Si])\n" " - fepv ([Fe][Si])\n" " - alpv ([Al][Al])\n" - " - hepv + hlpv ([Fef][Fef,Fels])\n" + " - hepv + hlpv ([Fef][Fef,Fefls])\n" " - fapv ([Fef][Al])\n" " - crpv ([Cr][Cr])\n" " The equilibrium spin state is calculated automatically.\n" @@ -537,6 +557,8 @@ def rfloat(x, m=1.0): docstring += "Endmembers (and site species distributions) are given in the order:\n" for i in range(prm["n_mbrs"]): docstring += f'- {prm["mbr_names"][i]} ({prm["mbr_site_formulae"][i]})\n' + if key == "mw": + docstring += "The entropy from the first site in magnetite is not counted.\n" docstring += '"""' print( diff --git a/burnman/eos/bukowinski_electronic.py b/burnman/eos/bukowinski_electronic.py index cbe6fd12..825fa63d 100644 --- a/burnman/eos/bukowinski_electronic.py +++ b/burnman/eos/bukowinski_electronic.py @@ -12,33 +12,52 @@ """ -def helmholtz(temperature, volume, V_0, bel_0, gel): - return -0.5 * bel_0 * np.power(volume / V_0, gel) * temperature * temperature +def helmholtz(temperature, volume, T_0, V_0, bel_0, gel): + return ( + -0.5 + * bel_0 + * np.power(volume / V_0, gel) + * (temperature * temperature - T_0 * T_0) + ) -def pressure(temperature, volume, V_0, bel_0, gel): +def pressure(temperature, volume, T_0, V_0, bel_0, gel): + """ + P = -dF/dV + """ return ( 0.5 * gel * bel_0 * np.power(volume / V_0, gel) - * temperature - * temperature + * (temperature * temperature - T_0 * T_0) / volume ) -def KToverV(temperature, volume, V_0, bel_0, gel): - return -(gel - 1.0) * pressure(temperature, volume, V_0, bel_0, gel) / volume +def KToverV(temperature, volume, T_0, V_0, bel_0, gel): + """ + KT = -V dP/dV + """ + return -(gel - 1.0) * pressure(temperature, volume, T_0, V_0, bel_0, gel) / volume def entropy(temperature, volume, V_0, bel_0, gel): + """ + S = -dF/dT + """ return bel_0 * np.power(volume / V_0, gel) * temperature def CVoverT(volume, V_0, bel_0, gel): + """ + CV = T dS/dT + """ return bel_0 * np.power(volume / V_0, gel) def aKT(temperature, volume, V_0, bel_0, gel): + """ + aKT = dP/dT + """ return gel * bel_0 * np.power(volume / V_0, gel) * temperature / volume diff --git a/burnman/eos/slb.py b/burnman/eos/slb.py index 7ebb6224..09438419 100644 --- a/burnman/eos/slb.py +++ b/burnman/eos/slb.py @@ -66,7 +66,14 @@ def _delta_pressure( nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41 gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f) - Pel = 0.5 * gel * bel_0 * np.power(x / V_0, gel) * temperature * temperature / x + Pel = ( + 0.5 + * gel + * bel_0 + * np.power(x / V_0, gel) + * (temperature * temperature - T_0 * T_0) + / x + ) return ( (1.0 / 3.0) @@ -272,10 +279,14 @@ def pressure(self, temperature, volume, params): ) + gr * ( E_th - E_th_ref ) / volume # EQ 21 - if self.conductive: P += el.pressure( - temperature, volume, params["V_0"], params["bel_0"], params["gel"] + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], ) return P @@ -309,7 +320,12 @@ def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params): if self.conductive: K += volume * el.KToverV( - temperature, volume, params["V_0"], params["bel_0"], params["gel"] + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], ) return K @@ -408,7 +424,12 @@ def helmholtz_free_energy(self, pressure, temperature, volume, params): if self.conductive: F += el.helmholtz( - temperature, volume, params["V_0"], params["bel_0"], params["gel"] + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], ) return F diff --git a/burnman/minerals/SLB_2024.py b/burnman/minerals/SLB_2024.py index acaa408a..68b67d9a 100644 --- a/burnman/minerals/SLB_2024.py +++ b/burnman/minerals/SLB_2024.py @@ -14,6 +14,8 @@ from ..eos.slb import SLB3 +from ..eos import debye +from ..eos import birch_murnaghan as bm from ..classes.mineral import Mineral from ..classes.solution import Solution, RelaxedSolution from ..classes.solutionmodel import ( @@ -1898,8 +1900,15 @@ def __init__(self): class SLB3Stishovite(SLB3): def shear_modulus(self, pressure, temperature, volume, params): - G_bare = SLB3.shear_modulus(self, pressure, temperature, volume, params) + """ + Returns the shear modulus. :math:`[Pa]` + SLB2024 first calculate the shear modulus at T_0 using + the third order Birch-Murnaghan equation of state, + then apply the softening relevant to the temperature + and pressure of interest, and then apply the usual thermal term. + """ + G_bare = bm.shear_modulus_third_order(volume, params) Pc = 51.6e9 + 11.1 * (temperature - 300.0) * 1.0e6 # Eq B4 PcQ = ( Pc + 50.7e9 @@ -1917,7 +1926,14 @@ def shear_modulus(self, pressure, temperature, volume, params): G = 0.5 * (13.0 / 15.0 * G_bare + 2.0 / 15.0 * Cs) + 0.5 / ( 13.0 / 15.0 / G_bare + 2.0 / 15.0 / Cs ) # Eq B5 - return G + + debye_T = self._debye_temperature(params["V_0"] / volume, params) + eta_s = self._isotropic_eta_s(params["V_0"] / volume, params) + + E_th = debye.thermal_energy(temperature, debye_T, params["n"]) + E_th_ref = debye.thermal_energy(params["T_0"], debye_T, params["n"]) + + return G - eta_s * (E_th - E_th_ref) / volume class st(Mineral): @@ -1990,18 +2006,7 @@ def __init__(self): } self.property_modifiers = [ - [ - "magnetic_chs", - { - "structural_parameter": 0.6, - "curie_temperature": [191.0, 0.0], - "magnetic_moment": [623.9214193067143, 0.0], - }, - ], - [ - "linear", - {"delta_E": 10138.21935759185, "delta_S": 53.5254, "delta_V": 0.0}, - ], + ["landau_slb_2022", {"Tc_0": 191.0, "S_D": 53.5254, "V_D": 0.0}] ] Mineral.__init__(self) @@ -2202,7 +2207,8 @@ def __init__(self, molar_fractions=None): - wu ([Fe]2[Fe]2) - wuls ([Fels]2[Fels]2) - anao ([Na]2[Al]2) - - mag ([Vac1/2Fe1/2]2[Fef]2) + - mag ([Fef]0[Fef]2) + The entropy from the first site in magnetite is not counted. """ self.name = "ferropericlase" self.solution_model = AsymmetricRegularSolution( @@ -2211,7 +2217,7 @@ def __init__(self, molar_fractions=None): [wu(), "[Fe]2[Fe]2"], [wuls(), "[Fels]2[Fels]2"], [anao(), "[Na]2[Al]2"], - [mag(), "[Vac1/2Fe1/2]2[Fef]2"], + [mag(), "[Fef]0[Fef]2"], ], alphas=[1.0, 1.0, 1.0, 1.0, 0.08293], energy_interaction=[ @@ -2238,7 +2244,8 @@ def __init__(self): - pe ([Mg]2[Mg]2) - wu and wuls ([Fe,Fels]2[Fe,Fels]2) - anao ([Na]2[Al]2) - - mag ([Vac1/2Fe1/2]2[Fef]2) + - mag ([Fe]0[Fef]2) + The entropy from the first site in magnetite is not counted. The equilibrium spin state is calculated automatically. """ solution = ferropericlase() @@ -2408,7 +2415,7 @@ def __init__(self): - mgpv ([Mg][Si]) - fepv ([Fe][Si]) - alpv ([Al][Al]) - - hepv + hlpv ([Fef][Fef,Fels]) + - hepv + hlpv ([Fef][Fef,Fefls]) - fapv ([Fef][Al]) - crpv ([Cr][Cr]) The equilibrium spin state is calculated automatically. diff --git a/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png b/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png new file mode 100644 index 00000000..e63059a7 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_fO2.png b/misc/benchmarks/figures/SLB_2024_fper_fO2.png new file mode 100644 index 00000000..7137ce98 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_fO2.png differ diff --git a/misc/benchmarks/slb_2024_Fe_endmembers.dat b/misc/benchmarks/slb_2024_Fe_endmembers.dat new file mode 100644 index 00000000..c268b7cb --- /dev/null +++ b/misc/benchmarks/slb_2024_Fe_endmembers.dat @@ -0,0 +1,13 @@ + # id spec Pi Ti Enthalpy Entropy Cp Gibbs Volume Smix Sconf CV + 1 wu 0.00 833.00 -787.6607 460.0267 207.0322 -1170.8630 50.0110 6.8091 -0.0000 197.0969 414.6426 + 2 mag 0.00 833.00 -863.5328 345.8858 288.7717 -1151.6557 45.2057 18.1363 -0.0000 282.8666 518.6847 + 3 wu 0.00 833.00 -787.6607 460.0267 207.0322 -1170.8630 50.0110 11.5263 -0.0000 197.0969 414.6426 + 4 mag 0.00 833.00 -863.5328 345.8858 288.7717 -1151.6557 45.2057 11.5263 -0.0000 282.8666 518.6847 + 5 hem 0.00 833.00 -646.7135 221.0119 164.2131 -830.8164 30.7971 0.0000 -0.0000 158.3723 636.3544 + 6 hepv 0.00 833.00 -618.1526 232.5549 126.4491 -811.8709 30.0608 5.7632 -0.0000 121.2145 630.6739 + 7 hlpv 0.00 833.00 -515.1353 204.1707 125.4698 -685.2095 27.9927 5.7632 -0.0000 119.9336 739.8142 + 8 hppv 0.00 833.00 -533.6867 227.6978 128.9032 -723.3590 28.2999 0.0000 -0.0000 120.9283 656.4890 + 9 fea 0.00 833.00 25.0948 59.0890 39.5161 -24.1263 7.2463 0.0000 -0.0000 37.7495 383.9275 + 10 feg 0.00 833.00 30.4729 64.4276 30.2506 -23.1952 7.1594 0.0000 -0.0000 28.0479 272.5673 + 11 fee 0.00 833.00 27.4641 58.1046 31.3715 -20.9371 6.9617 0.0000 -0.0000 28.3096 356.9375 + 12 hmag 0.00 833.00 -840.9691 343.3957 291.6978 -1127.0177 42.5359 0.0000 -0.0000 282.7694 526.2834 diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py index ba6bf694..ae6fd7fc 100644 --- a/misc/benchmarks/slb_2024_benchmarks.py +++ b/misc/benchmarks/slb_2024_benchmarks.py @@ -7,18 +7,61 @@ from burnman.tools.unitcell import molar_volume_from_unit_cell_volume from burnman.minerals import SLB_2024 from burnman.minerals.HP_2011_fluids import O2 - +from burnman.eos.property_modifiers import calculate_property_modifications import numpy as np import matplotlib.pyplot as plt import matplotlib.image as mpimg +def check_iron_bearing_mineral_property_values(): + print("Checking iron bearing mineral properties...") + m_names = np.genfromtxt("slb_2024_Fe_endmembers.dat", dtype=str)[:, 1] + minerals = [eval("SLB_2024." + m)() for m in m_names] + d = np.genfromtxt("slb_2024_Fe_endmembers.dat", dtype=float)[:, 2:] + for i, m in enumerate(minerals): + P, T, _, entropy, Cp, gibbs, volume, _, _, _, _ = d[i] + m.set_state(P * 1.0e4, T) + delta_gibbs = gibbs * 1.0e3 - m.gibbs + if np.abs(delta_gibbs) > 1.0: + print(m_names[i]) + print("delta Gibbs:", delta_gibbs) + print("delta S:", entropy - m.S) + print("delta Cp:", Cp - m.C_p) + print("delta V:", volume * 1.0e-6 - m.V) + else: + rel_error = max( + [ + (entropy - m.S) / entropy, + (Cp - m.C_p) / Cp, + (volume - m.V * 1.0e6) / volume, + ] + ) + print( + f"{m_names[i]} in agreement with HeFESTo (max rel error {rel_error*100.:.1g} %)" + ) + + def check_bcc_iron_consistency(): + print("") fe_bcc = SLB_2024.alpha_bcc_iron() - check_eos_consistency(fe_bcc, verbose=True) + assert check_eos_consistency(fe_bcc, verbose=True) + + +def check_fper_entropy(): + print("\nChecking fper configurational entropy...") + fper = SLB_2024.ferropericlase() + print("Endmember entropies (should all be zero):") + print(np.abs(fper.solution_model.endmember_configurational_entropies)) + + fper.set_composition([0.0, 0.5, 0.0, 0.0, 0.5]) + Sxs = -2 * gas_constant * np.log(0.5) + print( + f"Equimolar wu-mag: {fper.excess_entropy:.5f} J/K/mol (should be {Sxs:.5f} J/K/mol)" + ) def check_fig_1_fper_relaxed(): + print("\nChecking Figure 1...") fper = SLB_2024.ferropericlase_relaxed() c = molar_volume_from_unit_cell_volume(1.0, SLB_2024.periclase().params["Z"]) @@ -49,6 +92,13 @@ def check_fig_1_fper_relaxed(): ax[0].set_ylim(50.0, 80.0) ax[1].set_ylim(100.0, 700.0) + fper.set_composition([0.5, 0.5, 0.0, 0.0]) + fper.set_state(70.0e9, 300.0) + print("Mg2Fe2O4 ferropericlase at 70 GPa, 300 K:") + fstr = ", ".join([f"{f:.2f}" for f in fper.molar_fractions]) + print(f"Molar proportions: {fstr}") + print(f"Volume: {fper.V*1.e6:.2f} cm^3/mol") + fper = SLB_2024.ferropericlase() assemblage = simplify_composite_with_composition( Composite([fper]), {"Mg": 0.5, "Fe": 0.5, "O": 1.0} @@ -115,38 +165,104 @@ def check_fig_1_fper_relaxed(): def check_fig_3_fcc_ferric_fper(): + print("\nChecking Figure 3...") + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 3 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)] + + fig1 = mpimg.imread("figures/SLB_2024_Fe_O_phase_diagram.png") + fig2 = mpimg.imread("figures/SLB_2024_fper_fO2.png") + + ax[0].imshow(fig1, extent=[-1, -0.75, 700.0, 1800.0], aspect="auto") + ax[1].imshow(fig2, extent=[0.0, 1.0, 0.0, 0.12], aspect="auto") + + # Subplot a + bcc = SLB_2024.alpha_bcc_iron() fcc = SLB_2024.gamma_fcc_iron() fper = SLB_2024.ferropericlase() + mag = SLB_2024.smag() + + assemblage = simplify_composite_with_composition( + Composite([fper]), {"Fe": 1.0, "O": 1.1} + ) + fper = assemblage.phases[0] - x_MgOs = np.linspace(0.001, 0.999, 101) - x_Fe3oversumFes = np.empty_like(x_MgOs) + composition = {"Fe": 1.0, "O": 1.0} + assemblage = Composite([bcc, fper, mag], [0.5, 0.49, 0.01]) + fper.set_composition([0.05, 0.0, 0.95]) + equality_constraints = [["P", 1.0e5], ["phase_fraction", (mag, 0.0)]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + T_a_wu_mag = sol[0].assemblage.temperature + print(f"BCC-fper-mag triple point at 1 bar: {T_a_wu_mag:.2f} K") + + assemblage = Composite([bcc, fcc], [0.5, 0.5]) + equality_constraints = [["P", 1.0e5], ["phase_fraction", (bcc, 0.0)]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + T_a_g = sol[0].assemblage.temperature + + assemblage = Composite([bcc, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_wu_mag, T_a_g, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + assemblage = Composite([fcc, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_g, 1800.0, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + assemblage = Composite([mag, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_wu_mag, 1800.0, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + # Subplot b + fcc = SLB_2024.gamma_fcc_iron() + fper = SLB_2024.ferropericlase() - fig = plt.figure(figsize=(10, 5)) - fig.suptitle("Figure 3 (Stixrude and Lithgow-Bertelloni, 2024)") - ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + x_MgOs = np.linspace(0.01, 0.99, 33) + x_Mgs = np.empty_like(x_MgOs) * np.nan + x_Fe3oversumFes = np.empty_like(x_MgOs) * np.nan + composition = {"Mg": 0.01, "Fe": 2.0, "O": 1.0} + assemblage = Composite([fper, fcc], [0.5, 0.5]) + x_MgO = 0.01 + fper.set_composition([x_MgO, 0.95 * (1.0 - x_MgO), 0.0, 0.0, 0.05 * (1.0 - x_MgO)]) + assemblage = simplify_composite_with_composition(assemblage, composition) for equality_constraints in [ [["P", 1.0e5], ["T", 1473.0]], [["P", 21.0e9], ["T", 1673.0]], ]: 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], [0.9, 0.1]) - assemblage = simplify_composite_with_composition(assemblage, composition) - + composition = {"Mg": x_MgO, "Fe": 2.0, "O": 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) + sol = equilibrate( + composition, assemblage, equality_constraints, tol=1.0e-5 + ) + if sol[0].success: + f = assemblage.phases[0].formula + x_Mgs[i] = f["Mg"] / 4.0 + 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) + ax[1].plot(x_Mgs, x_Fe3oversumFes, linestyle=":", linewidth=3.0) plt.show() def check_fig_6a_iron_Cp_V(): + print("\nChecking Figure 6a...") fe_bcc = SLB_2024.alpha_bcc_iron() fe_fcc = SLB_2024.gamma_fcc_iron() @@ -166,17 +282,22 @@ def check_fig_6a_iron_Cp_V(): ax[0].plot(temperatures, Cp, linestyle=":", linewidth=3.0) ax[1].plot(temperatures, V * 1.0e6, linestyle=":", linewidth=3.0) + fe_bcc.set_state(1.0e5, 1000.0) + print( + f"BCC iron Cp, V at 1 bar, 1000 K: {fe_bcc.C_p:.2f} J/K/mol, {fe_bcc.V*1.e6:.2f} cm^3/mol" + ) ax[0].set_ylim(0.0, 70.0) ax[1].set_ylim(6.9, 7.8) plt.show() def check_fig_6c_iron_phase_diagram(): + print("\nChecking Figure 6c...") fe_bcc = SLB_2024.alpha_bcc_iron() fe_fcc = SLB_2024.gamma_fcc_iron() fe_hcp = SLB_2024.epsilon_hcp_iron() - fig = plt.figure(figsize=(10, 5)) + fig = plt.figure(figsize=(8, 5)) fig.suptitle("Figure 6c (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] fig1 = mpimg.imread("figures/SLB_2024_Fe_phase_diagram.png") @@ -197,7 +318,7 @@ def check_fig_6c_iron_phase_diagram(): T_invariant = irons.temperature print( - f"The bcc-fcc-hcp invariant is at {P_invariant/1.e9:.3f} GPa, {T_invariant:.3f} K" + f"The bcc-fcc-hcp invariant is at {P_invariant/1.e9:.2f} GPa, {T_invariant:.1f} K" ) for Tmin, Tmax, m1, m2 in [ @@ -226,6 +347,7 @@ def check_fig_6c_iron_phase_diagram(): def check_fig_7_fO2(): + print("\nChecking Figure 7...") hem = SLB_2024.hematite() mag = SLB_2024.smag() hmag = SLB_2024.high_pressure_magnetit() @@ -247,7 +369,7 @@ def check_fig_7_fO2(): O2_gas.set_state(1.0e5, 298.15) f = O2_gas.S * 298.15 - fig = plt.figure(figsize=(10, 5)) + fig = plt.figure(figsize=(8, 5)) fig.suptitle("Figure 7 (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] @@ -288,7 +410,7 @@ def check_fig_7_fO2(): ["phase_fraction", [a.phases[0], 0.0]], ["T", 1500.0], ] - equilibrate(a.phases[0].formula, a, equality_constraints) + equilibrate(a.phases[0].formula, a, equality_constraints, tol=1.0e-5) P_bounds[i][1] = a.pressure P_bounds[i + 1][0] = a.pressure @@ -304,7 +426,7 @@ def check_fig_7_fO2(): a8 = Composite([Fe_fcc, wu], [0.5, 0.5]) a9 = Composite([Fe_hcp, wu], [0.5, 0.5]) - fig = plt.figure(figsize=(10, 5)) + fig = plt.figure(figsize=(8, 5)) ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] fig1 = mpimg.imread("figures/SLB_2024_Fe_O_fO2.png") @@ -312,16 +434,21 @@ def check_fig_7_fO2(): for i, assemblage in enumerate([a1, a2, a3, a4, a5, a6, a7, a8, a9]): pressures = np.linspace(P_bounds[i][0], P_bounds[i][1], 11) - logfO2 = np.empty_like(pressures) + logfO2 = np.empty_like(pressures) * np.nan T = 1500.0 for i, P in enumerate(pressures): assemblage.set_state(P, T) - 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) + try: + sol = equilibrate( + assemblage.formula, assemblage, [["P", P], ["T", T]], tol=1.0e-7 + )[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 ax[0].plot(pressures / 1.0e9, logfO2, linestyle=":", linewidth=3.0) ax[0].set_ylim(-14, 22) @@ -329,12 +456,13 @@ def check_fig_7_fO2(): def check_fig_a1_fe2o3(): + print("\nChecking Figure A1...") hem = SLB_2024.hematite() hppv = SLB_2024.hppv() bdg = SLB_2024.bridgmanite_relaxed() bdg.set_composition([0.0, 0.0, 0.0, 1.0, 0.0, 0.0]) - fig = plt.figure(figsize=(10, 5)) + fig = plt.figure(figsize=(8, 5)) fig.suptitle("Figure A1 (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] @@ -351,13 +479,16 @@ def check_fig_a1_fe2o3(): volumes = bdg.evaluate(["V"], pressures, temperatures)[0] ax[0].plot(pressures / 1.0e9, volumes * 1.0e6, linestyle=":", linewidth=3.0) + bdg.set_state(60.0e9, 300.0) + print(f"(Mg0.5Fe0.5)(Fe0.5Si0.5)O3 V at 60 GPa, 300 K: {bdg.V*1.e6:.3f} cm^3/mol") plt.show() def check_fig_a4_stishovite(): + print("\nChecking Figure A4...") stv = SLB_2024.stishovite() - fig = plt.figure(figsize=(10, 5)) + fig = plt.figure(figsize=(8, 5)) fig.suptitle("Figure A4 (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] @@ -369,11 +500,16 @@ def check_fig_a4_stishovite(): temperatures = pressures * 0.0 + T Vs = stv.evaluate(["shear_wave_velocity"], pressures, temperatures)[0] ax[0].plot(pressures / 1.0e9, Vs / 1.0e3, linestyle=":", linewidth=3.0) + + stv.set_state(1.0e10, 2000.0) + print(f"Stv Vs at 10 GPa, 2000 K: {stv.shear_wave_velocity/1.e3:.3f} km/s") plt.show() if __name__ == "__main__": + check_iron_bearing_mineral_property_values() check_bcc_iron_consistency() + check_fper_entropy() check_fig_1_fper_relaxed() check_fig_3_fcc_ferric_fper() check_fig_6a_iron_Cp_V() diff --git a/misc/ref/slb_2024_benchmarks.py.out b/misc/ref/slb_2024_benchmarks.py.out index 1b400a15..4a846b15 100644 --- a/misc/ref/slb_2024_benchmarks.py.out +++ b/misc/ref/slb_2024_benchmarks.py.out @@ -1,3 +1,17 @@ +Checking iron bearing mineral properties... +wu in agreement with HeFESTo (max rel error 0.0001 %) +mag in agreement with HeFESTo (max rel error 0.0001 %) +wu in agreement with HeFESTo (max rel error 0.0001 %) +mag in agreement with HeFESTo (max rel error 0.0001 %) +hem in agreement with HeFESTo (max rel error 0.0001 %) +hepv in agreement with HeFESTo (max rel error 0.0001 %) +hlpv in agreement with HeFESTo (max rel error 0.0001 %) +hppv in agreement with HeFESTo (max rel error 0.0001 %) +fea in agreement with HeFESTo (max rel error 0.0002 %) +feg in agreement with HeFESTo (max rel error 0.0006 %) +fee in agreement with HeFESTo (max rel error 0.0003 %) +hmag in agreement with HeFESTo (max rel error 0.0001 %) + Checking EoS consistency for 'burnman.minerals.SLB_2024.fea' Expressions within tolerance of 0.000100 G = F + PV : True @@ -15,13 +29,39 @@ Vphi = np.sqrt(K_S/rho) : True Vp = np.sqrt((K_S + 4G/3)/rho) : True Vs = np.sqrt(G_S/rho) : True All EoS consistency constraints satisfied for 'burnman.minerals.SLB_2024.fea' -The bcc-fcc-hcp invariant is at 11.250 GPa, 827.042 K + +Checking fper configurational entropy... +Endmember entropies (should all be zero): +[0. 0. 0. 0. 0.] +Equimolar wu-mag: 11.52629 J/K/mol (should be 11.52629 J/K/mol) + +Checking Figure 1... +Mg2Fe2O4 ferropericlase at 70 GPa, 300 K: +Molar proportions: 0.50, 0.24, 0.26, 0.00, 0.00 +Volume: 35.49 cm^3/mol + +Checking Figure 3... +BCC-fper-mag triple point at 1 bar: 825.98 K + +Checking Figure 6a... +BCC iron Cp, V at 1 bar, 1000 K: 52.29 J/K/mol, 7.30 cm^3/mol + +Checking Figure 6c... +The bcc-fcc-hcp invariant is at 11.27 GPa, 825.7 K + +Checking Figure 7... fO2 values at 1500 K: ['Fayalite', 'gamma_(fcc)_Iron', 'Quartz'] -log10fO2 = -12.063 +log10fO2 = -12.076 ['Fayalite', 'alpha_(bcc)_Iron', 'Quartz'] -log10fO2 = -12.063 +log10fO2 = -12.077 ['Fayalite', 'Magnetite', 'Quartz'] log10fO2 = -8.233 ['Magnetite', 'Hematite', 'Hematite'] log10fO2 = -2.635 + +Checking Figure A1... +(Mg0.5Fe0.5)(Fe0.5Si0.5)O3 V at 60 GPa, 300 K: 21.634 cm^3/mol + +Checking Figure A4... +Stv Vs at 10 GPa, 2000 K: 6.421 km/s diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index 223214f9..9476236b 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -810,6 +810,36 @@ def test_subregular_model_ternary_hessian_multicomponent_change(self): self.assertArraysAlmostEqual(H0.dot(df), dGdx2 - dGdx1) + def test_slb_ferropericlase_partial_gibbs_multicomponent_change(self): + ss = burnman.minerals.SLB_2024.ferropericlase() + f0 = [0.15, 0.25, 0.35, 0.05, 0.2] + ss.set_composition(f0) + ss.set_state(1.0e5, 300.0) + mu0 = ss.partial_gibbs + + df = np.array([2.0e-6, -1.5e-6, -0.5e-6, 0.25e-6, -0.25e-6]) + ss.set_composition(f0 - df / 2.0) + G1 = ss.gibbs + ss.set_composition(f0 + df / 2.0) + G2 = ss.gibbs + + self.assertAlmostEqual(mu0.dot(df), G2 - G1) + + def test_slb_ferropericlase_hessian_multicomponent_change(self): + ss = burnman.minerals.SLB_2024.ferropericlase() + f0 = [0.15, 0.25, 0.35, 0.05, 0.2] + ss.set_composition(f0) + ss.set_state(1.0e5, 300.0) + H0 = ss.gibbs_hessian + + df = np.array([2.0e-6, -1.5e-6, -0.5e-6, 0.25e-6, -0.25e-6]) + ss.set_composition(f0 - df / 2.0) + dGdx1 = ss.partial_gibbs + ss.set_composition(f0 + df / 2.0) + dGdx2 = ss.partial_gibbs + + self.assertArraysAlmostEqual(H0.dot(df), dGdx2 - dGdx1) + def test_polynomial_model_partial_entropy_multicomponent_change(self): ss = two_site_ss_polynomial_ternary() f0 = np.array([0.25, 0.35, 0.4])