-
Notifications
You must be signed in to change notification settings - Fork 41
/
example_optimize_pv.py
105 lines (83 loc) · 3.44 KB
/
example_optimize_pv.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2015 by the BurnMan team, released under the GNU
# GPL v2 or later.
"""
example_optimize_pv
-------------------
Vary the amount perovskite vs. ferropericlase and compute the error in the
seismic data against PREM. For more extensive comments on this setup, see tutorial/step_2.py
*Uses:*
* :doc:`mineral_database`
* :class:`burnman.composite.Composite`
* :class:`burnman.seismic.PREM`
* :func:`burnman.geotherm.brown_shankland`
* :func:`burnman.material.Material.evaluate`
* :func:`burnman.tools.compare_l2`
*Demonstrates:*
* compare errors between models
* loops over models
"""
from __future__ import absolute_import
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
import burnman_path # adds the local burnman directory to the path
import burnman
from burnman import minerals
assert burnman_path # silence pyflakes warning
if __name__ == "__main__":
# Define reference model and depth to evaluate
seismic_model = burnman.seismic.PREM()
number_of_points = 20
depths = np.linspace(700e3, 2800e3, number_of_points)
seis_p, seis_rho, seis_vp, seis_vs, seis_vphi, seis_K, seis_G = seismic_model.evaluate(
['pressure', 'density', 'v_p', 'v_s', 'v_phi', 'K', 'G'], depths)
# Define geotherm
temperature = burnman.geotherm.brown_shankland(depths)
# Define solid solutions
perovskite = minerals.SLB_2011.mg_fe_perovskite()
# Set molar_fraction of fe_perovskite and al_perovskite:
perovskite.set_composition([0.94, 0.06, 0.])
ferropericlase = minerals.SLB_2011.ferropericlase()
# Set molar_fraction of MgO and FeO:
ferropericlase.set_composition([0.8, 0.2])
def material_error(amount_perovskite):
# Define composite using the values
rock = burnman.Composite([perovskite, ferropericlase],
[amount_perovskite, 1.0 - amount_perovskite])
# Compute velocities
mat_rho, mat_vp, mat_vs, mat_vphi, mat_K, mat_G = \
rock.evaluate(
['density', 'v_p', 'v_s', 'v_phi', 'K_S', 'G'], seis_p, temperature)
print("Calculations are done for:")
rock.debug_print()
# Calculate errors
[vs_err, vphi_err, rho_err, K_err, G_err] = \
burnman.tools.compare_l2(depths, [mat_vs, mat_vphi, mat_rho, mat_K, mat_G], [
seis_vs, seis_vphi, seis_rho, seis_K, seis_G])
# Normalize errors
vs_err = vs_err / np.mean(seis_vs) ** 2.
vphi_err = vphi_err / np.mean(seis_vphi) ** 2.
rho_err = rho_err / np.mean(seis_rho) ** 2.
K_err = K_err / np.mean(seis_K) ** 2.
G_err = G_err / np.mean(seis_G) ** 2.
return vs_err, vphi_err, rho_err, K_err, G_err
# Run through fractions of perovskite
xx = np.linspace(0.0, 1.0, 40)
errs = np.array([material_error(x) for x in xx])
# Plot results
yy_vs = errs[:, 0]
yy_vphi = errs[:, 1]
yy_rho = errs[:, 2]
yy_K = errs[:, 3]
yy_G = errs[:, 4]
plt.plot(xx, yy_vs, "r-x", label=("vs error"))
plt.plot(xx, yy_vphi, "b-x", label=("vphi error"))
plt.plot(xx, yy_rho, "m-x", label=("rho error"))
plt.plot(xx, yy_K, "g-x", label=("K error"))
plt.plot(xx, yy_G, "y-x", label=("G error"))
plt.yscale('log')
plt.xlabel('% Perovskite')
plt.ylabel('Error')
plt.legend()
plt.show()