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, os.path.pardir, 'python'))
Then load everything we need from sz_problems and also set our default plotting preferences.
import utils
from sz_problems.sz_params import default_params, allsz_params
from sz_problems.sz_slab import create_slab
from sz_problems.sz_geometry import create_sz_geometry
from sz_problems.sz_steady_isoviscous import SteadyIsoSubductionProblem
from sz_problems.sz_steady_dislcreep import SteadyDislSubductionProblem
import pathlib
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 = 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, extra_width,
coast_distance, lc_depth, uc_depth)
# set up a subduction zone problem
sz = SteadyIsoSubductionProblem(geom, A=A, Vs=Vs, sztype=sztype, qs=qs)
# solve it using a steady state assumption and an isoviscous rheology
sz.solve()
# evaluate the diagnostics
diag = sz.get_diagnostics()
return diag
resscales = [4.0, 2.0, 1.0]
diagnostics_case1 = []
for resscale in resscales:
diagnostics_case1.append((resscale, solve_benchmark_case1(resscale)))
T_ndof = 8871
T_(200,-100) = 517.05 deg C
slab_diag_length = 111.80
T_slab = 451.76 deg C
wedge_diag_area = 5500.00
T_wedge = 926.53 deg C
V_rms,w = 34.64 mm/yr
T_ndof = 32749
T_(200,-100) = 516.93 deg C
slab_diag_length = 111.80
T_slab = 451.68 deg C
wedge_diag_area = 5500.00
T_wedge = 926.26 deg C
V_rms,w = 34.64 mm/yr
T_ndof = 126633
T_(200,-100) = 516.85 deg C
slab_diag_length = 111.80
T_slab = 451.63 deg C
wedge_diag_area = 5500.00
T_wedge = 926.14 deg C
V_rms,w = 34.64 mm/yr
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.values()))
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 |
Which we can test against our values.
values_wvk_case1 = [
{'resscale': 2, 'T_ndof': 21403, 'T_{200,-100}': 517.17, 'Tbar_s': 451.83, 'Tbar_w': 926.62, 'Vrmsw': 34.64},
{'resscale': 1, 'T_ndof': 83935, 'T_{200,-100}': 516.95, 'Tbar_s': 451.71, 'Tbar_w': 926.33, 'Vrmsw': 34.64},
{'resscale': 0.5, 'T_ndof': 332307, 'T_{200,-100}': 516.86, 'Tbar_s': 451.63, 'Tbar_w': 926.15, 'Vrmsw': 34.64},
]
print('')
print('{:<12} {:<12}'.format('', 'error'))
for key, val in diagnostics_case1[-1][-1].items():
if key != "T_ndof":
err = abs(values_wvk_case1[1][key]-val)/val
print('{:<12} {:<12.4g}'.format(key, err))
assert(err < 1.e-3) # check error is less than 0.1%
error
T_{200,-100} 0.0001886
Tbar_s 0.0001768
Tbar_w 0.0002022
Vrmsw 6.596e-05
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 = 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, extra_width,
coast_distance, lc_depth, uc_depth)
# set up a subduction zone problem
sz = SteadyDislSubductionProblem(geom, A=A, Vs=Vs, sztype=sztype, qs=qs)
# solve it using a steady state assumption and a dislocation creep rheology
sz.solve()
# evaluate the diagnostics
diag = sz.get_diagnostics()
return diag
resscales = [4.0, 2.0, 1.0]
diagnostics_case2 = []
for resscale in resscales:
diagnostics_case2.append((resscale, solve_benchmark_case2(resscale)))
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.23905e-06
T_ndof = 8867
T_(200,-100) = 679.07 deg C
slab_diag_length = 111.80
T_slab = 569.41 deg C
wedge_diag_area = 5500.00
T_wedge = 938.74 deg C
V_rms,w = 41.31 mm/yr
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
slab_diag_length = 111.80
T_slab = 573.07 deg C
wedge_diag_area = 5500.00
T_wedge = 935.15 deg C
V_rms,w = 40.83 mm/yr
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
slab_diag_length = 111.80
T_slab = 572.36 deg C
wedge_diag_area = 5500.00
T_wedge = 936.62 deg C
V_rms,w = 40.78 mm/yr
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.values()))
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 |
Which we can test against our values.
values_wvk_case2 = [
{'resscale': 2, 'T_ndof': 21403, 'T_{200,-100}': 683.05, 'Tbar_s': 571.58, 'Tbar_w': 936.65, 'Vrmsw': 40.89},
{'resscale': 1, 'T_ndof': 83935, 'T_{200,-100}': 682.87, 'Tbar_s': 572.23, 'Tbar_w': 936.11, 'Vrmsw': 40.78},
{'resscale': 0.5, 'T_ndof': 332307, 'T_{200,-100}': 682.80, 'Tbar_s': 572.05, 'Tbar_w': 937.37, 'Vrmsw': 40.77},
]
print('')
print('{:<12} {:<12}'.format('', 'error'))
for key, val in diagnostics_case2[-1][-1].items():
if key != "T_ndof":
err = abs(values_wvk_case2[1][key]-val)/val
print('{:<12} {:<12.4g}'.format(key, err))
assert(err < 1.e-3) # check error is less than 0.1%
error
T_{200,-100} 0.0001233
Tbar_s 0.000219
Tbar_w 0.0005452
Vrmsw 0.0001139
Finish up#
Convert this notebook to a python module (saving first and ignoring markdown cells and those tagged as “main” or “ipy”).
from ipylab import JupyterFrontEnd
app = JupyterFrontEnd()
app.commands.execute('docmanager:save')
!jupyter nbconvert --TagRemovePreprocessor.enabled=True --TagRemovePreprocessor.remove_cell_tags="['main', 'ipy']" --TemplateExporter.exclude_markdown=True --TemplateExporter.exclude_input_prompt=True --TemplateExporter.exclude_output_prompt=True --NbConvertApp.export_format=script --ClearOutputPreprocessor.enabled=True --FilesWriter.build_directory=../../python/sz_problems --NbConvertApp.output_base=sz_benchmark 3.4b_sz_benchmark.ipynb
[NbConvertApp] Converting notebook 3.4b_sz_benchmark.ipynb to script
[NbConvertApp] Writing 2518 bytes to ../../python/sz_problems/sz_benchmark.py