Steady-State Subduction Zone Setup#
Themes and variations - global suite geometries#
In this notebook we will try using more realistic geometries from the global suite.
Preamble#
Let’s start by adding the path to the modules in the python folder to the system path (so we can find the our custom modules).
import sys, os
basedir = ''
if "__file__" in globals(): basedir = os.path.dirname(__file__)
sys.path.insert(0, os.path.join(basedir, os.path.pardir, os.path.pardir, 'python'))
Then load everything we need from sz_problems and other modules.
import fenics_sz.utils
from fenics_sz.sz_problems.sz_params import default_params, allsz_params
from fenics_sz.sz_problems.sz_slab import create_slab
from fenics_sz.sz_problems.sz_geometry import create_sz_geometry
from fenics_sz.sz_problems.sz_steady_isoviscous import SteadyIsoSubductionProblem
from fenics_sz.sz_problems.sz_steady_dislcreep import SteadyDislSubductionProblem
import pathlib
output_folder = pathlib.Path(os.path.join(basedir, "output"))
output_folder.mkdir(exist_ok=True, parents=True)
Alaska Peninsula (dislocation creep, low res)#
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.
resscale_ak = 5.0
szdict_ak = allsz_params['01_Alaska_Peninsula']
slab_ak = create_slab(szdict_ak['xs'], szdict_ak['ys'], resscale_ak, szdict_ak['lc_depth'])
geom_ak = create_sz_geometry(slab_ak, resscale_ak, szdict_ak['sztype'], szdict_ak['io_depth'], szdict_ak['extra_width'],
szdict_ak['coast_distance'], szdict_ak['lc_depth'], szdict_ak['uc_depth'])
sz_ak = SteadyDislSubductionProblem(geom_ak, **szdict_ak)
print("\nSolving steady state flow with isoviscous rheology...")
sz_ak.solve()
Solving steady state flow with isoviscous rheology...
Iteration Residual Relative Residual
------------------------------------------
0 7125.26 1
1 811.298 0.113862
2 108.535 0.0152324
3 60.4568 0.00848486
4 26.631 0.00373754
5 10.465 0.00146871
6 4.44722 0.000624148
7 2.16681 0.000304102
8 1.1928 0.000167404
9 0.719893 0.000101034
10 0.467766 6.5649e-05
11 0.3237 4.54299e-05
12 0.236329 3.31677e-05
13 0.180352 2.53116e-05
14 0.142387 1.99834e-05
15 0.115173 1.6164e-05
16 0.0946823 1.32883e-05
17 0.0786273 1.1035e-05
18 0.06567 9.21651e-06
19 0.0549949 7.7183e-06
20 0.0460837 6.46764e-06
21 0.0385896 5.41588e-06
22 0.0322656 4.52834e-06
plotter_ak = fenics_sz.utils.plot.plot_scalar(sz_ak.T_i, scale=sz_ak.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})
fenics_sz.utils.plot.plot_vector_glyphs(sz_ak.vw_i, plotter=plotter_ak, gather=True, factor=0.1, color='k', scale=fenics_sz.utils.mps_to_mmpyr(sz_ak.v0))
fenics_sz.utils.plot.plot_vector_glyphs(sz_ak.vs_i, plotter=plotter_ak, gather=True, factor=0.1, color='k', scale=fenics_sz.utils.mps_to_mmpyr(sz_ak.v0))
sz_ak.geom.pyvistaplot(plotter=plotter_ak, color='green', width=2)
cdpt = sz_ak.geom.slab_spline.findpoint('Slab::FullCouplingDepth')
fenics_sz.utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter_ak, render_points_as_spheres=True, point_size=10.0, color='green')
fenics_sz.utils.plot.plot_show(plotter_ak)
fenics_sz.utils.plot.plot_save(plotter_ak, output_folder / "sz_steady_tests_ak_solution.png")
eta_ak = sz_ak.project_dislocationcreep_viscosity()
plotter_eta_ak = fenics_sz.utils.plot.plot_scalar(eta_ak, scale=sz_ak.eta0, gather=True, log_scale=True, show_edges=True, scalar_bar_args={'title': 'Viscosity (Pa s)', 'bold':True})
sz_ak.geom.pyvistaplot(plotter=plotter_eta_ak, color='green', width=2)
cdpt = sz_ak.geom.slab_spline.findpoint('Slab::FullCouplingDepth')
fenics_sz.utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter_eta_ak, render_points_as_spheres=True, point_size=10.0, color='green')
fenics_sz.utils.plot.plot_show(plotter_eta_ak)
fenics_sz.utils.plot.plot_save(plotter_eta_ak, output_folder / "sz_steady_tests_ak_eta.png")
# clean up
del plotter_eta_ak
del sz_ak
del geom_ak
del slab_ak
N Antilles (dislocation creep, low res)#
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.
resscale_ant = 5.0
szdict_ant = allsz_params['19_N_Antilles']
slab_ant = create_slab(szdict_ant['xs'], szdict_ant['ys'], resscale_ant, szdict_ant['lc_depth'])
geom_ant = create_sz_geometry(slab_ant, resscale_ant, szdict_ant['sztype'], szdict_ant['io_depth'], szdict_ant['extra_width'],
szdict_ant['coast_distance'], szdict_ant['lc_depth'], szdict_ant['uc_depth'])
sz_ant = SteadyDislSubductionProblem(geom_ant, **szdict_ant)
print("\nSolving steady state flow with isoviscous rheology...")
sz_ant.solve()
Solving steady state flow with isoviscous rheology...
Iteration Residual Relative Residual
------------------------------------------
0 4821.23 1
1 2813.69 0.583604
2 233.101 0.0483489
3 88.591 0.0183752
4 38.3363 0.00795156
5 19.2591 0.00399464
6 10.9435 0.00226986
7 6.64963 0.00137924
8 4.35547 0.000903393
9 3.09337 0.000641614
10 2.33995 0.000485342
11 1.83278 0.000380147
12 1.45825 0.000302465
13 1.16782 0.000242225
14 0.940906 0.000195159
15 0.757093 0.000157033
16 0.604835 0.000125452
17 0.479977 9.95548e-05
18 0.379908 7.87989e-05
19 0.299543 6.213e-05
20 0.235121 4.87677e-05
21 0.183763 3.81154e-05
22 0.143067 2.96743e-05
23 0.110997 2.30225e-05
24 0.0858494 1.78065e-05
25 0.066217 1.37344e-05
26 0.0509492 1.05677e-05
27 0.0391163 8.11333e-06
28 0.0299733 6.21693e-06
29 0.0229277 4.75557e-06
plotter_ant = fenics_sz.utils.plot.plot_scalar(sz_ant.T_i, scale=sz_ant.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})
fenics_sz.utils.plot.plot_vector_glyphs(sz_ant.vw_i, plotter=plotter_ant, gather=True, factor=0.25, color='k', scale=fenics_sz.utils.mps_to_mmpyr(sz_ant.v0))
fenics_sz.utils.plot.plot_vector_glyphs(sz_ant.vs_i, plotter=plotter_ant, gather=True, factor=0.25, color='k', scale=fenics_sz.utils.mps_to_mmpyr(sz_ant.v0))
sz_ant.geom.pyvistaplot(plotter=plotter_ant, color='green', width=2)
cdpt = sz_ant.geom.slab_spline.findpoint('Slab::FullCouplingDepth')
fenics_sz.utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter_ant, render_points_as_spheres=True, point_size=10.0, color='green')
fenics_sz.utils.plot.plot_show(plotter_ant)
fenics_sz.utils.plot.plot_save(plotter_ant, output_folder / "sz_steady_tests_ant_solution.png")
eta_ant = sz_ant.project_dislocationcreep_viscosity()
plotter_eta_ant = fenics_sz.utils.plot.plot_scalar(eta_ant, scale=sz_ant.eta0, gather=True, log_scale=True, show_edges=True, scalar_bar_args={'title': 'Viscosity (Pa s)', 'bold':True})
sz_ant.geom.pyvistaplot(plotter=plotter_eta_ant, color='green', width=2)
cdpt = sz_ant.geom.slab_spline.findpoint('Slab::FullCouplingDepth')
fenics_sz.utils.plot.plot_points([[cdpt.x, cdpt.y, 0.0]], plotter=plotter_eta_ant, render_points_as_spheres=True, point_size=10.0, color='green')
fenics_sz.utils.plot.plot_show(plotter_eta_ant)
fenics_sz.utils.plot.plot_save(plotter_eta_ant, output_folder / "sz_steady_tests_ant_eta.png")
# clean up
del plotter_eta_ant
del sz_ant
del geom_ant
del slab_ant
Having played with some test examples using a steady state solver we will next formally examine the published benchmark solutions and test how our implementation performs.