Skip to content

Commit 9074046

Browse files
bobmyhillbobmyhill
committed
new anharmonicity example
1 parent ed383fb commit 9074046

File tree

2 files changed

+200
-0
lines changed

2 files changed

+200
-0
lines changed

‎examples/example_modular_mgd_with_anharmonicity.py‎

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
2+
# the Earth and Planetary Sciences
3+
# Copyright (C) 2012 - 2025 by the BurnMan team, released under the GNU
4+
# GPL v2 or later.
5+
6+
7+
"""
8+
9+
example_modular_mgd_with_anharmonicity
10+
--------------------------------------
11+
12+
This example shows how to create a mineral object using
13+
the Modular Mie-Grüneisen-Debye With Anharmonicity equation of state.
14+
15+
This equation of state is derived from the Modular Mie-Grüneisen-Debye equation of state,
16+
whose use was demonstrated in example_modular_mgd.py. Like that equation of state,
17+
the modified version demonstrated here allows users to pick and choose
18+
different options for two key aspects of the mineral's behavior:
19+
20+
1. The equation of state (EOS) used to describe the mineral's thermodynamic properties at the reference temperature.
21+
2. The model used to describe the mineral's Debye temperature as a function of volume.
22+
23+
In the Modular Mie-Grüneisen-Debye equation of state, anharmonicity is incorporated
24+
via extra terms in the Helmholtz energy expression; in other words, the model is
25+
somewhat removed from the underlying quasiharmonic form
26+
(although it shares the same volume dependence on the reference Debye temperature).
27+
In the modified equation of state described in this example, anharmonic effects are
28+
more tightly coupled to the basic quasiharmonic framework, and use the thermodynamic
29+
properties of that framework to describe the mineral's anharmonic behavior. The
30+
model used here is based on the empirical ansatz introduced by Wu and Wentzcovitch
31+
(2009; https://doi.org/10.1103/PhysRevB.79.104304), simplified somewhat to avoid
32+
in internal V(P(V, T), T=0) inversion to calculate the anharmonic contributions.
33+
This has the dual benefits of reducing computational complexity and improving
34+
the shape of the thermodynamic property functions as a function of temperature,
35+
especially at high volumes. The simplification still requires only one parameter,
36+
`c_anh`, which controls the strength of the anharmonic effects, but the value of that
37+
parameter may need to be adjusted relative to the model proposed by Wu and Wentzcovitch.
38+
39+
*Uses:*
40+
41+
* :class:`burnman.Mineral`
42+
43+
44+
*Demonstrates:*
45+
46+
* How to create a mineral object using the Modular Mie-Grüneisen-Debye equation of state.
47+
48+
"""
49+
# First, a couple of standard imports
50+
import numpy as np
51+
import matplotlib.pyplot as plt
52+
53+
# Now some specific BurnMan imports
54+
# The Mineral class is the standard class for creating mineral objects
55+
from burnman import Mineral
56+
57+
# For this example, we will use the SLB_2011 periclase mineral
58+
# as a starting point
59+
from burnman.minerals.SLB_2011 import periclase as mineral
60+
61+
# This helper function creates an Equation of State object
62+
from burnman.eos.helper import create
63+
64+
# These imports contain the different models available for the EOS
65+
from burnman.eos import debye_temperature_models
66+
67+
# Create the SLB mineral object
68+
m_SLB = mineral()
69+
70+
# Now create a new mineral object using the modular MGD with
71+
# anharmonicity equation of state.
72+
params = m_SLB.params.copy()
73+
params["equation_of_state"] = "modular_mgd_with_anharmonicity"
74+
75+
# In these two lines, we set the reference EOS and Debye temperature model
76+
# to be an exact clone of the original SLB object
77+
params["reference_eos"] = create("bm3")
78+
params["debye_temperature_model"] = debye_temperature_models.SLB()
79+
80+
# Now we just need to set the anharmonicity parameter
81+
# Positive values reduce the isochoric heat capacity, as expected from
82+
# classic anharmonicity theory.
83+
# Wu and Wentzcovitch (2009) proposed a value of c=0.1 for periclase,
84+
# but as mentioned above, the model implemented here is somewhat simplified
85+
# and c_anh needs to be adjusted downwards to avoid excessive depression of
86+
# the isochoric heat capacity, Grueneisen parameter, and bulk modulus.
87+
# Note also that the quasiharmonic reference state is somewhat different to
88+
# the one used by Wu and Wentzcovitch (2009); in particular, the thermal
89+
# expansion in the SLB model is significantly larger at high temperatures:
90+
# at 0 GPa and 3000 K alpha ~ 8e-5 / K in Wu and Wentzcovitch (2009),
91+
# but in the SLB model, alpha ~ 21.6e-5 / K at the same conditions.
92+
params["c_anh"] = 0.04
93+
94+
# Create the new mineral object with the property modifiers from the SLB model
95+
m = Mineral(params, property_modifiers=m_SLB.property_modifiers)
96+
97+
# And we're done! Now we can find the properties of the new mineral object
98+
# using set_state or evaluate.
99+
m.set_state(5.0e9, 4000.0)
100+
101+
# Print some properties to standard output
102+
print(
103+
f"Properties of the new {m.name} object at "
104+
f"{m.pressure/1.e9:.1f} GPa and {m.temperature:.1f} K:"
105+
)
106+
print(f"Volume: {m.V * 1.e6:.1f} cm^3/mol")
107+
print(f"Entropy: {m.S:.1f} J/K/mol")
108+
print(f"Isochoric heat capacity: {m.molar_heat_capacity_v:.1f} J/K/mol")
109+
print(f"Isobaric heat capacity: {m.molar_heat_capacity_p:.1f} J/K/mol")
110+
print(f"Thermal expansion coefficient: {m.alpha*1.e6:.1f} * 1/MK")
111+
print(f"Grueneisen parameter: {m.gr:.2f}")
112+
print(f"Isothermal bulk modulus: {m.K_T / 1.e9:.1f} GPa")
113+
print(f"Isentropic bulk modulus: {m.K_S / 1.e9:.1f} GPa")
114+
115+
# Create the figure and axes for the plots
116+
fig = plt.figure(figsize=(12, 6))
117+
ax = [fig.add_subplot(2, 4, i) for i in range(1, 8)]
118+
119+
# Create plots for the mineral properties as a function of temperature
120+
# at different pressures.
121+
temperatures = np.linspace(10.0, 3000.0, 101)
122+
for P_GPa in [0, 10, 30, 60, 100]:
123+
P = P_GPa * 1.0e9
124+
pressures = temperatures * 0.0 + P
125+
126+
volumes_SLB, Cv_SLB, Cp_SLB, alpha_SLB, gr_SLB, K_T_SLB = m_SLB.evaluate(
127+
[
128+
"V",
129+
"molar_heat_capacity_v",
130+
"molar_heat_capacity_p",
131+
"alpha",
132+
"gr",
133+
"K_T",
134+
],
135+
pressures,
136+
temperatures,
137+
)
138+
ln = ax[0].plot(temperatures, volumes_SLB / m_SLB.params["V_0"], linestyle=":")
139+
ax[1].plot(temperatures, Cv_SLB, linestyle=":", color=ln[0].get_color())
140+
ax[2].plot(temperatures, Cp_SLB, linestyle=":", color=ln[0].get_color())
141+
ax[3].plot(temperatures, alpha_SLB * 1.0e5, linestyle=":", color=ln[0].get_color())
142+
ax[4].plot(temperatures, gr_SLB, linestyle=":", color=ln[0].get_color())
143+
ax[5].plot(temperatures, K_T_SLB / 1.0e9, linestyle=":", color=ln[0].get_color())
144+
ax[6].plot(
145+
temperatures,
146+
alpha_SLB * K_T_SLB / 1.0e6,
147+
linestyle=":",
148+
color=ln[0].get_color(),
149+
)
150+
151+
volumes, Cv, Cp, alpha, gr, K_T = m.evaluate(
152+
[
153+
"V",
154+
"molar_heat_capacity_v",
155+
"molar_heat_capacity_p",
156+
"alpha",
157+
"gr",
158+
"K_T",
159+
],
160+
pressures,
161+
temperatures,
162+
)
163+
ax[0].plot(temperatures, volumes / m.params["V_0"], c=ln[0].get_color())
164+
ax[1].plot(temperatures, Cv, c=ln[0].get_color(), label=f"{P/1.e9:.0f} GPa")
165+
ax[2].plot(temperatures, Cp, c=ln[0].get_color())
166+
ax[3].plot(temperatures, alpha * 1.0e5, c=ln[0].get_color())
167+
ax[4].plot(temperatures, gr, c=ln[0].get_color())
168+
ax[5].plot(temperatures, K_T / 1.0e9, c=ln[0].get_color())
169+
ax[6].plot(temperatures, alpha * K_T / 1.0e6, c=ln[0].get_color())
170+
171+
ax[1].legend()
172+
173+
for i in range(6):
174+
ax[i].set_xlabel("Temperature (K)")
175+
ax[i].set_xlim(0.0, 3000.0)
176+
ax[1].set_ylim(0.0, 80.0)
177+
ax[2].set_ylim(0.0, 80.0)
178+
ax[3].set_ylim(0.0, 10.0)
179+
ax[4].set_ylim(1.0, 2.0)
180+
ax[5].set_ylim(100.0, 300.0)
181+
182+
ax[0].set_ylabel("$V/V_0$ (cm³/mol)")
183+
ax[1].set_ylabel("$C_V$ (J/mol/K)")
184+
ax[2].set_ylabel("$C_P$ (J/mol/K)")
185+
ax[3].set_ylabel("$\\alpha$ (10$^{{-5}}$/K)")
186+
ax[4].set_ylabel("$\\gamma$")
187+
ax[5].set_ylabel("$K_T$ (GPa)")
188+
ax[6].set_ylabel("$\\alpha K_T$ (MPa/K)")
189+
190+
fig.set_layout_engine("tight")
191+
plt.show()

‎misc/ref/example_modular_mgd_with_anharmonicity.py.out‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Properties of the new Periclase object at 5.0 GPa and 4000.0 K:
2+
Volume: 13.1 cm^3/mol
3+
Entropy: 158.8 J/K/mol
4+
Isochoric heat capacity: 47.3 J/K/mol
5+
Isobaric heat capacity: 71.7 J/K/mol
6+
Thermal expansion coefficient: 82.3 * 1/MK
7+
Grueneisen parameter: 1.56
8+
Isothermal bulk modulus: 68.5 GPa
9+
Isentropic bulk modulus: 103.6 GPa

0 commit comments

Comments
 (0)