Geodynamics - Group: Seismology ~ Wiki: Tutorial 3

The following instructions assume that you have installed SPECFEM3D and familiarized yourself with how you will run the package based on your computer configuration, as detailed in the “ SPECFEM3D User manual ” (Chapter 2 provides installation help). Additionally, we will make use of an external, hexahedral mesher CUBIT. Please make sure you have these packages installed on your system.

The example is distributed with the package under the examples/ directory. However, you might need to edit these example scripts slightly to launch them on your system.

Mount St. Helens

This is a step-by-step tutorial how to create a mesh for a region around Mount St. Helens, export it into a SPECFEM3D file format and run the mesh partitioning and database generation.

Change to the example directory

[]$ cd SPECFEM3D/examples/Mount_StHelens

Please also see the README file in this example directory. An example topography file “ptopo.mean.utm” is provided as well, in the following we explain in detail how to construct such a file.

Downloading topography data

You can get SRTM 90m Digital Elevation Data for a region of interest at: http://srtm.csi.cgiar.org

For this example, we choose Mount St.Helens as region of interest. Mount St. Helens is located at: 46.197 N 122.186 W This region is contained in a downloadable file srtm_12_03.zip, which you will have to download from http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp using the region of interest:

Latitude min: 45 N  max: 50 N
Longitude min: 125 W max: 120 W

unzipping the file

[Mount_StHelens]$ unzip srtm_12_03.zip

leads to:

..
srtm_12_03.tif
..

Converting Geotif topography data

To convert the Geotif-file into an longitude/latitude/elevation format, you can use the package FWTools at: http://fwtools.maptools.org/

Install the package and use their gdal2xyz executable to extract the tif file into xyz format:

[Mount_StHelens]$ FWTools-2.0.6/bin_safe/gdal2xyz.py srtm_12_03.tif > srtm_12_03.xyz

the newly created file srtm_12_03.xyz has now the format:

#longitude #latitude #elevation (m)

the file size is ~ 963 MB.

Extracting the region of interest

To further extract and manipulate the topography data, you can use the package GMT at: http://gmt.soest.hawaii.edu/

For our purpose, the region of interest will be:

\ -R-122.3/-122.1/46.1/46.3//

that is a region of ~23 km x 23 km extent.

Using the blockmean executable from the GMT package, we extract and interpolate the topography data for the detailed region, using an interpolated grid spacing of 0.006 degrees (~ 700 m):

[Mount_StHelens]$ blockmean srtm_12_03.xyz -R-122.3/-122.1/46.1/46.3 -I0.006/0.006 > ptopo.mean.xyz

This will create a new file ptopo.mean.xyz with a longitude/latitude/elevation format of the region of interest

Converting to Cartesian coordinates

Since the CUBIT mesh will need Cartesian coordinates, we convert the topography file from longitude/latitude/elevation to X/Y/Z coordinates using a UTM projection. Mount St.Helens lies in the UTM zone: 10 (T).

Use the script “convert_lonlat2utm.pl” provided in this example folder:

[Mount_StHelens]$ ./convert_lonlat2utm.pl ptopo.mean.xyz 10 > ptopo.mean.utm

to create a new file ptopo.mean.utm which will have a file format:

#UTM_X (m) #UTM_Y (m) #Z (m)

Meshing

In the following, we will run a python script within CUBIT to create the needed mesh files for the partitioner.

if not done so yet, change your working directory to the example folder:
```
[]$ cd SPECFEM3D/examples/Mount_StHelens
```
start the graphical user interface (GUI) of CUBIT:
```
[Mount_StHelens]$ claro
```
run the python file “ read_topo.py ”:
```
use the (Play Journal File) button and open the file 'read_topo.py'
```
The script creates a topography surface topo.stl in STL file format as well as topo.cub based on file ptopo.mean.utm.

run the python file “mesh_mount.py”:

use the (Play Journal File) button and open the file 'mesh_mount.py'

In case the script executed successfully, it should create you a single volume with a refined mesh for the top layer, a triplication layer and a coarser mesh at the bottom. [600px] check that the script created a new folder “ MESH/ ” containing the files:

[Mount_StHelens]$ ls -1 MESH/
absorbing_surface_file_bottom
absorbing_surface_file_xmax
absorbing_surface_file_xmin
absorbing_surface_file_ymax
absorbing_surface_file_ymin
free_surface_file
materials_file
mesh_file
nodes_coords_file
nummaterial_velocity_file

you can also check if the export was successful by examining the output in the Command line window of CUBIT.

Setting up example folder for simulations

We will set up the example folder for simulation runs: All the steps and following decomposition, database generation and solver run are put in a process.sh bash script file in the example folder. 1. In case you can run parallel programs on your desktop (needs an MPI installation), you can simply run the script:

[Mount_StHelens]$ ./process.sh

to do the setup and following steps for you. The simulation will run on 4 CPU cores by default. Please modify and adapt the script to your needs.

Note, this example should take about 15 minutes to finish the simulation.

2. check the output file output_solver.txt in the output directory in_out_files/OUTPUT_FILES/ to see if the forward simulation was successfully finishing. the seismograms will look like this, using gnuplot:

gnuplot> plot "X10.DB.HXZ.semd" w l,"X20.DB.HXZ.semd" w l,"X30.DB.HXZ.semd" w l,"X40.DB.HXZ.semd" w l

Created on 02 Mar 2022, Last modified on 02 Mar 2022