48 N Honshu

Authors: Cameron Seebeck, Cian Wilson

Time-dependent 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, 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_tdep_dislcreep import TDDislSubductionProblem
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 and target Courant number cfl of the model.

Resolution

By default the resolution (both spatial and temporal) is low to allow for a quick runtime and smaller website size. If sufficient computational resources are available set a lower resscale and a lower cfl to get higher spatial and temporal resolutions respectively. This is necessary to get results with sufficient accuracy for scientific interpretation.

name = "48_N_Honshu"
resscale = 3.0
cfl      = 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))

48_N_Honshu:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       280
extra_width          22
lc_depth             40
io_depth             195
uc_depth             15
dirname              48_N_Honshu
As                   40.0
qs                   0.065
A                    100.0
sztype               continental
Vs                   81.6
xs                   [0, 40.4, 166.0, 248.0, 265.8, 269.9, 298.0, 425.0, 538.0]
ys                   [-6, -15.0, -40.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 TDDislSubductionProblem class.

TDDislSubductionProblem?

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 TDDislSubductionProblem problem class using the dictionary of parameters.

sz = TDDislSubductionProblem(geom, **szdict)

Solve

Solve using a dislocation creep rheology.

# Select the timestep based on the approximate target Courant number
dt = cfl*resscale/szdict['Vs']
# Reduce the timestep to get an integer number of timesteps
dt = szdict['As']/np.ceil(szdict['As']/dt)
sz.solve(szdict['As'], dt, theta=0.5, rtol=1.e-1, verbosity=1)

Entering timeloop with 363 steps (dt = 0.110193 Myr, final time = 40 Myr)

Finished timeloop after 363 steps (final time = 40 Myr)

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 / "{}_td_solution_resscale_{:.2f}_cfl_{:.2f}.png".format(name, resscale, cfl,))

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

filename = output_folder / "{}_td_solution_resscale_{:.2f}_cfl_{:.2f}.bp".format(name, resscale, cfl,)
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/48_N_Honshu_td_solution_resscale_3.00_cfl_3.00.bp/ (stored 0%)
  adding: output/48_N_Honshu_td_solution_resscale_3.00_cfl_3.00.bp/mmd.0 (deflated 64%)
  adding: output/48_N_Honshu_td_solution_resscale_3.00_cfl_3.00.bp/profiling.json (deflated 62%)
  adding: output/48_N_Honshu_td_solution_resscale_3.00_cfl_3.00.bp/md.idx (deflated 59%)
  adding: output/48_N_Honshu_td_solution_resscale_3.00_cfl_3.00.bp/md.0 (deflated 75%)
  adding: output/48_N_Honshu_td_solution_resscale_3.00_cfl_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_td', szdict['dirname']+'_minres_2.00_cfl_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.7526821847454038