Skip to content

Commit 70fda93

Browse files
bobmyhillbobmyhill
committed
added example
1 parent df720a9 commit 70fda93

File tree

3 files changed

+392
-1
lines changed

3 files changed

+392
-1
lines changed

‎examples/example_modular_mgd.py‎

Lines changed: 383 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,383 @@
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
10+
-------------------
11+
12+
This example shows how to create a mineral object using
13+
the Modular Mie-Grüneisen-Debye equation of state.
14+
15+
This equation of state allows users to pick and choose
16+
different options for two key aspects of the mineral's behavior:
17+
18+
1. The equation of state (EOS) used to describe the mineral's thermodynamic properties at the reference temperature.
19+
2. The model used to describe the mineral's Debye temperature as a function of volume.
20+
21+
and a further two options if the user wants to account for anharmonic effects:
22+
23+
1. The model used to describe the mineral's anharmonicity as a function of temperature.
24+
2. The model used to describe the prefactor as a function of volume.
25+
26+
*Uses:*
27+
28+
* :class:`burnman.Mineral`
29+
30+
31+
*Demonstrates:*
32+
33+
* How to create a mineral object using the Modular Mie-Grüneisen-Debye equation of state.
34+
35+
"""
36+
# First, a couple of standard imports
37+
import numpy as np
38+
import matplotlib.pyplot as plt
39+
40+
# Now some specific BurnMan imports
41+
# The Mineral class is the standard class for creating mineral objects
42+
from burnman import Mineral
43+
44+
# For this example, we will use the SLB_2011 fayalite mineral
45+
# as a starting point
46+
from burnman.minerals.SLB_2011 import fayalite as mineral
47+
48+
# This helper function creates an Equation of State object
49+
from burnman.eos.helper import create
50+
51+
# These imports contain the different models available for the EOS
52+
from burnman.eos import debye_temperature_models
53+
from burnman.eos import anharmonic_prefactor_models
54+
from burnman.eos import anharmonic_thermal_models
55+
56+
# Create the SLB mineral object
57+
m_SLB = mineral()
58+
59+
# Now create a new mineral object using the modular MGD EOS
60+
params = m_SLB.params.copy()
61+
params["equation_of_state"] = "modular_mgd"
62+
63+
# In these two lines, we set the reference EOS and Debye temperature model
64+
# to be an exact clone of the original SLB object
65+
params["reference_eos"] = create("bm3")
66+
params["debye_temperature_model"] = debye_temperature_models.SLB()
67+
68+
# We are also able to add an anharmonicity model
69+
# that is not a part of the SLB model. This model has two parts.
70+
# The first is a volume-dependent prefactor model,
71+
# which is here chosen to follow a power law with parameters
72+
# a_anh and m_anh: i.e, A(V) = a_anh * (V/V_0)^m_anh
73+
params["anharmonic_prefactor_model"] = anharmonic_prefactor_models.PowerLaw()
74+
params["a_anh"] = 1.0
75+
params["m_anh"] = 10.0
76+
77+
# The second part is a temperature-dependent anharmonic model,
78+
# which is here chosen so that the modification to the isochoric heat capacity
79+
# is based on the CDF log-normal distribution with parameters mu_anh and sigma_anh.
80+
params["anharmonic_thermal_model"] = anharmonic_thermal_models.LogNormal()
81+
params["mu_anh"] = 0.0
82+
params["sigma_anh"] = 1.0
83+
84+
# Create the new mineral object with the property modifiers from the SLB model
85+
m = Mineral(params, property_modifiers=m_SLB.property_modifiers)
86+
87+
# And we're done! Now we can find the properties of the new mineral object
88+
# using set_state or evaluate.
89+
m.set_state(5.0e9, 4000.0)
90+
91+
# Print some properties to standard output
92+
print(
93+
"Properties of the new mineral object at "
94+
f"{m.pressure/1.e9:.1f} GPa and {m.temperature:.1f} K:"
95+
)
96+
print(f"Volume: {m.V * 1.e6:.1f} cm^3/mol")
97+
print(f"Entropy: {m.S:.1f} J/K/mol")
98+
print(f"Isochoric heat capacity: {m.molar_heat_capacity_v:.1f} J/K/mol")
99+
print(f"Isobaric heat capacity: {m.molar_heat_capacity_p:.1f} J/K/mol")
100+
print(f"Thermal expansion coefficient: {m.alpha*1.e6:.1f} * 1/MK")
101+
print(f"Grueneisen parameter: {m.gr:.2f}")
102+
print(f"Isothermal bulk modulus: {m.K_T / 1.e9:.1f} GPa")
103+
print(f"Isentropic bulk modulus: {m.K_S / 1.e9:.1f} GPa")
104+
105+
# Create the figure and axes for the plots
106+
fig = plt.figure(figsize=(12, 6))
107+
ax = [fig.add_subplot(2, 4, i) for i in range(1, 8)]
108+
109+
# Create plots for the mineral properties as a function of temperature
110+
# at different pressures.
111+
temperatures = np.linspace(10.0, 5000.0, 101)
112+
for P in np.linspace(0.0, 30.0e9, 7):
113+
pressures = temperatures * 0.0 + P
114+
try:
115+
volumes_SLB, Cv_SLB, Cp_SLB, alpha_SLB, gr_SLB, K_T_SLB = m_SLB.evaluate(
116+
[
117+
"V",
118+
"molar_heat_capacity_v",
119+
"molar_heat_capacity_p",
120+
"alpha",
121+
"gr",
122+
"K_T",
123+
],
124+
pressures,
125+
temperatures,
126+
)
127+
ln = ax[0].plot(temperatures, volumes_SLB / m_SLB.params["V_0"], linestyle=":")
128+
ax[1].plot(temperatures, Cv_SLB, linestyle=":", color=ln[0].get_color())
129+
ax[2].plot(temperatures, Cp_SLB, linestyle=":", color=ln[0].get_color())
130+
ax[3].plot(temperatures, alpha_SLB, linestyle=":", color=ln[0].get_color())
131+
ax[4].plot(temperatures, gr_SLB, linestyle=":", color=ln[0].get_color())
132+
ax[5].plot(
133+
temperatures, K_T_SLB / 1.0e9, linestyle=":", color=ln[0].get_color()
134+
)
135+
ax[6].plot(
136+
temperatures,
137+
alpha_SLB * K_T_SLB / 1.0e6,
138+
linestyle=":",
139+
color=ln[0].get_color(),
140+
)
141+
except Exception:
142+
temperatures_lower = np.linspace(10.0, 2695.0 + P / 7.0e6, 101)
143+
volumes_SLB, Cv_SLB, Cp_SLB, alpha_SLB, gr_SLB, K_T_SLB = m_SLB.evaluate(
144+
[
145+
"V",
146+
"molar_heat_capacity_v",
147+
"molar_heat_capacity_p",
148+
"alpha",
149+
"gr",
150+
"K_T",
151+
],
152+
pressures,
153+
temperatures_lower,
154+
)
155+
ln = ax[0].plot(
156+
temperatures_lower, volumes_SLB / m_SLB.params["V_0"], linestyle=":"
157+
)
158+
ax[1].plot(temperatures_lower, Cv_SLB, linestyle=":", color=ln[0].get_color())
159+
ax[2].plot(temperatures_lower, Cp_SLB, linestyle=":", color=ln[0].get_color())
160+
ax[3].plot(
161+
temperatures_lower, alpha_SLB, linestyle=":", color=ln[0].get_color()
162+
)
163+
ax[4].plot(temperatures_lower, gr_SLB, linestyle=":", color=ln[0].get_color())
164+
ax[5].plot(
165+
temperatures_lower, K_T_SLB / 1.0e9, linestyle=":", color=ln[0].get_color()
166+
)
167+
ax[6].plot(
168+
temperatures_lower,
169+
alpha_SLB * K_T_SLB / 1.0e6,
170+
linestyle=":",
171+
color=ln[0].get_color(),
172+
)
173+
try:
174+
volumes, Cv, Cp, alpha, gr, K_T = m.evaluate(
175+
[
176+
"V",
177+
"molar_heat_capacity_v",
178+
"molar_heat_capacity_p",
179+
"alpha",
180+
"gr",
181+
"K_T",
182+
],
183+
pressures,
184+
temperatures,
185+
)
186+
ax[0].plot(temperatures, volumes / m.params["V_0"], c=ln[0].get_color())
187+
ax[1].plot(temperatures, Cv, c=ln[0].get_color(), label=f"{P/1.e9:.0f} GPa")
188+
ax[2].plot(temperatures, Cp, c=ln[0].get_color())
189+
ax[3].plot(temperatures, alpha, c=ln[0].get_color())
190+
ax[4].plot(temperatures, gr, c=ln[0].get_color())
191+
ax[5].plot(temperatures, K_T / 1.0e9, c=ln[0].get_color())
192+
ax[6].plot(temperatures, alpha * K_T / 1.0e6, c=ln[0].get_color())
193+
except Exception:
194+
temperatures_lower = np.linspace(10.0, 2695.0 + P / 1.0e7, 101)
195+
volumes, Cv, Cp, alpha, gr, K_T = m.evaluate(
196+
[
197+
"V",
198+
"molar_heat_capacity_v",
199+
"molar_heat_capacity_p",
200+
"alpha",
201+
"gr",
202+
"K_T",
203+
],
204+
pressures,
205+
temperatures_lower,
206+
)
207+
ax[0].plot(temperatures_lower, volumes / m.params["V_0"], c=ln[0].get_color())
208+
ax[1].plot(
209+
temperatures_lower, Cv, c=ln[0].get_color(), label=f"{P/1.e9:.0f} GPa"
210+
)
211+
ax[2].plot(temperatures_lower, Cp, c=ln[0].get_color())
212+
ax[3].plot(temperatures_lower, alpha, c=ln[0].get_color())
213+
ax[4].plot(temperatures_lower, gr, c=ln[0].get_color())
214+
ax[5].plot(temperatures_lower, K_T / 1.0e9, c=ln[0].get_color())
215+
ax[6].plot(temperatures_lower, alpha * K_T / 1.0e6, c=ln[0].get_color())
216+
217+
ax[1].legend()
218+
219+
for i in range(6):
220+
ax[i].set_xlabel("Temperature (K)")
221+
ax[i].set_xlim(0.0, 5000.0)
222+
ax[1].set_ylim(0.0, 180.0)
223+
ax[2].set_ylim(0.0, 250.0)
224+
ax[3].set_ylim(0.0, 3.0e-4)
225+
226+
ax[0].set_ylabel("$V/V_0$ (cm³/mol)")
227+
ax[1].set_ylabel("$C_V$ (J/mol/K)")
228+
ax[2].set_ylabel("$C_P$ (J/mol/K)")
229+
ax[3].set_ylabel("$\\alpha$ (1/K)")
230+
ax[4].set_ylabel("$\\gamma$")
231+
ax[5].set_ylabel("$K_T$ (GPa)")
232+
ax[6].set_ylabel("$\\alpha K_T$ (MPa/K)")
233+
234+
fig.set_layout_engine("tight")
235+
plt.show()
236+
237+
238+
# Now do exactly the same thing again, but using evaluate_with_volumes
239+
fig = plt.figure(figsize=(12, 6))
240+
ax = [fig.add_subplot(2, 4, i) for i in range(1, 8)]
241+
242+
temperatures = np.linspace(10.0, 5000.0, 101)
243+
for Vrels in np.linspace(1.2, 0.8, 5):
244+
volumes = temperatures * 0.0 + Vrels * m_SLB.params["V_0"]
245+
246+
try:
247+
pressures_SLB, Cv_SLB, Cp_SLB, alpha_SLB, gr_SLB, K_T_SLB = (
248+
m_SLB.evaluate_with_volumes(
249+
[
250+
"pressure",
251+
"molar_heat_capacity_v",
252+
"molar_heat_capacity_p",
253+
"alpha",
254+
"gr",
255+
"K_T",
256+
],
257+
volumes,
258+
temperatures,
259+
)
260+
)
261+
ln = ax[0].plot(temperatures, pressures_SLB / 1.0e9, linestyle=":")
262+
ax[1].plot(temperatures, Cv_SLB, linestyle=":", color=ln[0].get_color())
263+
ax[2].plot(temperatures, Cp_SLB, linestyle=":", color=ln[0].get_color())
264+
ax[3].plot(temperatures, alpha_SLB, linestyle=":", color=ln[0].get_color())
265+
ax[4].plot(temperatures, gr_SLB, linestyle=":", color=ln[0].get_color())
266+
ax[5].plot(
267+
temperatures, K_T_SLB / 1.0e9, linestyle=":", color=ln[0].get_color()
268+
)
269+
ax[6].plot(
270+
temperatures,
271+
alpha_SLB * K_T_SLB / 1.0e6,
272+
linestyle=":",
273+
color=ln[0].get_color(),
274+
)
275+
except Exception:
276+
temperatures_lower = np.linspace(10.0, 2695.0 + (1.0 - Vrels) * 1000.0, 101)
277+
pressures_SLB, Cv_SLB, Cp_SLB, alpha_SLB, gr_SLB, K_T_SLB = (
278+
m_SLB.evaluate_with_volumes(
279+
[
280+
"pressure",
281+
"molar_heat_capacity_v",
282+
"molar_heat_capacity_p",
283+
"alpha",
284+
"gr",
285+
"K_T",
286+
],
287+
volumes,
288+
temperatures_lower,
289+
)
290+
)
291+
ln = ax[0].plot(temperatures_lower, pressures_SLB / 1.0e9, linestyle=":")
292+
ax[1].plot(temperatures_lower, Cv_SLB, linestyle=":", color=ln[0].get_color())
293+
ax[2].plot(temperatures_lower, Cp_SLB, linestyle=":", color=ln[0].get_color())
294+
ax[3].plot(
295+
temperatures_lower, alpha_SLB, linestyle=":", color=ln[0].get_color()
296+
)
297+
ax[4].plot(temperatures_lower, gr_SLB, linestyle=":", color=ln[0].get_color())
298+
ax[5].plot(
299+
temperatures_lower, K_T_SLB / 1.0e9, linestyle=":", color=ln[0].get_color()
300+
)
301+
ax[6].plot(
302+
temperatures_lower,
303+
alpha_SLB * K_T_SLB / 1.0e6,
304+
linestyle=":",
305+
color=ln[0].get_color(),
306+
)
307+
try:
308+
pressures, Cv, Cp, alpha, gr, K_T = m.evaluate_with_volumes(
309+
[
310+
"pressure",
311+
"molar_heat_capacity_v",
312+
"molar_heat_capacity_p",
313+
"alpha",
314+
"gr",
315+
"K_T",
316+
],
317+
volumes,
318+
temperatures,
319+
)
320+
ax[0].plot(temperatures, pressures / 1.0e9, c=ln[0].get_color(), linestyle="--")
321+
ax[1].plot(
322+
temperatures,
323+
Cv,
324+
c=ln[0].get_color(),
325+
linestyle="--",
326+
label=f"$V_{{rel}}$ = {Vrels:.1f}",
327+
)
328+
ax[2].plot(temperatures, Cp, c=ln[0].get_color(), linestyle="--")
329+
ax[3].plot(temperatures, alpha, c=ln[0].get_color(), linestyle="--")
330+
ax[4].plot(temperatures, gr, c=ln[0].get_color(), linestyle="--")
331+
ax[5].plot(temperatures, K_T / 1.0e9, c=ln[0].get_color(), linestyle="--")
332+
ax[6].plot(
333+
temperatures, alpha * K_T / 1.0e6, c=ln[0].get_color(), linestyle="--"
334+
)
335+
except Exception:
336+
temperatures_lower = np.linspace(10.0, 2695.0 + (1.0 - Vrels) * 1000.0, 101)
337+
pressures, Cv, Cp, alpha, gr, K_T = m.evaluate_with_volumes(
338+
[
339+
"pressure",
340+
"molar_heat_capacity_v",
341+
"molar_heat_capacity_p",
342+
"alpha",
343+
"gr",
344+
"K_T",
345+
],
346+
volumes,
347+
temperatures_lower,
348+
)
349+
ax[0].plot(
350+
temperatures_lower, pressures / 1.0e9, c=ln[0].get_color(), linestyle="--"
351+
)
352+
ax[1].plot(
353+
temperatures_lower,
354+
Cv,
355+
c=ln[0].get_color(),
356+
linestyle="--",
357+
label=f"$V_{{rel}}$ = {Vrels:.1f}",
358+
)
359+
ax[2].plot(temperatures_lower, Cp, c=ln[0].get_color(), linestyle="--")
360+
ax[3].plot(temperatures_lower, alpha, c=ln[0].get_color(), linestyle="--")
361+
ax[4].plot(temperatures_lower, gr, c=ln[0].get_color(), linestyle="--")
362+
ax[5].plot(temperatures_lower, K_T / 1.0e9, c=ln[0].get_color(), linestyle="--")
363+
ax[6].plot(
364+
temperatures_lower, alpha * K_T / 1.0e6, c=ln[0].get_color(), linestyle="--"
365+
)
366+
367+
ax[1].legend()
368+
369+
for i in range(6):
370+
ax[i].set_xlabel("Temperature (K)")
371+
ax[i].set_xlim(0.0, 5000.0)
372+
ax[2].set_ylim(0.0, 400.0)
373+
ax[3].set_ylim(0.0, 3.0e-4)
374+
375+
ax[0].set_ylabel("$P$ (GPa)")
376+
ax[1].set_ylabel("$C_V$ (J/mol/K)")
377+
ax[2].set_ylabel("$C_P$ (J/mol/K)")
378+
ax[3].set_ylabel("$\\alpha$ (1/K)")
379+
ax[4].set_ylabel("$\\gamma$")
380+
ax[5].set_ylabel("$K_T$ (GPa)")
381+
ax[6].set_ylabel("$\\alpha K_T$ (MPa/K)")
382+
fig.set_layout_engine("tight")
383+
plt.show()

‎misc/ref/example_modular_mgd.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 mineral object at 5.0 GPa and 4000.0 K:
2+
Volume: 51.8 cm^3/mol
3+
Entropy: 605.2 J/K/mol
4+
Isochoric heat capacity: 171.4 J/K/mol
5+
Isobaric heat capacity: 247.8 J/K/mol
6+
Thermal expansion coefficient: 87.3 * 1/MK
7+
Grueneisen parameter: 1.28
8+
Isothermal bulk modulus: 48.4 GPa
9+
Isentropic bulk modulus: 69.9 GPa

0 commit comments

Comments
 (0)