Add Gmsh Mesh Generation for subduction-3d Example

developer/format_source.py

@@ -219,7 +219,7 @@ def _parse_command_line(self):
         parser.add_argument("--python", action="store", dest="python",
                             default=None, help="[None, 'all', FILE, PATTERN]")
         parser.add_argument("--config-dir", action="store", dest="config_dir",
-                            default="doc/developer", help="Directory containing config files for formatters.")
+                            default="developer", help="Directory containing config files for formatters.")
         parser.add_argument("--quiet", action="store_false",
                             dest="show_progress", default=True)
         parser.add_argument("--debug", action="store_true",

examples/subduction-3d/download_topography.sh

@@ -0,0 +1,13 @@
+mkdir -p topography
+cd topography
+
+#Download
+curl -L -o 50_120.zip https://dds.cr.usgs.gov/download/eyJpZCI6ODEzODA3MzE3LCJjb250YWN0SWQiOjI3NTc1ODA4fQ==/
+curl -L -o 50_150.zip https://dds.cr.usgs.gov/download/eyJpZCI6ODEzODA4MTEwLCJjb250YWN0SWQiOjI3NTc1ODA4fQ==/
+curl -L -o 30_150.zip https://dds.cr.usgs.gov/download/eyJpZCI6ODEzNzY2NDk0LCJjb250YWN0SWQiOjI3NTc1ODA4fQ==/
+curl -L -o 30_120.zip https://dds.cr.usgs.gov/download/eyJpZCI6ODEzNzQ0ODMwLCJjb250YWN0SWQiOjI3NTc1ODA4fQ==/
+
+unzip 50_120.zip -d tifs_50_120
+unzip 50_150.zip -d tifs_50_150
+unzip 30_150.zip -d tifs_30_150
+unzip 30_120.zip -d tifs_30_120

examples/subduction-3d/generate_localdem.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+"""Python script to extract elevation data in region around Cascadia Subduction Zone using GDAL Python API."""
+from osgeo import gdal
+
+FILENAME_DEMS = [
+    "topography/tifs_30_150/30n150w_20101117_gmted_med300.tif",
+    "topography/tifs_50_150/50n150w_20101117_gmted_med300.tif",
+    "topography/tifs_30_120/30n120w_20101117_gmted_med300.tif",
+    "topography/tifs_50_120/50n120w_20101117_gmted_med300.tif",
+]
+FILENAME_DEMSUB = "topography.nc"
+BBOX = (-130.195718069, 51.5863472565, -115.821471093, 39.0650222389)
+RESOLUTION = 0.2  # degrees
+
+vrt_path = "/vsimem/merged.vrt"
+gdal.BuildVRT(vrt_path, FILENAME_DEMS)
+
+gdal.UseExceptions()
+options = gdal.TranslateOptions(
+    projWin=BBOX,
+    xRes=RESOLUTION,
+    yRes=RESOLUTION,
+    creationOptions=["COMPRESS=DEFLATE"],
+    resampleAlg="cubic")
+gdal.Translate(FILENAME_DEMSUB, vrt_path, options=options)
+
+gdal.Unlink(vrt_path)

examples/subduction-3d/pylithapp.cfg

@@ -25,6 +25,7 @@ pylithapp = 1
 timedependent = 1
 solution = 1
 meshiocubit = 1
+meshiopetsc = 1
 isotropiclinearelasticity = 1
 dirichlettimedependent = 1
 faultcohesivekin = 1

examples/subduction-3d/step01_axialdisp_gmsh.cfg

@@ -0,0 +1,33 @@
+[pylithapp.metadata]
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = Axial compression in the east-west direction using Dirichlet boundary conditions.
+keywords = [axial compression]
+arguments = [step01_axialdisp_gmsh.cfg, step01_axialdisp.cfg, mat_elastic.cfg]
+version = 2.0.0
+
+features = [
+    Static simulation,
+    pylith.materials.IsotropicLinearElasticity,
+    spatialdata.spatialdb.SimpleDB
+    ]
+
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step01_axialdisp_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+bc_xneg.label_value = 11
+bc_xpos.label_value = 13
+bc_yneg.label_value = 10
+bc_ypos.label_value = 12
+bc_zneg.label_value = 14
+
+# End of file

examples/subduction-3d/step02_coseismic_gmsh.cfg

@@ -0,0 +1,42 @@
+[pylithapp.metadata]
+# Materials
+#   The 'wedge' and 'crust' blocks use the the linear, isotropic
+#   elastic material model, whereas the 'slab' and 'mantle' use the
+#   linear Maxwell viscoelastic material model (mat_viscoelastic.cfg).
+#
+# Boundary Conditions
+#   Roller boundary conditions (pylithapp.cfg).
+#
+# Fault
+#   Uniform oblique slip on the slab_top fault rupture patch.
+#
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = Postseismic viscoelastic relaxation from coseismic slip on a patch in the middle of the top of the slab.
+keywords = [prescribed slip, coseismic slip]
+arguments = [step02_coseismic_gmsh.cfg, step02_coseismic.cfg, mat_viscoelastic.cfg]
+version = 2.0.0
+
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step02_coseismic_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.problem.interfaces.fault_slabtop]
+label_value = 24
+edge_value = 25
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+bc_xneg.label_value = 11
+bc_xpos.label_value = 13
+bc_yneg.label_value = 10
+bc_ypos.label_value = 12
+bc_zneg.label_value = 14
+
+# End of file

