21 Scotia

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 = "21_Scotia"
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))

21_Scotia:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       135
extra_width          23
lc_depth             30
io_depth             166
dirname              21_Scotia
As                   20.0
Ac                   22.0
A                    59.1
sztype               oceanic
Vs                   60.8
xs                   [0, 34.8, 75.5, 134.0, 141.2, 142.9, 153.0, 176.0, 217.0]
ys                   [-6, -15.0, -30.0, -70.0, -80.0, -82.5, -100.0, -150.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           3481.17      1           
1           407.837      0.117155    

2           201.429      0.0578626   
3           75.4074      0.0216615   

4           30.5725      0.00878224  
5           15.5055      0.0044541   
6           9.67147      0.00277822  

7           6.64093      0.00190767  
8           4.72215      0.00135648  

9           3.41029      0.00097964  
10          2.49192      0.00071583  
11          1.8427       0.000529333 

12          1.37948      0.000396269 
13          1.04519      0.000300243 

14          0.800958     0.000230083 
15          0.620034     0.000178111 

16          0.484022     0.00013904  
17          0.380377     0.000109267 
18          0.300462     8.63108e-05 

19          0.238198     6.84246e-05 
20          0.189339     5.43896e-05 

21          0.150798     4.33183e-05 
22          0.120274     3.455e-05   

23          0.0960295    2.75854e-05 
24          0.0767343    2.20427e-05 

25          0.0613591    1.7626e-05  
26          0.0490957    1.41032e-05 

27          0.0393081    1.12916e-05 
28          0.0314927    9.0466e-06  
29          0.0252496    7.2532e-06  

30          0.0202602    5.81993e-06 
31          0.0162705    4.67385e-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/21_Scotia_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/21_Scotia_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/21_Scotia_ss_solution_resscale_3.00.bp/profiling.json (deflated 62%)
  adding: output/21_Scotia_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/21_Scotia_ss_solution_resscale_3.00.bp/md.0 (deflated 75%)
  adding: output/21_Scotia_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.619307611645316