Subduction Zone Benchmark#
Authors: Kidus Teshome, Cian Wilson
Convergence testing#
Preamble#
Let’s start by adding the path to the modules in the python folder to the system path (so we can find the our custom modules).
import sys, os
basedir = ''
if "__file__" in globals(): basedir = os.path.dirname(__file__)
sys.path.append(os.path.join(basedir, os.path.pardir, 'python'))
Then load everything we need from sz_problem and also set our default plotting preferences.
import utils
from sz_base import default_params, allsz_params
from sz_slab import create_slab
from sz_geometry import create_sz_geometry
from sz_problem import SubductionProblem, plot_slab_temperatures
import pyvista as pv
if __name__ == "__main__" and "__file__" in globals():
pv.OFF_SCREEN = True
import pathlib
if __name__ == "__main__":
output_folder = pathlib.Path(os.path.join(basedir, "output"))
output_folder.mkdir(exist_ok=True, parents=True)
Benchmark case 1#
def solve_benchmark_case1(resscale):
xs = [0.0, 140.0, 240.0, 400.0]
ys = [0.0, -70.0, -120.0, -200.0]
lc_depth = 40
uc_depth = 15
coast_distance = 0
extra_width = 0
sztype = 'continental'
io_depth_1 = 139
A = 100.0 # age of subducting slab (Myr)
qs = 0.065 # surface heat flux (W/m^2)
Vs = 100.0 # slab speed (mm/yr)
# create the slab
slab = create_slab(xs, ys, resscale, lc_depth)
# construct the geometry
geom = create_sz_geometry(slab, resscale, sztype, io_depth_1, extra_width,
coast_distance, lc_depth, uc_depth)
# set up a subduction zone problem
sz = SubductionProblem(geom, A=A, Vs=Vs, sztype=sztype, qs=qs)
# solve it using a steady state assumption and an isoviscous rheology
sz.solve_steadystate_isoviscous()
return sz
if __name__ == "__main__":
resscales = [4.0, 2.0, 1.0]
diagnostics_case1 = []
for resscale in resscales:
sz = solve_benchmark_case1(resscale)
diagnostics_case1.append((resscale, sz.get_diagnostics()))
T_ndof = 8871
T_(200,-100) = 517.05 deg C
T_slab = 451.76 deg C
T_wedge = 926.53 deg C
V_rms,w = 34.64 mm/yr
Warning : Gmsh has aleady been initialized
T_ndof = 32749
T_(200,-100) = 516.93 deg C
T_slab = 451.68 deg C
T_wedge = 926.26 deg C
V_rms,w = 34.64 mm/yr
Warning : Gmsh has aleady been initialized
T_ndof = 126633
T_(200,-100) = 516.85 deg C
T_slab = 451.63 deg C
T_wedge = 926.14 deg C
V_rms,w = 34.64 mm/yr
if __name__ == "__main__":
print('')
print('{:<12} {:<12} {:<12} {:<12} {:<12} {:<12}'.format('resscale', 'T_ndof', 'T_{200,-100}', 'Tbar_s', 'Tbar_w', 'Vrmsw'))
for resscale, diag in diagnostics_case1:
print('{:<12.4g} {:<12d} {:<12.4f} {:<12.4f} {:<12.4f} {:<12.4f}'.format(resscale, *diag))
resscale T_ndof T_{200,-100} Tbar_s Tbar_w Vrmsw
4 8871 517.0494 451.7642 926.5325 34.6383
2 32749 516.9332 451.6764 926.2604 34.6379
1 126633 516.8525 451.6301 926.1427 34.6377
For comparison here are the values reported for case 1 using TerraFERMA in Wilson & van Keken, 2023:
|
\(T_{\text{ndof}} \) |
\(T_{(200,-100)}^*\) |
\(\overline{T}_s^*\) |
\( \overline{T}_w^* \) |
\(V_{\text{rms},w}^*\) |
|---|---|---|---|---|---|
2.0 |
21403 |
517.17 |
451.83 |
926.62 |
34.64 |
1.0 |
83935 |
516.95 |
451.71 |
926.33 |
34.64 |
0.5 |
332307 |
516.86 |
451.63 |
926.15 |
34.64 |
Benchmark case 2#
def solve_benchmark_case2(resscale):
xs = [0.0, 140.0, 240.0, 400.0]
ys = [0.0, -70.0, -120.0, -200.0]
lc_depth = 40
uc_depth = 15
coast_distance = 0
extra_width = 0
sztype = 'continental'
io_depth_2 = 154
A = 100.0 # age of subducting slab (Myr)
qs = 0.065 # surface heat flux (W/m^2)
Vs = 100.0 # slab speed (mm/yr)
# create the slab
slab = create_slab(xs, ys, resscale, lc_depth)
# construct the geometry
geom = create_sz_geometry(slab, resscale, sztype, io_depth_2, extra_width,
coast_distance, lc_depth, uc_depth)
# set up a subduction zone problem
sz = SubductionProblem(geom, A=A, Vs=Vs, sztype=sztype, qs=qs)
# solve it using a steady state assumption and a dislocation creep rheology
sz.solve_steadystate_dislocationcreep()
return sz
if __name__ == "__main__":
resscales = [4.0, 2.0, 1.0]
diagnostics_case2 = []
for resscale in resscales:
sz = solve_benchmark_case2(resscale)
diagnostics_case2.append((resscale, sz.get_diagnostics()))
Warning : Gmsh has aleady been initialized
Iteration Residual Relative Residual
------------------------------------------
0 14370.5 1
1 1493.84 0.103952
2 313.238 0.0217973
3 114.859 0.00799269
4 50.178 0.00349174
5 25.488 0.00177363
6 14.7551 0.00102677
7 9.26217 0.000644528
8 6.18334 0.000430281
9 4.32205 0.000300759
10 3.12156 0.00021722
11 2.31066 0.000160792
12 1.74374 0.000121342
13 1.3363 9.29895e-05
14 1.03696 7.21592e-05
15 0.813199 5.65882e-05
16 0.643615 4.47873e-05
17 0.513608 3.57405e-05
18 0.412917 2.87337e-05
19 0.334197 2.32558e-05
20 0.272119 1.8936e-05
21 0.22277 1.55019e-05
22 0.183246 1.27516e-05
23 0.151373 1.05336e-05
24 0.125503 8.73341e-06
25 0.104384 7.26377e-06
26 0.08705 6.05756e-06
27 0.0727553 5.06283e-06
28 0.0609171 4.23904e-06
T_ndof = 8867
T_(200,-100) = 679.07 deg C
T_slab = 569.41 deg C
T_wedge = 938.74 deg C
V_rms,w = 41.31 mm/yr
Warning : Gmsh has aleady been initialized
Iteration Residual Relative Residual
------------------------------------------
0 12652.9 1
1 2593.5 0.204972
2 354.193 0.027993
3 126.582 0.0100041
4 56.6955 0.00448081
5 31.9665 0.0025264
6 20.3208 0.00160601
7 13.862 0.00109555
8 9.80868 0.00077521
9 7.05601 0.000557657
10 5.15867 0.000407705
11 3.82963 0.000302667
12 2.88071 0.000227671
13 2.18946 0.000173039
14 1.67703 0.000132541
15 1.29174 0.00010209
16 0.998644 7.89258e-05
17 0.773531 6.11344e-05
18 0.599318 4.73659e-05
19 0.463747 3.66513e-05
20 0.35788 2.82843e-05
21 0.275088 2.1741e-05
22 0.210376 1.66267e-05
23 0.159927 1.26395e-05
24 0.120782 9.54576e-06
25 0.0906273 7.16254e-06
26 0.0676413 5.3459e-06
27 0.0503836 3.98197e-06
T_ndof = 32725
T_(200,-100) = 683.72 deg C
T_slab = 573.07 deg C
T_wedge = 935.15 deg C
V_rms,w = 40.83 mm/yr
Warning : Gmsh has aleady been initialized
Iteration Residual Relative Residual
------------------------------------------
0 9605.5 1
1 4470.14 0.465373
2 475.321 0.0494843
3 211.402 0.0220084
4 107.367 0.0111777
5 67.8616 0.00706486
6 47.536 0.00494884
7 34.7909 0.00362198
8 25.8324 0.00268934
9 19.2772 0.00200689
10 14.42 0.00150122
11 10.808 0.00112519
12 8.11782 0.000845122
13 6.11086 0.000636184
14 4.61053 0.000479988
15 3.48603 0.000362921
16 2.64058 0.000274903
17 2.00269 0.000208494
18 1.51971 0.000158213
19 1.15284 0.000120019
20 0.873397 9.09268e-05
21 0.660128 6.8724e-05
22 0.497184 5.17603e-05
23 0.372685 3.87991e-05
24 0.277676 2.8908e-05
25 0.205363 2.13797e-05
26 0.150566 1.56749e-05
27 0.10932 1.1381e-05
28 0.0785894 8.18171e-06
29 0.0560496 5.83516e-06
30 0.0399335 4.15736e-06
T_ndof = 127428
T_(200,-100) = 682.95 deg C
T_slab = 572.36 deg C
T_wedge = 936.62 deg C
V_rms,w = 40.78 mm/yr
if __name__ == "__main__":
print('')
print('{:<12} {:<12} {:<12} {:<12} {:<12} {:<12}'.format('resscale', 'T_ndof', 'T_{200,-100}', 'Tbar_s', 'Tbar_w', 'Vrmsw'))
for resscale, diag in diagnostics_case2:
print('{:<12.4g} {:<12d} {:<12.4f} {:<12.4f} {:<12.4f} {:<12.4f}'.format(resscale, *diag))
resscale T_ndof T_{200,-100} Tbar_s Tbar_w Vrmsw
4 8867 679.0696 569.4063 938.7441 41.3075
2 32725 683.7207 573.0692 935.1450 40.8269
1 127428 682.9542 572.3554 936.6206 40.7846
For comparison here are the values reported for case 2 using TerraFERMA in Wilson & van Keken, 2023:
|
\(T_{\text{ndof}} \) |
\(T_{(200,-100)}^*\) |
\(\overline{T}_s^*\) |
\( \overline{T}_w^* \) |
\(V_{\text{rms},w}^*\) |
|---|---|---|---|---|---|
2.0 |
21403 |
683.05 |
571.58 |
936.65 |
40.89 |
1.0 |
83935 |
682.87 |
572.23 |
936.11 |
40.78 |
0.5 |
332307 |
682.80 |
572.05 |
937.37 |
40.77 |
Finish up#
Convert this notebook to a python script (making sure to save first)
if __name__ == "__main__" and "__file__" not in globals():
from ipylab import JupyterFrontEnd
app = JupyterFrontEnd()
app.commands.execute('docmanager:save')
!jupyter nbconvert --NbConvertApp.export_format=script --ClearOutputPreprocessor.enabled=True sz_benchmark.ipynb
[NbConvertApp] Converting notebook sz_benchmark.ipynb to script
[NbConvertApp] Writing 5440 bytes to sz_benchmark.py