27 Java

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 = "27_Java"
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))

27_Java:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       270
extra_width          19
lc_depth             33.5
io_depth             184
uc_depth             15
dirname              27_Java
As                   40.0
qs                   0.065
A                    100.0
sztype               continental
Vs                   67.8
xs                   [0, 59.2, 142.0, 217.0, 229.5, 232.4, 251.0, 337.0, 371.0]
ys                   [-6, -15.0, -33.5, -70.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           11685.2      1           

1           4594.82      0.393216    
2           242.827      0.0207807   

3           99.1234      0.0084828   
4           44.034       0.00376835  

5           18.7917      0.00160816  
6           9.55297      0.000817525 

7           6.08937      0.000521117 
8           4.34096      0.000371491 

9           3.23025      0.000276438 
10          2.43387      0.000208286 

11          1.83614      0.000157133 
12          1.38193      0.000118263 

13          1.03734      8.8774e-05  
14          0.777357     6.65247e-05 

15          0.582113     4.98161e-05 
16          0.435886     3.73023e-05 

17          0.326545     2.79451e-05 
18          0.244871     2.09556e-05 

19          0.183907     1.57384e-05 
20          0.138418     1.18455e-05 

21          0.104475     8.94077e-06 
22          0.0791375    6.77244e-06 

23          0.0602085    5.15253e-06 
24          0.0460483    3.94073e-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/27_Java_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/27_Java_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/27_Java_ss_solution_resscale_3.00.bp/profiling.json (deflated 62%)
  adding: output/27_Java_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/27_Java_ss_solution_resscale_3.00.bp/md.0 (deflated 76%)
  adding: output/27_Java_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 = 0.8778995633121633