-
Notifications
You must be signed in to change notification settings - Fork 41
/
example_averaging.py
140 lines (107 loc) · 4.73 KB
/
example_averaging.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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
# 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_averaging
-----------------
This example shows the effect of different averaging schemes. Currently four
averaging schemes are available:
1. Voight-Reuss-Hill
2. Voight averaging
3. Reuss averaging
4. Hashin-Shtrikman averaging
See :cite:`Watt1976` Journal of Geophysics and Space Physics for explanations
of each averaging scheme.
*Specifically uses:*
* :class:`burnman.averaging_schemes.VoigtReussHill`
* :class:`burnman.averaging_schemes.Voigt`
* :class:`burnman.averaging_schemes.Reuss`
* :class:`burnman.averaging_schemes.HashinShtrikmanUpper`
* :class:`burnman.averaging_schemes.HashinShtrikmanLower`
*Demonstrates:*
* implemented averaging schemes
"""
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__":
# Create a rock out of MgSiO3 perovskite and MgO periclase
amount_perovskite = 0.6
rock = burnman.Composite([minerals.SLB_2011.mg_perovskite(),
minerals.SLB_2011.periclase()],
[amount_perovskite, 1.0 - amount_perovskite])
perovskitite = minerals.SLB_2011.mg_perovskite()
perovskitite.name = 'Mg Perovskite'
periclasite = minerals.SLB_2011.periclase()
periclasite.name = 'Periclase'
# seismic model for comparison:
# pick from .prem() .slow() .fast() (see burnman/seismic.py)
seismic_model = burnman.seismic.PREM()
# set on how many depth slices the computations should be done
number_of_points = 20
# we will do our computation and comparison at the following depth values:
depths = np.linspace(700e3, 2800e3, number_of_points)
# alternatively, we could use the values where prem is defined:
# depths = seismic_model.internal_depth_list(mindepth=700.e3,
# maxdepth=2800.e3)
pressures, seis_rho, seis_vp, seis_vs, seis_vphi = seismic_model.evaluate(
['pressure', 'density', 'v_p', 'v_s', 'v_phi'], depths)
temperatures = burnman.geotherm.brown_shankland(depths)
print("Calculations are done for:")
rock.debug_print()
# calculate the seismic velocities of the rock using a whole battery of
# averaging schemes:
# Voigt Reuss Hill / Voigt / Reuss averaging
rock_v = rock.copy()
rock_v.set_averaging_scheme(burnman.averaging_schemes.Voigt())
rock_v.name = 'Voigt'
rock_r = rock.copy()
rock_r.set_averaging_scheme(burnman.averaging_schemes.Reuss())
rock_r.name = 'Reuss'
rock_vrh = rock.copy()
rock_vrh.set_averaging_scheme(burnman.averaging_schemes.VoigtReussHill())
rock_vrh.name = 'Voigt-Reuss-Hill'
# Upper/lower bound for Hashin-Shtrikman averaging
rock_hsu = rock.copy()
rock_hsu.set_averaging_scheme(burnman.averaging_schemes.HashinShtrikmanUpper())
rock_hsu.name = 'Hashin-Shtrikman upper'
rock_hsl = rock.copy()
rock_hsl.set_averaging_scheme(burnman.averaging_schemes.HashinShtrikmanLower())
rock_hsl.name = 'Hashin-Shtrikman lower'
rocks = [rock_v, rock_r, rock_vrh, rock_hsu, rock_hsl,
perovskitite, periclasite]
markers = ['^', 'v', 'x', 'x', 'x', 'x', 'x']
colors = ['c', 'k', 'b', 'r', 'r', 'y', 'g']
v_ss = [rock.evaluate(['shear_wave_velocity'],
pressures, temperatures)[0] for rock in rocks]
# PLOTTING
# plot vs
fig = plt.figure(figsize=(10, 7))
for i, rock in enumerate(rocks):
plt.plot(pressures / 1.e9, v_ss[i] / 1.e3, color=colors[i],
linestyle='-', marker=markers[i], markersize=4,
label=rock.name)
plt.xlim(min(pressures) / 1.e9, max(pressures) / 1.e9)
plt.legend(loc='upper left', prop={'size': 11}, frameon=False)
plt.xlabel('pressure (GPa)')
plt.ylabel('Vs (km/s)')
v_s_norms = [(v_s - v_ss[-1]) / (v_ss[-2] - v_ss[-1]) for v_s in v_ss]
ax = fig.add_axes([0.58, 0.18, 0.3, 0.3])
for i, rock in enumerate(rocks):
plt.plot(pressures / 1.e9, v_s_norms[i], color=colors[i],
linestyle='-', marker=markers[i], markersize=4,
label=rock.name)
ax.tick_params(labelsize=10)
plt.title("normalized by mixture endmembers", fontsize=10)
plt.xlim(min(pressures) / 1.e9, max(pressures) / 1.e9)
plt.ylim(-0.005, 1.005)
plt.xlabel('pressure (GPa)', fontsize=10)
plt.ylabel('normalized Vs', fontsize=10)
plt.savefig("output_figures/example_averaging_normalized.png")
plt.show()