40 S Marianas

Authors: Kidus Teshome, Cian Wilson

Steady state implementation

Preamble

Set some path information.

import sys, os
basedir = ''
if "__file__" in globals(): basedir = os.path.dirname(__file__)
sys.path.append(os.path.join(basedir, os.path.pardir))
sys.path.append(os.path.join(basedir, os.path.pardir, os.path.pardir, 'python'))

Loading everything we need from sz_problem and also set our default plotting and output preferences.

import utils
from sz_base import allsz_params
from sz_slab import create_slab, plot_slab
from sz_geometry import create_sz_geometry
from sz_problem import SubductionProblem
import numpy as np
import dolfinx as df
import pyvista as pv
if __name__ == "__main__" and "__file__" in globals():
    pv.OFF_SCREEN = True
import pathlib
output_folder = pathlib.Path(os.path.join(basedir, "output"))
output_folder.mkdir(exist_ok=True, parents=True)
import hashlib
import zipfile
import requests

Parameters

We first select the name and resolution scale, resscale of the model.

Resolution

By default the resolution is low to allow for a quick runtime and smaller website size. If sufficient computational resources are available set a lower resscale to get higher resolutions and results with sufficient accuracy.

name = "40_S_Marianas"
resscale = 3.0

Then load the remaining parameters from the global suite.

szdict = allsz_params[name]
print("{}:".format(name))
print("{:<20} {:<10}".format('Key','Value'))
print("-"*85)
for k, v in allsz_params[name].items():
    if v is not None and k not in ['z0', 'z15']: print("{:<20} {}".format(k, v))

40_S_Marianas:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       255
extra_width          21
lc_depth             30
io_depth             172
dirname              40_S_Marianas
As                   40.0
Ac                   65.0
A                    100.0
sztype               oceanic
Vs                   50.0
xs                   [0, 30.3, 78.7, 184.0, 205.4, 209.5, 230.0, 277.3, 314.0]
ys                   [-6, -15.0, -30.0, -70.0, -80.0, -82.5, -100.0, -175.0, -240.0]

Any of these can be modified in the dictionary.

Several additional parameters can be modified, for details see the documentation for the SubductionProblem class.

if __name__ == "__main__" and "__file__" not in globals():
    SubductionProblem?

The if __name__ == "__main__" and "__file__" not in globals(): logic above is only necessary to make sure that this only runs in the Jupyter notebook version of this code and not the python version. It is not generally necessary when getting the docstring of a function or class in Jupyter.

Setup

Setup a slab.

slab = create_slab(szdict['xs'], szdict['ys'], resscale, szdict['lc_depth'])
_ = plot_slab(slab)

Create the subduction zome geometry around the slab.

geom = create_sz_geometry(slab, resscale, szdict['sztype'], szdict['io_depth'], szdict['extra_width'], 
                             szdict['coast_distance'], szdict['lc_depth'], szdict['uc_depth'])
_ = geom.plot()

Finally, declare the SubductionZone problem class using the dictionary of parameters.

sz = SubductionProblem(geom, **szdict)

Solve

Solve using a dislocation creep rheology and assuming a steady state.

sz.solve_steadystate_dislocationcreep()

Iteration   Residual     Relative Residual
------------------------------------------
0           6958.54      1           

1           1822.13      0.261855    
2           305.525      0.0439064   

3           118.717      0.0170606   
4           46.703       0.00671162  

5           22.8324      0.00328121  
6           13.3765      0.00192232  

7           8.6081       0.00123706  
8           5.81486      0.000835645 

9           4.06661      0.000584406 
10          2.89632      0.000416225 

11          2.09177      0.000300605 
12          1.53048      0.000219943 

13          1.13621      0.000163283 
14          0.858115     0.000123318 

15          0.660855     9.49704e-05 
16          0.519458     7.46505e-05 

17          0.416244     5.98178e-05 
18          0.338964     4.8712e-05  

19          0.279344     4.01441e-05 
20          0.231978     3.33371e-05 

21          0.193422     2.77963e-05 
22          0.16149      2.32075e-05 

23          0.134764     1.93667e-05 
24          0.112276     1.61349e-05 

25          0.0933249    1.34116e-05 
26          0.0773705    1.11188e-05 

27          0.0639696    9.19297e-06 
28          0.0527473    7.58023e-06 

