-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsensitivity_test.py
237 lines (222 loc) · 9.45 KB
/
sensitivity_test.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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
"""
Modified version of the run_model version to
run sensitivity tests and evaluate the stability
of our models.
=============================================
Requirements:
*
*
=============================================
VUW
May 2019
Author: Konstantinos Michailos
"""
import numpy as np
from scipy.optimize import minimize
import os
from main_functions.functions import *
# Define fixed parameters
# Year in seconds
year = 365.25*24*60*60
# Material properties
# average crustal density of 2700 kg m-3
# J/kgK Volumetric heat capacity = density * specific H capacity
C = 2700.*790.
# W/mK Thermal conductivity
k = 3.2
# W/m^3s Volumetric heat productivity - typical value for granite/
# heat production
H = 3.e-6
# Sea level; top of exhuming block
z0 = 0.
# Temperature at the top of exhuming block
T0 = 13.
# Regularization parameter (Q = alpha * Q_therm + Q_seis)
alpha = 1.0
# Decollement depth; bottom of exhuming block
z1 = -35000.
# Temperature at the bottom of exhuming block
T1 = [400., 450., 500., 550., 600., 650., 700.]
# Choose which model to run
# Steady state exhumation OR initial geotherm
# assumption = '_steady_state'
assumption = '_initial_state'
# Define working directory
work_dir = os.getcwd()
# Define path of thermochronological data
Tc_path = (work_dir + '/dataset/tc_data/')
# Define path of seismicity data
npzfilespath = (work_dir + '/dataset/npz_files/')
# Read them observations
observations = read_observations(Tc_path, npzfilespath)
seis_obs = observations[1]
thermochron_obs = observations[0]
# Initial guesses for adjustable parameters
# Temperature at the brittle-ductile depth (initial guess)
T_BDT = [200., 250., 300., 350., 400., 450., 500.]
# First line of grid points
vel_A1 = 2.0
vel_A2 = 2.0
vel_A3 = 2.
vel_A4 = 4.
vel_A5 = 8.
vel_A6 = 6.5
vel_A7 = 4.
vel_A8 = 2.
vel_A9 = 1.5
vel_A10 = 1.
vel_A11 = 1.
# Second line of grid points
vel_B1 = 1.5
vel_B2 = 1.5
vel_B3 = 1.5
vel_B4 = 3.
vel_B5 = 6.
vel_B6 = 6.
vel_B7 = 4.
vel_B8 = 2.
vel_B9 = 2.
vel_B10 = 1.
vel_B11 = 1.
# Third line of grid points
vel_C1 = 1.5
vel_C2 = 1.5
vel_C3 = 1.5
vel_C4 = 2.
vel_C5 = 4.
vel_C6 = 4.
vel_C7 = 3.
vel_C8 = 2.
vel_C9 = 2.
vel_C10 = 1.
vel_C11 = 1.
initial_uplifts = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
vel_A1, vel_A2, vel_A3, vel_A4, vel_A5, vel_A6,
vel_A7, vel_A8, vel_A9, vel_A10, vel_A11,
vel_B1, vel_B2, vel_B3, vel_B4, vel_B5, vel_B6,
vel_B7, vel_B8, vel_B9, vel_B10, vel_B11,
vel_C1, vel_C2, vel_C3, vel_C4, vel_C5, vel_C6,
vel_C7, vel_C8, vel_C9, vel_C10, vel_C11,
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
for i, T1_ in enumerate(T1):
for k, T_BDT_ in enumerate(T_BDT):
if T_BDT_ < T1_:
# Set model's run name to be used
# to create a folder in the outputs directory
run_name = ('Tbdt_' + str(int(T_BDT_)) + '_T1_' +
str(int(T1_)) + '_' + str(int(z1/(-1000))) + assumption)
print(run_name)
# Parameters for the minimization process
# List of adjustable parameters
# (note that z1 and T1 are in the adjustable but are used as fixed)
p = [vel_A1, vel_A2, vel_A3, vel_A4, vel_A5,
vel_A6, vel_A7, vel_A8, vel_A9,
vel_A10, vel_A11,
vel_B1, vel_B2, vel_B3, vel_B4, vel_B5,
vel_B6, vel_B7, vel_B8, vel_B9,
vel_B10, vel_B11,
vel_C1, vel_C2, vel_C3, vel_C4, vel_C5,
vel_C6, vel_C7, vel_C8, vel_C9,
vel_C10, vel_C11, z1, T1_, T_BDT_]
# List of fixed parameters
fixed_par = [T1_, z1, thermochron_obs, seis_obs, alpha, work_dir, run_name]
pfix = np.array(fixed_par)
x0 = np.array(p)
# Make directories for outputs of model
output_path = os.path.join(work_dir, 'outputs', run_name)
make_dirs(output_path)
log_dir = os.path.join(output_path, 'log_files')
make_dirs(log_dir)
plot_dir = os.path.join(output_path, 'plots')
make_dirs(plot_dir)
# If these directories exist delete the files within em
filenames = [log_dir + '/depths.csv',
log_dir + '/ages.csv',
output_path + '/mod_uplifts.txt',
log_dir + '/thermo_res.csv',
log_dir + '/seis_res.csv',
log_dir + '/temps.txt',
log_dir + '/logfile_SWISS.txt']
for filepath in filenames:
try:
os.remove(filepath)
except Exception as e:
print(e)
print("Error while deleting file ", filepath, ", file doesn't exist!")
def printx(Xi):
"""Function used to output the values of the parameters as the iterations
progress.
"""
global Nfeval
global fout
# print(Nfeval, Q(Xi))
fout.write('{0: 3.1f}, {1: 3.1f}, {2: 3.1f},{3: 3.1f}, {4: 3.1f}, {5: 3.1f},\
{6: 3.1f}, {7: 3.1f}, {8: 3.1f},{9: 3.1f},\
{10: 3.1f}, {11: 3.1f}, {12: 3.1f},{13: 3.1f}, {14: 3.1f}, {15: 3.1f},\
{16: 3.1f}, {17: 3.1f}, {18: 3.1f},{19: 3.1f},{20: 3.1f}, {21: 3.1f}, {22: 3.1f},{23: 3.1f}, {24: 3.1f}, {25: 3.1f},\
{26: 3.1f}, {27: 3.1f}, {28: 3.1f},{29: 3.1f}, {30: 3.1f}, {31: 3.1f},\
{32: 3.1f}, {33: 3.1f}, {34:4d}, {35: 5.2f}'
.format(Xi[0], Xi[1], Xi[2],
Xi[3], Xi[4], Xi[5],
Xi[6], Xi[7], Xi[8],
Xi[9], Xi[10], Xi[11],
Xi[12], Xi[13], Xi[14],
Xi[15], Xi[16], Xi[17],
Xi[18], Xi[19], Xi[20],
Xi[21], Xi[22], Xi[23],
Xi[24], Xi[25], Xi[26],
Xi[27], Xi[28], Xi[29],
Xi[30], Xi[31], Xi[32], Xi[33],
Nfeval, Q(Xi, pfix)) + '\n')
Nfeval += 1
Nfeval = 1
fout = open(log_dir + '/logfile_SWISS.txt', 'w')
# Choosing which Q function to use from functions.py
if assumption == '_steady_state':
Q = Q_steady_state
else:
Q = Q_initial_state
# Run minimization process
res = minimize(Q, x0, args=(pfix), method='powell',
callback=printx, options={'disp': True, 'maxiter': 200})
fout.close()
# Read minimization process results and set the final modelled
# exhumation rates on the model's grid points
final_uplifts = define_grid(res.x[0], res.x[1], res.x[2], res.x[3], res.x[4],
res.x[5], res.x[6], res.x[7], res.x[8], res.x[9],
res.x[10], res.x[11], res.x[12], res.x[13],
res.x[14], res.x[15], res.x[16], res.x[17],
res.x[18], res.x[19], res.x[20], res.x[21],
res.x[22], res.x[23], res.x[24], res.x[25],
res.x[26], res.x[27], res.x[28], res.x[29],
res.x[30], res.x[31], res.x[32])
# Write output for with final exhumation rates
for i, j in enumerate(final_uplifts[0]):
print(final_uplifts[1][i])
with open(output_path + '/mod_uplifts.txt', 'a') as of:
of.write('{} {} {}\n'.
format(final_uplifts[1][i],
final_uplifts[0][i],
round(final_uplifts[2][i], 2)))
# changes in uplift rates
for i, j in enumerate(final_uplifts[0]):
print(final_uplifts[1][i])
print(len(final_uplifts[2]), len(initial_uplifts))
dif = round(final_uplifts[2][i], 2)- initial_uplifts[i]
with open(output_path + '/UPLIFT_CHANGES.txt', 'a') as of:
of.write('{}, {}, {}, {}, {}\n'.
format(final_uplifts[1][i],
final_uplifts[0][i],
round(final_uplifts[2][i], 2),
initial_uplifts[i],
dif))
# Make plots
# Read residuals from the last iteration of the minimization
# and create lists that will be inputs for the maps
seis_res = sort_residuals('seis_res.csv', log_dir)
ther_res = sort_residuals('thermo_res.csv', log_dir)
plot3 = plot_res(seis_res, ther_res, run_name, plot_dir)
plots1 = plot_temps(final_uplifts, run_name, p, res, work_dir, plot_dir)
plots2 = plot_exhum(final_uplifts, run_name, plot_dir)
# This line creates the input for the gmt script cross_sections_temperat.sh
cross_section_prep = temp_cross_section(final_uplifts, p, work_dir)