14 NC Chile

14 NC Chile#

Authors: Kidus Teshome, Cian Wilson

Steady state implementation#

Preamble#

We are going to run this notebook in parallel using ipyparallel. To do this we launch a parallel engine with nprocs processes. Following this, all cells we wish to run in parallel need to start with the Jupyter magic %%px.

%%px

Cells without the %%px magic will run locally and will not have access to anything loaded on the parallel engine.

import ipyparallel as ipp
nprocs = 2
rc = ipp.Cluster(engine_launcher_class="mpi", n=nprocs).start_and_connect_sync()

Starting 2 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>

Set some path information.

%%px
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.

%%px
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
from mpi4py import MPI
comm = MPI.COMM_WORLD
my_rank = comm.rank

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.

%%px
name = "14_NC_Chile"
resscale = 3.0

Then load the remaining parameters from the global suite.

%%px
szdict = allsz_params[name]
if my_rank == 0:
    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))

[stdout:0] 14_NC_Chile:
Key                  Value     
-------------------------------------------------------------------------------------
coast_distance       120
extra_width          20
lc_depth             55
io_depth             190
uc_depth             15
dirname              14_NC_Chile
As                   40.0
qs                   0.07999999999999999
A                    42.8
sztype               continental
Vs                   77.4
xs                   [0, 27.9, 129.0, 171.0, 197.0, 203.2, 245.0, 364.0, 576.0]
ys                   [-6, -15.0, -55.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.

%%px
if my_rank == 0: help(SteadyDislSubductionProblem.__init__)

[stdout:0] Help on function __init__ in module sz_problems.sz_base:

__init__(self, geom, **kwargs)
    Initialize a BaseSubductionProblem.

    Arguments:
      * geom  - an instance of a subduction zone geometry

    Keyword Arguments:
     required:
      * A      - age of subducting slab (in Myr) [required]
      * Vs     - incoming slab speed (in mm/yr) [required]
      * sztype - type of subduction zone (either 'continental' or 'oceanic') [required]
      * Ac     - age of the over-riding plate (in Myr) [required if sztype is 'oceanic']
      * As     - age of subduction (in Myr) [required if sztype is 'oceanic']
      * qs     - surface heat flux (in W/m^2) [required if sztype is 'continental']

     optional:
      * Ts   - surface temperature (deg C, corresponds to non-dim)
      * Tm   - mantle temperature (deg C, corresponds to non-dim)
      * kc   - crustal thermal conductivity (non-dim) [only has an effect if sztype is 'continental']
      * km   - mantle thermal conductivity (non-dim)
      * rhoc - crustal density (non-dim) [only has an effect if sztype is 'continental']
      * rhom - mantle density (non-dim)
      * cp   - isobaric heat capacity (non-dim)
      * H1   - upper crustal volumetric heat production (non-dim) [only has an effect if sztype is 'continental']
      * H2   - lower crustal volumetric heat production (non-dim) [only has an effect if sztype is 'continental']

     optional (dislocation creep rheology):
      * etamax - maximum viscosity (Pas) [only relevant for dislocation creep rheologies]
      * nsigma - stress viscosity power law exponents (non-dim) [only relevant for dislocation creep rheologies]
      * Aeta   - pre-exponential viscosity constant (Pa s^(1/n)) [only relevant for dislocation creep rheologies]
      * E      - viscosity activation energy (J/mol) [only relevant for dislocation creep rheologies]

Setup#

Setup a slab.

%%px
slab = create_slab(szdict['xs'], szdict['ys'], resscale, szdict['lc_depth'])
if my_rank == 0: _ = plot_slab(slab)

[output:0]
../../../_images/09ac70c88b8cc9a6ba871988803676f6fb42451830a2bddf150e0e1070203b3d.png

Create the subduction zome geometry around the slab.

%%px
geom = create_sz_geometry(slab, resscale, szdict['sztype'], szdict['io_depth'], szdict['extra_width'], 
                             szdict['coast_distance'], szdict['lc_depth'], szdict['uc_depth'])
if my_rank == 0: _ = geom.plot()

[output:0]
../../../_images/76e95a9fc362e8b1f18351777c6b410692a4b0fbec75e58c0cc54c19f1d52da7.png

Finally, declare the SubductionZone problem class using the dictionary of parameters.

%%px
sz = SteadyDislSubductionProblem(geom, **szdict)

Solve#

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

%%px
sz.solve()

[stdout:0] Iteration   Residual     Relative Residual
------------------------------------------
0           2582.26      1           
1           302.649      0.117203    
2           94.121       0.0364491   
3           41.1085      0.0159196   
4           21.531       0.00833806  
5           12.5361      0.00485469  
6           7.43818      0.0028805   
7           4.47803      0.00173415  
8           2.74031      0.00106121  
9           1.70162      0.000658966 
10          1.07275      0.00041543  
11          0.69076      0.000267502 
12          0.459276     0.000177858 
13          0.320986     0.000124305 
14          0.239856     9.28863e-05 
15          0.192053     7.4374e-05  
16          0.162236     6.28272e-05 
17          0.141632     5.4848e-05  
18          0.125796     4.87157e-05 
19          0.112595     4.36033e-05 
20          0.101031     3.91251e-05 
21          0.0906324    3.50982e-05 
22          0.0811654    3.1432e-05  
23          0.0725126    2.80811e-05 
24          0.0646067    2.50195e-05 
25          0.0574       2.22286e-05 
26          0.0508529    1.96932e-05 
27          0.0449279    1.73987e-05 
28          0.039587     1.53304e-05 
29          0.0347915    1.34733e-05 
30          0.0305021    1.18122e-05 
31          0.0266793    1.03318e-05 
32          0.0232841    9.01697e-06 
33          0.0202787    7.85308e-06 
34          0.0176264    6.82595e-06 
35          0.0152924    5.92212e-06 
36          0.0132442    5.12892e-06 
37          0.0114512    4.43457e-06

Plot#

Plot the solution.

%%px
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))
geom.pyvistaplot(plotter=plotter, color='green', width=2)
cdpt = slab.findpoint('Slab::FullCouplingDepth')
utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter, render_points_as_spheres=True, point_size=10.0, color='green')
if my_rank == 0: 
    utils.plot.plot_show(plotter)
    utils.plot.plot_save(plotter, output_folder / "{}_ss_solution_resscale_{:.2f}.png".format(name, resscale))