examples/subduction-3d/step03_interseismic_gmsh.cfg

@@ -0,0 +1,52 @@
+[pylithapp.metadata]
+# Materials
+#   The 'wedge' and 'crust' blocks use the the linear, isotropic
+#   elastic material model, whereas the 'slab' and 'mantle' use the
+#   linear Maxwell viscoelastic material model (mat_viscoelastic.cfg).
+#
+# Boundary Conditions
+#   Roller boundary conditions (pylithapp.cfg).
+#
+# Fault
+#   Uniform creep on the bottom of the slab.
+#   Uniform creep on the deeper portion of the subduction interface.
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = Interseismic deformation, combining aseismic slip on the top and bottom of the subducting slab with viscoelastic relaxation in the mantle and deep part of the slab.
+keywords = [prescribed slip, aseismic creep]
+arguments = [step03_interseismic_gmsh.cfg, step03_interseismic.cfg, mat_viscoelastic.cfg]
+version = 2.0.0
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step03_interseismic_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.problem.interfaces.fault_slabbot]
+label_value = 21
+edge_value = 23
+
+[pylithapp.problem.interfaces.fault_slabtop]
+label_value = 20
+edge_value = 22
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+# boundary_xneg_noslab
+bc_xneg.label_value = 16
+bc_xpos.label_value = 13
+# boundary_yneg_noslab
+bc_yneg.label_value = 17
+# boundary_ypos_noslab
+bc_ypos.label_value = 18
+bc_zneg.label_value = 14
+
+# End of file

examples/subduction-3d/step04_eqcycle_gmsh.cfg

@@ -0,0 +1,55 @@
+[pylithapp.metadata]
+# Materials
+#   The 'wedge' and 'crust' blocks use the the linear, isotropic
+#   elastic material model, whereas the 'slab' and 'mantle' use the
+#   linear Maxwell viscoelastic material model (mat_viscoelastic.cfg).
+#
+# Boundary Conditions
+#   Roller boundary conditions (pylithapp.cfg).
+#
+# Fault
+#   Uniform creep on the bottom of the slab.
+#   Uniform creep on the deeper portion of the subduction interface.
+#   + Coseismic slip on the shallow portion of the subduction interface.
+#   + Coseismic slip on the splay fault.
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = Coseismic and interseismic deformation, combining earthquake ruptures and aseismic slip on the top and bottom of the subducting slab with viscoelastic relaxation in the mantle and deep part of the slab.
+keywords = [prescribed slip, aseismic creep, coseismic slip]
+arguments = [step04_eqcycle_gmsh.cfg, step04_eqcycle.cfg, mat_viscoelastic.cfg]
+version = 2.0.0
+
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step04_eqcycle_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.problem.interfaces.fault_slabbot]
+label_value = 21
+edge_value = 23
+
+[pylithapp.problem.interfaces.fault_slabtop]
+label_value = 20
+edge_value = 22
+
+[pylithapp.problem.interfaces.splay]
+label_value = 30
+edge_value = 31
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+# boundary_xneg_noslab
+bc_xneg.label_value = 16
+bc_xpos.label_value = 13
+# boundary_yneg_noslab
+bc_yneg.label_value = 17
+# boundary_ypos_noslab
+bc_ypos.label_value = 18
+bc_zneg.label_value = 14
+
+# End of file

