36 Kermadec

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 = "36_Kermadec"
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))

36_Kermadec:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       185
extra_width          18.5
lc_depth             30
io_depth             170
dirname              36_Kermadec
As                   20.0
Ac                   45.0
A                    100.0
sztype               oceanic
Vs                   64.6
xs                   [0, 24.2, 60.4, 96.0, 117.7, 119.2, 129.0, 180.0, 201.5]
ys                   [-6, -15.0, -30.0, -50.0, -80.0, -82.5, -100.0, -200.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           9642.89      1           
1           873.086      0.0905419   

2           258.092      0.0267649   
3           156.868      0.0162677   
4           79.1667      0.00820984  

5           38.6126      0.00400425  
6           21.3806      0.00221724  
7           13.4904      0.001399    

8           9.22918      0.000957096 
9           6.60939      0.000685415 
10          4.87466      0.000505518 

11          3.67343      0.000380947 
12          2.8142       0.000291842 
13          2.18345      0.000226431 

14          1.71025      0.000177359 
15          1.34876      0.000139871 
16          1.06863      0.00011082  

17          0.849186     8.80634e-05 
18          0.67595      7.00983e-05 
19          0.538477     5.58419e-05 

20          0.429026     4.44914e-05 
21          0.341732     3.54388e-05 
22          0.272063     2.82139e-05 

23          0.216461     2.24477e-05 
24          0.172105     1.78478e-05 
25          0.136744     1.41808e-05 

26          0.108577     1.12598e-05 
27          0.08616      8.93508e-06 
28          0.068335     7.08656e-06 

29          0.0541727    5.61789e-06 
30          0.0429291    4.45189e-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/36_Kermadec_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/36_Kermadec_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/36_Kermadec_ss_solution_resscale_3.00.bp/profiling.json (deflated 62%)
  adding: output/36_Kermadec_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/36_Kermadec_ss_solution_resscale_3.00.bp/md.0 (deflated 76%)
  adding: output/36_Kermadec_ss_solution_resscale_3.00.bp/data.0 (deflated 63%)

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 = 1.4672748644251876