/
performance.py
97 lines (79 loc) · 3.1 KB
/
performance.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
from __future__ import absolute_import
from __future__ import print_function
# 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.
import numpy as np
import burnman_path # adds the local burnman directory to the path
import burnman
from burnman import minerals
import timeit
assert burnman_path # silence pyflakes warning
if True:
def test_set_state(runs=10000):
g = burnman.minerals.SLB_2011.garnet()
g.set_composition([0.1, 0.2, 0.4, 0.2, 0.1])
x = 0.0
for x in range(0,runs):
g.set_state(10.e9+x, 1500.)
x += g.activities
print(x)
test_set_state(1)
if False:
import cProfile,pstats,StringIO
pr = cProfile.Profile()
pr.enable()
test_set_state()
pr.disable()
pr.dump_stats("test_set_state.stats")
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print (s.getvalue())
if True:
def test_grueneisen():
for aa in range(0, 500):
a = burnman.eos.SLB3()
m = minerals.SLB_2011.fayalite()
x = 0
for i in range(0, 10000):
x += a.grueneisen_parameter(1e9, 1e4, 1.1, m.params)
return x
test_grueneisen()
# r=timeit.timeit('__main__.test_grueneisen()', setup="import __main__", number=3)
# print("test_grueneisen", r)
if True:
seismic_model = burnman.seismic.PREM()
number_of_points = 10 # set on how many depth slices the computations should be done
depths = np.linspace(1000e3, 2500e3, 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)
temperature = burnman.geotherm.brown_shankland(depths)
def calc_velocities(a, b, c):
amount_perovskite = a
pv = minerals.SLB_2011.mg_fe_perovskite([b, 1.0 - b, 0.0])
fp = minerals.SLB_2011.ferropericlase([c, 1.0 - c])
rock = burnman.Composite(
[pv, fp], [amount_perovskite, 1.0 - amount_perovskite])
mat_rho, mat_vp, mat_vs = rock.evaluate(
['density', 'v_phi', 'v_s'], seis_p, temperature)
return mat_vp, mat_vs, mat_rho
def error(a, b, c):
mat_vp, mat_vs, mat_rho = calc_velocities(a, b, c)
vs_err = burnman.tools.math.l2(depths, mat_vs, seis_vs) / 1e9
vp_err = burnman.tools.math.l2(depths, mat_vp, seis_vp) / 1e9
den_err = burnman.tools.math.l2(depths, mat_rho, seis_rho) / 1e9
return vs_err + vp_err + den_err
def run():
n = 5
m = 1e10
for a in np.linspace(0.0, 1.0, n):
for b in np.linspace(0.0, 1.0, n):
for c in np.linspace(0.0, 1.0, n):
m = min(m, error(a, b, c))
print(m)
return m
run()
# r=timeit.timeit('__main__.run()', setup="import __main__", number=3)
# print("test_seismic_velocities", r)