05 Mexico

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, 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_problems.sz_params import allsz_params
from sz_problems.sz_slab import create_slab, plot_slab
from sz_problems.sz_geometry import create_sz_geometry
from sz_problems.sz_steady_dislcreep import SteadyDislSubductionProblem
import numpy as np
import dolfinx as df
import pyvista as pv
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 = "05_Mexico"
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))

05_Mexico:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       80
extra_width          20
lc_depth             45
io_depth             172
uc_depth             15
dirname              05_Mexico
As                   40.0
qs                   0.065
A                    10
sztype               continental
Vs                   47.0
xs                   [0, 69.5, 139.0, 268.0, 292.0, 334.9, 337.0, 350.0, 416.0]
ys                   [-6, -15.0, -35.0, -45.0, -50.0, -80.0, -82.5, -100.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 SteadyDislSubductionProblem class.

SteadyDislSubductionProblem?

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 = SteadyDislSubductionProblem(geom, **szdict)

Solve

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

sz.solve()

Iteration   Residual     Relative Residual
------------------------------------------
0           781.594      1           

1           72.9536      0.0933394   
2           13.997       0.0179083   

3           9.01409      0.011533    
4           4.9923       0.00638732  

5           2.743        0.0035095   
6           1.60259      0.00205041  

7           0.96979      0.00124079  
8           0.608755     0.000778864 

9           0.393506     0.000503466 
10          0.260941     0.000333857 

11          0.177086     0.00022657  
12          0.122866     0.0001572   

13          0.0870112    0.000111325 
14          0.062717     8.02424e-05 

15          0.0458583    5.86728e-05 
16          0.0339082    4.33834e-05 

17          0.0252878    3.23541e-05 
18          0.0189833    2.4288e-05  

19          0.0143234    1.83259e-05 
20          0.010851     1.38832e-05 

21          0.00824733   1.05519e-05 
22          0.00628541   8.04179e-06 

23          0.00480127   6.14292e-06 
24          0.00367488   4.70178e-06 

Plot

Plot the solution.

plotter = pv.Plotter()
utils.plot.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.plot_vector_glyphs(sz.vw_i, plotter=plotter, gather=True, factor=0.1, color='k', scale=utils.mps_to_mmpyr(sz.v0))
utils.plot.plot_vector_glyphs(sz.vs_i, plotter=plotter, gather=True, factor=0.1, color='k', scale=utils.mps_to_mmpyr(sz.v0))
utils.plot.plot_geometry(geom, plotter=plotter, color='green', width=2)
utils.plot.plot_couplingdepth(slab, plotter=plotter, render_points_as_spheres=True, point_size=10.0, color='green')
utils.plot.plot_show(plotter)
utils.plot.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
zipfilename = filename.with_suffix(".zip")
!zip -r $zipfilename $filename

  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/profiling.json (deflated 62%)
  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/md.0 (deflated 76%)
  adding: output/05_Mexico_ss_solution_resscale_3.00.bp/data.0

 (deflated 64%)

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, 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(basedir, os.path.pardir, os.path.pardir, os.path.pardir, 'data')
with zipfile.ZipFile(zipfilename, 'r') as z:
    z.extract(tffilename, path=tffilepath)

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

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

diffgrid = utils.plot.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.plot_geometry(geom, plotter=plotter_diff, color='green', width=2)
utils.plot.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.16470986759196654