[output:0]

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

%%px
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

[stdout:1]   adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/data.0 (deflated 66%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/md.0 (deflated 85%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/profiling.json (deflated 72%)

[stdout:0]   adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/ (stored 0%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/data.0 (deflated 66%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/md.idx (deflated 59%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/mmd.0 (deflated 64%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/md.0 (deflated 85%)
  adding: output/14_NC_Chile_ss_solution_resscale_3.00.bp/profiling.json (deflated 72%)

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.

%%px
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"))
# only one process should download the data
if my_rank == 0:
    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)
# wait until rank 0 has downloaded the data
comm.barrier()
assert hashlib.md5(open(zipfilename, 'rb').read()).hexdigest() == 'a8eca6220f9bee091e41a680d502fe0d'

%%px
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')
# one process should extract the file
if my_rank == 0:
    with zipfile.ZipFile(zipfilename, 'r') as z:
        z.extract(tffilename, path=tffilepath)
# other processes should wait
comm.barrier()

%%px
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)

%%px
# first gather the data onto rank 0
diffgrid_g = utils.plot.pyvista_grids(diffgrid.cells, diffgrid.celltypes, diffgrid.points, comm, gather=True)
T_g = comm.gather(diffgrid.point_data['T'], root=0)
for r, grid in enumerate(diffgrid_g):
    grid.point_data['T'] = T_g[r]
    grid.set_active_scalars('T')
    grid.clean(tolerance=1.e-2)
diffgrid_g

# then plot it
plotter_diff = pv.Plotter()
clim = None
for grid in diffgrid_g: 
    plotter_diff.add_mesh(grid, cmap='coolwarm', clim=clim, scalar_bar_args={'title': 'Temperature Difference (deg C)', 'bold':True})
geom.pyvistaplot(plotter=plotter_diff, color='green', width=2)
cdpt = slab.findpoint('Slab::FullCouplingDepth')
#utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter_diff, render_points_as_spheres=True, point_size=10.0, color='green')
plotter_diff.enable_parallel_projection()
plotter_diff.view_xy()
if my_rank == 0: plotter_diff.show()

[output:0]

%%px
integrated_data = diffgrid.integrate_data()
totalarea = comm.allreduce(integrated_data['Area'][0], op=MPI.SUM)
error = comm.allreduce(integrated_data['T'][0], op=MPI.SUM)/totalarea
if my_rank == 0: print("Average error = {}".format(error,))
assert np.abs(error) < 5

[stdout:0] Average error = 0.3757773115207858
Contents