examples/subduction-3d/step06_slowslip_gmsh.cfg

@@ -0,0 +1,58 @@
+[pylithapp.metadata]
+# Materials
+#   All four blocks use the linear, isotropic elastic material model
+#   (mat_elastic.cfg).
+#
+# Boundary Conditions
+#   Roller boundary conditions (pylithapp.cfg).
+#
+# Fault
+#   Uniform oblique slip with a duration of one month on the slab
+#   fault rupture patch.
+#
+# Output
+#   We output the displacements for the slip impulses on the ground
+#   surface ang at the fake cGNSS locations defined by cgnss_stations.txt.
+#
+# BEFORE RUNNING THE SIMULATION, run the Python script that generates
+# the spatial database with the final slip distribution and the
+# temporal database with the slip time function.
+#
+#   ./utils/generate_slowslip.py
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = Slow slip on a patch in the middle of the top of the slab.
+keywords = [prescribed slip, slow slip, generate spatial database]
+arguments = [step06_slowslip_gmsh.cfg, step06_slowslip.cfg, mat_elastic.cfg]
+version = 2.0.0
+pylith_version = [>3.0]
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step06_slowslip_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.problem.interfaces.fault_slabtop]
+# patch
+label_value = 24
+edge_value = 25
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+bc_xneg.label_value = 11
+bc_xpos.label_value = 13
+bc_yneg.label_value = 10
+bc_ypos.label_value = 12
+bc_zneg.label_value = 14
+
+
+
+# End of file

examples/subduction-3d/step07a_leftlateral_gmsh.cfg

@@ -0,0 +1,51 @@
+[pylithapp.metadata]
+# Materials
+#   All four blocks use the linear, isotropic elastic material model
+#   (mat_elastic.cfg).
+#
+# Boundary Conditions
+#   Roller boundary conditions (pylithapp.cfg).
+#
+# Fault
+#   Slip impulses for left-lateral slip on the suduction interface (slab top).
+#
+# Output
+#   We output the displacements for the slip impulses on the ground
+#   surface at at the fake cGNSS locations defined by cgnss_stations.txt.
+base = [pylithapp.cfg, mat_elastic.cfg]
+description = "Generate static Green's functions for left-lateral slip on the subduction interface."
+arguments = [step07a_leftlateral_gmsh.cfg, step07a_leftlateral.cfg, mat_elastic.cfg]
+version = 2.0.0
+pylith_version = [>3.0]
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step07a_leftlateral_gmsh
+
+[pylithapp.mesh_generator]
+reader = pylith.meshio.MeshIOPetsc
+reader.filename = mesh_tet.msh
+
+[pylithapp.greensfns]
+label_value = 24
+
+[pylithapp.problem.interfaces.fault_slabtop]
+# patch
+label_value = 24
+edge_value = 25
+
+[pylithapp.problem.solution_observers.groundsurf]
+label_value = 15
+
+[pylithapp.problem.bc]
+bc_xneg.label_value = 11
+bc_xpos.label_value = 13
+bc_yneg.label_value = 10
+bc_ypos.label_value = 12
+bc_zneg.label_value = 14
+
+# End of file