Home

User`s guide - California Institute of Technology

1. 00 00000 ee eee 2 4 Fault interface conditions 0 000 eee ee 24 1 Linear slip laws aaa Bi a Ona ie eA ee n ee 2 4 2 Normal stress response 1 a 2 4 3 O 4 8 4 orar O tA es De eI Ate oe et 3 Mesh generation 3 1 General guidelines i sars ea eR Oe Re wee a as 3 2 Meshing features included in the solver 0 0 00002 eae 3 3 Meshing with the MESH2D Matlab utilities 3 4 Generating a mesh with EMC2 e e 3 4 1 The mesh generator EMC2 020200000 AZ Notations a ae Rodd at Gene A Pe ER 3 43 Basic step by step sise a eaa o e ioa eaa aa A a a SAAN Further tipsesans inga aa ta ee in pe ye A i 3 5 Importing a mesh from CUBIT EZ4U or other mesh generation software OO WOM ONNA H oro 11 12 12 13 13 14 15 15 16 16 16 17 19 19 20 20 21 21 22 22 25 25 CONTENTS 4 4 The solver SEM2D 27 4 1 Abovt the method roscas he ae ES a a od nee A 27 42 Basic usape ow sac co be aoe A By eS 28 4 3 General format of the input fle 2 000 28 4 4 Input Blocks Reference Guide oa aaa a 02000002 eee 32 4 5 Verifying the settings and running a simulation 54 4 6 Outputs their visualization and manipulation 59 4 6 1 Spectral element grid e 59 4 6 2 Source time function e 59 4 6 3 OMApshotS ss ias dbo wig a aS ER ek a a Ae 59 464 Seis MOSTAMS erro a de A
2. wave propagation in anisotropic TTI media Dewangan et al 2006 fault reflections from fluid infiltrated faults Haney et al 2007 benchmark for anisotropic wave propagation De la Puente et al 2007 dynamic earthquake rupture with rate and state friction Kaneko et al 2008 benchmark for dynamic earthquake rupture simulation De la Puente et al 2009 non linear wave propagation in damaged rocks Lyakhovsky et al 2009 modeling marine seismic profiles Roberts et al 2009 surface wave propagation in applied geophysics Vignoli and Cassiani 2010 Vignoli et al 2011 Vignoli et al 2012 Boaga et al 2012 1 4 Download and updates 7 wave propagation around a prototype nuclear waste storage tunnel Smith and Snieder 2010 earthquake dynamic rupture with off fault plasticity Harris and et al 2011 Gabriel et al 2012a Xu et al 2012a Xu et al 2012b earthquake rupture in heterogeneous media Huang and Ampuero 2011 Huang and Ampuero 2012 benchmark for wave propagation in heterogeneus media O Brien and Bean 2011 dynamic rupture model of the 2011 Tohoku earthquake Huang et al 2012 dynamic rupture model of the 2012 off Sumatra earthquake Meng and Ampuero 2012 earthquake rupture with velocity and state dependent friction Gabriel et al 2012b 1 4 Download and updates SEM2DPACK is hosted by SourceForge at http sourceforge net pro
3. field D Axis of the seismogram plot irepr D Snapshot Outputs Timestep of first snapshot output itl 0 Number of timesteps between snapshots itd 100000 Save results in PS file or not ps F Save grid triangulation for GMT gmt F Save results in AVS file or not avs F Save results in Visual3 file or not visual3 F Save results in binary file or not bin F Selected fields Displacement F Velocity s T Acceleration F Strain F Stress E Selected components for PostScript snapshots x ee gee Ve SS F Y T Ye ae O F Amplitude F FI I II I I IO II II IE ok x 11 A Set ab sae en a ed on phase FI I OO AR RA I IO I oe uorgepus e SUIUUNA pue s3ur33os 94 SUIAJLIOA GP Thursday March 06 2008 info 1 2 LS Printed by Jean Paul Ampuero solver Time step secs 17 621E 03 Number of time steps 1987 Total duration secs 35 013E 00 Courant number 300 000E 03 STABILITY CFL number 300 000E 03 Defining work arrays for elasticity OK Initializing kinematic fields OK Max displ 0 000E 00 Max veloc 0 000E 00 Building the mass matrix OK Defining boundary conditions OK Initializing receivers eceivers Receivers have been relocated to the nearest GLL node Receiver x requested z requested x obtained z obtained distance 1 0 000E 00 0 000E 00 0 000E 00 0 000E 00 0 000E 00 2 5 000E 00 0 00
4. Vilotte 1998 Elastic wave propagation in an irregularly layered medium Soil Dyn Earthquake Eng 18 1 11 18 Vignoli G and G Cassiani 2010 Identification of lateral discontinuities via multi offset phase analysis of surface wave data Geoph Prosp 58 3 389 413 doi 10 1111 5 1365 2478 2009 00838 x Vignoli G G Cassiani M Rossi R Deiana J Boaga and P Fabbri 2012 Geophysical characterization of a small pre Alpine catchment J App Geophys 80 32 42 Vignoli G C Strobbia G Cassiani and P Vermeer 2011 Statistical multioffset phase analysis for surface wave processing in laterally varying media Geophys 76 2 U1 U11 doi 10 1190 1 3542076 Xu S Y Ben Zion and J P Ampuero 2012a Properties of inelastic yielding zones gen erated by in plane dynamic ruptures I Model description and basic results Geophys J Int submitted Xu S Y Ben Zion and J P Ampuero 2012b Properties of inelastic yielding zones generated by in plane dynamic ruptures II Detailed parameter space study Geophys J Int submitted
5. see stf_gen 90 e spatial distributions see distribution _general f90 The source files listed above contain step by step instructions just follow the comments start ing by Chapter 6 Frequently Asked Questions 6 1 SEM2D Segmentation fault This problem is often related to a small stack size in your computer settings In your Linux shell do ulimit s unlimited under bash or limit stacksize unlimited under csh Place this command in your startup files login bashrc or cshrc Instabilities on very distorted elements Very distorted elements with very small or very large angles are usual close to wedges of sedimentary basins fault branching points etc In general distorted elements are less stable than square elements spurious motions with exponentially increasing amplitude might appear in their vicinity In most cases these instabilities can be suppressed by reducing the Courant input parameter There is currently no simple recipe to determine the maximum value of this parameter so trial and error is required 6 2 EMC2 I can t get rid of a few triangles Obtaining a quality quad mesh is not always a trivial task Trial and error and experience is needed This can be by far the most time consuming stage of modeling First make sure that the total number of element edges along the perimeter of each mesh domain is even This is a necessary topological condition to generate a quad only mesh 6 2 EMC2 68 Whe
6. weak form as described in any textbook on the finite element method Two types of 2D problems are solved Plane strain Also known in seismology as P SV and in fracture mechanics as inplane mode or mode IT It is assumed that ug 0 and 0 0x3 0 Hence 13 23 33 0 and there are two degrees of freedom per node uz and uz Antiplane shear Also known in seismology as SH and in fracture mechanics as antiplane mode or mode III It is assumed that uy ug 0 and 0 0x3 0 Hence only 13 and 23 are non zero and there is one degree of freedom per node uy 2 2 Material rheologies We describe here the constitutive equations implemented in SEM2DPACK relating stress o strain eij and internal variables 2 2 1 Linear elasticity Linear isotropic elasticity Stress and strain are linearly related by Hooke s law oj cijxieij Where C 1 is the tensor of elastic moduli In particular for isotropic elasticity Oij AEkk ij 2UEij 2 3 where and y are Lam s first and second parameters respectively In 2D plane strain the only relevant stress components are 011 022 and 012 The intermediate stress 933 although not null does not enter in the equations of motion The S and P wave speeds are cg y p p and cp y 2u p respectively In 2D antiplane shear only the stress components 013 and 033 are relevant and only S waves are generated Linear anisotropic elasticity Transverse
7. Ae ou x 61 4 6 5 Faulti0UGpUtS Va a a A a 61 A 6 67 Stress lutea a5 a a ar eee te Tee ae e 63 AOA Energies e aa da nee a pee Ala le we a ed a ee ee Ate 63 4 6 8 Matlab utilities 22 4426 ce 9 da ea a 63 5 Adding features to SEM2D notes for advanced users 65 5 1 Overview of the code architecture ee 65 5 2 Accessible areas of the code 2 2 ee 65 6 Frequently Asked Questions 67 Ord AS BIND es hers fog a Side nape hat Se ab tay toe go Te head Senile sated Ua abet ache 67 0 24 EMOL raae eet ca Aa A ae A a ei ed eee a 67 Chapter 1 Introduction 1 1 Overview The SEM2DPACK package is a set of software tools for the simulation and analysis of 2D wave propagation and dynamic fracture with emphasis on computational seismology and earthquake dynamics The core of the package is SEM2D a solver for the 2D elastic wave equations and dynamic earthquake rupture based on the Spectral Element Method SEM with explicit time stepping Chapter 2 of this User s Guide summarizes the range of problems that can be solved with SEM2DPACK Section 4 1 provides some background on the SEM The essential properties of the method are its high order accuracy affordable at competitive computational cost its geometrical flexibility to treat realistic complicated crustal structures and its natural treatment of mixed boundary conditions such as fault friction SEM2DPACK provides tools for each step of the general flow of a simulation proj
8. H gt SPP TOUS MSEt GING S Ss SAS SR amp SNAP_DEF itd 100000 ps F bin F Thursday March 06 2008 1 1 Figure 4 1 Input file Par inp for an elementary example in EXAMPLES TestSH a boxed region with a structured mesh 4 3 General format of the input file 31 Printed by Jean Paul Ampuero Mar 06 08 10 49 Par inp Page 1 1 Some general parameters amp GENERAL lexec 0 Ngll 6 fmax 1 5 ndof 1 Title Palos Grandes NS meshed with EMC2 Verbose 1111 ItInfo 1000 Build the mesh amp MESH_DEF Method EMC2 amp MESH_EMC2 File NSO3qb ftq Elastic material parameters amp MATERIAL tag 1 kind ELAST amp MAT_ELASTIC rho 1800 d0 cp 850 d0 cs 450 d0 amp MATERIAL tag 2 kind ELAST amp MAT_ELASTIC rho 2100 d0 cp 1800 d0 cs 650 d0 amp MATERIAL tag 3 kind ELAST amp MAT_ELASTIC rho 2400 d0 cp 2300 d0 cs 850 d0 amp MATERIAL tag 4 kind ELAST amp MAT_ELASTIC rho 2600 d0 cp 3800 d0 cs 2200 d0 amp MAT_ELASTIC rho 2500 d0 cp 5000 d0 cs 2900 d0 pos Boundary conditions Sa 2shsass SSeS amp BC_DEF Tag 2 Kind ABSORB amp BC_ABSORB Stacey F amp BC_DEF Tag 3 Kind ABSORB amp BC_ABSORB Stacey F let_wave T amp BC_DEF Tag 4 Kind ABSORB amp BC_ABSORB Stacey F Ho Tame scheme
9. Settings 352 2222532 ES amp TIME TotalTime 25 d0 Courant 0 55d0 amp TIME_NEWMARK alpha 1 d0 beta 0 d0 gamma 0 5d0 to SOURCES O E amp SRC_DEF stf RICKER Mechanism WAVE Coord 1160000 d0 2000 d0 STF_RICKER f0 1 d0 Onset 1 5d0 Ampli 1 d0 amp SRC_FORCE Angle 90 amp SRC_WAVE Angle 30 phase S RECEIVE LS RS SAS SS SA RRE SS receivers located at the surface by giving a very large vertical position locating them at the nearest computational node AtNode true is the default amp REC_LINE Number 31 First 1163068 0d0 1 d3 Last 1159697 36d0 1 d3 Isamp 10 fo Sa SS aaa ee ee ee amp SNAP_DEF itd 100000 fields V components x itd 3500 amp SNAP_PS Mesh T Vectors F Color T Interpol T DisplayPts 7 ScaleField 0 2d0 Thursday March 06 2008 1 1 Figure 4 2 Input file Par inp for a more realistic example a sedimentary basin with an unstructured mesh generated by EMC2 Available in EXAMPLES UsingEMC2 4 4 Input Blocks Reference Guide 32 4 4 Input Blocks Reference Guide NAME BC_ABSORB GROUP BOUNDARY_CONDITION PURPOSE Absorbing boundary SYNTAX amp BC_ABSORB stacey let_wave stacey log F Apply Stacey absorbing conditions for P SV Higher order than Clayton Engquist the default let_wave log T Allow incident waves across this boundary if mechanism WAVE in amp SRC_DEF NOTE Only i
10. ampli fO dble 0d0 Amplitude dble 0d0 Frequency STF_RICKER SOURCE TIME FUNCTIONS PURPOSE SYNTAX The Ricker wavelet is the second derivative of a gaussian amp STF_RICKER ampli fO onset real 1 Signed amplitude of the central peak real gt 0 0 Fundamental frequency Hz distribution it has a peak at fO and an exponential decay at high frequency The cut off high frequency is usually taken as fmax 2 5 x f0 real gt 1 f0 0 Delay time secs with respect to the peak value The spectrum has a peak at f0 and decays exponentially at high frequencies Beyond 2 5 f0 there is little energy this is a recommended value for fmax onset gt 1 f0 is needed to avoid a strong jump at t 0 which can cause numerical oscillations Ignore if using incident waves STF_TAB SOURCE TIME FUNCTIONS PURPOSE Source time function spline interpolated from values in a file 4 4 Input Blocks Reference Guide 51 SYNTAX file NOTE NAME GROUP PURPOSE SYNTAX ampli onset pari pari pari pari NAME PURPOSE SYNTAX kind Dt Courant NbSteps amp STF_TAB file string stf tab ASCII file containing the source time function two columns time and value Time can be irregularly sampled and must increase monotonically assumes value t lt min time value min time and value t gt max time value max time STF_USER SOURCE TIME FUNCTIONS A templ
11. code is written in FORTRAN 95 Resources available for programmers include A ToDo file included with SEM2DPACK contains a wish list that ranges from basic functionalities to complex code re engineering Chapter 5 gives some guidelines for programmers A Developers Forum to discuss the implementation of new features is available at http sourceforge net forum forum php forum_id 635737 1 10 License This software is freely available for academic research purposes If you use this software in writing scientific papers include proper attributions to its author Jean Paul Ampuero This program is free software you can redistribute it and or modify it under the terms of the GNU General Public License as published by the Free Software Foundation either version 2 of the License or at your option any later version This program is distributed in the hope that it will be useful but WITHOUT ANY WAR RANTY without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License for more details 1 10 License 10 You should have received a copy of the GNU General Public License along with this program if not write to the Free Software Foundation Inc 59 Temple Place Suite 330 Boston MA 02111 1307 USA Chapter 2 Physical background This chapter summarizes the physical assumptions and notations in SEM2DPACK Footnotes provide reference to the input arguments de
12. conditions 17 Version with a regularization slip scale E Ve o _ _ i a o 2 28 where V is slip rate and V and L are constitutive parameters 2 4 3 Friction Slip weakening friction Slip occurs when the fault shear stress reaches the shear strength 7 yo or T po if the Prakash Clifton law is assumed The friction coefficient y is a function of the cumulated slip D given by one of the following laws Linear slip weakening law pb max pd Hs A 2 29 De Exponential slip weakening law H Hs Hs Ma exp D De 2 30 Fast rate and state dependent friction Friction with fast power law velocity weakening at fast slip speed is a first order proxy for physical weakening processes that operate on natural fault zones at coseismic slip veloci ties A rate and state dependent friction law with fast velocity weakening is implemented in SEM2DPACK similar to that adopted e g by Ampuero and Ben Zion 2008 The friction coefficient depends on slip velocity V and a state variable 0 0 V 2 z ton 2 31 ce ea E ca The state variable has units of slip and is governed by the evolution equation V 0V De 2 32 The friction law is defined by the following constitutive parameters us is the static fric tion coefficient a and b are positive coefficients of a direct effect and an evolution effect respectively V is a characteristic velocity scale and
13. is encouraged to inspect the Matlab s function POST sem2d_snapshot_read m to find more about the data formats and their manipulation Snapshots can also be exported as PostScript files xx_XXX_sem2d ps These can be merged into an animated GIF movie file movie gif by the script POST movie csh and displayed by xanim movie gif or animate movie gif An animated GIF can also be created by the Matlab function POST sem2d_snapshot_movie m 4 6 Outputs their visualization and manipulation 60 TOWA YIM poysout SN SOPue p soed popon mowa Tenoedg 117 WOVdATINAS HER k Max 2 252E 00 Time 7 912E 00 s xX Time step 5000 Figure 4 4 Sample snapshot from EXAMPLES UsingEMC2 an obliquely incident SH plane wave impinging on a sedimentary basin The unstructured mesh of spectral elements is plotted on background 4 6 Outputs their visualization and manipulation 61 4 6 4 Seismograms The seismograms are stored using the SEP format a simple binary block of single precision floats The components of the vector field velocity by default are stored in separate files U _sem2d dat where is x or z in P SV and y in SH The seismograms header is in the file SeisHeader_sem2d hdr Its second line contains the sampling timestep DT the number of samples NSAMP and the number of
14. stations NSTA The stations coordinates XSTA and ZSTA are listed from the third line to the end of file With this notations U _sem2d dat contains a NSAMPXNSTA single precision matrix You can view the seismograms using any tool that is able to read the SEP format which is the case of almost all the softwares able to deal with seismic data sem2dsolve generates scripts for the XSU Seismic Unix visualization tool Xline_sem2d csh displays all seismograms together on screen PSline_sem2d csh plots all seismograms on PostScript files U Poly_sem2d ps Xtrace_sem2d csh prompts the user for a trace number between 1 and NSTA and then displays this particular trace on screen PStrace_sem2d csh does the same as Xtrace but exports the traces as PostScript files U TraceXXX_sem2d ps where XXX is the number of that particular trace The program post_seis exe performs similar basic manipulation and plotting through gnuplot of the seismograms Its interactive menu is self explanatory It is usually called in side a script as in POST seis_b2a csh converts all seismograms to ASCII or POST seis_plot csh plots all seismograms together an example is shown in Figure 4 5 The script POST sample_seis m shows how to manipulate and plot seismogram data in Mat lab It uses the functions POST sem2d_read_seis m and POST plot_seis m 4 6 5 Fault outputs Fault data from dynamic rupture simulations is stored in three files where XX is the
15. the cumulative stress glut tensors defined respectively as si t fone dz dz 4 2 f f loy t Spayefa t de de 4 3 Ww a ane A I where e is the plastic strain e the elastic strain the absolute stress and the tensor of elastic moduli of the undamaged medium To enable this feature set COMPUTE_STRESS_GLUT true in file SRC constants f90 then re compile the code The stress glut output is exported in the file stress_glut_sem2d tab in 7 columns time sf Sba Sha SH Shay S40 4 6 7 Energies The solver can export the cumulative plastic energy the kinetic energy and the total change of elastic energy defined respectively by BP t i ot E 0 de dede 4 4 EX 1 2 I J mie ands 4 5 E t I I Ule t Ulee 0 dx dz 4 6 where U is the elastic potential and summation over subindices is implied To enable this feature set COMPUTE_ENERGIES true in file SRC constants f90 then re compile the code The energy output is exported in the file energy_sem2d tab in 4 columns time EP EF and Ee 4 6 8 Matlab utilities A range of functions and sample scripts for Matlab are available to read manipulate and plot output data Add the directory POST to your Matlab path addpath For an overview of existing utilities type help POST 4 6 Outputs their visualization and manipulation 64 SEM2DPACK POST provides Matlab utilities for the manipulation and visualization of SEM2
16. to better see the boundaries The same procedure applies to define split node interfaces such as faults and cracks you must assign a different tag to each side of the fault Save your work in EMC2 format lt SAVE lt name The resulting file is name emc2_bd EDIT generates the mesh 1 Switch to the edit mesh tool lt EDIT_MESH lt Press ENTER 4 times A triangles mesh appears You must convert it to a quad mesh Convert the triangle mesh to a quad mesh lt QUADRANGULATE gt lt ALL gt You can smooth the mesh with lt REGULARIZE gt lt ALL gt The final mesh is displayed If there remain some triangles come back to the previous step and figure out how to modify the points per edge to help the mesher Some experience is needed here Renumber the mesh in order to optimize computations lt RENUMBER gt Define the boundary condition for the 4 corner nodes of the model these nodes belong to 2 external boundaries so they were given a reference number 0 MODIF_REF n corner click close to corner inside element Where n is the reference number of one of the 2 boundaries containing the corner node Zooming can be useful The same operation must be performed for the corner nodes of the subdomains belonging to an external boundary and for the the crack tip nodes However as a special case crack tip nodes must be assigned the 1 tag 3 5 Importing a mesh from CUBIT EZ4U or other mesh gene
17. 0E 00 5 000E 00 0 000E 00 0 000E 00 3 10 000E 00 0 000E 00 10 000E 00 0 000E 00 0 000E 00 4 15 000E 00 0 000E 00 15 000E 00 0 000E 00 0 000E 00 Mar 06 08 18 04 info Page 3 4 Mar 06 08 18 04 info Page 4 4 5 20 000E 00 0 000E 00 20 000E 00 0 000E 00 0 000E 00 6 25 000E 00 0 000E 00 25 000E 00 0 000E 00 0 000E 00 Defining the FEM mesh OK 7 30 000E 00 0 000E 00 30 000E 00 0 000E 00 0 000E 00 Saving node coordinates in file MeshNodesCoord_sem2d tab Saving element connectivity in file ElmtNodes_sem2d tab Maximum distance between asked and real 0 000E 00 Sipe Gs tex aL elements grid Sampling rate Hz 56 751E 00 Sampling timestep secs 17 621E 03 Total number of samples 1988 Numbering GLL points OK Number of receivers 7 Total number of GLL points 90601 OK Saving element node table in binary file ibool_sem2d dat OK Initializing sources Defining nodes coordinates OK Saving the grid coordinates coord in a text file OK Saving the grid coordinates coord in a binary file OK Sources have been relocated to the nearest GLL node Material properties Source x requested z requested x obtained z obtained distance Translating input model OK T 0 000E 00 0 000E 00 0 000E 00 0 000E 00 0 000E 00 Exporting model OK Maximum distance between requested and real 0 000E 00 OK Timestep 0 t 0 000E 00 vmax 0 000E 00 dmax 0 000E 00 Checking mesh OK FOI I
18. Blocks Reference Guide Al NOTE NAME GROUP SYNTAX eta ETAxDT NOTE NAME GROUP SYNTAX cp cs rho Some combinations of material kinds can be assigned to the same domain Any material type can be combined with KV for instance MATERIAL tag 1 kind ELAST KV followed by a amp MAT_ELAST block and a amp MAT_KV block sets an elastic material with Kelvin Voigt damping MAT_KV MATERIALS PURPOSE Sets material properties for a Kelvin Voigt viscous material Adds a damping term C v K eta v where eta is a viscous time This produces attenuation with frequency dependent quality factor Q 1 etax2x pixf Its main usage is for artificial damping of high frequency numerical artifacts generated by dynamic faults which requires a thin layer of Kelvin Voigt elements surrounding the fault with eta dt 0 1 to 0 3 and a layer thickness of 1 to 2 elements on each side of the fault amp MAT_KV eta ETAxDT amp MAT_KV etaH ETAxDT followed by a DIST_XXX input block dble 0d0 Viscosity coefficient log T If eta is given in units of dt timestep Kelvin Voigt viscosity modifies the stability of time integration The timestep or the Courant number must be set to a value smaller than usual The critical timestep for a Kelvin Voigt material integrated with the leapfrog time scheme is dtc_kv eta sqrt 1 dtc 2 eta 2 1 where dtc is the critical timestep for a purely e
19. DIST_SPLINE file dim file name none Name of the ASCII file containing the interpolation data one line per point two columns one line per point two columns position X or Z value dim int 1 Interpolate along X dim 1 or along Z dim 2 NAME GENERAL PURPOSE General parameters SYNTAX amp GENERAL iexec ngll fmax title verbose itInfo iexec int 0 Run level O just check 1 solve ngll int 9 Number of GLL nodes per edge on each spectral element polynomial order 1 Usually 5 to 9 fmax dble 1 d0 The code checks if this maximum frequency is well resolved by the mesh and issues a warning if not This parameter is not used in computations only for checking To improve the resolution for a given fmax you must increase ngll but you will have to use shorter timesteps or refine the mesh ndof int 2 Number of degrees of freedom per node 1 SH waves anti plane 2 P SV waves in plane title word none Title of the simulation verbose char 4 1101 Print progress information during each phase verbose 1 input phase verbose 2 initialization phase verbose 3 check phase verbose 4 solver phase Example 0001 is verbose only during solver itInfo int 100 Frequency in number of timesteps for printing progress information during the solver phase 4 4 Input Blocks Reference Guide 40 NAME MAT_DAMAGE GROUP MATERIALS PURPOSE Set material properties for the dama
20. DPACK simulations results Reading simulation data SEM2D_READ_SPECGRID SEM2D_SNAPSHOT_READ SEM2D_READ_SEIS SEM2D_READ_FAULT Data manipulation SEM2D_EXTRACT_POINT SEM2D_EXTRACT_LINE ARIAS_INTENSITY RESPONSE_SPECTRUM Data visualization SEM2D_PLOT_GRID SEM2D_SNAPSHOT_PLOT SEM2D_SNAPSHOT_GUI reads a spectral element grid reads snapshot data reads seismogram data reads fault data extracts field values at an arbitrary point extracts field values along a vertical or horizontal line computes Arias Intensity and Significant Duration computes response spectra peak dynamic response of single degree of freedom systems plots a spectral element grid plots snapshot data interactively plots snapshot data SEM2D_SNAPSHOT_MOVIE makes an animation of snapshot data PLOT_MODEL PLOT_SEIS PLOT_FRONTS SAMPLE_FAULT SAMPLE_SEIS Miscellaneous tools XCORRSHIFT SPECSHIFT SPECFILTER plots velocity and density model plots multiple seismograms space time plot of rupture front and process zone tail example of visualization of fault data example of visualization of seismogram data cross correlation time delay measurement Signal time shift by non integer lag via spectral domain zero phase Butterworth filter via spectral domain Chapter 5 Adding features to SEM2D notes for advanced users Sometimes you will need to add new capabilities to the SEM2DPACK solver by modifying the program The following notes are intend
21. De is a characteristic slip scale The steady state 6 0 friction coefficient NV V V 4In amp BC_DYNFLT_NOR this law is set by kind 3 and the two relevant parameters are V and L vstar in amp BC_DYNFLT_RSF Hf Hs a b 2 33 2 4 Fault interface conditions 18 weakens asymptotically as 1 V when V gt V if a lt b approaching its dynamic value ua fis a b over a relaxation timescale D V The value of the relaxation time D V tunes the weakening mechanism between two limit cases slip weakening and velocity weakening If D V is much longer than the typical time scale of fluctuation of the state variable 0 0 Equation 2 32 becomes 6 V implying that 6 is proportional to slip and that the evolution term of the friction coefficient is effectively slip weakening with characteristic slip weakening distance De Conversely if De Ve is short the relaxation to steady state is fast 0 De V V and the friction is effectively velocity weakening with characteristic velocity scale Vo Logarithmic rate and state friction Not implemented yet Dieterich and Ruina classical rate and state laws with aging or slip state evolution law Chapter 3 Mesh generation 3 1 General guidelines The Spectral Element Method SEM requires an initial decomposition of the space domain into quadrilateral elements a quad mesh Obtaining the best performance accuracy cost out of the SEM imposes constraints on the
22. II IOI IO IOI IO IO II IO II Rak Max mesh size 142 616E 03 Solver phase x Min mesh size 58 736E 03 FOI II IOI OOOO Ratio max min 2 428E 00 RESOLUTION nodes per min wavelength 8 000E 00 CPU TIME ESTIMATES in seconds for maximum frequency 1 250E 00 Hz CPU time for initialization 806 877E 03 minimum wavelength 1 600E 00 m CPU time per timestep 18 997E 03 Total solver CPU time 37 747E 00 Dump PostScript Resolution_sem2d ps OK mins 629 118E 03 Dump PostScript Stability_sem2d pS OK hours 34 10 485E 03 Time Timestep 1000 t 17 621E 00 vmax 91 661E 03 dmax 28 653E 03 CPU TIME INFORMATION in seconds CPU time for initialization 806 877E 03 CPU time per timestep 20 395E 03 Total solver CPU time 40 524E 00 mins 675 397E 03 hours 11 257E 03 Storing sismos data SEP format Program SPECFEM end 06 03 2008 Thursday March 06 2008 2 2 uorgepus e SUIUUNA pue s3ur33os 94 SUIAJLIOA Sp 8S 4 6 Outputs their visualization and manipulation 59 4 6 Outputs their visualization and manipulation In addition to the screen output described above sem2dsolve generates different files and scripts that allow the user to control the parameters of the simulation and to display the results All the outputs files follow the naming convention SomeName_sem2d xxx where xxx is one of the following e
23. OSE Define a boundary condition SYNTAX amp BC_DEF tag tags kind possibly followed by amp BC_kind blocks tag Lint none A number assigned to the boundary If you are using SEM2D built in structured mesher the conventions are 1 bottom 2 right 3 up 4 left If you are importing a mesh you must use the tags assigned to the boundaries during the mesh construction tags int 2 none Two tags are needed for split node interfaces faults and for periodic boundaries kind char 6 none Type of boundary condition The following are implemented gt DIRNEU ABSORB PERIOD LISFLT DYNFLT 4 4 Input Blocks Reference Guide 36 NOTE NAME GROUP PURPOSE SYNTAX Ktang Ctang Knorm Cnorm Most of the boundary conditions need additional data given in a BC_kind input block of the BOUNDARY_CONDITIONS group immediately following the BC_DEF block BC_LSF BOUNDARY_CONDITION Linear slip fault a displacement discontinuity interface where stress and slip are linearly related amp BC_LSF Ktang Ctang Knorm Cnorm dble Inf Tangential stiffness dble 0d0 Tangential compliance dble Inf Normal stiffness dblel 0d0 Normal compliance NOTE For each component You can set K _or_ C but _not_both_ If C 0d0 or K Inf then no discontinuity is allowed transparent If K 0d0 the fault is free stress boundary NAME DIST_GAUSSIAN GROUP DISTRIBUTIONS_2D PURPOSE Bell shaped Gaus
24. P MESH_DEF PURPOSE Define mesh parameters for one layer SYNTAX amp MESH_LAYER nz ztop ztopH tag followed by a DIST_XXXX block if ztopH is set nz int none Number of elements in layer along Z direction ztop dble none Only for layers with flat top surface vertical position of top boundary ztopH string none Name of the type of spatial distribution to generate an irregular non flat top boundary In general it is one of the 1D distribution available through a DIST_XXXX block ztopH LINEAR or ztopH SPLINE etc There are two methods to generate a curve with a smooth normal typically to guarantee smooth boundary conditions on curved faults The first method is based on quadratic splines and sometimes produces degenerated elements ztopH QSPLINE followed by a amp QC_SPLINE block The second method is based on cubic splines and is more robust ztopH CSPLINE followed by a amp QC_SPLINE block tag int none Material tag If not given a tag is automatically assigned to the layer sequentially numbered from top to bottom top layer tag 1 NOTE If ztopH LINEAR the mesh uses linearly deformed Q4 elements otherwise it uses quadratically deformed Q9 elements 4 4 Input Blocks Reference Guide 45 NAME QC_SPLINE GROUP MESH_LAYER PURPOSE Define the boundary of a layer using quadratic or cubic splines and enforcing smooth continuous normal between elements for instance to
25. Q4 elements 4 control nodes are supported For a smoother description of boundaries Q9 would be desirable 3 5 Importing a mesh from CUBIT EZ4U or other mesh generation software The most convenient way of doing this is by writing a Matlab function that reads the external mesh file performs mesh manipulations if required and exports a mesh file with 3 5 Importing a mesh from CUBIT EZ4U or other mesh generation software 26 the MESH2D format Section 3 3 Two examples of this are provided in POST mesh2d The functions READ_DCM and READ_INP read a mesh in respectively the DCM format of EZ4U http www lacan upc es ez4u htm and in the INP format of ABAQUS exported by CUBIT Chapter 4 The solver SEM2D 4 1 About the method Given a crustal model meshed with quadrilateral elements and a set of material properties sources receivers and boundary conditions SEM2D solves the elastic wave equation applying the Spectral Element Method SEM for the space discretization and a second order explicit scheme for the time discretization The range of physical problems solved by SEM2D material constitutive equations and boundary conditions is described in more detail in Chapter 2 The SEM introduced by Patera 1984 in Computational Fluid Dynamics can be seen as a domain decomposition version of the Spectral Method or as a high order version of the Finite Element Method It inherits from its parent methods the accuracy spectral convergence
26. RIAL followed by a amp MAT Material block amp BC_DEF one for each boundary condition each followed by a amp BC_Kind block amp TIME amp SRC_DEF followed by amp STF_SourceTimeFunction and amp SRC_Mechanism blocks amp REC_LINE 4 3 General format of the input file 30 Printed by Jean Paul Ampuero Mar 06 08 10 48 Par inp Page 1 1 Parameter file for SEM2DPACK 2 0 Some general parameters amp GENERAL iexec 1 ngll 6 fmax 1 25d0 ndof 1 title Test SH verbose 1111 ItInfo 1000 SS Builds themes hws Sa Se a SS amp MESH_DEF method CARTESIAN amp MESH_CART xlim 0 d0 30 d0 zlim 0 d0 30 d0 nelem 60 60 Elastic material parameters amp MATERIAL tag 1 kind ELAST amp MAT_ELASTIC rho 1 d0 cp 1 7321d0 cs 1 d0 aoa Boundary conditions amp BC_DEF tag 2 kind ABSORB amp BC_ABSORB stacey F Rh BC_DEF tag 3 kind ABSORB BC_ABSORB stacey F QR tes Tame scheme settings 2322222222 ssaesses amp TIME TotalTime 35 d0 courant 0 3d0 to SOU CSS SaaS RO A eS amp SRC_DEF stf RICKER coord 0 d0 0 d0 mechanism FORCE amp STF_RICKER f0 0 5d0 onset 3 d0 ampli 0 25d0 amp SRC_FORCE angle 0d0 Receivers 2329 992 SS amp REC_LINE number 7 field D first 0 d0 0 d0 last 30d0 0 d0 isamp 1
27. SEM2DPACK A Spectral Element Method tool for 2D wave propagation and earthquake source dynamics User s Guide Version 2 3 8 August 2012 Jean Paul Ampuero California Institute of Technology Seismological Laboratory 1200 E California Blvd MC 252 21 Pasadena CA 91125 2100 USA E mail ampuero gps caltech edu Web http web gps caltech edu ampuero Phone 626 395 6958 Fax 626 564 0715 2003 2012 Jean Paul Ampuero Contents 1 Introduction Tale Overviews 008 ed A eo woe ae Poe a A eh oe eS 1 2 History and credits yur Gin Goede he hy Bets Se ote ot ee Hae 1 3 Application examples 02000 eee ee ee 1 4 Download and Updates 0 020000 2 eee ee eee 1 5 Requirements 5 anra fied A a a ee 1 6 Installation 24 gods aoc A GA ee i ae Ee he BE LT Documentation a i Bo eS oe Ge ee PP ee oh ee ds 1 8 SUPPO vii ic Eo te eee Be eaten g ky Bick ip eb gr ang edb aa ot cae 1 94 Contributions e 2 3 2 gh dab A UY Als ae wii eps den Seg D ele LAOS w 25 42 eee a Pa SENN Bi e teed Geko eee es isn aia dE 2 Physical background 2 1 General assumptions and conventions 0 0 000002 ee ee 22 Material theologies ici a Roe ee ee a 2 2 Linear elasticity gt pozos ia a ade oe 2 2 2 Linear visco elasticity 2 2 ee 2 2 3 Coulomb plasticity and visco plasticity 0 4 2 2 4 Continuum damage 0 000 a 2 3 Boundary conditions 0 0 a 2 3 1 Absorbing boundaries
28. Ti pcs Ut 2 21 Tn pcp n 2 22 The formulation above is for P SV mode In SH mode the ABC is Ty peg ty 2In amp MAT_DMG the input argument R is defined as R oC 2 4 Fault interface conditions 16 Stacey ABC The second order accurate ABC introduced by Stacey 1988 under the name P3 is Our QU ee Bo yes fence n 2 23 t Sa cp cg T 2 23 dun Our Be eee E ee ks 2 24 n TA cp cs T 2 24 Its formulation as a mixed boundary condition is T pcs ty pes 2cs cp T 2 25 Tt u Tn pcp n peg 2cs cp TA 2 26 Tt This ABC is only implemented in P SV mode 2 4 Fault interface conditions 2 4 1 Linear slip law Represents a compliant fault zone with elastic contact See Haney et al 2007 2 4 2 Normal stress response Unilateral contact No interpenetration during contact free stress during opening Modified Prakash Clifton regularization Regularization of the normal stress response as required for bimaterial rupture problems is implemented following Rubin and Ampuero 2007 The frictional strength is proportional to a modified normal stress o related to the real fault normal stress by either of the following evolution laws Version with a regularization time scale a 0 2 27 where T is a constitutive parameter 3In amp BC_DYNFLT_NOR this law is set by kind 2 and the relevant parameter is T 2 4 Fault interface
29. Useful to set a damping layer near the fault If ezflt 0 a fault is assumed at the bottom boundary int 1 vertical size number of elements of near fault layer log F Same as ezflt 1 Obsolete will be deprecated following tags are automatically assigned to the boundaries 1 Bottom Right Top Left Fault bottom side Fault top side an FF WD MESH_CART_DOMAIN PURPOSE Define a subdomain within a structured meshed box SYNTAX tag ex ez amp MESH_CART_DOMAIN tag ex ez int none Tag number assigned to this domain int 2 mone Horizontal index of the first and last elements The leftmost element column has horizontal index 1 int 2 mone Vertical index of the first and last elements The bottom element row has vertical index 1 4 4 Input Blocks Reference Guide 43 NOTE If you ignore this input block a single domain tag 1 will span the whole box NAME MESH_EMC2 GROUP MESH_DEF PURPOSE Imports a mesh from INRIA s EMC2 mesh generator in FTQ format SYNTAX amp MESH_EMC2 file file name none Name of the FTQ file including suffix NAME MESH_DEF PURPOSE Selects a method to import generate a mesh SYNTAX amp MESH_DEF method followed by a amp MESH_method input block method name none Meshing method name CARTESIAN LAYERED EMC2 MESH2D NAME MESH_LAYERED GROUP MESH_DEF PURPOSE Structured mesh for layered medium with surface and interfac
30. abriel A A J P Ampuero L A Dalguer and P M Mai 2012b The transition of dynamic rupture styles in elastic media under velocity weakening friction J Geophys Res submitted Hamiel Y Y F Liu V Lyakhovsky Y Ben Zion and D Lockner 2004 A viscoelas tic damage model with applications to stable and unstable fracturing Geophys J Int 159 3 1155 1165 BIBLIOGRAPHY 70 Haney M R Snieder J P Ampuero and R Hofmann 2007 Spectral element mod elling of fault plane reflections arising from fluid pressure distributions Geophys J Int 170 2 933 951 Harris R A and et al 2011 Verifying a computational method for predicting extreme ground motion Seis Res Let 82 5 638 644 doi 10 1785 gssrl 82 5 638 Huang Y and J P Ampuero 2011 Pulse like ruptures induced by low velocity fault zones J Geophys Res 116 B12307 doi 10 1029 2011JB008684 Huang Y and J P Ampuero 2012 The competitive effects on slip pulse properties of low velocity fault zones and strong velocity weakening friction Geophys Res Let submitted Huang Y L Meng and J P Ampuero 2012 A dynamic model of the frequency dependent rupture process of the 2011 Tohoku Oki earthquake EPS in press doi 10 5047 eps 2012 05 011 Kaneko Y N Lapusta and J P Ampuero 2008 Spectral element modeling of earth quake rupture on rate and state faults Effects of velocity strengthening friction at shallow dept
31. ackground only for vector plots at none P P velocity model 25 S velocity model T domains 4 4 Input Blocks Reference Guide 47 isubsamp boundaries symbols numbers legend color ScaleField Interpol DisplayPts NAME REC_LINE lint 3 Subsampling of the GLL nodes for the output of velocity model The default samples every 3 GLL points log T Colors every tagged boundary log T Plots symbols for sources and receivers log F Plots the element numbers log T Writes legends log T Color output dblel 0d0 Fixed amplitude scale saturation convenient for comparing snapshots and making movies The default scales each snapshot by its maximum amplitude log F Interpolate field on a regular subgrid inside each element log 3 Size of interpolation subgrid inside each element is DisplayPts DisplayPts The default plots at vertices mid edges and element center PURPOSE Defines a line of receivers SYNTAX If single receiver line amp REC_LINE number first last AtNode isamp field irepr If receiver locations from file amp REC_LINE file AtNode isamp field irepr number int 0 Number of stations in the line first dble 2 Receivers can be located along a line last file AtNode isamp field irepr NOTE this is the position x z of the first receiver dble 2 Position x z of the last receiver other receivers will be located with regular spacing between First and Las
32. anisotropy with vertical symmetry axis VTI is implemented for 2D P SV Ko matitsch et al 2000 The stress strain constitutive relation for P SV in Voigt notation Ora ci cg 0 Err Ozz c13 c33 0 Ezz 2 4 Oxz 0 0 655 ZExz where the cj are elastic moduli ln the GENERAL input block plane strain is selected by ndof 2 and antiplane shear by ndof 1 2 2 Material rheologies 13 2 2 2 Linear visco elasticity Generalized Maxwell material Not implemented yet Kelvin Voigt material Kelvin Voigt damping can be combined with any of the other constitutive equations by re placing the elastic strain e by e y 0e 0t where y is a viscosity timescale The resulting quality factor Q is frequency dependent Q f 27nf This rheology is not approriate to model crustal attenuation with constant Q unless the source has a narrow frequency band and y is selected to achieve a given Q value at the dominant frequency of the source The main application of Kelvin Voigt viscosity is the artificial damping of high frequency numerical artifacts generated by dynamic faults Dynamic source simulations using meth ods that discretize the bulk such as finite difference finite element and spectral element methods are prone to high frequency numerical noise when the size of the process zone is not well resolved Efficient damping is typically achieved by a thin layer of Kelvin Voigt elements surrounding the fault with thickness of 1 or 2 elem
33. ar grid and values at grid points The format of this ASCII file is Line 1 ncol nx nz x0 z0 dx dz ncol int number of data columns nx nz 2 int number of nodes along x and z x0 z0 2 dble bottom left corner dx dz 2 dble spacing along x and z Line 2 to nx nz 1 ncol dble values at grid points listed from left to right x0 to x0 nx dx then from bottom to top z0 to zO nzx dx int 1 Column of the file to be read The same file can contain values for ncol different properties e g rho vp vs but each DIST_HETE1 block will read only one Even if the original model domain has an irregular shape the regular grid where input values are defined must be rectangular and large enough to contain the whole model domain The regular grid possibly contains buffer areas with dummy values These dummy values should be assigned carefully not random nor zero because SEM2D might use them during nearest neighbor interpolation DIST_LINEAR DISTRIBUTIONS_1D Piecewise linear 1D distribution along X or Z amp DIST_LINEAR n dim length 4 4 Input Blocks Reference Guide 38 dim file length NAME GROUP PURPOSE SYNTAX xn zn NAME GROUP PURPOSE SYNTAX num ref followed immediately by the interpolation data one line per point two columns position X or Z value or amp DIST_LINEAR file dim length and the interpolation data is read from a two column file int 0 Number of points
34. ark parameter If beta 0 the scheme is fully explicit the update of displacement depends only on the last value of acceleration otherwise it is a single predictor corrector scheme gamma dble 0 5d0 Second Newmark parameter Second order requires gamma 1 2 NAME TIME_HHTA GROUP TIME SCHEMES PURPOSE Explicit HHT alpha time integration scheme second order SYNTAX amp TIME_HHTA alpha rho alpha dble 0 5d0 Parameter in the HHT alpha method Values in 0 1 Defined here as 1 HHT s original definition of alpha When alpha 1 it reduces to second order explicit Newmark beta 0 gamma 0 5 rho dble 0 5d0 Minimum damping factor for high frequencies Values in 0 5 1 Rho 1 is non dissipative NOTE We consider only second order schemes for which alphatgamma 3 2 If alpha lt 1 Newmark s beta is related to the HHT parameters by beta 1 alpha rho 2 rho 1 1 alpha 1 rho 3 If alpha 1 we set rho 1 beta 0 gamma 0 5 NOTE Dissipative schemes rho lt 1 require slightly smaller Courant number 0 56 for rho 0 5 compared to 0 6 for rho 1 NOTE This is an explicit version of the HHT alpha scheme of H M Hilber T J R Hughes and R L Taylor 1977 Improved numerical dissipation for time integration algorithms in structural dynamics Earthquake Engineering and Structural Dynamics 5 283 292 implemented with a slightly different definition of alpha 1 original Its properties can be derived from the EG a
35. ate for user supplied source time function File stf_user f90 must be modified by the user to fit special needs amp STF_USER ampli onset pari par2 iparl ipar2 dble 1 Amplitude dble 0 Delay time secs dble 0 Example parameter dble 0 Example parameter int 0 Example parameter int 0 Example parameter TIME Defines time integration scheme amp TIME kind Dt or Courant NbSteps or TotalTime Possibly followed by a TIME_XXXX block char 10 leapfrog Type of scheme gt newmark Explicit Newmark gt HHT alpha Explicit HHT alpha leapfrog Central difference gt symp_PV Position Verlet gt symp_PFR Position Forest Ruth 4th order gt symp_PEFRL Extended PFR 4th order dble none Timestep in seconds dble 0 5d0 the maximum value of the Courant Friedrichs Lewy stability number CFL defined as CFL Dt wave_velocity dx where dx is the distance between GLL nodes Tipically CFL lt 0 5 int none Total number of timesteps TotalTime int none Total duration in seconds 4 4 Input Blocks Reference Guide 52 NOTE The leap frog scheme is recommended for dynamic faults It is equivalent to the default Newmark scheme beta 0 gamma 1 2 However it is faster and requires less memory NAME TIME_NEWMARK GROUP TIME SCHEMES PURPOSE Explicit Newmark time integration scheme SYNTAX amp TIME_NEWMARK gamma beta beta dble 0d0 First Newm
36. boundary tag of the first side of the fault tags 1 of the BC_SWFFLT input block F1tXX_sem2d hdr contains the information needed to read the other fault data files Its format line by line is NPTS NDAT NSAMP DELT name of parameters Value of parameters above Name of fields exported in F1tXX_sem2d dat separated by XPTS ZPTS name of coordinate axis a F WwW N Fe from here to the end of file a two column table of coordinates of the output fault nodes Seismic Unix is freely available from the Colorado School of Mines at http timna mines edu cwpcodes 4 6 Outputs their visualization and manipulation 62 3500 r r r 3000 E E 2500 E Ne J K E 1500 E Q 1000 F E 500 E VO 0 500 fi fi fi fi 0 5 10 15 20 25 Time secs Figure 4 5 Sample seismograms from EXAMPLES UsingEMC2 generated with POST seis_plot csh F1tXX_sem2d dat contains the space time distribution of fault data such as slip slip rate stress and strength Every DELT seconds a block of fault data values is written The total number of blocks is NSAMP Each block has NDAT lines one per fault data field and NPTS columns one per fault node Stresses are relative to their initial values F1tXX_init_sem2d tab contains the spatial distribution o
37. can afford Estimates of total CPU time are printed at the end of check mode Details about memory usage can be found in MemoryInfo_sem2d txt 4 5 Verifying the settings and running a simulation 55 Stability 7 7 7 r 7 350 3000 300 2000 1000 250 0 200 1000 2000 150 3000 aol 4000 l i i 50 1 164 1 162 1 16 1 158 1 156 1 154 x 10 0 SS i 1 0 5 0 0 5 1 Resolution 7 7 7 7 400 3000 350 2000 300 1000 250 0 1000 200f 2000 150 3000 100 4000 50 F 1 164 1 162 1 16 1 158 1 156 1 154 x 10 0 5 0 6 0 4 0 2 0 0 2 0 4 0 6 Figure 4 3 Checking the quality of a mesh with PRE ViewMeshQuality m for the example in EXAMPLES UsingEMC2 The balance of the stability and resolution properties of the mesh can be analyzed logarithmic stability index top and logarithmic resolution index bottom Histograms of these indices in number of elements are shown on the right 4 5 Verifying the settings and running a simulation 56 The quality of the mesh can be inspected with the Matlab script PRE ViewMeshQuality m which produces plots like Figure 4 3 The proper balance of the mesh with respect to the following two criteria can be analyzed e Stability criterion related to the largest stable timestep The stability of each element is quantified by S min Az cp We also define a stability index as SI log S median S where the median value is
38. cation x z of the source m file name none Name of file containing source parameters The file format is ASCII with one line per source and 2 3 or 4 columns per line 1 X position in m 2 Z position in m 3 time delay in seconds 4 relative amplitude If column 4 is absent amplitude 1 If columns 3 and 4 are absent delay O and amplitude 1 NAME SRC_DOUBLE_COUPLE GROUP SOURCE MECHANISM PURPOSE Define a double couple source SYNTAX amp SRC_DOUBLE_COUPLE dip dip dble 90 Dip angle in degrees clockwise from the positive X direction 4 4 Input Blocks Reference Guide 49 NOTE NOTE NOTE Sign convention if the source amplitude is positive the right block moves up positive Z direction in PSV and forward positive Y direction in SH The source time function gives the cumulative seismic moment Mo t NOT the seismic moment rate The seismic moment Mo must be rescaled because a 2D point source is equivalent to a 3D line source A proper scaling is obtained by dividing the original 3D moment by the characteristic size of the rupture area in the off plane dimension An approximate scaling for a fault area with aspect ratio close to unity is Mo_2D Mo_3D dtau 2 3 dtau where dtau is the stress drop typically a few MPa NAME GROUP PURPOSE SYNTAX Mxx Mxz Myx Myz SRC_MOMENT SOURCE MECHANISM Define a moment tensor source amp SRC_MOMENT Mxx Mxz M
39. ch 1997 under the direction of Prof Jean Pierre Vilotte at the Institut de Physique du Globe de Paris IPGP The elastic anisotropic solver and several significant improvements were added by D Komatitsch later as part of a research contract with DIA Consultants Further functionalities were added by Jean Paul Ampuero as part of a Ph D thesis Ampuero 2002 also directed by Prof Vilotte at IPGP Most of these additional features were motivated by an ECOS NORD FONACYT research project for the study of the seismic response of the valley of Caracas Venezuela That became the version 1 0 of the SEM2DPACK released in April 2002 For version 2 0 most of the solver was rewritten in preparation for the implementation of higher level functionalities Developments for the simulation of earthquake dynamics Am puero 2002 were included in the main branch of SEM2DPACK in October 2003 version 2 2 Spontaneous rupture along multiple non planar faults can be currently modelled with a range of friction laws Non linear inelastic materials were introduced in March 2008 version 2 3 Damage and visco plastic rheologies are included especially for the modeling of earthquake rupture with off fault dissipation 1 3 Application examples SEM2DPACK has been utilized in a variety of applications in Earth sciences and has con tributed to more than 20 publications dynamic rupture on non planar faults and seismic wave radiation Madariaga et al 2006
40. d state dependent friction Some friction types can be combined E g to set the friction coefficient to the minimum of SWF and TWF set friction SWF TWF cohesion dble 0d0 part of the strength not proportional to normal stress opening log T Allow fault opening instead of tensile normal stress Tn dble 0d0 Initial normal traction positive tensile Tt dble 0d0 Initial tangent traction positive antiplane y gt 0 Sxx dble 0d0 Initial stress sigma_xx Sxy dble 0d0 Initial stress sigma_xy Sxz dble 0d0 Initial stress sigma_xz Syz dble 0d0 Initial stress sigma_yz Szz dble 0d0 Initial stress sigma_zz otd dble 0 d0 Time lag between outputs in seconds Internally adjusted to the nearest multiple of the timestep Its value can be found in the output file F1tXX_sem2d hdr The default internally resets otd timestep oti dble 0 d0 Time of first output in seconds Internally adjusted to the nearest multiple of the timestep Its value can be found in the output file F1tXX_sem2d hdr oxi int 3 1 huge 1 First last node and stride for output The default resets oxi 2 last fault node osides log F Export displacement and velocities on each side of the fault NOTE the initial stress can be set as a stress tensor Sxx etc as initial tractions on the fault plane Tn and Tt or as the sum of both NOTE we recommend to use dynamic faults with the leapfrog time scheme and a layer of Kelvi
41. e not been tested 1 6 Installation 8 1 6 Installation 1 Uncompress and expand the SEM2DPACK package tar xvfz sem2dpack tgz 2 Go to the source directory cd SEM2DPACK SRC 3 Edit the Makefile according to your FORTRAN 95 compiler following the instructions therein 4 Modify the optimization parameters declared and described in SRC constant f90 5 Compile make 6 Move to the SEM2DPACK POST directory edit the Makefile and compile On normal termination you should end up with a set of executable files among which sem2dsolve in home yourhome bin 1 7 Documentation Documentation is available through the following resources This User s Manual The EXAMPLES directory contains several examples some have a README file The pre processing and post processing tools for Matlab are documented through Mat lab s help For instance help mesh2d provides an overview of the MESH2D utilities and help mesh2d_wedge provides detailed documentation for the wedge meshing function The ToDo file contains a list of known issues 1 8 Support Support for users of SEM2DPACK is available through a tracking system at http sourceforge net tracker group_id 182742 Three separate tracker lists deal with the following aspects Feature Requests requests for implementation of new features Support Requests questions related to the usage of SEM2DPACK Bugs bug reports Before submitting an issue make sure
42. e topography SYNTAX amp MESH_LAYERED xlim zmin nx file nlayer ezflt fztag xlim dble 2 none X limits of the box min and max zmin dble none bottom Z limit of the box nx int 1 Number of elements along the X direction Not needed if ztopH QSPLINE in a amp MESH_LAYER block file string Only for flat layers name of ASCII file containing layer parameters one line per layer listed from top to bottom 3 columns per line 1 vertical position of top boundary 2 number of elements along Z direction 3 material tag nlayer int mone Number of layers If a file name is not given the layer parameters must be given immediately after the amp MESH_LAYERED block by nlayer amp MESH_LAYER input blocks one for each layer listed from top to bottom 4 4 Input Blocks Reference Guide 44 ezflt int 0 introduce a fault between the ezflt th and the ezf1t 1 th element rows numbered from bottom to top If ezflt 0 default no fault is introduced If ezflt 1 a horizontal fault is introduced at near the middle of the box ezflt is reset to int nelem 2 2 fztag int 0 tag for elements near the fault Useful to set a damping layer near the fault fznz Lint 1 vertical size of near fault layer half thickness in number of elements NOTE the following tags are automatically assigned to the boundaries 1 Bottom Right Top Left Fault lower side Fault upper side aon F WN NAME MESH_LAYER GROU
43. ect 1 Mesh generation partition the domain into deformed quadrilateral elements Whereas no general mesh generation code is included SEM2DPACK contains basic meshing util ities for structured and semi structured grids and can import unstructured quadrilateral meshes generated externally These features are described in Chapter 3 2 Mesh quality verification check the accuracy stability and computational cost applying tools described in Section 4 5 Return to previous step if needed 3 Numerical simulation run the SEM2D solver Chapter 4 explains its usage 4 Post processing visualization and analysis of the output A number of post processing and graphic tools are included as decribed in Section 4 6 Outputs are in the form of binary data files ASCII data files and PostScript figure files Scripts are provided for graphic display and analysis on Seismic Unix Gnuplot and Matlab We recommend the usage of the Matlab functions included see Section 4 6 8 This is a research code constantly under development and provided as is and therefore it should not be considered by the user as a 100 bug free software package We welcome 1 2 History and credits 6 comments suggestions feature requests bug reports see Section 1 8 and contributions to the code itself see Section 1 9 1 2 History and credits The main part of the elastic isotropic solver was written by Dimitri Komatitsch as part of his Ph D thesis Komatits
44. ed to guide you through this process We will not give here a comprehensive description of the code architecture only enough details to get you started in performing safely the most usual and evident modifications 5 1 Overview of the code architecture in progresss This code uses a mixture of procedural imperative and object oriented paradigms Histori cally it evolved from a purely procedural code Extensive use of modularity Object Oriented Programming OOP features principles applied in this code encapsula tion classes static polymorphism These are not applied everywhere in the code for different reasons reusage of legacy code performance difficulty related to the limits of Fortra 90 or sections of code yet to be updated Added cost of structures containing pointer components the possibility of pointer aliasing prevents more agressive compiler optimizations and adds overhead for safety checks 5 2 Accessible areas of the code Some areas of the code have been written in such a way that a moderately experienced Fortran 95 programmer with a limited understanding of the code architecture can introduce new fea 5 2 Accessible areas of the code 66 tures without breaking the whole system This is achieved through modularity encapsulation and templates The modifications that are currently accessible are e boundary conditions see bc_gen 90 e material rheology see mat_gen 90 e source time functions
45. ents on each side of the fault and n At faut 0 1 to 0 3 where Atfauiz is the critical time step size based on the size of the spectral elements along the fault not necessarily equal to the critical time step over the whole mesh The value of At gyi can be obtained with the Matlab function PRE critical_timestep nm 2 2 3 Coulomb plasticity and visco plasticity Perfect plasticity Perfect plasticity with a Coulomb yield function is implemented for 2D plane strain as in Andrews 2005 The total strain is the sum of an elastic and a plastic contribution e e e The plastic strain is assumed to be purely deviatoric 0 Plastic yield occurs when the maximum shear stress over all orientations Tmax o2 ale xz Oy 14 2 5 reaches the yield strength Y c cos ose ozz sin 2 2 6 where c is the cohesion and is the internal friction angle 2 2 Material rheologies 14 Visco plasticity In classical Duvaut Lions visco plasticity the visco plastic strain rate is proportional to the excess of stress over the yield strength p 1 Tij o E 2 7 Ekl on Ts 2 7 Tmax where T is the visco plastic relaxation time 1 x a 2 is the ramp function and Tij Tij 30 dij is the deviatoric part of the stress tensor Visco plasticity is often employed as a regularization of plasticity to avoid or delay the occur rence of strain localization features such as shear bands that are mesh de
46. esh for a single quadrilateral domain possibly with curved sub horizontal boundaries and curved sub horizontal layer interfaces The domain can be cut in the horizontal direction by a single fault possibly curved or kinked For further details see the Reference Guide for the input blocks MESH_CART and MESH_LAYERED in Section 4 4 3 3 Meshing with the MESH2D Matlab utilities A number of Matlab functions for 2D meshing are provided in POST mesh2d These can generate structured meshes for quadrilateral domains with curved boundaries and merge several such meshes to generate a more complicated globally unstructured mesh Functions for manipulating visualizing and exporting these meshes are included Here is an overview of available tools SEM2DPACK PRE mesh2d provides Matlab utilities for the generation manipulation and visualization of structured 2D quadrilateral meshes and unstructured compositions of structured meshes The mesh database is stored in a structure described in MESH2D_TFI s help Mesh generation MESH2D_TFI generates a structured mesh by transfinite interpolation 3 4 Generating a mesh with EMC2 21 MESH2D_QUAD generates a structured mesh for a quadrilateral domain MESH2D_CIRC_HOLE generates a mesh for a square domain with a circular hole MESH2D_WEDGE generates a mesh for a triangular wedge domain MESH2D_EXO mesh for a vertical fault MESH2D_EX1 mesh for a shallow layer over half space with dipping fault Mesh mani
47. f initial shear stress initial normal stress and initial friction 3 columns F1tXX_potency_sem2d tab contains time series of seismic potency and potency rate The seismic potency tensor p is defined by the following integral along the fault 1 Pij F n Auj njAui dx 4 1 fault 2 where Au is slip and n is the local unit vector normal to the fault The file contains one line per timestep In SH ndof 1 each line has 4 columns 2 components of potency pig and p23 and 2 components of potency rate p13 and p23 In P SV ndof 2 each line has 3 components of potency p11 p22 and pi2 and 3 components of potency rate P11 P22 and p12 Some tools are available to manipulate the data in F1tXX_sem2d dat The script F1tXX_sem2d csh shows how to extract ASCII time series of different fields at given locations on the fault using Seismic Unix tools The actual number of columns is NPTS 2 Fortran adds a one word tag at the front and end of each record 4 6 Outputs their visualization and manipulation 63 The program post_fault exe performs basic manipulations of the fault data including conversion to an ASCII file readable by gnuplot Its interactive menu is self explanatory The script POST sample_fault m and function POST sem2d_read_fault m show how to manipulate and plot fault data in Matlab 4 6 6 Stress glut For damage and plastic materials the solver can export the plastic and damage components of
48. f the user prefers to sacrifice accuracy by using a non adapted structured mesh i e a logically cartesian mesh where the element faces do not necessarily follow the material interfaces the basic built in meshing capabilities of the solver SEM2D described in Section 3 2 are sufficient 3 2 Meshing features included in the solver 20 2 Moderately complicated meshes can be generated with the included Matlab tools de scribed in Section 3 3 3 Adapted meshes for more complicated geological models must be generated with some external software As an illustration the usage of the freely available 2D mesh genera tion software EMC2 is described in Section 3 4 SEM2DPACK provides only basic meshing capabilities and does not include an unstructured mesh generator for complicated realistic geological models This chapter describes how to achieve that task with an external software EMC2 Generating a high quality unstructured quad mesh can be a time consuming task Let s note that for wave propagation problems without dynamic faults if the acceptable accuracy is low or large computational resources are available to work with a very fine mesh a structured mesh in which the element faces do not necessarily follow the material interfaces can be generated with the basic built in meshing capabilities of SEM2D 3 2 Meshing features included in the solver The solver itself has very limited meshing capabilities It can only generate a structured m
49. frog Number of time steps e NbSteps will be set later Time step increment Dt will be set later Courant number Courant 0 30 Total simulation duration TotalTime 35 000E 00 Materia Properties Number of materials 1 Material index tag 1 Material type A kind Elastic P wave velocity lt sa w a 0p 1 732E 00 S wave velocity 08 1 000E 00 Mass density rho 1 000E 00 Poisson s ratio 7 First Lame parameter Lambda Second Lame parameter Mu Bulk modulus K Young s modulus E Tondi ti I w o undary ons Boundary tag Boundary condition Type of absorbing boundary Allow incident wave Ane Boundary tag Boundary condition Type of absorbing boundary Allow incident wave et urces X position meters Y position meters Source time function Fundamental frequency Time delay s e Multiplying factor Source Type If P SV Hz ceivers coord 1 coord 2 let_wave let_wave We E 0 onset E counterclockwise angle A up 250 021E 03 000E 00 000E 00 667E 00 500E 00 Nere 2 ABSORB Clayton Engquist T 3 ABSORB Clayton Engquist T 0 000E 00 0 000E 00 Ricker 500 000E 03 3 000E 00 250 000E 03 Collocated Force 0 00 Number of receivers y ate a ak number 7 Subsampling for seismograms recording isamp 1 Field recorded 5
50. g a mesh with EMC2 23 STEP II a with the mouse POINT mouse click in figure window b by coordinates POINT xy pt x y This is the safest way to get really vertical and horizontal boundaries needed for the absorbing conditions in SPECFEM90 You probably need to get the coordinates of an existing reference point POINT lt QUERY lt point click on point c you can also reload another point file 12 Delete points POINT lt DESTRUCT lt point click on point Now you must define the geometry of the domains These macro blocks are intended to be internally meshed by deformed quadrilaterals Their geometry follows the geometry of the geological model one domain per material Each domain must be bounded by segments or splines Segments SEGMENT point 2 click extreme point Splines SPLINE point click point You will see the spline evolve as you click points PREPARE defines the properties of the discrete spectral element mesh 1 Switch to the preparing mesh tool lt PREP MESH lt Define domains with rock n DOMAIN REF n any click inside domain You will see the domains edges get colored and the domains get numbered with n At any moment you can decide to show or not the domain decomposition To hide the domain decomposition gt REFRESH gt Show the domain decomposition SHOW ALL Remove a domain definition REMOVE DOMAINE any cl
51. ge rheology of Lyakhovsky Ben Zion and Agnon J Geophys Res 1997 and Hamiel et al Geophys J Int 2004 SYNTAX amp MAT_DAMAGE cp cs rho phi alpha Cd R e0 ep cp dble 0d0 P wave velocity cs dble 0d0 S wave velocity rho dble 0d0 density phi dble 0d0 internal friction angle alpha dble 0d0 initial value of damage variable Cd dble 0d0 damage evolution coefficient R dble 0d0 damage related plasticity coefficient Cv normalized by the inverse of the intact shear modulus e0 dble 3 040 initial total strain 11 22 and 12 ep dble 3 0d0 initial plastic strain 11 22 and 12 NAME MAT_ELASTIC GROUP MATERIALS PURPOSE Set material properties for a linear elastic medium SYNTAX For isotropic material amp MAT_ELASTIC rho rhoH cplcpH cslcsH followed by DIST_XXXX blocks for arguments with suffix H if present in the same order as listed above For transverse anisotropy with vertical symmetry axis amp MAT_ELASTIC rho c11 c13 c33 c55 cp dble 0d0 P wave velocity cs dble 0d0 S wave velocity rho dble 0d0 density c11 c13 c33 c55 dblel 0d0 anisotropic elastic moduli NAME MATERIAL PURPOSE Define the material type of a tagged domain SYNTAX amp MATERIAL tag kind followed by one or two MAT_XXXX input blocks tag int none Number identifying a mesh domain kind name 2 ELAST Material types ELAST DMG PLAST KV 4 4 Input
52. guarantee smooth boundary conditions on curved faults SYNTAX amp QC_SPLINE file file string Name of ASCII file containing information of all the element vertex nodes lying on the boundary curve One line per node ordered by increasing x 3 columns per line 1 x position 2 z position 3 derivative dz dx of the curve at the node All QC_SPLINE curves in a mesh must have the same number of nodes The parameter nx in amp MESH_LAYERED is automatically reset nx number of nodes in QC_SPLINE 1 NAME MESH_MESH2D GROUP MESH_DEF PURPOSE Imports a mesh in mesh2d format as defined by the PRE mesh2d mesh generator tools for Matlab SYNTAX amp MESH_MESH2D file file name none Name of the MESH2D file including suffix The format of this file is NEL NPEL NNOD NBC 1 line with 4 integers nb of elements nodes per element total nb of nodes nb of boundaries NID X Y NNOD lines one per node with 1 integer and 2 reals node id x y EID NODES TAG NEL lines one per element with NPEL 2 integers element id NPEL node ids tag BCTAG NBEL 2 integers boundary tag nb of boundary elements BID EID EDGE repeat for each of NBEL lines one per boundary element 3 integers the NBC boundaries boundary element id bulk element id edge id 4 4 Input Blocks Reference Guide 46 NAME SNAP_DEF GROUP SNAPSHOT_OUTPUTS PURPOSE Set preferences for exporting snapshots SYNTAX am
53. hs J Geophys Res 113 B9 B09317 doi 10 1029 2007JB005553 Komatitsch D 1997 M thodes spectrales et l ments spectraux pour l quation de V lastodynamique 2D et 3D en milieu h t rog ne Ph D thesis Institut de Physique du Globe de Paris Paris Komatitsch D C Barnes and J Tromp 2000 Simulation of anisotropic wave propagation based upon a spectral element method Geophys 65 4 1251 1260 doi 10 1190 1 1444816 Komatitsch D and J Tromp 1999 Introduction to the spectral element method for 3 D seismic wave propagation Geophys J Int 139 806 822 Komatitsch D and J P Vilotte 1998 The spectral element method an efficient tool to simulate the seismic response of 2D and 3D geological structures Bull Seis Soc Am 88 368 392 Komatitsch D J P Vilotte R Vai and F J S nchez Sesma 1999 The Spectral Element method for elastic wave equations application to 2D and 3D seismic problems Int J Num Meth Engng 45 9 1139 1164 Lyakhovsky V Y Ben Zion and A Agnon 1997 Distributed damage faulting and friction J Geophys Res 102 B12 27635 27649 Lyakhovsky V Y Hamiel J P Ampuero and Y Ben Zion 2009 Non linear dam age rheology and wave resonance in rocks Geophys J Int 178 2 910 920 doi 10 1111 3 1365 246X 2009 04205 x Madariaga R J P Ampuero and M Adda Bedia 2006 Seismic radiation from simple models of earthquakes In A McGarr R Abercro
54. ick inside domain WARNING corrections to the domain decomposition are sometimes displayed only after refreshing the figure window Now you must define the subdivision of each domain in quadrilateral finite elements Define the number n of elements on each edge NB INTERVAL n any click edge You will see the intermediate points appear The number of intervals n is mainly dictated by the resolution criterion elements should be smaller than the smallest wavelength you want to propagate Moreover a domain can 3 4 Generating a mesh with EMC2 24 STEP III be quadrangulated only if the total number of intervals along its perimeter is even the sum of all n along its boundaries However a quality mesh is not always guaranteed and you need to proceed by trial and error emc2 allows you to jump back and forth between the different steps of the meshing procedure Finally you must define the external boundaries of the modelled region which will have a special treatment You must associate a tag a number to each absorbing boundary No convention is assumed but you should remember those tags later when setting the boundary conditions in SEM2D It is also useful to assign a tag to the free surface boundary that will be eventually used by SEM2D to locate the receivers or sources Define a boundary with index n LINE REF n any click edge Of course each boundary can be composed of many domain edges Refresh the display
55. ist each item is followed by two informations within brackets The first bracketed information is the type of the argument double precision dble integer int logical log single character char fixed length word e g char 6 is a 6 characters word arbitrary length word name or vectors e g int 2 is a two element integer vector The second bracketed information is the default value of the argument Some arguments are optional or when absent they are automatically assigned the default values Some arguments have a second version with a suffix H that allows to set values that are spatially non uniform The H version of the argument must be set to the name of any of the input blocks of the DISTRIBUTIONS group The appropriate DIST_xxxx block must follow immediately For example to set the argument eta to a Gaussian distribution 4 3 General format of the input file 29 amp MAT_KV etaH GAUSSIAN amp DIST_GAUSSIAN length 1d6 100d0 ampli 0 1d0 Arguments that accept an H version are indicated in Section 4 4 When more than one H version argument is present the amp DIST_xxxx blocks must appear in the same order as in the argument list of Section 4 4 In the next section Input Block Reference Guide you should get acquainted with the syntax of the blocks you are most likely to use The mandatory or more important input blocks are amp GENERAL amp MESH_DEF followed by a MESH_Method block amp MATE
56. jects sem2d All versions of the code can be downloaded from the package repository at http sourceforge net projects sem2d files sem2dpack Taking full advantage of the convenient features offered by SourceForge subscribe to new release announcements submit and track bug reports requires a SourceForge net account which can be created at http sourceforge net account registration SEM2DPACK is updated regularly typically every few months To receive email notifications about new releases you must sign up for the Update Notifications in the project s main page scroll down a bit http sourceforge net projects sem2d 1 5 Requirements Compiling the solver code requires the make utility and a Fortran 95 compiler The code is being developed with the Intel compiler for Linux It works properly with the Intel compiler starting with version 8 0 046_pe047 1 so make sure you have a recent version of ifort Other compilers are not being tested on a regular basis so please report any related problems The solver runs under the Linux operating system In particular input output file name conventions are specific to Linux Other operating systems have not been tested Pre processing and post processing tools including graphic visualization are provided for Seismic Unix Gnuplot GMT and Matlab The included Matlab tools are by far the most complete so a Matlab license is highly recommended Matlab clone softwares hav
57. lastic medium eta 0 In terms of the normalized viscosity if ETAxDT T dtc_kv dtc sqrt 1 2 eta MAT_PLASTIC MATERIALS PURPOSE Set material properties for elasto plastic material with Mohr Coulomb yield criterion and non dilatant null volumetric plastic strain amp MAT_PLASTIC cp cs rho phi coh Tv e0 dble 040 P wave velocity dble 040 S wave velocity dble 0d0 density 4 4 Input Blocks Reference Guide 42 phi coh Tv e0 NAME GROUP PURPOSE SYNTAX xlim zlim nelem ezflt fztag fznz FaultX NOTE the NAME dble 040 internal friction angle dble 0d0 cohesion dble 0d0 visco plastic relaxation time dble 3 0d0 initial total strain 11 22 and 12 MESH_CART MESH_DEF Rectangular box with structured mesh amp MESH_CART xlim zlim nelem ezflt fztag FaultX dble 2 none X limits of the box min and max dble 2 none Z limits of the box min and max int 2 none Number of elements along each direction int 0 introduce a horizontal fault between the ezflt th and the ezflt 1 th element rows Rows are numbered from bottom to top starting at ezflt 1 If ezflt 0 default no fault is introduced inside the box for symmetric problems a fault can still be set at an external boundary If ezflt 1 a fault is introduced at near the middle of the box ezflt is reset to int nelem 2 2 int 0 fault zone tag for elements close to the fault
58. lpha scheme of 4 4 Input Blocks Reference Guide 53 G M Hulbert and J Chung 1996 Explicit time integration algorithms for structural dynamics with optimal numerical dissipation Comp Methods Appl Mech Engrg 137 175 188 by setting alpha_m 0 and alpha 1 alpha_f 4 5 Verifying the settings and running a simulation 54 4 5 Verifying the settings and running a simulation Once the code has been successfully compiled the simulation can be started by typing sem2dsolve from your working directory which contains the file Par inp The computa tions can be run in background and the screen output saved in a file e g info by typing sem2dsolve gt info amp A typical screen output of SEM2D corresponding to the first example is shown on the following pages The parameters of the simulation and some verification information are reported there in a self explanatory form You are advised to do a first run with iexec 0 in the GENERAL input block and check all these informations prior to the real simulation You should always verify the following e Stability the CFL stability number should be smaller than 0 55 0 60 for second order time schemes and much smaller for highly deformed meshes see Section 6 1 on Instabilities in very distorted elements This number is defined at each computational node as CFL cp AN Az where At is the timestep cp the P wave velocity and Az the local grid spacing Note that Az is usually much s
59. maller than the element size h Ng11 times smaller because SEM internally subdivides each element onto a non regular grid of Ng11xNg11 nodes clustered near the element edges Gauss Lobatto Legendre nodes If the computation is unstable the maximum displacement printed every ItInfo time steps increases exponentially with time Stability can be controlled by decreasing Dt or Courant in Par inp e Resolution the number of nodes per shortest wavelength Amin should be larger than 4 5 5 The minimum wavelength is defined as Amin min cg fmaz where cg is the S wave velocity and fmax the highest frequency you would like to resolve e g the maximum frequency at which the source spectrum has significant power for a Ricker wavelet fimax 2 5 x fo For an element of size h and polynomial order p Ng11 1 the number of nodes per wavelength G is P Amin a Typical symptoms of poor resolution are ringing and dispersion of the higher frequencies However in heterogeneous media these spurious effects might be hard to distinguish from a physically complex wavefield so mesh resolution must be checked beforehand If resolution is too low the mesh might be refined by increasing Ng11 in Par inp p refinement or by generating a denser mesh h refinement If you were using EMC2 as a mesh generator the script PRE href csh can be useful for h refinement G Cost the total CPU time an memory required for the simulation are as much as you
60. mbie H Kanamori and G di Toro Eds Earthquakes Radiated Energy and the Physics of Earthquake Faulting Volume 170 of Geophysical Monograph pp 223 236 Am Geophys Union Meng L and J P Ampuero 2012 Slow rupture and weakly pressure sensitive strength enables compressional branching Dynamic rupture simulations of the 2012 off Sumatra earthquake Geophys Res Lett submitted BIBLIOGRAPHY 71 O Brien G and C Bean 2011 An irregular lattice method for elastic wave propagation Geophys J Int 187 3 1699 1707 doi 10 1111 j 1365 246X 2011 05229 x Patera A 1984 A spectral element method for fluid dynamics laminar flow in a channel expansion J Comp Phys 54 468 488 Roberts A W R S White and P A F Christie 2009 Imaging igneous rocks on the North Atlantic rifted continental margin Geophys J Int 179 2 1024 1038 doi 10 1111 j 1365 246X 2009 04306 x Rubin A M and J P Ampuero 2007 Aftershock asymmetry on a bimaterial interface J Geophys Res 112 B05307 doi 10 1029 2006JB004337 Smith S and R Snieder 2010 Seismic modeling and analysis of a prototype heated nuclear waste storage tunnel Yucca Mountain Nevada Geophys 75 1 T1 T8 doi 10 1190 1 3273868 Stacey R 1988 Improved transparent boundary formulations for the elastic wave equa tion Bull Seism Soc Am 78 6 2089 2097 Vai R J M Castillo Covarrubias F J S nchez Sesma D Komatitsch and J P
61. mesh design The interfaces between different materials at which sharp contrasts of material properties occur should preferably coincide with faces of the elements This is sometimes called an adapted mesh and is the only way to preserve spectral accuracy at material interfaces Fault planes across which displacement discontinuities occur must coincide with element faces Faults are implemented with a split node formulation Elements can be deformed but extremely small and extremely large angles between faces of a same element must be avoided This would penalize both accuracy and stability The linear size of the elements must be small enough so that each element contains enough computational nodes per minimum wavelength and each fault boundary element contains enough nodes per rupture process zone Unnecessarily small elements should be avoided they penalize the stability of the method Generating high quality quad meshes for complicated geological models is not yet a fully automated process and can be very time consuming Iterations between mesh generation and mesh quality check are sometimes required The last two constraints above are addressed more quantitatively in Section 4 5 Mesh quality assessment tools are also presented in Section 4 5 The remainder of this chapter describes three possible ways to generate quad meshes by order of complexity 1 If the geometrical structure of the model is simple or i
62. mplemented for vertical and horizontal boundaries NAME BC_DIRNEU GROUP BOUNDARY_CONDITION PURPOSE Dirichlet null displacement and or Neumann null or time dependent traction boundary conditions on vertical or horizontal boundaries SYNTAX amp BC_DIRNEU h v hsrc vsrc possibly followed by one or two STF_XXXX blocks h char N Boundary condition on the horizontal component v char N Boundary condition on the vertical component N Neumann 7D Dirichlet hsrc name none Name of the source time function for a time dependent horizontal traction RICKER TAB USER etc see STF_XXXX input blocks vsrc name none Same for the vertical component NAME BC_DYNFLT GROUP BOUNDARY_CONDITION DYNAMIC_FAULT PURPOSE Dynamic fault with friction SYNTAX amp BC_DYNFLT friction cohesion cohesionH opening Tn TnH Tt TtH 4 4 Input Blocks Reference Guide 33 Sxx SxxH Sxy SxyH Sxz SxzH Syz SyzH Szz SzzH oti otd oxi osides followed in order by 1 amp DIST_XXX blocks from the DISTRIBUTIONS group for arguments with suffix H if present in the order listed above 2 amp BC_DYNFLT_SWF amp BC_DYNFLT_TWF or amp BC_DYNFLT_RSF block s if absent default values are used 3 amp BC_DYNFLT_NOR block if absent default values are used friction name 2 SWF Friction law type SWF slip weakening friction TWF time weakening friction RSF rate an
63. n amp BC_DYNFLT_SWF Dc DcH MuS MuSH MuD MuDH healing followed by amp DIST_XXX blocks from the DISTRIBUTIONS group for arguments with suffix H if present in the order listed above Lint 1 Type of slip weakening function 1 linear 2 exponential dble 0 5d0 Critical slip dble 0 6d0 Static friction coefficient 4 4 Input Blocks Reference Guide 35 MuD dble 0 5d0 Dynamic friction coefficient healing log F Instantaneous healing upon slip arrest Healing is currently valid only with the leapfrog time scheme NAME BC_DYNFLT_TWF GROUP DYNAMIC_FAULT PURPOSE Time weakening friction for dynamic faults with prescribed rupture speed SYNTAX amp BC_DYNFLT_TWF kind MuS MuD Mu0 X Z V L T kind int 1 Type of time weakening history 1 expansion at constant speed V up to time T 2 expansion at decreasing speed then contraction as in Andrews and Ben Zion JGR 1997 eqs 2 and 3 MuS dble 0 6d0 Static friction coefficient MuD dble 0 5d0 Dynamic friction coefficient Mu0 dble 0 6d0 Friction coefficient at the hypocenter at time 0 X Z dble 0d0 Position of hypocenter V dble 1d3 Rupture propagation speed initial speed if kind 2 L dble 1d0 Size of weakening zone T dble huge Total duration NOTE Time weakening is usually applied as an artificial nucleation procedure The maximum size of the nucleation region is 2 V T if kind 1 V T 2 if kind 2 NAME BC_DEF PURP
64. n Voigt damping material near the fault NAME BC_DYNFLT_NOR GROUP DYNAMIC_FAULT 4 4 Input Blocks Reference Guide 34 PURPOSE SYNTAX kind NAME GROUP kind NAME GROUP kind MuS Normal stress response for dynamic faults amp BC_DYNFLT_NOR kind V L T int 1 Type of normal stress response O shear strength is independent of normal stress the cohesive strength is set as the product of friction coefficient and initial normal stress 1 Coulomb 2 Prakash Clifton with regularizing time scale 3 Prakash Clifton with regularizing length scale dble 1d0 Regularization time scale if kind 2 dble 1d0 Characteristic velocity if kind 3 dble 1d0 Regularization length scale if kind 3 BC_DYNFLT_RSF DYNAMIC_FAULT PURPOSE SYNTAX Velocity and state dependent friction amp BC_DYNFLT_RSF kind Dc DcH Mus MusH a aH b bH Vstar VstarH followed by amp DIST_XXX blocks from the DISTRIBUTIONS group for arguments with suffix H if present in the order listed above int 1 Type of rate and state friction law 1 strong velocity weakening at high speed as in Ampuero and Ben Zion 2008 dble 0 5d0 Critical slip dble 0 6d0 Static friction coefficient dble 0 01d0 Direct effect coefficient dble 0 02d0 Evolution effect coefficient dble 1d0 Characteristic or reference slip velocity BC_DYNFLT_SWF DYNAMIC_FAULT PURPOSE SYNTAX Slip weakening frictio
65. n the geometry seems too complicated for quad meshing you should consider simplifying the geometry especially those details that are much smaller than the dominant wavelength If the above fails or does not apply you have to help the mesher The recommended procedure in EMC2 is 1 Divide your original mesh into simple domains in such a way that most domains have exactly four sides possibly curved and the remaining non four sided domains are as small as possible 2 Generate a structured quad mesh a regular grid inside each four sided domain with the QUADRANGULATE tool of the PREP_MESH mode as described in section 5 2 13 of EMC2 s manual note that this is not the same as the lt QUADRANGULATE gt button in the EDIT_MESH mode 3 Proceed as usual triangulation followed by quadrangulation inside the remaining non four sided domains If these are small enough EMC2 should not have problems doing a correct tri to quad meshing Bibliography Ampuero J P 2002 Etude physique et num rique de la nucl ation des s ismes Ph D thesis Universit Paris 7 Denis Diderot Paris Ampuero J P and Y Ben Zion 2008 Cracks pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity weakening friction Geophys J Int 173 2 674 692 doi 10 1111 j 1365 246X 2008 03736 x Andrews D 1999 Test of two methods for faulting in finite difference calculations Bull Seis Soc Am 89 931 937 A
66. ndrews D J 2005 Rupture dynamics with energy loss outside the slip zone J Geophys Res 110 B1 B01307 doi 10 1029 2004JB003191 Ben Zion Y and V Lyakhovsky 2006 Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology Geophys J Int 165 1 197 210 Boaga J S Renzi G Vignoli R Deiana and G Cassiani 2012 From surface wave inversion to seismic site response prediction Beyond the 1d approach Soil Dyn Earthq Engng 36 38 51 Clayton R and B Engquist 1977 Absorbing boundary conditions for acoustic and elastic wave equations Bull Seis Soc Am 67 6 1529 1540 De la Puente J J P Ampuero and M K ser 2009 Dynamic rupture modeling on unstructured meshes using a discontinuous Galerkin method J Geophys Res 114 B10302 doi 10 1029 2008JB006271 De la Puente J M K ser M Dumbser and H Igel 2007 An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes IV Anisotropy Geophys J Int 169 3 1210 1228 doi 10 1111 3 1365 246X 2007 03381 x Dewangan P I Tsvankin M Batzle K Van Wijk and M Haney 2006 PS wave move out inversion for tilted TI media A physical modeling study Geophys 71 4 D135 D143 doi 10 1190 1 2212274 Gabriel A A J P Ampuero L A Dalguer and P M Mai 2012a Source properties of dynamic rupture pulses with off fault plasticity Geophys Res Let submitted G
67. om the sediment mesh Similar information is plotted by gv Stability sem2d ps and gv Resolution sem2d ps The indices in these files are however not logarithmic and are not normalized by the median 3In future releases of SEM2DPACK this penalty on computational efficiency will be reduced by non conformal meshing with mortar elements by timestep subcycling or by implicit explicit timestep partitioning Printed by Jean Paul Ampuero Mar 06 08 18 04 info Page 1 4 Mar 06 08 18 04 info Page 2 4 SPECFEM Program 06 03 2008 Time KAKAK I I RI I II A I I Ie 1 NP u E phase Ed ORO I I RI II A A Parameter o Execution mode P iexec solve Number of nodes per edge ng11 6 Number of d o f per node ndof 1 Highest frequency to be resolved fmax 1 250E 00 Print progress information during input phase verbose 1 T initialization phase verbose 2 T checking phase verbose 3 T solver phase 4 verbose 4 T Frequency for solver progress information itInfo 1000 Mesh Generation Method method CARTESIAN Minimum X xlim 1 0 000E 00 Maximum X xlim 2 3 000E 01 Minimum Z zlim 1 0 000E 00 Maximum Z zlim 2 3 000E 01 Number of elements along x nelem 1 60 Number of elements along Z nelem 2 60 Cut by horizontal fault faultx F Time integration Scheme e kind leap
68. p SNAP_DEF iti itd fields components bin visual3 avs ps gmt Followed by a amp SNAP_PS block if ps T itl int 0 Time step of first snapshot output itd Lint 100 Number of timesteps between snapshots fields char V fields to export in snapshots the prefix of the output file names is given in parenthesis pe displacements dx dy dz da V velocity vx vy vz va A acceleration ax ay az aa E strain e11 e22 e12 e23 e13 AS stress s11 522 812 533 e13 e23 d divergence rate dvx dx dvz dz HE curl rate dvx dz dvz dx components char ya components for PostScript outputs in P SV x z and or a amplitude y is ignored in SH y only Other values are ignored ps log T PostScript see amp SNAP_PS input block gmt log F output triangulation file grid_sem2d gmt to be used in pscontour T of the General Mapping Tool GMT avs log F AVS only for D V and A fields visual3 log F Visual3 only for D V and A fields bin log T binary NOTE E and S fields are exported only as binary NAME SNAP_PS GROUP SNAPSHOT_OUTPUTS PURPOSE Preferences for PostScript snapshots SYNTAX amp SNAP_PS vectors mesh background color isubsamp boundaries symbols numbers legend ScaleField Interpol DisplayPts vectors log F Plots a vectorial field with arrows mesh log F Plots the mesh on background background char Filled b
69. pendent For that particular application T is typically set to the average P wave traveltime across a few grid points i e a few times the average spacing between GLL nodes divided by the P wave speed 2 2 4 Continuum damage The continuum damage formulation by Lyakhovsky et al 1997 including damage related plasticity as introduced by Hamiel et al 2004 is implemented with modifications for 2D plane strain The first and second invariants of the 2D elastic strain tensor are defined as I ej and In ejej respectively A strain invariant ratio is defined as 1 YT2 The following non linear stress strain relation is assumed Lyakhovsky et al 1997 eq 12 Tij A y 105 2u VE Fy 2 8 where y is an additional elastic modulus The elastic moduli depend on a scalar damage variable 0 lt a lt 1 through Lyakhovsky et al 1997 eq 19 A Ao 2 9 u uo Yr o a 2 10 Y Wa 2 11 where Ay and uo are Lam s parameters for the intact material a 0 The parameter o is the threshold value of the strain invariant ratio at the onset of damage It is related to the internal friction angle in a cohensionless Mohr Coulomb yield criterion by the 2D plane strain version of Lyakhovsky et al 1997 eq 37 2 AAA 2 12 o 1 A0 po 1 sin o The scaling factor y is chosen such that convexity is lost at 1 when amp It is derived from the 2D plane strain version of Lyakho
70. pulation MESH2D_ROTATE rotates the node coordinates MESH2D_TRANSLATE translates the node coordinates MESH2D_MERGE merges several meshes into a single mesh MESH2D_WRITE writes a 2D mesh database file mesh2d MESH2D_READ reads a 2D mesh database from a mesh2d file Reading mesh files from other mesh generation software READ_DCM reads a 2D mesh in the DCM format of EZ4U http www lacan upc es ez4u htm READ_INP reads a 2D mesh in the INP format of ABAQUS exported by CUBIT Mesh visualization MESH2D_PLOT plots a 2D mesh Miscellaneous tools SAMPLE_SEGMENTS generates points that regularly sample multiple segments of a line The functions MESH2D_TFI and MESH2D_MERGE are the core tools The script MESH2D_EXAMPLE1 is a good starting point The syntax of the mesh database file mesh2d is described in Section 4 4 3 4 Generating a mesh with EMC2 3 4 1 The mesh generator EMC2 EMC2 is one of the few public domain 2D mesh generation softwares that includes quadri lateral elements and a Graphical User Interface Its C code sources and executables can be freely downloaded from http www ann jussieu fr hecht ftp emc2 We show here an example featuring the most useful functionalities of EMC2 For further details you should refer to the complete documentation of the EMC2 package Before starting you must prepare files containing in 2 column data the coordinates X Z of all the points needed to define the geometry of the model topog
71. raphy sediment bottom 3 4 Generating a mesh with EMC2 22 Once installed you can run EMC2 by typing emc2 3 4 2 Notations The following notations are assumed in the next section e XXX click XXX on top menu bar xxx click xxx on bottom menu bar lt XXX lt gt XXX gt click XXX on left menu bar click XXX on right menu bar xxx enter xxx from keyboard or from the calculator in the right panel xxx type xxx in bottom prompt xxx perform action xxx xxx do xxx as many times as needed n xxx do xxx n times 3 4 3 Basic step by step A typical EMC2 session has three steps STEP I CONSTRUCT defines the geometry of the model 1 4 Switch to the construction tool lt CONSTRUCTION lt Load the points POINT xy file palosgrandes dat You must give the full path to your points file the root directory being the one where you launched emc2 Reset the figure window to fit all points gt SHOW ALL gt The original data has some geometrical features that are too complex to be meshed by quadrilaterals for instance the corners at the N and S ends of the basin you may want to smooth out these features You also need to define the extreme boundaries of the region to be modelled N S and bottom absorbing boundaries and some additional points on the free surface outside the basin You must modify the data set add and delete points Add new points 3 4 Generatin
72. ration software 25 5 Export the mesh lt SAVE lt Two questions are asked in the bottom prompt e Format of the file you must select ftg e Prefix name for the file name The resulting file name will be name ftq 3 4 4 Further tips Whenever possible it is better to mesh a domain with a structured mesh a deformed cartesian grid This can be done with QUADRANGULATE during the PREPARE step See our FAQ for further details To load an existing project in the construction tool or in the preparation mesh tool lt RESTORE lt name EMC2 will look for the file name emc2_bd Beware the project loaded will replace the actual project if any there is no superposition BUG WARNING 13 07 01 the Sun release of EMC2 has a bug with the reference indices in the ftq format This bug is fixed in the 2 12c version If you work on a Sun station download the most recent version of the sources rather than the executable and compile it yourself To densify h refinement an existing mesh use the script SEM2DPACK POST href csh It edits the emc2_bd file You can then restore it in EMC2 and save it in ftq format To create a fault in EDIT MESH mode a Crack an existing edge CRACK segment b Give a reference number to each side of the fault MODIF_REF n segment c Give the tag 1 to crack tip nodes MODIF_REF 1 corner click close to crack tip node inside element Note that only
73. scribed in Chapter 4 2 1 General assumptions and conventions The coordinate sytem is Cartesian rectangular SEM2DPACK works in the two dimensional x z plane where z is the horizontal coordinate with positive direction pointing to the right and z is the vertical coordinate with positive direction pointing upwards The coordinates x y z will be also denoted as 11 12 13 This notations carry also for subscripts For instance the k th component of displacement is denoted as uz with k 1 2 3 or with k 2X Y 2 The reference frame is Eulerian Infinitesimal strain is assumed The symmetric infinitesi mal strain tensor e is defined as 1 Ou Ou Se A aE seed 2 1 ij 2 E e Material density is deonted p x z The displacements and stresses relative to an initial equilibrium configuration are denoted uz xw z t and o4 x z t respectively External forces sources are denoted f x z t SEM2DPACK solves the following equations of motion to obtain the relative displacements uz 2 z t lu e 00 OB On tH 2 where summation over repeated indices is assumed The initial conditions are uz 0 and Ou Ot 0 Stresses are related to strain and possibly to other internal variables by constitutive equations described in Section 2 2 The governing equations are supplemented by boundary conditions described in Section 2 3 SEM2DPACK actually solves the governing 2 2 Material rheologies 12 equations in variational
74. sian 2D distribution SYNTAX amp DIST_GAUSSIAN centered_at length offset ampli order centered_at dble 2 0 0 Coordinates of the center point length dble 2 1 Characteristic lengths on each axis offset dble 0 Background level ampli dable 1 Amplitude from background order int 1 Exponent NAME DIST_GRADIENT GROUP DISTRIBUTIONS_2D PURPOSE Constant gradient 2D distribution SYNTAX amp DIST_GRADIENT file valref grad angle file name none Name of the file containing the coordinates of the points defining the reference line It is an ASCII file with 2 columns per line 1 X position in m and 4 4 Input Blocks Reference Guide 37 valref grad angle NOTE NAME GROUP PURPOSE SYNTAX file col NOTE NOTE NAME GROUP PURPOSE SYNTAX 2 Z position in m dble none Value along the reference line dble gt 0 none Positive gradient valref_units meter dble none Angle degrees between the vertical down and the grad direction Anticlockwise convention grad points down if 0 right if 90 Make sure the angle and ref line are compatible The code will abort if the ref line is too short some points of the domain cannot be projected to ref line in the angle direction DIST_HETE1 DISTRIBUTIONS_2D Linear interpolation of values from a regular 2D grid amp DIST_HETE1 file col name none Name of the file containing the definition of the regul
75. t name none Station positions can instead be read from an ASCII file with 2 columns X and Z in meters log T Relocate the stations at the nearest GLL node int 1 Sampling stride in number of timesteps Note that for stability reasons the timestep can be very small char V The field in the seismogram D displacement V velocity A acceleration char D Abscissa for the seismic multitrace plot X Horizontal position Z Depth 7D Distance to the first station to locate receivers at the free surface set their vertical position 4 4 Input Blocks Reference Guide 48 above the free surface and AtNode T NAME SRC_FORCE GROUP SOURCE MECHANISM PURPOSE Point force source SYNTAX amp SRC_FORCE angle angle dble 0d0 For P SV the angle of the applied force in degrees counterclockwise from Z UP e g 90 points left 180 points down For SH angle is ignored and the SRC_FORCE block is not required NAME SRC_DEF PURPOSE Define the sources SYNTAX amp SRC_DEF stf mechanism coord amp SRC_DEF stf mechanism file followed by one SOURCE TIME FUNCTION block STF_XXXX and one SOURCE MECHANISM block SRC_XXXX stf name none Name of the source time function gt RICKER TAB HARMONIC BRUNE or USER mechanism name none Name of the source mechanism FORCE EXPLOSION DOUBLE_COUPLE MOMENT or WAVE coord dble 2 huge Lo
76. taken over the whole mesh Red elements small SI are relatively unstable and require small timesteps At Because At is constant over the whole mesh and the computational cost is inversely proportional to At these red ele ments penalize the computational efficiency The mesh should be redesigned to increase their size as much as possible while keeping them small enough to resolve the shortest wavelength see next e Resolution criterion related to the number of nodes per shortest wavelength The resolution of each element is quantified by R min cs h We also define a resolution index as RI log R median R where the median value is taken over the whole mesh Red elements small RI have relatively poor resolution in their vicinity the maximum frequency resolvable by the mesh is limited The mesh should be redesigned to decrease their size as much as possible Conversely elements with very high RI blue are smaller than required and might increase the computational cost To minimize the CPU and memory cost of a simulation an ideal mesh design should minimize the spread of the two indices above by aiming at a ratio of element size to wave velocity h c as uniform as possible across the whole mesh However in some cases a poorly balanced mesh is inevitable in the example of Figure 4 3 the worst elements are near the edges of the sedimentary basin at a sharp velocity contrast Small element sizes on the rock side are inherited fr
77. that 1 you have read the documentation see Section 1 7 including the Frequently Asked Questions Chapter 6 Suggestions on how to improve the documentation are also welcome 1 9 Contributions 9 2 you are running the most recent version of SEM2DPACK Your issue might have been already fixed in a more recent version 3 you understand the changes listed in SEM2DPACK s ChangeLog file especially changes in the format of the input files 4 your problem has not been treated in previous submissions You can browse the tracker message titles or search for keywords A new submission must include the input files needed to reproduce your problem Par inp ftq mesh2d etc You will receive email notifications of any update of your submitted item until it is closed If the item is declared Pending you are expected to reply to the last message of the developer within two weeks otherwise the item will be closed For more instructions see http sourceforge net support getsupport php group_id 182742 1 9 Contributions Contributions to SEM2DPACK by experienced programmers are always welcome and en couraged Although the code is stable for typical applications in computational seismology and earthquake dynamics there is still a number of missing features Their implementa tion could make SEM2DPACK interesting for a broader audience in mechanical engineering geotechnical engineering applied geophysics and beyond The solver
78. the geometrical flexibility and the natural implementation of mixed boundary conditions Introductory texts to the SEM can be found at www math lsa umich edu karni m501 boyd pdf chapter draft by J P Boyd at www mate tue nl people vosse docs vosse96b pdf a tutorial exposition of the SEM and its connection to other methods by F N van de Vosse and P D Minev and at www siam org siamnews 01 04 spectral pdf a perspective paper Details about the elastodynamic algorithm and study of some of its properties are presented by Komatitsch 1997 Komatitsch and Vilotte 1998 Komatitsch et al 1999 Komatitsch and Tromp 1999 and Vai et al 1998 The implementation of fault dynamics is similar to that in FEM with the traction at split nodes method explained by Andrews 1999 More details can be found in the author s Ph D dissertation Ampuero 2002 in Gaetano Festa s Ph D dissertation and in Kaneko et al 2008 A more accesible tutorial code SBIEMLAB written in Matlab can be downloaded from the author s website at web gps caltech edu ampuero software html lyeb gps caltech edu ampuero publications html people na infn it festa 4 2 Basic usage flow 28 4 2 Basic usage flow In general a simulation requires the following steps p Prepare the input file Par inp Section 4 3 and Section 4 4 Run the solver in check mode iexec 0 in the GENERAL input block of Par inp sem2dsol
79. to be interpolated int 1 Interpolate along X dim 1 or along Z dim 2 name none Name of the ASCII file containing the data dble 0 Smoothing length for sliding average window No smoothing if length 0 DIST_ORDERO DISTRIBUTIONS_2D Blockwise constant 2D distribution amp DIST_ORDERO xn zn x 1 x xn 1 z 1 z zn 1 v 1 1 v xn 1 v 1 zn v xn zn int none Number of zones along X int none Number of zones along Z dble xn 1 none Boundaries of X zones first zone X lt x 1 second zone x 1 lt X lt x 2 last zone x xn 1 lt X dble zn 1 mone Boundaries of Z zones dble xn zn none Values inside each zone DIST_PWCONR DISTRIBUTIONS_2D Piecewise constant radial 2D distribution This distribution defines a set of annular zones centered at an arbitrary reference point and assigns constant values within each zone amp DIST_PWCONR num ref FILI r num 1 v 1 vQ v num 1 v num int none Number of annular zones including inner and exterior dble 2 0d0 0d0 Reference point center of radial zones dble num 1 none External radius of zones 4 4 Input Blocks Reference Guide 39 first zone R lt r 1 second r 1 lt R lt r 2 last r num 1 lt R v dble num none Value inside each zone NAME DIST_SPLINE GROUP DISTRIBUTIONS_1D PURPOSE Spline interpolated 1D distribution along X or Z SYNTAX amp
80. ve gt info amp Verify the resolution stability estimated CPU cost and memory cost Section 4 5 If needed go back to step 1 and modify Par inp Section 4 5 else proceed to next step Run the solver in production mode iexec 1 sem2dsolve Plot and manipulate the solver results Section 4 6 N Doe Ww Full details are given in the following sections 4 3 General format of the input file The input file must be called Par inp Its typical structure is illustrated by two examples in Figure 4 1 and Figure 4 2 Most of the file is made of standard FORTRAN 90 NAMELIST in put blocks Each block gives input for a specific aspect of the simulation material properties sources receivers boundary conditions etc The general syntax of a NAMELIST block can be found in any FORTRAN 90 textbook In summary a block named STUFF with possible input arguments a b and c must be given as amp STUFF a b c where are user input values Line breaks and comments preceded by are allowed within an input block The complete Reference Guide of the input blocks is presented in Section 4 4 For each block the documentation includes its name possibly the name of a group of blocks to which it belongs its purpose its syntax the list of its arguments with their description and some important notes In the syntax description a vertical bar between two arguments means one or the other In the argument l
81. vsky et al 1997 eq 15 Y p VP 2109 213 2 3 Boundary conditions 15 where q 2u0 2 0 2 amp 2 14 p amp q o 2 15 The evolution equation for the damage variable is Lyakhovsky et al 1997 eq 20 _ J Cable o if gt EN 0 otherwise 220 No healing is assumed below o The evolution of the plastic strain e j 1S driven by the damage variable a Hamiel et al 2004 eq 9 e Ti C 4 if amp gt 0 ij 0 otherwise 2 17 1 where Tij oi 30h 6 is the deviatoric part of the stress tensor The parameter C is of order 1 19 and is related to the seismic coupling coefficient 0 lt x lt 1 by Ben Zion and Lyakhovsky 2006 Cy 2 18 2 3 Boundary conditions 2 3 1 Absorbing boundaries Two approximate absorbing boundary conditions ABC are implemented to model the out wards propagation of waves at the boundaries of the computational domain Both conditions are of paraxial type Their performance is appropriate at normal incidence but degrades at grazing incidence Clayton Engquist ABC In the local coordinate frame t n related to the tangential t and outgoing normal n directions to the boundary the first order accurate ABC proposed by Clayton and Engquist 1977 reads Ou in cg ai 2 19 Oun i ora 2 20 The implementation is based on an equivalent formulation as a mixed boundary condition that relates tractions T to displacements u
82. xtensions tab for ASCII data files txt for other text files dat for binary data files etc This makes it easy to clean a working directory with a single command like rm f _sem2d 4 6 1 Spectral element grid As explained in the previous section sem2dsolve generates two PostScript files for mesh quality checking purposes Stability_sem2d ps and Resolution_sem2d ps The relevant information is contained in the files Stability_sem2d tab and Resolution_sem2d tab and can also be inspected with the Matlab script PRE ViewMeshQuality m 4 6 2 Source time function sem2dsolve generates a file called SourcesTime_sem2d tab containing the source time func tion sampled at the same rate as the receivers It is important to verify that the spectrum of the source has little power at those high frequencies that are not well resolved by the mesh those that correspond to less than 5 nodes per wavelength If this is not the case you must be very cautious in the interpretation of the seismograms in the high frequency range or low pass filter the results 4 6 3 Snapshots sem2dsolve generates snapshots at a constant interval defined in number of solver timesteps by the input parameter itd of the SNAP_DEF input block An example is shown in Fig ure 4 4 Requested fields are exported in binary data files called xx_XXX_sem2d dat where xx is the field code defined in the documentation of the PLOTS input block and XXX is the 3 digit snapshot number The user
83. zx Mzz amp SRC_MOMENT Myx Myz Mzx Mzz dble 0 Tensor components for PSV dble 0 Tensor components for SH NAME GROUP PURPOSE SYNTAX angle phase NOTE SRC_WAVE SOURCE MECHANISM Incident plane wave through the absorbing boundaries amp SRC_WAVE angle phase dable 0d0 Incidence angle in degrees within 180 180 counterclockwise from the positive Z up direction to the wave vector direction Exs incidence from below if angle in 90 90 normal incidence from below if angle 0 from bottom right if angle 45 from bottom left if angle 45 char S S or P only needed in PSV ignored in SH Incident waves enter through the absorbing boundaries An incident wave is applied on every absorbing boundary unless let_wave F in the respective BC_ABSO block Incident waves are not implemented for Stacey absorbing boundaries 4 4 Input Blocks Reference Guide 50 NAME GROUP SYNTAX NAME GROUP ampli fO NAME GROUP ampli fO onset NOTE NOTE NAME GROUP STF_BRUNE SOURCE TIME FUNCTIONS PURPOSE Brune 1970 s model with omega squared spectral fall off stf t amplix 1 1 2 pi fc t exp 2 pi fc t amp STF_BRUNE ampli fc dble 1d0 Amplitude usually the seismic moment dble 1d0 Corner frequency Hz STF_HARMONIC SOURCE TIME FUNCTIONS PURPOSE SYNTAX Harmonic source time function f t ampli sin 2 pi t f0 amp STF_HARMONIC

User`s guide - California Institute of Technology

Contents

Download Pdf Manuals

Related Search

Related Contents