Skip to content

Commit 5c8b8ff

Browse files
bobmyhillbobmyhill
committed
added new VTp fitting example
1 parent cf45710 commit 5c8b8ff

File tree

2 files changed

+272
-0
lines changed

2 files changed

+272
-0
lines changed

‎examples/example_fit_eos_with_independent_volumes.py‎

Lines changed: 235 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,235 @@
1+
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
2+
# the Earth and Planetary Sciences
3+
# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU
4+
# GPL v2 or later.
5+
6+
7+
"""
8+
example_fit_eos_with_independent_volumes
9+
----------------------------------------
10+
11+
This example demonstrates BurnMan's functionality to fit data to
12+
an EoS of the user's choice using volumes as independent variables.
13+
14+
The first example deals with simple P(V, T) fitting.
15+
The first example deals with E(V, T) fitting.
16+
17+
teaches:
18+
- least squares fitting
19+
20+
"""
21+
import numpy as np
22+
import matplotlib.pyplot as plt
23+
24+
import burnman
25+
26+
from burnman.utils.unitcell import molar_volume_from_unit_cell_volume
27+
from burnman.utils.misc import attribute_function
28+
from burnman.utils.misc import pretty_print_values
29+
30+
if __name__ == "__main__":
31+
print(
32+
"Least squares equation of state fitting with volume as "
33+
"an independent variable.\n"
34+
)
35+
36+
print("1) Fitting to stishovite room temperature P(V) data\n")
37+
# Let's fit an EoS to stishovite data from Andrault et al. (2003)
38+
PV = np.array(
39+
[
40+
[0.0001, 46.5126, 0.0061],
41+
[1.168, 46.3429, 0.0053],
42+
[2.299, 46.1756, 0.0043],
43+
[3.137, 46.0550, 0.0051],
44+
[4.252, 45.8969, 0.0045],
45+
[5.037, 45.7902, 0.0053],
46+
[5.851, 45.6721, 0.0038],
47+
[6.613, 45.5715, 0.0050],
48+
[7.504, 45.4536, 0.0041],
49+
[8.264, 45.3609, 0.0056],
50+
[9.635, 45.1885, 0.0042],
51+
[11.69, 44.947, 0.002],
52+
[17.67, 44.264, 0.002],
53+
[22.38, 43.776, 0.003],
54+
[29.38, 43.073, 0.009],
55+
[37.71, 42.278, 0.008],
56+
[46.03, 41.544, 0.017],
57+
[52.73, 40.999, 0.009],
58+
[26.32, 43.164, 0.006],
59+
[30.98, 42.772, 0.005],
60+
[34.21, 42.407, 0.003],
61+
[38.45, 42.093, 0.004],
62+
[43.37, 41.610, 0.004],
63+
[47.49, 41.280, 0.007],
64+
]
65+
)
66+
67+
Z = 2.0 # number of formula units per unit cell in stishovite
68+
VTP = np.array(
69+
[
70+
molar_volume_from_unit_cell_volume(PV[:, 1], Z),
71+
298.15 * np.ones_like(PV[:, 0]),
72+
PV[:, 0] * 1.0e9,
73+
]
74+
).T
75+
76+
# Here, we assume that the pressure uncertainties are equal to 3% of the total pressure,
77+
# that the temperature uncertainties are negligible, and take the unit cell volume
78+
# uncertainties from the paper.
79+
# We also assume that the uncertainties in pressure and volume are uncorrelated.
80+
nul = np.zeros_like(VTP[:, 0])
81+
VTP_covariances = np.array(
82+
[
83+
[molar_volume_from_unit_cell_volume(PV[:, 2], Z), nul, nul],
84+
[nul, nul, nul],
85+
[nul, nul, 0.03 * VTP[:, 2]],
86+
]
87+
).T
88+
VTP_covariances = np.power(VTP_covariances, 2.0)
89+
90+
# Here's where we fit the data
91+
# The mineral parameters are automatically updated during fitting
92+
stv = burnman.minerals.HP_2011_ds62.stv()
93+
params = ["V_0", "K_0", "Kprime_0"]
94+
fitted_eos = burnman.eos_fitting.fit_VTP_data(
95+
stv, params, VTP, VTP_covariances, verbose=False
96+
)
97+
98+
# Print the optimized parameters
99+
print("Optimized equation of state for stishovite (HP):")
100+
pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)
101+
print("")
102+
103+
# Create a corner plot of the covariances
104+
fig = burnman.nonlinear_fitting.corner_plot(
105+
fitted_eos.popt, fitted_eos.pcov, params
106+
)
107+
plt.show()
108+
109+
# Finally, let's plot our equation of state
110+
T = 298.15
111+
pressures = np.linspace(1.0e5, 60.0e9, 101)
112+
volumes = np.empty_like(pressures)
113+
114+
VTPs = np.empty((len(pressures), 3))
115+
for i, P in enumerate(pressures):
116+
stv.set_state(P, T)
117+
VTPs[i] = [stv.V, stv.temperature, stv.pressure]
118+
119+
# Plot the 95% confidence and prediction bands
120+
cp_bands = burnman.nonlinear_fitting.confidence_prediction_bands(
121+
model=fitted_eos,
122+
x_array=VTPs,
123+
confidence_interval=0.95,
124+
f=attribute_function(stv, "P", using_volume=True),
125+
flag="P",
126+
)
127+
128+
plt.fill_between(
129+
VTPs[:, 0] * 1.0e6,
130+
cp_bands[0] / 1.0e9,
131+
cp_bands[1] / 1.0e9,
132+
alpha=0.3,
133+
linewidth=0.0,
134+
color="r",
135+
label="95% confidence bands",
136+
)
137+
plt.fill_between(
138+
VTPs[:, 0] * 1.0e6,
139+
cp_bands[2] / 1.0e9,
140+
cp_bands[3] / 1.0e9,
141+
alpha=0.3,
142+
linewidth=0.0,
143+
color="b",
144+
label="95% prediction bands",
145+
)
146+
147+
plt.plot(
148+
VTPs[:, 0] * 1.0e6, VTPs[:, 2] / 1.0e9, label="Optimized fit for stishovite"
149+
)
150+
plt.errorbar(
151+
VTP[:, 0] * 1.0e6,
152+
VTP[:, 2] / 1.0e9,
153+
xerr=np.sqrt(VTP_covariances.T[0][0]) * 1.0e6,
154+
yerr=np.sqrt(VTP_covariances.T[2][2]) / 1.0e9,
155+
linestyle="None",
156+
marker="o",
157+
label="Andrault et al. (2003)",
158+
)
159+
160+
plt.xlabel("Volume (cm^3/mol)")
161+
plt.ylabel("Pressure (GPa)")
162+
plt.legend(loc="upper right")
163+
plt.title("Stishovite EoS (room temperature)")
164+
plt.show()
165+
166+
# Now let's fit the same parameters for the SLB EoS.
167+
stv_SLB = burnman.minerals.SLB_2011.stishovite()
168+
169+
# Removing the property modifiers makes fitting the SLB EoS
170+
# a lot faster as we avoid an iterative solve
171+
stv_SLB.property_modifiers = []
172+
173+
fitted_eos_SLB = burnman.eos_fitting.fit_VTP_data(
174+
stv_SLB, params, VTP, VTP_covariances, verbose=False
175+
)
176+
# Print the optimized parameters
177+
print("Optimized equation of state for stishovite (SLB):")
178+
pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)
179+
print("")
180+
181+
# 2) Now let's fit some F(V, T) data
182+
print("2) Fitting to synthetic F(V, T) data\n")
183+
184+
# Create the mineral instance that we'll optimise
185+
per_SLB = burnman.minerals.SLB_2011.periclase()
186+
187+
V_0 = per_SLB.params["V_0"]
188+
volumes = np.linspace(0.7 * V_0, V_0, 4)
189+
temperatures = np.linspace(500.0, 2000.0, 4)
190+
VV, TT = np.meshgrid(volumes, temperatures)
191+
V_data = VV.flatten()
192+
T_data = TT.flatten()
193+
F_data = per_SLB.evaluate_with_volumes(["helmholtz"], V_data, T_data)[0]
194+
195+
# Add some random noise
196+
np.random.seed(100)
197+
F_stdev = 10.0
198+
F_data = F_data + np.random.normal(scale=F_stdev, size=F_data.shape)
199+
200+
# Collect all the objects we need to supply as arguments to the
201+
# fitting function
202+
VTF_data = np.array([V_data, T_data, F_data]).T
203+
nul = 0.0 * V_data
204+
Fcov = F_stdev * F_stdev * np.ones_like(nul)
205+
VTF_covariances = np.array([[nul, nul, nul], [nul, nul, nul], [nul, nul, Fcov]]).T
206+
flags = ["helmholtz"] * len(F_data)
207+
fit_params = ["V_0", "K_0", "Kprime_0", "grueneisen_0", "q_0", "Debye_0", "F_0"]
208+
209+
print("Actual equation of state:")
210+
for p in fit_params:
211+
print(f"{p}: {per_SLB.params[p]:.4e}")
212+
print("")
213+
214+
# Let's change some parameters in our equation of state
215+
# so that fitting is not trivial
216+
per_SLB.params["F_0"] = per_SLB.params["F_0"] + 2300.0
217+
per_SLB.params["V_0"] = per_SLB.params["V_0"] * 0.9
218+
per_SLB.params["K_0"] = per_SLB.params["K_0"] * 0.9
219+
220+
# And now we can fit our data!
221+
fitted_eos = burnman.eos_fitting.fit_VTp_data(
222+
mineral=per_SLB,
223+
flags=flags,
224+
fit_params=fit_params,
225+
data=VTF_data,
226+
data_covariances=VTF_covariances,
227+
verbose=False,
228+
)
229+
230+
# Print the optimized parameters
231+
print("Optimized equation of state:")
232+
pretty_print_values(fitted_eos.popt, fitted_eos.pcov, fitted_eos.fit_params)
233+
print("\nGoodness of fit:")
234+
print(fitted_eos.goodness_of_fit)
235+
print("")

