Skip to content

Commit

Permalink
added benchmarks, stv softening
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed May 16, 2024
1 parent d518750 commit c263178
Show file tree
Hide file tree
Showing 6 changed files with 296 additions and 135 deletions.
103 changes: 96 additions & 7 deletions burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,45 @@ def rfloat(x, m=1.0):


# Hard-coded into HeFESTo
low_spin_endmembers = ["wuls", "hlpv", "lppv"]
"""
Ferropericlase:
Crystallography of wustites shows that the locations of the tetrahedral site
and the octahedral vacancies are highly correlated, with the tetrahedral
site being surrounded by a tetrahedral array of octahedral vacancies
(Gavarri et al. 1979 ).Because of this highly structured arrangement,
we ignore the contribution to the ideal entropy of mixing from the
tetrahedral site and the vacant octahedral sites.
...the formula of mag : [Fe]Fe2O4, with the square brackets indicating
that the first site does not participate in mixing. Thus, the ideal entropy
of mixing of an equimolar solution of wu and mag is 2Rln2 per 4 oxygens,
where R is the gas constant.
Spinel:
Pure magnetite is known to have an inverse arrangement of cations at
ambient conditions, rather than normal, as we have assumed.
Bridgmanite:
hepv and hlpv where only the cation on the
smaller B site undergoes the high spin to low spin transition...
"""

hardcoded_site_occupancies = {
"wu": "[Fe]2[Fe]2",
"wuls": "[Fels]2[Fels]2",
"mag": "[Vac1/2Fe1/2]2[Fef]2",
"smag": "[Fe][Fef][Fef]",
"hepv": "[Fef][Fef]",
"hlpv": "[Fef][Fefls]",
"fapv": "[Fef][Al]",
"hppv": "[Fef][Fef]",
"hmag": "[Fe][Fef][Fef]",
"acm": "[Na][Fef][Si]2",
"andr": "[Ca]3[Fef][Fef]",
"hem": "[Fef][Fef]",
"lppv": "[Fefls][Fefls]",
}

magnetic_structural_parameters = {"fea": 0.4, "mag": 0.4, "wu": 0.6}

hefesto_path = "HeFESTo_parameters_010123"
Expand Down Expand Up @@ -379,17 +417,66 @@ def rfloat(x, m=1.0):
'"""\n\n'
"from __future__ import absolute_import\n"
"\n\n"
"from ..eos.slb import SLB3\n"
"from ..classes.mineral import Mineral\n"
"from ..classes.solution import Solution\n"
"from ..classes.solution import Solution, RelaxedSolution\n"
"from ..classes.solutionmodel import (\n"
" IdealSolution,\n"
" SymmetricRegularSolution,\n"
" AsymmetricRegularSolution,\n"
")\n\n"
)


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"
" 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"
" PminusPc = pressure - Pc\n"
" if (PminusPc < 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"
"\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"
)

fper_relaxed_class = (
"class ferropericlase_relaxed(RelaxedSolution):\n"
" def __init__(self):\n"
' """RelaxedSolution model for ferropericlase (mw).\n'
" Endmembers (and site species distributions) are given in the order:\n"
" - 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"
" The equilibrium spin state is calculated automatically.\n"
' """\n'
" solution = ferropericlase()\n"
" vrel = [[0.0, 0.5, -0.5, 0.0, 0.0]]\n"
" vunrel = [\n"
" [1.0, 0.0, 0.0, 0.0, 0.0],\n"
" [0.0, 0.5, 0.5, 0.0, 0.0],\n"
" [0.0, 0.0, 0.0, 1.0, 0.0],\n"
" [0.0, 0.0, 0.0, 0.0, 1.0],\n"
" ]\n"
" RelaxedSolution.__init__(self, solution, vrel, vunrel)\n"
)

