-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathREADME
executable file
·128 lines (96 loc) · 7.55 KB
/
README
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
This code is used for creating a 3D mesh (Box or semi-sphere models, see schematic diagram and jpg files in the work directory) for Specfem3D.
The Box model creates four blocks by cutting a long cube with free surface, fault and several planes. The free surface and fault could be planar or curved surface, by setting Interface and Topography as False or True in ./create_box_mesh.py. Blocks 2 and 4 are the seismogenic layer, bounded by upper and lower seismogenic planes. Please setup the parameters of seismogenic layer, Upper_cutoff and Lower_cutoff, in the script create_box_mesh.py. If Upper_cutoff is a value larger than or equal 0, then the fault cuts through the free surface (i.e., block 1 is removed). The horizontal boundaries of the box domain is aligned with latitude and longitude.
The semi-sphere model creates three blocks. We firstly cut a cylinder with free surface and fault interface. The semi-sphere is created by cutting a sphere with the free surfce, then subtract block 1 and 2 from the semi-sphere to obtain block 3. Finally we merge all blocks together. Block 1 and 2 have a finite seismogenic layer and the fault can cut through this seismogenic layer at any orientation. Please setup Lower_cutoff of seismogenic layer in the script create_semisphere_mesh.py.
~~~\/\/\/\/\~~~/\/\/\/\/~~~~
/ /|
/ Free surface / |
/ / /|
~~~\/\/\/\/\~~~/\/\/\/\/~~~/ / |
| 1 |/ |
---------------------------/ |
| / | |
| /Fault | /|
| 2 \ 4 | / |
| / |/ |
---------/-----------------/ /
| | /
| 3 |/
---------------------------/
Box model
~~~~~~~\/\/\/\~~~/\~~~~~/\/\/\/~~~~~~~~
* | / | *
* |1 / | *
* | / 2 | *
* |/ | *
* -------- *
* *
* 3 *
* *
*******
Semi-sphere model
Huihui Weng, Now associate professor at Nanjing University, contact me by weng@nju.edu.cn or qfkq7850@mail.ustc.edu.cn. Jul. 16 2024
Postdoc at Geoazur UCA, Oct. 31 2018
Added semi-sphere model and revised the bugs in the box model. Dec. 7 2018
Make the fault can cut through the free surface in the box model. Dec. 17 2018
Thank Jean-Paul Ampuero for his helpful suggestion.
Software and package requirements:
CUBIT (or Trelis)
GEOCUBIT (a package provided in the Specfem3D code in the directory: specfem3d-master/CUBIT_GEOCUBIT/)
Python and Python packages: numpy and scipy
It is better to set up the path of CUBIT and GEOCUBIT from ~/.bashrc or ~/.bash_profile. In this case, you can import these two Python packages easily.
Commands for testing the necessary python packages:
> python
Python 2.7.13 (default, Jan 03 2017, 17:41:54) [GCC] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
>>> import scipy
>>> import cubit
>>> import geocubitlib
If not, you need to provide the lib path in the scripts ./create_box_mesh.py and ./python_scripts/playback_jou.py by the line: sys.path.append('PATH to LIBS').
Step-by-step procedure
1. Do you need curved free surface or fault?
If yes, please set up the ranges and reference point of the model in geographical coordinates in Surface/Create_surface.py. Note that the range of the longitude and latitude shall be larger than the actual model range. It is a good choice to set the epicenter as the reference point, i.e., the original point in Cartesian coordinates.
For creating a curved free surface, you need to setup 4 parameters in Surface/Create_surface.py, in addition to the geographical coordinates:
Sur_GRD_data Use a grid-format data or a xyz-format data
Sur_input_data Path of the data
grid_size The sample interval
sigma It is the parameter (sigma) of gaussian_filter in Python. Larger values mean stronger smoothing. The guidelines for gaussian_filter in Python: https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
For creating a curved fault, you need to setup 5 parameters in Interface/Create_interface.py:
Int_GRD_data Use a grid-format data or a xyz-format data
Int_input_data Path of the data
Int_inc The sample interval
sigma It is the parameter (sigma) of gaussian_filter in Python. Larger values mean stronger smoothing. The guidelines for gaussian_filter in Python: https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html
The global topography data can be downloaded from NOAA https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/ . Put the downloaded data "ETOPO1_Ice_g_gmt4.grd.gz" in the directory ./Surface/data and unzip it by the command "gunzip -k ./Surface/data/ETOPO1_Ice_g_gmt4.grd.gz" (or using other ways to unzip a ~.gz file). The downloaded data file is large (377MB), which contains the global topography data. If you need a higher or lower resolution topograhpy data, please put it in ./Surface/data and set up its path to the new data in Surface/Create_surface.py.
The fault interface data can be downloaded from USGS https://earthquake.usgs.gov/data/slab/models.php. Put the downloaded data in the directory ./Interface/data
Run the scripts Surface/Create_surface.py and Interface/Create_interface.py to create the .sat files (CUBIT-format file) by typing:
cd ./Surface
python Create_surface.py
cd ..
cd ./Interface
python Create_interface.py
cd ..
You can check the quality of mesh by the created *.jou script (CUBIT script) in ./Surface/CUBIT_jou/ and ./Interface/CUBIT_jou/ under the GUI of CUBIT (for debug). It is recommanded to change different sample interval and smooth factor to create a lot surface for later usage. The created free surface and fault are saved as *.sat in ../output/
2. If you need to use curved free surface and fault, set Interface = True or Topography = True, then provide the created file by
Int_name = work_dir + "/output/interface_sigma_3_inc_6.sat"
Top_name = work_dir + "/output/surface_sigma_0_inc_8.sat"
For planar free surface, set Topography = False
For planar fault, set Interface = False, and you need to set the strike, dip and depth of the fault. Firstly, the python script creates a Zplane, then rotates about Y axis with angle Dip, rotates about Z axis with angle Strike, and finally moves along Z direction with distance Dep
For example:
Strike = 90
Dip = 60
Dep = -10
Set up other parameters in the first section of create_box_mesh.py and create_semisphere_mesh.py (there are comments), such as dimension of model, materials properties, mesh scheme, mesh size, and output name.
For mesh scheme, you need to set up the grid size, mesh scheme (thex or map), element type, or mesh refinement. Note that we cannot use "map" mesh scheme for semi-sphere model.
grid_size = 1
mesh_scheme = "thex"
#mesh_scheme = "map"
element_type = "HEX8"
#element_type = "HEX27"
fault_refine_numsplit = 0
fault_refine_depth = 5
3. Run the script for the box model by type
python create_box_mesh.py
and for the semi-sphere model by type
python create_semisphere_mesh.py
The created files are seved in the directory ./output/model_name
If you have any question or comment, please feel free to contact me.