29          0.0433807    6.23417e-06 
30          0.0355894    5.1145e-06  

31          0.0291302    4.18626e-06 

Plot

Plot the solution.

plotter = pv.Plotter()
utils.plot_scalar(sz.T_i, plotter=plotter, scale=sz.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})
utils.plot_vector_glyphs(sz.vw_i, plotter=plotter, gather=True, factor=0.1, color='k', scale=utils.mps_to_mmpyr(sz.v0))
utils.plot_vector_glyphs(sz.vs_i, plotter=plotter, gather=True, factor=0.1, color='k', scale=utils.mps_to_mmpyr(sz.v0))
utils.plot_geometry(geom, plotter=plotter, color='green', width=2)
utils.plot_couplingdepth(slab, plotter=plotter, render_points_as_spheres=True, point_size=10.0, color='green')
utils.plot_show(plotter)
utils.plot_save(plotter, output_folder / "{}_ss_solution_resscale_{:.2f}.png".format(name, resscale))

Save it to disk so that it can be examined with other visualization software (e.g. Paraview).

filename = output_folder / "{}_ss_solution_resscale_{:.2f}.bp".format(name, resscale)
with df.io.VTXWriter(sz.mesh.comm, filename, [sz.T_i, sz.vs_i, sz.vw_i]) as vtx:
    vtx.write(0.0)
# zip the .bp folder so that it can be downloaded from Jupyter lab
if __name__ == "__main__" and "__file__" not in globals():
    zipfilename = filename.with_suffix(".zip")
    !zip -r $zipfilename $filename

  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/profiling.json (deflated 60%)
  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/data.0

 (deflated 64%)
  adding: output/40_S_Marianas_ss_solution_resscale_3.00.bp/md.0 (deflated 76%)

Comparison

Compare to the published result from Wilson & van Keken, PEPS, 2023 (II) and van Keken & Wilson, PEPS, 2023 (III). The original models used in these papers are also available as open-source repositories on github and zenodo.

First download the minimal necessary data from zenodo and check it is the right version.

zipfilename = pathlib.Path(os.path.join(basedir, os.path.pardir, os.path.pardir, "data", "vankeken_wilson_peps_2023_TF_lowres_minimal.zip"))
if not zipfilename.is_file():
    zipfileurl = 'https://zenodo.org/records/13234021/files/vankeken_wilson_peps_2023_TF_lowres_minimal.zip'
    r = requests.get(zipfileurl, allow_redirects=True)
    open(zipfilename, 'wb').write(r.content)
assert hashlib.md5(open(zipfilename, 'rb').read()).hexdigest() == 'a8eca6220f9bee091e41a680d502fe0d'

tffilename = os.path.join('vankeken_wilson_peps_2023_TF_lowres_minimal', 'sz_suite_ss', szdict['dirname']+'_minres_2.00.vtu')
tffilepath = os.path.join(os.pardir, os.pardir, 'data')
with zipfile.ZipFile(zipfilename, 'r') as z:
    z.extract(tffilename, path=tffilepath)

fxgrid = utils.grids_scalar(sz.T_i)[0]

tfgrid = pv.get_reader(os.path.join(tffilepath, tffilename)).read()

diffgrid = utils.pv_diff(fxgrid, tfgrid, field_name_map={'T':'Temperature::PotentialTemperature'}, pass_point_data=True)

diffgrid.set_active_scalars('T')
plotter_diff = pv.Plotter()
clim = None
plotter_diff.add_mesh(diffgrid, cmap='coolwarm', clim=clim, scalar_bar_args={'title': 'Temperature Difference (deg C)', 'bold':True})
utils.plot_geometry(geom, plotter=plotter_diff, color='green', width=2)
utils.plot_couplingdepth(slab, plotter=plotter_diff, render_points_as_spheres=True, point_size=5.0, color='green')
plotter_diff.enable_parallel_projection()
plotter_diff.view_xy()
plotter_diff.show()

integrated_data = diffgrid.integrate_data()
error = integrated_data['T'][0]/integrated_data['Area'][0]
print("Average error = {}".format(error,))
assert np.abs(error) < 5

Average error = 0.550617734830479

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 40_S_Marianas.ipynb

[NbConvertApp] Converting notebook 40_S_Marianas.ipynb to script

[NbConvertApp] Writing 7034 bytes to 40_S_Marianas.py