‎misc/ref/example_fit_eos_with_independent_volumes.py.out‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
Least squares equation of state fitting with volume as an independent variable.
2+
3+
1) Fitting to stishovite room temperature P(V) data
4+
5+
Optimized equation of state for stishovite (HP):
6+
V_0: (1.4005 +/- 0.0001) x 1e-05
7+
K_0: (3.12 +/- 0.03) x 1e+11
8+
Kprime_0: (4.4 +/- 0.2) x 1e+00
9+
10+
Optimized equation of state for stishovite (SLB):
11+
V_0: (1.4005 +/- 0.0001) x 1e-05
12+
K_0: (3.12 +/- 0.03) x 1e+11
13+
Kprime_0: (4.4 +/- 0.2) x 1e+00
14+
15+
2) Fitting to synthetic F(V, T) data
16+
17+
Actual equation of state:
18+
V_0: 1.1244e-05
19+
K_0: 1.6138e+11
20+
Kprime_0: 3.8405e+00
21+
grueneisen_0: 1.3613e+00
22+
q_0: 1.7217e+00
23+
Debye_0: 7.6710e+02
24+
F_0: -5.6944e+05
25+
26+
Optimized equation of state:
27+
V_0: (1.1245 +/- 0.0001) x 1e-05
28+
K_0: (1.613 +/- 0.002) x 1e+11
29+
Kprime_0: (3.841 +/- 0.004) x 1e+00
30+
grueneisen_0: (1.360 +/- 0.003) x 1e+00
31+
q_0: (1.71 +/- 0.01) x 1e+00
32+
Debye_0: (7.671 +/- 0.001) x 1e+02
33+
F_0: (-5.69450 +/- 0.00008) x 1e+05
34+
35+
Goodness of fit:
36+
0.673607559420123
37+

0 commit comments

Comments
 (0)