Skip to content

Commit

Permalink
more benchmarking
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed May 25, 2024
1 parent 8feba02 commit c73a4e0
Show file tree
Hide file tree
Showing 11 changed files with 375 additions and 81 deletions.
10 changes: 8 additions & 2 deletions burnman/classes/solutionmodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
52 changes: 37 additions & 15 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 @@ -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]",
Expand All @@ -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"
Expand Down Expand Up @@ -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"
Expand All @@ -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 = (
Expand All @@ -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"
Expand All @@ -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"
Expand Down Expand Up @@ -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(
Expand Down
33 changes: 26 additions & 7 deletions burnman/eos/bukowinski_electronic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
31 changes: 26 additions & 5 deletions burnman/eos/slb.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down
43 changes: 25 additions & 18 deletions burnman/minerals/SLB_2024.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -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
Expand All @@ -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):
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand All @@ -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=[
Expand All @@ -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()
Expand Down Expand Up @@ -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.
Expand Down
Binary file added misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.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_fO2.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.
13 changes: 13 additions & 0 deletions misc/benchmarks/slb_2024_Fe_endmembers.dat
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit c73a4e0

Please sign in to comment.