-
Notifications
You must be signed in to change notification settings - Fork 41
/
example_premite_isothermal.py
130 lines (98 loc) · 4.2 KB
/
example_premite_isothermal.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
# 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.
"""
This example is under construction.
requires:
teaches:
"""
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
import pymc
assert burnman_path # silence pyflakes warning
seismic_model = burnman.seismic.PREM() # pick a seismic model, see seismiic.py for details
number_of_points = 20 # set on how many depth slices the computations should be done
depths = np.linspace(750.e3, 2890.e3, number_of_points)
seis_p, seis_rho, seis_vp, seis_vs, seis_vphi = seismic_model.evaluate(
['pressure', 'density', 'v_p', 'v_s', 'v_phi'], depths)
print("preparations done")
def calc_velocities(ref_rho, K_0, K_prime, G_0, G_prime):
rock = burnman.Mineral()
rock.params['V_0'] = 10.e-6
rock.params['molar_mass'] = ref_rho * rock.params['V_0']
rock.params['K_0'] = K_0
rock.params['Kprime_0'] = K_prime
rock.params['G_0'] = G_0
rock.params['Gprime_0'] = G_prime
rock.set_method('bm3')
temperature = np.empty_like(seis_p)
mat_rho, mat_vphi, mat_vs = \
rock.evaluate(['density', 'v_phi', 'v_s'], seis_p, temperature)
return mat_rho, mat_vphi, mat_vs
def error(ref_rho, K_0, K_prime, G_0, G_prime):
rho, vphi, vs = calc_velocities(ref_rho, K_0, K_prime, G_0, G_prime)
vphi_chi = burnman.tools.math.chi_factor(vphi, seis_vphi)
vs_chi = burnman.tools.math.chi_factor(vs, seis_vs)
rho_chi = burnman.tools.math.chi_factor(rho, seis_rho)
return rho_chi + vphi_chi + vs_chi
if __name__ == "__main__":
# Priors on unknown parameters:
ref_rho = pymc.Uniform('ref_rho', lower=3300., upper=4500.)
K_0 = pymc.Uniform('K_0', lower=200.e9, upper=300.e9)
K_prime = pymc.Uniform('Kprime_0', lower=3., upper=6.)
G_0 = pymc.Uniform('G_0', lower=50.e9, upper=250.e9)
G_prime = pymc.Uniform('Gprime_0', lower=0., upper=3.)
minerr = 1e100
@pymc.deterministic
def theta(p1=ref_rho, p2=K_0, p3=K_prime, p4=G_0, p5=G_prime):
global minerr
if (p1 < 0 or p2 < 0 or p3 < 0 or p4 < 0 or p5 < 0):
return 1e30
try:
e = error(p1, p2, p3, p4, p5)
if (e < minerr):
minerr = e
print("best fit", e, "values:",
p1, p2 / 1.e9, p3, p4 / 1.e9, p5)
return e
except ValueError:
return 1e20
sig = 1e-4
misfit = pymc.Normal('d', mu=theta, tau=1.0 / (
sig * sig), value=0, observed=True, trace=True)
model = dict(ref_rho=ref_rho, K_0=K_0, K_prime=K_prime,
G_0=G_0, G_prime=G_prime, misfit=misfit)
things = ['ref_rho', 'K_0', 'K_prime', 'G_0', 'G_prime']
S = pymc.MAP(model)
S.fit(method='fmin')
rho, vphi, vs = calc_velocities(
S.ref_rho.value, S.K_0.value, S.K_prime.value, S.G_0.value, S.G_prime.value)
plt.subplot(2, 2, 1)
plt.plot(seis_p / 1.e9, vs / 1000., color='r', linestyle='-',
marker='^', markerfacecolor='r', markersize=4)
plt.plot(seis_p / 1.e9, seis_vs / 1000., color='k',
linestyle='-', marker='v', markerfacecolor='k', markersize=4)
plt.ylim([4, 8])
plt.title("Vs (km/s)")
plt.subplot(2, 2, 2)
plt.plot(seis_p / 1.e9, vphi / 1000., color='r', linestyle='-',
marker='^', markerfacecolor='r', markersize=4)
plt.plot(seis_p / 1.e9, seis_vphi / 1000., color='k',
linestyle='-', marker='v', markerfacecolor='k', markersize=4)
plt.ylim([7, 12])
plt.title("Vphi (km/s)")
# plot density
plt.subplot(2, 2, 3)
plt.plot(seis_p / 1.e9, rho / 1000., color='r', linestyle='-',
marker='^', markerfacecolor='r', markersize=4, label='model 1')
plt.plot(seis_p / 1.e9, seis_rho / 1000., color='k', linestyle='-',
marker='v', markerfacecolor='k', markersize=4, label='ref')
plt.title("density ($\cdot 10^3$ kg/m$^3$)")
plt.legend(loc='upper left')
plt.ylim([3, 7])
plt.show()
print("done")