[cig-commits] r18653 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . src src/devel
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Mon Jun 20 00:40:10 PDT 2011
Author: percygalvez
Date: 2011-06-20 00:40:10 -0700 (Mon, 20 Jun 2011)
New Revision: 18653
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
new_database for splay faults working with tapper
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2011-06-20 07:40:10 UTC (rev 18653)
@@ -0,0 +1,245 @@
+This package contains a version of the 3D spectral element code SPECFEM3D-SESAME (Komatitsch et al)
+modified by Galvez, Ampuero and Nissen-Meyer to model dynamic earthquake on non-planar faults
+with slip-weakening friction. The main modifications are two new modules (fault_object.f90 and
+fault_solver.f90). We also include some examples and Matlab function for post-processing and
+visualization.
+
+This is a preliminary version, still under development and testing. The format of inputs, outputs
+and features are subject to change.
+
+Details about the original package are described in its manual (manual_SPECFEM3D.pdf).
+The current README file describes how to install and run our modified version the code.
+
+This document assumes that the root directory containing the code is ~/SPECFEM3D.
+
+
+
+
+I. INSTALLATION
+----------------
+
+1. Copy the start-up and source files to the root directory:
+
+ >> cd ~/SPECFEM3D
+ >> cp START_UP/* src/* .
+
+
+2. Modify compiler options if needed:
+ The file flags.guess contains default values for compiler flags, which in general work fine.
+ To force your preferred compiler settings you must edit the variables FLAGS_CHECK and
+ FLAGS_NO_CHECK in file flags.guess.
+
+
+3. Configure the package, for instance:
+
+ >> ./configure FC=ifort MPIFC=mpif90
+
+ Several settings can be tuned for your system, see Chapter 2 of the manual for more details.
+ The code runs by default in single precision. If you prefer to run in double precision do instead:
+
+ >> ./configure --enable-double-precision FC=ifort MPIFC=mpif90
+
+
+4. Compile the package:
+
+ >> make
+
+ This creates the following executables:
+ xgenerate_dabases (database generation)
+ xspecfem3d (solver)
+ xcombine_AVS_DX
+ xconvolve_source_timefunction.
+
+
+5. Install SCOTCH:
+
+ >> cd decompose_mesh_SCOTCH
+ >> tar -xvf scotch.5.1.7.tar
+ Follow the instructions in the scotch_5.1.7/INSTALL.txt file.
+
+
+6. Compile the SCOTCH-to-SPECFEM3D interface program, also in the directory decompose_mesh_SCOTCH:
+
+ Edit the Makefile:
+ + the F90 variable to match your compiler.
+ + the SCOTCH_LIB_PATH to match your SCOTCH library path.
+ >> make
+
+ This creates the executable xdecompose_mesh_SCOTCH in directory decompose_mesh_SCOTCH.
+
+
+7. Obtain and install the mesh generation software CUBIT (cubit.sandia.gov).
+ If you don't have a CUBIT license yet, you can still run the examples in this package using
+ the mesh files we provide.
+
+
+
+
+II. RUNNING A SIMULATION
+-------------------------
+
+1. Create a mesh with CUBIT:
+ (If you have not installed CUBIT yet skip this step and use the mesh files included
+ with the examples.)
+
+ Move to the directory CUBIT, open the cubit GUI, run a CUBIT script:
+ >> cd ~/SPECFEM3D/CUBIT
+ >> cubit
+ In CUBIT's menu "Tools", select "Play Journal File" and select a script file,
+ for instance EXAMPLES/tpv5/fault_tpv5.py
+
+ This creates several mesh files in directory CUBIT/MESH/:
+ absorbing_surface_file_* (5 files)
+ free_surface_file
+ materials_file
+ mesh_file
+ nodes_coords_file
+ nummaterial_velocity_file
+
+ The CUBIT graphics window should show the mesh (e.g. EXAMPLES/tpv5/planar_fault.png)
+
+
+2. Partition the mesh with the domain decomposition software SCOTCH.
+
+ >> cd ~/SPECFEM3D/decompose_mesh_SCOTH
+ >> ./xdecompose_mesh_SCOTCH 'nproc' ../CUBIT/MESH/ ../DATABASES_MPI/
+
+ where 'nproc' is the number of processors that you will use in the simulation.
+ The second argument (../CUBIT/MESH/ in this example) is the path to the directory
+ containing the mesh files that were generated by CUBIT.
+ This creates one mesh file proc000***_Database per processor in directory DATABASES_MPI.
+ These files can be large, so you might want to place DATABASES_MPI in a scratch disk
+ and specify the path accordingly when executing xdecompose_mesh_SCOTCH.
+
+
+3. Edit the file DATA/Par_file:
+
+ LOCAL_PATH = should be the path to DATABASES_MPI
+ NPROC = number of processors. The same as the number of partitions in SCOTCH (step 2).
+
+
+4. Generate databases:
+
+ >> cd ~/SPECFEM3D
+ >> mpirun -np nproc ./xgenerate_databases
+ or
+ >> ./run/run.xdatabases
+
+ This creates several binary mesh files for each processor (proc000***.bin)
+ in directory DATABASES_MPI.
+
+
+5. Run the solver:
+
+ >> mpirun -np nproc ./xspecfem3D
+ or
+ >>./run/run.xspecfem3d
+
+
+
+
+III. EXAMPLES
+--------------
+
+The package includes a few examples, the SCEC benchmark problems:
+ + TPV5, a planar vertical strike-slip fault
+ + TPV14 and TPV15, vertical strike-slip fault system with a fault branch
+
+To run the examples first replace the contents of directory DATA/* by one of the EXAMPLES/tpv*/DATA,
+for instance:
+
+ >> cp -r EXAMPLES/tpv5/DATA ./
+
+Then follow all the steps in section II above.
+
+
+
+
+IV. INPUT FILES
+----------------
+
+DATA/Par_file See SPECFEM3D manual page 17.
+ A first version of this file is generated by ./configure.
+
+DATA/STATIONS List of stations outside the fault (see manual page 23).
+
+DATA/FAULT/Par_file_faults.in contains parameters of the fault:
+ Line 1: Number of faults (NF)
+ Lines 2 to NF+1: Two columns: domain tag on side #1 of the fault
+ and domain tag on side #2 of the fault. These domain tags are assigned
+ to blocks of elements on each side of the fault during mesh generation in CUBIT.
+ #3 eta (Kelvin Voigt factor).
+ Line NF+2: Number of time steps between updates of the time series outputs at selected
+ fault points (see DATA/FAULT/FAULT_STATIONS.in), usually a large number (100s or 1000s).
+ Note that the sampling rate of the time series is usually much higher.
+ Line NF+3: Number of time steps between fault snapshot outputs (quantities at every fault
+ point exported at selected times), usually a large number (100s or 1000s).
+ Line NF+4: Slip velocity threshold below which frictional healing is set (friction coefficient
+ is reset to its static value). If this value is negative healing is disabled.
+ Line NF+5: Slip velocity threshold to define the rupture front. Only used for outputs.
+ Line NF+6: Fault parameters id.
+ Line NF+7: Initial fault stress for fault #1:
+ S1 = along-strike shear
+ S2 = along-dip shear
+ S3 = normal stress (negative in compresion)
+ See figure 1 below.
+ Line NF+8: Static and dynamic friction coefficients.
+ Repeat NF+6 to NF+8 for each fault.
+
+DATA/FAULT/FAULT_STATIONS.in Stations in the fault plane.
+ Line 1: number of stations.
+ Line 2 to end: 5 columns: X, Y, Z (-depth), station name, fault-id
+ The fault-id identifies the fault that contains the station.
+ It is the index of appearance in the faults list after line 2 of Par_file_faults.in
+
+Several files are generated automatically by xgenerate_databases in directory DATABASES_MPI
+and do not need to be modified by the user.
+
+
+
+
+IV. OUTPUT FILES
+-----------------
+
+Several output files are saved in ~/SPECFEM3D/OUTPUT_FILES:
+
+1. Stations on the fault plane. Their locations and names are given in
+ DATA/FAULT/FAULT_STATIONS.in. Their output files are named after the station.
+
+# Column #1 = Time (s)
+# Column #2 = horizontal right-lateral slip (m)
+# Column #3 = horizontal right-lateral slip rate (m/s)
+# Column #4 = horizontal right-lateral shear stress (MPa)
+# Column #5 = vertical up-dip slip (m)
+# Column #6 = vertical up-dip slip rate (m/s)
+# Column #7 = vertical up-dip shear stress (MPa)
+# Column #8 = normal stress (MPa)
+
+Fault coordinates system : (strike,dip,normal)
+
+ .(normal)
+ .
+ .
+ -------------> (strike)
+ | |
+ | Fault |
+ |---------
+ .
+ . (dip)
+ Figure 1.
+
+2. Stations in the bulk, outside the fault plane.
+ Their locations and names are given in DATA/STATIONS.
+ Output format: see manual page 51.
+
+3. Rupture time files are named Rupture_time_FAULT-id. Their format is 4 columns:
+ X, Y, Z and rupture time.
+
+4. The face FAULT plane reference used in the computatinos is the first face defined by the user.
+
+
+V. POST-PROCESSING AND VISUALIZATION
+-------------------------------------
+
+Some Matlab functions for post-processing and visualization are included in directory
+Post-processing. The functions are internally documented (see their matlab help).
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-06-19 18:25:06 UTC (rev 18652)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-06-20 07:40:10 UTC (rev 18653)
@@ -376,11 +376,11 @@
integer, intent(in) :: iin,n
real(kind=CUSTOM_REAL) :: b(size(a))
- character(len=10) :: shape
- real(kind=CUSTOM_REAL) :: val, xc, yc, zc, r, l, lx,ly,lz
+ character(len=20) :: shape
+ real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
integer :: i
- NAMELIST / DIST2D / shape, val, xc, yc, zc, r, l, lx,ly,lz
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
if (n==0) return
@@ -394,6 +394,8 @@
lx = 0e0_CUSTOM_REAL
ly = 0e0_CUSTOM_REAL
lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
read(iin,DIST2D)
select case(shape)
case ('circle')
@@ -403,18 +405,23 @@
case ('square')
b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL)
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
case ('rectangle')
b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL)
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
case default
stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
end select
-! a =a + b*val
-!Percy , assigning straight values of each patch .
- where (b /= 0) a = b*val
+ where (b /= 0) a = b
enddo
end subroutine init_2d_distribution
@@ -517,6 +524,8 @@
! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+! TO DO :
+! To avoid division by cero (temporal)
tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
t2 = T(2,:)/tnorm
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2011-06-19 18:25:06 UTC (rev 18652)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2011-06-20 07:40:10 UTC (rev 18653)
@@ -376,11 +376,11 @@
integer, intent(in) :: iin,n
real(kind=CUSTOM_REAL) :: b(size(a))
- character(len=10) :: shape
- real(kind=CUSTOM_REAL) :: val, xc, yc, zc, r, l, lx,ly,lz
+ character(len=20) :: shape
+ real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
integer :: i
- NAMELIST / DIST2D / shape, val, xc, yc, zc, r, l, lx,ly,lz
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
if (n==0) return
@@ -394,6 +394,8 @@
lx = 0e0_CUSTOM_REAL
ly = 0e0_CUSTOM_REAL
lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
read(iin,DIST2D)
select case(shape)
case ('circle')
@@ -403,18 +405,23 @@
case ('square')
b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL)
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
case ('rectangle')
b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL)
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
case default
stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
end select
-! a =a + b*val
-!Percy , assigning straight values of each patch .
- where (b /= 0) a = b*val
+ where (b /= 0) a = b
enddo
end subroutine init_2d_distribution
@@ -517,6 +524,8 @@
! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+! TO DO :
+! To avoid division by cero (temporal)
tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
t2 = T(2,:)/tnorm
More information about the CIG-COMMITS
mailing list