print('"""\n' "ENDMEMBERS\n" '"""\n')
for key, mbr_prm in sorted(mbr_params.items()):
if key == "st":
print(st_class)
mbr_prm = mbr_prm.__repr__().replace("'slb3'", "SLB3Stishovite()")
mbr_prm = mbr_prm.replace("slb3", "fault")
print(
f"\nclass {key}(Mineral):\n"
" def __init__(self):\n"
Expand All @@ -405,13 +492,12 @@ def rfloat(x, m=1.0):

print('"""\n' "SOLUTIONS\n" '"""\n')


for key, prm in sorted(sol_params.items()):
for i in range(prm["n_mbrs"]):
if prm["mbr_names"][i] in low_spin_endmembers:
prm["mbr_site_formulae"][i] = prm["mbr_site_formulae"][i].replace(
"Fe", "Fels"
)
if prm["mbr_names"][i] in hardcoded_site_occupancies:
prm["mbr_site_formulae"][i] = hardcoded_site_occupancies[
prm["mbr_names"][i]
]

docstring = '"""'
docstring += f'{prm["solution_type"]} model for {solution_aliases[key]} ({key}).\n'
Expand Down Expand Up @@ -445,6 +531,9 @@ def rfloat(x, m=1.0):
print(" )\n")
print(" Solution.__init__(self, molar_fractions=molar_fractions)\n")

if key == "mw":
print(fper_relaxed_class)


print('\n"""\n' "ENDMEMBER ALIASES\n" '"""\n')

Expand Down
90 changes: 68 additions & 22 deletions burnman/minerals/SLB_2024.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,9 @@
from __future__ import absolute_import


from ..eos.slb import SLB3
from ..classes.mineral import Mineral
from ..classes.solution import Solution
from ..classes.solution import Solution, RelaxedSolution
from ..classes.solutionmodel import (
IdealSolution,
SymmetricRegularSolution,
Expand Down Expand Up @@ -1867,12 +1868,36 @@ def __init__(self):
Mineral.__init__(self)


class SLB3Stishovite(SLB3):
def shear_modulus(self, pressure, temperature, volume, params):
G_bare = SLB3.shear_modulus(self, pressure, temperature, volume, params)

Pc = 51.6e9 + 11.1 * (temperature - 300.0) * 1.0e6 # Eq B4
PcQ = (
Pc + 50.7e9
) # From the HeFESTo code in stishtran.f. PcQ is T-independent in paper.
Cs0 = 128.0e9 * PcQ / Pc
Delta = 0.0 # 100.e9 in the HeFESTo code in stishtran.f, 0. in paper.
PminusPc = pressure - Pc
if PminusPc < 0.0:
Cs = Cs0 * PminusPc / (pressure - PcQ) # Eq B2
else:
Ast = (Cs0 - Delta) * (PcQ - Pc)
Cs = Cs0 - Ast / (PcQ - pressure + 3.0 * PminusPc) # Eq B3

# VRH-like average of the bare shear modulus with softened (C11-C12)/2.
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


class st(Mineral):
def __init__(self):
self.params = {
"name": "Stishovite",
"formula": {"Si": 1.0, "O": 2.0},
"equation_of_state": "slb3",
"equation_of_state": SLB3Stishovite(),
"F_0": -817124.2,
"V_0": 1.4017e-05,
"K_0": 305839630000.0,
Expand Down Expand Up @@ -2008,7 +2033,7 @@ def __init__(self, molar_fractions=None):
- mgcf ([Mg][Al][Al])
- fecf ([Fe][Al][Al])
- nacf ([Na][Al][Si])
- hmag ([Fe][Fe][Fe])
- hmag ([Fe][Fef][Fef])
- crcf ([Mg][Cr][Cr])
"""
self.name = "calcium_ferrite_structured_phase"
Expand All @@ -2017,7 +2042,7 @@ def __init__(self, molar_fractions=None):
[mgcf(), "[Mg][Al][Al]"],
[fecf(), "[Fe][Al][Al]"],
[nacf(), "[Na][Al][Si]"],
[hmag(), "[Fe][Fe][Fe]"],
[hmag(), "[Fe][Fef][Fef]"],
[crcf(), "[Mg][Cr][Cr]"],
],
alphas=[1.0, 1.0, 3.97647, 1.0, 1.0],
Expand All @@ -2041,7 +2066,7 @@ def __init__(self, molar_fractions=None):
- cen ([Mg][Mg][Si]2)
- cats ([Ca][Al][Si1/2Al1/2]2)
- jd ([Na][Al][Si]2)
- acm ([Na][Fe][Si]2)
- acm ([Na][Fef][Si]2)
"""
self.name = "clinopyroxene"
self.solution_model = AsymmetricRegularSolution(
Expand All @@ -2051,7 +2076,7 @@ def __init__(self, molar_fractions=None):
[cen(), "[Mg][Mg][Si]2"],
[cats(), "[Ca][Al][Si1/2Al1/2]2"],
[jd(), "[Na][Al][Si]2"],
[acm(), "[Na][Fe][Si]2"],
[acm(), "[Na][Fef][Si]2"],
],
alphas=[1.0, 1.0, 1.0, 3.5, 1.0, 1.0],
energy_interaction=[
Expand All @@ -2075,7 +2100,7 @@ def __init__(self, molar_fractions=None):
- gr ([Ca]3[Al][Al])
- mgmj ([Mg]3[Mg][Si])
- namj ([Na2/3Mg1/3]3[Si][Si])
- andr ([Ca]3[Fe][Fe])
- andr ([Ca]3[Fef][Fef])
- knor ([Mg]3[Cr][Cr])
"""
self.name = "garnet"
Expand All @@ -2086,7 +2111,7 @@ def __init__(self, molar_fractions=None):
[gr(), "[Ca]3[Al][Al]"],
[mgmj(), "[Mg]3[Mg][Si]"],
[namj(), "[Na2/3Mg1/3]3[Si][Si]"],
[andr(), "[Ca]3[Fe][Fe]"],
[andr(), "[Ca]3[Fef][Fef]"],
[knor(), "[Mg]3[Cr][Cr]"],
],
energy_interaction=[
Expand Down Expand Up @@ -2117,7 +2142,7 @@ def __init__(self, molar_fractions=None):
- mgil ([Mg][Si])
- feil ([Fe][Si])
- co ([Al][Al])
- hem ([Fe][Fe])
- hem ([Fef][Fef])
- esk ([Cr][Cr])
"""
self.name = "ilmenite"
Expand All @@ -2126,7 +2151,7 @@ def __init__(self, molar_fractions=None):
[mgil(), "[Mg][Si]"],
[feil(), "[Fe][Si]"],
[co(), "[Al][Al]"],
[hem(), "[Fe][Fe]"],
[hem(), "[Fef][Fef]"],
[esk(), "[Cr][Cr]"],
],
alphas=[1.0, 1.0, 1.0, 1.3, 1.0],
Expand All @@ -2149,7 +2174,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[Fe]2)
- mag ([Vac1/2Fe1/2]2[Fef]2)
"""
self.name = "ferropericlase"
self.solution_model = AsymmetricRegularSolution(
Expand All @@ -2158,7 +2183,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[Fe]2"],
[mag(), "[Vac1/2Fe1/2]2[Fef]2"],
],
alphas=[1.0, 1.0, 1.0, 1.0, 0.08293],
energy_interaction=[
Expand All @@ -2178,6 +2203,27 @@ def __init__(self, molar_fractions=None):
Solution.__init__(self, molar_fractions=molar_fractions)


class ferropericlase_relaxed(RelaxedSolution):
def __init__(self):
"""RelaxedSolution model for ferropericlase (mw).
Endmembers (and site species distributions) are given in the order:
- 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)
The equilibrium spin state is calculated automatically.
"""
solution = ferropericlase()
vrel = [[0.0, 0.5, -0.5, 0.0, 0.0]]
vunrel = [
[1.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.5, 0.5, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0],
]
RelaxedSolution.__init__(self, solution, vrel, vunrel)


class new_aluminous_phase(Solution):
def __init__(self, molar_fractions=None):
"""SymmetricRegularSolution model for new_aluminous_phase (nal).
Expand Down Expand Up @@ -2267,7 +2313,7 @@ def __init__(self, molar_fractions=None):
- mppv ([Mg][Si])
- fppv ([Fe][Si])
- appv ([Al][Al])
- hppv ([Fe][Fe])
- hppv ([Fef][Fef])
- cppv ([Cr][Cr])
"""
self.name = "post_perovskite"
Expand All @@ -2276,7 +2322,7 @@ def __init__(self, molar_fractions=None):
[mppv(), "[Mg][Si]"],
[fppv(), "[Fe][Si]"],
[appv(), "[Al][Al]"],
[hppv(), "[Fe][Fe]"],
[hppv(), "[Fef][Fef]"],
[cppv(), "[Cr][Cr]"],
],
energy_interaction=[
Expand All @@ -2297,9 +2343,9 @@ def __init__(self, molar_fractions=None):
- mgpv ([Mg][Si])
- fepv ([Fe][Si])
- alpv ([Al][Al])
- hepv ([Fe][Fe])
- hlpv ([Fels][Fels])
- fapv ([Fe][Al])
- hepv ([Fef][Fef])
- hlpv ([Fef][Fefls])
- fapv ([Fef][Al])
- crpv ([Cr][Cr])
"""
self.name = "bridgmanite"
Expand All @@ -2308,9 +2354,9 @@ def __init__(self, molar_fractions=None):
[mgpv(), "[Mg][Si]"],
[fepv(), "[Fe][Si]"],
[alpv(), "[Al][Al]"],
[hepv(), "[Fe][Fe]"],
[hlpv(), "[Fels][Fels]"],
[fapv(), "[Fe][Al]"],
[hepv(), "[Fef][Fef]"],
[hlpv(), "[Fef][Fefls]"],
[fapv(), "[Fef][Al]"],
[crpv(), "[Cr][Cr]"],
],
energy_interaction=[
Expand Down Expand Up @@ -2351,15 +2397,15 @@ def __init__(self, molar_fractions=None):
Endmembers (and site species distributions) are given in the order:
- sp ([Mg][Al][Al])
- hc ([Fe][Al][Al])
- smag ([Fe][Fe][Fe])
- smag ([Fe][Fef][Fef])
- picr ([Mg][Cr][Cr])
"""
self.name = "mg_fe_aluminous_spinel"
self.solution_model = SymmetricRegularSolution(
endmembers=[
[sp(), "[Mg][Al][Al]"],
[hc(), "[Fe][Al][Al]"],
[smag(), "[Fe][Fe][Fe]"],
[smag(), "[Fe][Fef][Fef]"],
[picr(), "[Mg][Cr][Cr]"],
],
energy_interaction=[
Expand Down
Binary file added misc/benchmarks/figures/SLB_2024_fper_K_T.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added misc/benchmarks/figures/SLB_2024_fper_V.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added misc/benchmarks/figures/SLB_2024_stv_VS.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading

0 comments on commit c263178

Please sign in to comment.