Swpc Guide E
User Manual: Pdf
Open the PDF directly: View PDF .
Page Count: 52
Download | |
Open PDF In Browser | View PDF |
OpenSWPC: An Open-source Seismic Wave Propagation Code User’s Guide Takuto Maeda Version 4.0-E (2017-09-21) Contents 1 Set Up 1.1 System Requirements . . . . . . . . . . 1.2 Code directory tree . . . . . . . . . . . 1.3 Compilation and Execution . . . . . . 1.3.1 make . . . . . . . . . . . . . . . 1.3.2 Specifying Compiler Options . 1.3.3 More about the NetCDF library 1.3.4 Preparing the Dataset . . . . . 1.3.5 On Embedding Parameters . . 1.3.6 Execution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 4 4 4 5 5 6 6 2 Parameter Settings 2.1 Notation of the Parameter File . . . . . . . . . 2.2 An Example Parameter File . . . . . . . . . . 2.3 Controlling Parameters . . . . . . . . . . . . . 2.4 Coordinate System and Parallel Computation 2.4.1 Staggered Grid . . . . . . . . . . . . . 2.5 Viscoelastic Bodies . . . . . . . . . . . . . . . 2.6 Simulation Data Output . . . . . . . . . . . . 2.6.1 Output Datafile Format . . . . . . . . 2.6.2 Snapshot Data Output . . . . . . . . . 2.6.3 Seismic Waveform Output . . . . . . . 2.6.4 Output Filename Conventions . . . . 2.7 Velocity Model . . . . . . . . . . . . . . . . . . 2.7.1 Choice of Velocity Model Type . . . . 2.7.2 Small-Scale Random Inhomogeneity . 2.8 Earthquake Source Specification . . . . . . . 2.8.1 Moment Rate Function . . . . . . . . . 2.8.2 Moment Tensor Source . . . . . . . . . 2.8.3 Body Force Mode . . . . . . . . . . . . 2.8.4 Plane Wave Mode . . . . . . . . . . . . 2.9 Absorbing Boundary Conditions . . . . . . . 2.10 Checkpointing and Restarting . . . . . . . . . 2.11 Reciprocity Mode . . . . . . . . . . . . . . . . 2.12 About Two-Dimensional Codes . . . . . . . . 2.13 Other Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 8 12 12 13 15 18 18 19 20 22 23 23 25 27 27 28 30 31 32 34 35 37 37 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Related Tools 3.1 Snapshot data handling . . . . . . . . . . . 3.1.1 read snp.x . . . . . . . . . . . . . . 3.1.2 diff snp.x . . . . . . . . . . . . . . 3.2 Supporting Parameter Settings . . . . . . . 3.2.1 fdmcond.x . . . . . . . . . . . . . . . 3.2.2 mapregion.x . . . . . . . . . . . . . 3.2.3 mapregion.gmt4, mapregion.gmt5 3.3 Velocity Structure . . . . . . . . . . . . . . . 3.3.1 qmodel tau.x . . . . . . . . . . . . . 3.3.2 grdsnp.x . . . . . . . . . . . . . . . . 3.4 Generation of Random Media . . . . . . . . 3.4.1 gen rmed3d.x . . . . . . . . . . . . . 3.4.2 gen rmed2d.x . . . . . . . . . . . . . 3.5 Miscellaneous Tools . . . . . . . . . . . . . . 3.5.1 timvis scripts . . . . . . . . . . . . . 3.5.2 Geographic Coordinate Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 39 39 40 40 40 41 42 42 42 42 42 42 44 44 44 44 4 Additional Materials 4.1 Hints for Parameter Settings . . . . . . . . . . . . 4.2 Hints for Modifying the Code . . . . . . . . . . . 4.2.1 Defining Your own Velocity Model . . . . 4.2.2 Defining Your own Source Time Function 4.2.3 Appending New Control Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 45 45 46 46 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments 48 Revision History 49 Copyright & License 50 2 Chapter 1 Set Up 1.1 System Requirements Executing OpenSWPC requires a Fortran compiler that can handle (at least a part of) the Fortran 2003 standard and an MPI library. The program can be run on a single CPU or CPU core without parallelization; however, the MPI library is still required. In addition, the NetCDF library, compiled by the same Fortran compiler, is recommended to use the direct input/output of the NetCDF-formatted files. The source code of SEISM almost strictly follows the language standard of Fortran2003. As an exception, system calls (the system() subroutine) are used. Note that this extension is supported by most available Fortran compilers. OpenSWPC uses stream I/O, which is part of the Fortran2003 standards. This functionality may not be implemented with older compilers. This code was developed in the following environment: • Apple OSX El Capitan • GNU gfortran 6.1.0 • OpenMPI 1.10.2 In addition, the following computers were confirmed to work with OpenSWPC: • EIC computer, ERI/UTokyo (ver. 2015; SGI Altix; intel fortran) • JAMSTEC Earth Simulator (NEC SX-ACE; NEC compiler) • AICS K computer (Fujitsu compiler) • Nagoya University (FX10/Fx100; Fujitsu compiler) • Linux Cent OS 6.6 (gfortran 4.9.2 & mpich) • Linux Ubuntu 16.04LTS (gfortran 5.4 & OpenMPI) 1.2 |-|-|-\-- Code directory tree doc : manuals bin : executable binaries (*.x) example : example input files src 3 arch name mac-intel mac-gfortran eic fx es3 ubuntu-gfortran |-|-|-|-\-- 1.3 1.3.1 swpc_3d swpc_psv swpc_sh tools shared Table 1.1: arch options for various environments. target NetCDF location Mac OSX + Intel Compiler + Open MPI ${HOME}/local Mac OSX + gfortran + Open MPI /usr/local EIC (ERI, UTokyo) with the Intel Compiler ${HOME}/local Fujitsu FX10, FX100 and the K-computer ${HOME}/xlocal The Earth Simulator Provided by the system Ubuntu 16.04LTS + gfortran + Open MPI Installation by apt : : : : : 3D problem 2D P-SV problem 2D SH problem Miscellaneous utility codes Modules commonly used in the above programs Compilation and Execution make The directories src/swpc 3d, src/swpc psv, src/swpc sh, and src/tools contain makefile. Execute the make command in each directory to generate the executable binaries. An executable file (with a *.x extension) will be stored in the bin directory. 1.3.2 Specifying Compiler Options In the makefiles, the following variables must be specified according to the environment: FC compiler name FFLAGS compiler option NCFLAG NetCDF flag NCLIB location of the NetCDF library directory NCINC location of the NetCDF include directory NETCDF linker option for NetCDF If NCFLAG = -D NETCDF is specified, make will try to compile OpenSWPC with NetCDF. A set of the above variables under different computer environments is defined in src/shared/makefile.arch and src/shared/makefile-tools.arch. The former is for the compilation of FDM codes, and the latter is for the compilation of misc tools. The user can specify the arch option in make as shown below: 1 make arch=eic debug=true The list of pre-defined architecture (arch) options is described in Table 1.1. 4 Figure 1.1: Area of JIVSM and eJIVSM. The colored area in the map is where the original JIVSM is defined. eJIVSM is extended to the gray-shaded area via an extrapolation. The surrounding graphs show the depth sections along the lines on the map of the model. 1.3.3 More about the NetCDF library The NetCDF library consists of the following items: libnetcdf.* NetCDF library file (static) libnetcdff.* NetCDF Fortran library file (only NetCDF version 4 or later) netcdf.mod Fortran module information file The extension of the library files may be *.a (static library) or *.so (dynamic library), depending on the installation. All these files are necessary for successful compilation with NetCDF. In particular, the netcdf.mod file must be created by the same Fortran compiler as OpenSWPC. If NetCDF is installed using packaging tools such as yum, apt, or homebrew, the use of gfortran is implicitly assumed. 1.3.4 Preparing the Dataset Subsurface Velocity Structure Model In OpenSWPC, a 3D inhomogeneous medium can be represented as a set of velocity discontinuities using NetCDF-formatted files (see Section 2.7 for details). As an example of the velocity structure beneath the Japanese Archipelago, an automatic model generation script for the Japan Integrated Velocity Structure Model (JIVSM; Figure 1.1; Koketsu et al. (2012)), developed and originally distributed by the Headquarters for Earthquake Research Promotion in Japan, is included. An extension of JIVSM (eJIVSM), which covers a wider area, is also provided. These velocity structure models contain the ground surface (topography and bathymetry), subsurface soil, Moho, and oceanic crust of the two subducting plates. To generate these models, GMT version 4 is required. If the user does not use this model, the following processing steps may not be necessary. First, download the original model files lp2012nankai-e str.zip and lp2012nankai-w str.zip and store them in dataset/vmodel. The URLs of these files can be found in the comments of the gen JIVSM.sh script. To generate eJIVSM, the ETOPO1 (ETOPO1 Bed g gmt4.grd) topography data are also necessary. Then, specify the Fortran compiler name to FC variables in the params.sh parameter file. The grid spacing (dlon, dlat) in the parameter file can be modified if necessary. Note that this spacing is not directly related to the grid width of the numerical simulations. At the time of the simulation, the OpenSWPC program automatically interpolates the velocity model data. After these steps, execute the generation script: 1 ./ gen_JIVSM .sh After a successful execution, 23 NetCDF-formatted files will be generated in the two model directories for JIVSM and eJIVSM. These files can be read and visualized in GMT by the grdimage or grd2xyz modules. The netcdf filename contains five integer numbers, 5 which correspond to mass density (in kg/m3 ), P wavespeed (m/s), S wavespeed (m/s), QP , and QS . They indicate the material information below the discontinuity defined in the file. List files of these NetCDF files (jivsm.lst and ejivsm.lst) for use in OpenSWPC will also be generated. Station Lists An example script to generate a station list file is stored in dataset/station/gen stlst hinet.sh. This script generates a formatted list of the high-sensitivity seismograph network Japan (Hinet) provided by the National Research Institute for Earth Science and Disaster Resilience (NIED). To use this script, first download the station csv list from the Hi-net website following the comments in the gen stlst hinet.sh file. Then, executing this bash script will result in the station list file for OpenSWPC. 1.3.5 On Embedding Parameters Although most of the behavior of OpenSWPC is controlled dynamically by the input parameter file, several parameters are embedded in the source code to achieve high-computational performance, as described below. These parameters are defined in src/swpc ??/m global.F90. If these parameters are modified, re-compilation is necessary. UC = 1E-15 (1E-12 for 2D codes) A number to convert the simulation results into the SI unit system. Modifications may be necessary to use a different unit system. MP = DP Precision of the finite-difference computation. By default (MP=DP), i.e., parts of the computation are performed in double precision, while other unnecessary variables are defined and calculated in single precision to save memory space and computation time. The user may change it to MP=SP to switch the entire computation to single precision, which will decrease the required memory up to 2/3 and allow faster computation speeds. However, in this case, a noisy seismic waveform might be observed, in particular, near the seismic source due to the overflow of floating point numbers. NM = 3 Number of generalized Zener viscoelastic bodies. If this number is larger than 1, it represents a nearly frequency-independent constant Q model in a specified frequency range. If this is set to zero, the simulation will be conducted with an elastic body without attenuation. Increasing this number enables the reproduction of a wider frequency range of constant Q; however, it may also result in a significant increase in the computational loads for 3D simulations. 1.3.6 Execution To run the program, the MPI program is necessary, such that 1 > mpirun -np ${NP} ./ bin/ swpc_3d .x -i ${ input} where ${NP} is the number of MPI processes and ${input} is the name of the input file. Note that the mpirun command may be different for different computational systems. 6 If the program runs properly, the following message will appear in the standard error output. The result may be slightly different for different programs (3D/P-SV/SH) or execution modes. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 -----------------------------------------------------------------------------SWPC_3D ( benchmark mode) -----------------------------------------------------------------------------Grid Size MPI Partitioning Total Memory Size Node Memory Size Stability Condition c Wavelength Condition r Minimum velocity Maximum velocity Maximum frequency : : : : : : : : : 384 x 384 x 384 4 x 6 12.705 [GiB] 0.529 [GiB] 0.980 (c <1) 7.000 (r >5 -10) 3.500 [km/s] 6.062 [km/s] 1.000 [Hz] -----------------------------------------------------------------------------it =0000050 , it =0000100 , it =0000150 , it =0000200 , 1.877 1.887 1.932 1.943 s/loop , s/loop , s/loop , s/loop , eta eta eta eta 000:29:43 , 000:28:18 , 000:27:22 , 000:25:54 , ( ( ( ( 5.00E -05 1.75E -05 1.02E -05 6.59E -06 5.00E -05 1.75E -05 1.02E -05 6.59E -06 4.96E -05 1.05E -05 5.41E -06 4.35E -06 ) ) ) ) it =0000950 , 1.986 s/loop , eta 000:01:39 , ( 4.89E -07 it =0001000 , 1.982 s/loop , eta 000:00:00 , ( 1.65E -07 4.89E -07 1.65E -07 1.81E -06 ) 1.54E -07 ) . . . -----------------------------------------------------------------------------Total time : 1982.348 s ------------------------------------------------------------------------------ The first part of the message contains information such as the estimated memory usage, stability condition, and wavelength condition. As shown in the above example, a stability condition of c < 1 is mandatory to execute; if the specified parameter violates this condition, the program aborts immediately. In addition, the wavelength condition (the ratio of spatial grid size and minimum wavelength) is recommended to r > 5−10. During the computation, the computation speed, remaining time (eta; estimated time of arrival), and maximum velocity amplitude of the components are shown. 7 Chapter 2 Parameter Settings 2.1 Notation of the Parameter File In the parameter files, one parameter is defined on each line in the following format. 1 variable_name = value The description of the values should follow Fortran notation. For example, logical (Boolean) values are denoted as .true. or .false.. Lines that do not contain an equal sign (=) will be neglected; in addition, lines starting with ! or # are regarded as comment lines and will be skipped. Comments can be followed by variable definitions, that is, comments can be written on the same line as the parameter definition. For example, the following parameter line will work without errors. 1 nx = 1024 ! number of grids The order of the parameter definition can be changed freely. If a parameter is not specified, OpenSWPC may use a pre-defined default variable in some cases. In such a case, the use of the default-parameter will be included in the standard error output. Note that there are parameters that must be defined explicitly. Multiple definitions of the same parameters in a single parameter file are not recommended; however, in this case, the first definition may be adopted. It is acceptable to leave blanks before and after the equal sign; however, it is not permitted to have a blank space between the minus character and succeeding numbers (e.g., ’- 35.0’ is not allowed). It is recommended to use quotation marks around string parameters. Without them, the directory path character (’/’) may be unexpectedly interpreted as a termination of the string parameter. 2.2 An Example Parameter File The following is a full set of example parameters. In the following sections, detailed descriptions of each parameter will be given. 1 2 3 4 5 6 !! ----------------------------------------------------------------------- !! !! !! SWPC input file !! !! ----------------------------------------------------------------------- !! 8 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 !! ----------------------------------------------------------------------- !! !! Control !! title odir ntdec_r = ’swpc ’ = ’./ out ’ = 50 !! exe title: used for output filenames !! output directory !! screen report timing (1/ cycle) !! ----------------------------------------------------------------------- !! !! Model /Grid Size and Area !! nproc_x nproc_y nx ny nz nt = = = = = = 4 6 384 384 384 1000 !! !! !! !! !! !! parallelization in x-dir parallelization in x-dir total grid number in x-dir total grid number in y-dir total grid number in z-dir time step number dx dy dz dt = = = = 0.5 0.5 0.5 0.02 !! !! !! !! grid grid grid time vcut = 1.5 !! minimum velocity !- smaller velocities will be increased xbeg ybeg zbeg tbeg = = = = -96.0 -96.0 -10.0 0.0 !! !! !! !! minimum in x-dir minimum in y-dir minimum in z-dir start time clon clat phi = 139.7604 = 35.7182 = 0.0 !! !! !! !- center longitude center latitude horizontal coordinate rotation measured clockwise from the north fq_min fq_max fq_ref = 0.02 = 2.00 = 1.0 !! minimum freq. for Q- constant model !! maximum freq. for Q- constant model !! ref. freq. for physical dispersion width in x-dir width in y-dir width in z-dir step width !! ----------------------------------------------------------------------- !! !! Snapshot Output !! snp_format = ’netcdf ’ !! snapshot format ( native or netcdf ) xy_ps %sw xz_ps %sw yz_ps %sw fs_ps %sw ob_ps %sw = = = = = . false. .true. . false. . false. .true. !! !! !! !! !! P&S P&S P&S P&S P&S xy_v%sw xz_v%sw yz_v%sw fs_v%sw ob_v%sw = = = = = . false. .true. . false. . false. .true. !! !! !! !! !! 3-comp. 3-comp. 3-comp. 3-comp. 3-comp. 9 amp. amp. amp. amp. amp. for for for for for xy section xz section yz section free surface ocean bottom velocity velocity velocity velocity velocity for for for for for xy section xz section yz section free surface ocean bottom 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 xy_u%sw xz_u%sw yz_u%sw fs_u%sw ob_u%sw = = = = = . false. .true. . false. . false. .true. z0_xy x0_yz y0_xz = = = ntdec_s = 5 !! time decimation of snapshot !- ( specify 1 for no decimation ) idec jdec kdec = 2 = 2 = 2 !! x- decimation for snapshot !! y- decimation for snapshot !! z- decimation for snapshot 7.0 0.0 0.0 !! !! !! !! !! 3-comp. 3-comp. 3-comp. 3-comp. 3-comp. disp. disp. disp. disp. disp. for for for for for xy section xz section yz section free surface ocean bottom !! depth for xy cross section !! x- value for yz cross section !! y- value for xz cross section !! ----------------------------------------------------------------------- !! !! Waveform Output !! sw_wav_v sw_wav_u ntdec_w st_format fn_stloc wav_format = = = = = = .true. !! velocity trace output at stations . false. !! displacement trace output at stations 5 !! time decimation of waveform output ’xy ’ !! station format : ’xy ’ or ’ll ’ ’./ example /stloc .xy ’ !! station location file ’sac ’ !! ’sac ’ or ’csf ’ !! ----------------------------------------------------------------------- !! !! Earthquake Source !! !! Moment tensor source format : !! xymoij / xym0dc / llm0ij / llm0dc / xymwij / xymwdc / llmwij / llmwdc !! Body force source format : !! xy or ll stf_format = ’xym0ij ’ !! Basis source time function !! ’boxcar ’ / ’triangle ’ / ’herrmann ’ / ’kupper ’ / ’cosine ’ / ’texp ’ stftype = ’kupper ’ fn_stf = "./ example / source .dat" !! Source grid file name !! source depth correction !! ’asis ’: use z value , ’bd{i}’: i-th boundary (i =0...9) sdep_fit = ’asis ’ !! --------------------------------------------------------------------- !! !! Body force source mode !! bf_mode = . false . !! --------------------------------------------------------------------- !! !! Plane wave source mode !! pw_mode = . false . !! plane wave input; neglects fn_stf pw_ztop = 100. !! top z- coordinate of the initial plane wave pw_zlen = 30. !! wavelength of the initial plane wave pw_ps = ’p’ !! ’p’ P-wave ’s’ S-wave 10 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 pw_strike pw_dip pw_rake = 0.0 = 0.0 = 0.0 !! strike direction of plane wave (deg .) !! dip of plane wave (deg .) !! rake of plane S-wave polarization (deg .) !! ----------------------------------------------------------------------- !! !! Absorbing Boundary Condition !! abc_type na stabilize_pml = ’pml ’ = 20 = .true. !! ’pml ’ or ’cerjan ’ !! absorbing layer thickness !! avoid low -v layer in PML region !! ----------------------------------------------------------------------- !! !! Velocity model !! vmodel_type is_ocean is_flatten = ’lhm ’ = .true. = . false. !! velocity model type ’uni ’/’grd ’/’lhm ’ !! topography z<0 is covered by ocean !! Force topography variation to zero !! --------------------------------------------------------------------- !! !! For uniform velocity model ’uni ’ !! vp0 = 5.0 !! P-wave velocity [km/s] vs0 = 3.0 !! S-wave velocity [km/s] rho0 = 2.7 !! mass density [g/cm ˆ3] qp0 = 200 !! Qp qs0 = 200 !! Qs topo0 = 0 !! topography location !! --------------------------------------------------------------------- !! !! For GMT grid file input ’grd ’ ( requires netcdf library ) !! dir_grd = ’${ DATASET }/ vmodel / ejivsm ’ !! directory for grd file fn_grdlst = ’./ example /grd.lst ’ !! grd file list node_grd = 0 !! input MPI node !! --------------------------------------------------------------------- !! !! For layered homogeneous medium model (’lhm ’) !! fn_lhm = ’example /lhm.dat ’ !! 1D velocity structure !! --------------------------------------------------------------------- !! !! For random medium models !! dir_rmed = ’./ in/’ !! location of random medium file fn_grdlst_rmed = ’./ example /grd.lst ’ !! grd file list rhomin = 1.0 !! minimum density threshold !! ----------------------------------------------------------------------- !! !! Checkpoint / Restart !! is_ckp = . false. !! perform checkpoint / restart ckpdir = ’./ out/ckp ’ !! output directory ckp_interval = 1000000 !! interval for checkpoint (1/ cycle ) ckp_time = 1000000. !! checkpoint time = .true. !! sequential output mode ckp_seq !! ----------------------------------------------------------------------- !! !! Reciprocity Green ’s Function Mode !! green_mode = . false. !! reciprocity Green ’s function mode 11 Figure 2.1: (a) Partitioning of the computational domain for MPI. (b) Schematic of the data exchange by the MPI protocol (modified from Maeda et al., 2013). 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 green_stnm green_cmp green_trise green_bforce green_maxdist green_fmt fn_glst = = = = = = = ’st01 ’ !! ’z’ !! 1.0 !! . false. !! 550. !! ’llz ’ !! ’example /green.lst ’ virtual station name from fn_stlst virtual source direction ’x’, ’y’, ’z’ rise time also calc. body force Green ’s function horizontal limit of source grid list file format : ’xyz ’ or ’llz ’ !! Green ’s function grid point list !! ----------------------------------------------------------------------- !! !! MISC !! stopwatch_mode benchmark_mode = .true. = .true. !! measure computation time at routines !! benchmark mode ipad jpad kpad = 0 = 0 = 0 !! memory padding size for tuning !! memory padding size for tuning !! memory padding size for tuning 2.3 Controlling Parameters title Title of the computation to be used for the output filename. odir Name of output directory. This is a relative directory path from the location of the program execution. If this directory does not exist at the time of run, OpenSWPC will automatically create it. ntdec r Number of Time-step DECimation factors for screen Reporting. The maximum amplitudes of the velocity components are reported in the standard error output every ntdec r steps. This screen output is generally used to confirm that the model is working correctly. A cycle that is too short (this parameter is too small) may slow down the computation. 2.4 Coordinate System and Parallel Computation For parallel computation, OpenSWPC performs 2D model partitioning for a 3D code (figure 2.1) and 1D partitioning for a 2D code, in the horizontal direction in both cases. The computation is performed in Cartesian coordinates. We adopt the computational coordinate system depicted in Figure 2.2. By default, the coordinate axes x, y, and z represent the north, east, and depth directions, respectively. They cover the region of xbeg ≤ x ≤ xend, ybeg ≤ y ≤ yend, and zbeg ≤ z ≤ zend. Note that the z-axis is defined as positive downward. Because the free surface is usually defined at z = 0, it is recommended to let zbeg be a negative value to include the free surface in the model. 12 Figure 2.2: Relation between computational coordinate and geographical coordinate systems. The volume is discretized into nx, ny, and nz grids with spatial grid widths of dx, dy, and dz, respectively, in each direction. The parameter file must provide definitions of xbeg, ybeg, and zbeg and nx, ny, and nz; other parameters (xend, yend, and zend) are automatically computed from them. The center of the Cartesian coordinate (x = 0, y = 0) corresponds to the center longitude (clon) and latitude (clat). The geographical coordinate is projected onto the Cartesian coordinate by the Gauss–Krüger transform as follows (see Figure 2.2): 1. First generate an evenly spaced grid in Cartesian coordinates from the input parameters phi and those related to the x, y coordinates. 2. Project the grid location onto the geographical coordinate by using the Gauss–Krüger transform with a center location of (clon, clat). 3. Obtain the medium parameter at the grid location via a bicubic interpolation of the input velocity structure model. If the specified area exceeds that of the input velocity model, the outermost value of the velocity structure is used for the extrapolation. 2.4.1 Staggered Grid OpenSWPC adopts the staggered-grid coordinate system shown in Figure 2.3. The unit volume shown in the figure is defined as a “voxel” at the grid indices (I,J,K). A grid location x belongs to the voxel number ⌈x − x ⌉ beg I= (2.1) , ∆x and if the voxel number I is given, its center coordinate location is ( ) 1 x = xbeg + I − ∆x, 2 (2.2) where ⌈·⌉ is a ceiling function and xbeg is the minimum value of the x-coordinate. Note that xbeg is set to belong to the voxel I=1. A voxel has a volume of xbeg + (I − 1)∆x < x ≤ xbeg + I∆x, ybeg + (J − 1)∆y < y ≤ ybeg + J∆y, (2.3) zbeg + (K − 1)∆z < z ≤ zbeg + K∆z, The normal stress tensor components are defined at the center of the voxel, the shear stress is defined on the edge, and velocity vector components are defined on its surface (Figure 2.3) . Medium parameters are defined at the center of the voxel at xbeg + (I − 1/2)∆x, ybeg + (J − 1/2)∆y, zbeg + (K − 1/2)∆z. 13 (2.4) Voxel (I,J,K) vz K∆z σyz σxz z σii vx vy σxy (K-1)∆z (J-1)∆y x y J∆y I∆x (I-1)∆x z x y Figure 2.3: Staggered grid layout in 3D space for the case of xbeg=ybeg=zbeg=0. If necessary, averaging will be performed between neighboring voxels. The spatial grid width, ∆x, ∆y, and ∆z, and the time step width, ∆t, must satisfy the stability condition. The stability condition in ND -dimensional space for the order of the finite difference method P is given by ∆t < 1 Vmax −1/2 P/2 −1 N D ∑ ∑ 1 Cp , 2 ∆x i=1 i (2.5) p=1 where Vmax is the maximum velocity of the medium, Cp are the coefficients of the finite difference formula, and ∆xi is the spatial grid width in the i-th direction. For the fourthorder formula of the finite difference method, which is used in the code, the coefficients are C1 = 9/8 and C2 = 1/24. For example, for the fourth-order finite difference with isotropic grid sizes (∆x = ∆y = ∆z = h) in three-dimensional space, the stability condition is reduced to ∆t < 6 7 √ Vmax 1 1 1 1 + + 2 2 2 ∆x ∆y ∆z = h 6h ≃ 0.495 . √ V max 7 3Vmax (2.6) This condition can be interpreted as “the distance that the seismic wave propagates within a single time step must be much smaller than the spatial grid width.” The numerical simulation will diverge immediately if this condition is not satisfied. In addition, the minimum wavelength of the simulated seismic waves should be much longer than the spatial grid width. If the wavelength becomes relatively small compared to this condition, a fictitious numerical dispersion will appear and result into inaccurate later phases. Usually, the wavelength is taken to be longer than 5–10 times the spatial grid width to avoid this effect. Therefore, the minimum velocity (usually the S-wave velocity) 14 in the velocity model should be selected carefully. One may specify a smaller spatial grid size to avoid this problem; however, in this case, the time-step size must also be shortened to satisfy the stability condition. nproc x, nproc y Number of partitions in the x- and y-directions (Figure 2.1). The total number of partitions will be nproc x × nproc y for the 3D case and nproc x for the 2D case. This total number of partitions must be equal to the number of processes given in mpirun. These numbers can be 1. If nproc x=nproc y=1, this will become a serial (non-parallel) computation in practice. nx, ny, nz Total number of spatial grids in each direction. nx and ny do not need to be multiples of nproc x and nproc y, respectively. dx, dy, dz Spatial grid width in each direction in units of km. The total computational size in the physical domain will be nx×dx, ny×dy, and nz×dz. The grid widths in different directions do not necessarily need to be equal. nt Number of time steps. dt Length of the time step in seconds. The total (physical) simulation time will be nt×dt. xbeg, ybeg, zbeg Minimum value of the coordinates. If specifications of xbeg or ybeg are omitted, they will automatically be set to xbeg = - nx × dx / 2 and ybeg = - ny × dy / 2. This setting is recommended to minimize distortion due to the map projection. The default value of zbeg is -30×dz. tbeg Starting time. Usually, it is set to zero. clon, clat Center longitude and latitude in degrees. The map projection will be performed with this location as a reference point. phi Horizontal rotation angle of the computational coordinate (see Figure 2.2). If phi=0, the x- and y-axes correspond to the north and east directions, respectively. Note that the output files (snapshot and waveform) will be rotated if this value is nonzero. 2.5 Viscoelastic Bodies OpenSWPC adopts the generalized Zener body (GZB) as a model of the viscoelastic body. It consists of several viscoelastic Zener bodies with different relaxation times to attain nearly constant Q in a wide frequency range. As a consequence, it accompanies the frequencydependent body wavespeed via physical dispersion (e.g., Aki and Richards, 2002); therefore, users should specify the reference frequency in which the velocity model is given. 15 ( τε1 τσ1 ) ( τε2 τσ2 ) ( τε3 τσ3 ) ( τεL τσL ) Figure 2.4: GZB. GZB consists of NM Zener bodies, as schematically shown in Figure 2.4. This viscoelastic body is described by the relaxation functions for an elastic moduli π ≡ λ + µ and µ, as ) NM ( εP ∑ τ 1 σ ψπ (t) = πR 1 − 1 − mσ e−t/τm H(t) NM τm m=1 (2.7) ( ) N M εS ∑ τ 1 σ ψµ (t) = µR 1 − 1 − mσ e−t/τm H(t), NM τm m=1 where τσm is the relaxation time of the m-th body, πR ≡ λR +2µR and µR are the relaxed moduli, εS and τεP m and τm are creep times of the P- and S-waves, respectively. The wide frequency range of constant Q is achieved by connecting Zener bodies with different relaxation times. In addition, the intrinsic attenuations of the P- and S-waves (QP and QS , respectively) can be defined independently by choosing different creep times between the elastic moduli π and µ. The constitutive equation between stress and strain (or particle velocity) is written as follows. ( ) σ̇ii (t) = ψ̇π (t) − 2ψ̇µ (t) ∗ ∂k vk (t) + 2ψ̇µ (t) ∗ ∂i vi (t) (Do not take the sum over i.) ( ) (2.8) σ̇ij (t) = ψ̇µ (t) ∗ ∂i v j (t) + ∂ j vi (t) . The convolution appearing in the constitutive equation can be avoided by defining the memory variables (Robertsson et al., 1994) and solving the auxiliary differential equations for them. We also adopt the τ-method of Blanch et al. (1994) to automatically choose the creep times that achieve a constant Q condition. 16 fq_min fq_max 0.02 0.01 Q−1 0.005 0.002 0.001 0.0005 0.0002 10−3 10−2 10−1 100 101 102 frequency [Hz] Figure 2.5: Example of frequency dependence of the intrinsic attenuation Q−1 for a GZB of NM=3. The thick solid line shows the attenuation of the entire model, while the dotted lines show the attenuation model for each constituent of the Zener body. The vertical lines show the specified minimum and maximum frequencies for the constant Q range. 17 fq min Minimum frequency of the constant-Q model. fq max Maximum frequency of the constant-Q model. fq ref Reference frequency at which the velocity model is given. The Q−1 value becomes nearly constant between the frequencies of fq min and fq max, as shown in Figure 2.5. Outside of the band, the attenuation becomes weaker with increasing/decreasing frequency. As shown in this figure, the nearly constant Q−1 is achieved using three different viscoelastic bodies. If one needs to make Q−1 constant over a wider frequency range, the hard-coded parameter NM should be increased. However, this leads to a significant increase in the memory usage and computational loads. The frequency dependence of Q−1 with the parameters specified above can be investigated using the program qmodel tau.x (see Section 3.3.1). 2.6 2.6.1 Simulation Data Output Output Datafile Format OpenSWPC can export two types of data: spatiotemporal snapshots and the seismic waveform at stations. For snapshot files, the user may choose from an originally defined binary format or a NetCDF file (recommended). The waveforms are usually exported in SAC format. The endian conversion is not performed at the time of the data output. However, the official libraries of NetCDF and SAC automatically detect the endian format and convert them if necessary. Therefore, users do not have to worry about the differences in endian formats between machines. There is a utility program to read the original-formatted data. Note that the binary format may have slight differences for different versions of OpenSWPC. Because the format change is tracked, backward compatibility is always assured. It is recommended to use the same version of the simulation code. For SAC files, the header components described in Table 2.1 are automatically set. The units of SAC files are nm/s for velocity and nm for displacement, following the standard of SAC. While the earthquake source may be represented by multiple point sources, the header always represents the source listed in the first line of the source input file. The snapshot file contains the header information listed in Table 2.2. These headers are commonly defined in either original format or NetCDF. For NetCDF, these headers are set as global attributes. The other headers are set following the COARDS Conventions 1 and the CF Convention 2 . Note that the horizontal directions of the snapshot and waveforms are same as the coordinate of computation. The x- and y-axes correspond to the north and east directions only if phi=0. For the waveform, this angle, phi, is stored in the cmpaz header. The vertical-component waveform is defined as positive upward. 1 http://ferret.wrc.noaa.gov/noaa_coop/coop_cdf_profile.html 2 http://cfconventions.org 18 Table 2.1: SAC headers automatically set by OpenSWPC. Header name Description kevnm title of the parameter file evlo,evla,evdp The location of the event (in degrees for horizontal, in m for depth) o Origin time of the event listed in the first line of the source list kzdate,kztime Date and time of the execution of the simulation code b tbeg of the parameter file delta ntdec w × dt mag The moment magnitude converted from the seismic moment Moment tensor (mxx , m yy , mzz , m yz , mxz , mxy ) of the first line of the source file user0,..,user5 user6,user7,user8 clon, clat, phi of the parameter file kstnm stnm of the parameter file stlo,stla,stdp Station location (in degrees for horizontal, in m for depth) kcmpnm Vx, Vy, Vz for velocities or Ux, Uy, Uz for displacements cmpinc,cmpaz Station directions according to the coordinate specification 7 for velocity, 7 for displacement idep MPI node #0 snap is=1 snap is=2 MPI node #1 snap is=3 snap is=4 MPI node #2 snap is=5 snap is=6 snap is=7 MPI node #3 snap is=8 snap is=9 i Figure 2.6: Schematic of the spatial decimation for the snapshot output. The vertical dotted lines show the borders of the MPI nodes. In this example, the data at the blue grids will be exported as the snapshot data. 2.6.2 Snapshot Data Output Spatiotemporal snapshot output may be created along cross sections of xy, yz, and xz profiles on the topography (fs) and/or on the bathymetry (ob). There are three types of snapshots: divergence and rotation of the velocity (ps), velocity (v), and displacement (u). The use of spatial and temporal decimations is recommended to reduce the I/O load and export data size. Decimation in time is specified by ntdec s starting from it=0 (before starting the computation). In space, the decimations are performed by factors of idec, jdec, and kdec, and OpenSWPC tries to export the center of the decimation window, as schematically shown in Figure 2.6. The numbers of exporting grids in each MPI node do not necessarily need to be the same for each node. The amplitudes of these snapshot points will be gathered to specific nodes (see Table 2.3) and exported as single files. snp format Datafile format of the snapshot files: "native" (original binary format) or "netcdf". 19 Var bintype codetype hdrver title exedate coordinate datatype ns1,ns2 beg1,beg2 ds1,ds2 dt na1,na2 nmed nsnp clon,clat v1,v2,v3 Table 2.2: Snapshot headers set by OpenSWPC. Type Description character(8) Fixed to “STREAMIO” character(8) “SWPC 3D”, “SWPC PV”, or “SWPC SH” depending on the code integer Header version character(80) title in the parameter file integer Date and time of the execution in POSIX time character(2) Snapshot cross section: ’xy’, ’xz’, ’yz’, ’fs’, or ’ob’ character(2) Data type: ’ps’, ’v2’, or ’v3’ integer Number of data samples along the first and second axes real Coordinate value at the first data point of the axes real Snapshot grid spacing real Time step width of the snapshot real Grid numbers of the absorbing boundary layer in the snapshot integer Number of stored medium parameters integer Number of snapshots per time step real clon, clat in the parameter file real Currently not being used Although the NetCDF file format is recommended for convenience in data handling, the use of this format may lead to a slight (∼ 10 %) increase in computation time. xy ps%sw, xz ps%sw, yz ps%sw, fs ps%sw, ob ps%sw Flags for exporting snapshot files of the PS files (.true. or .false.). If they are set to .true., the divergence and rotation vector of the particle velocity will be exported. xy v3%sw, xz v3%sw, yz v3%sw, fs v3%sw, ob v3%sw Flags for exporting snapshot files of the velocities. xy u3%sw, xz u3%sw, yz u3%sw, fs u3%sw, ob u3%sw Flags for exporting snapshot files of the displacements. z0 xy Depth (km) of the snapshot cross section. x0 yz X-coordinate value (km) of the snapshot cross section. y0 xz Y-coordinate value (km) of the snapshot cross section. ntdec s Temporal decimation factor of the snapshot output. Snapshots will be exported every ntdec s time steps. idec, jdec, kdec Spatial decimation factor of the snapshot output for the x, y, and z directions. 2.6.3 Seismic Waveform Output Seismic velocity and/or displacement records at specified stations can be obtained as SACformatted files by setting the parameters sw wav v and/or sw wav u to .true.. Displacement 20 Table 2.3: MPI node number for exporting snapshot files. Section Type Node yz PS 0 xz PS mod(1, nproc) xy PS mod(2, nproc) fs PS mod(3, nproc) ob PS mod(4, nproc) yz V mod(5, nproc) xz V mod(6, nproc) xy V mod(7, nproc) fs V mod(8, nproc) ob V mod(9, nproc) yz U mod(10, nproc) xz U mod(11, nproc) xy U mod(12, nproc) fs U mod(13, nproc) ob U mod(14, nproc) Table 2.4: Format of the station location file. Type ’xy’ ’ll’ Format x lon y lat z z name name zsw zsw records are calculated before the decimation, and therefore, they are expected to be more accurate than performing a numerical integration of the output velocity records. The traces are stored in the memory during the computation and are exported at the end. Station locations are given in Cartesian coordinates (xy) or geographical coordinates (ll), as in Table 2.4. In the station list, lines starting with # will be ignored. The depth of the station can be changed depending on the variable zsw in the station list, as shown in Table 2.5. This operation is important because the station near the free surface is occasionally located above the approximated ground surface in air due to the staircase approximation of the topography and bathymetry. Usually, it is recommended to set zsw=’obb’; this setting locates stations one-grid level below the ground surface (or seafloor). Multiple stations can be specified in the station list file. There is no fixed limit on the number of stations. The number of stations is automatically counted, and only the stations inside the computational region will be exported. sw wav v, sw wav u Output velocity (sw wav v) and displacement (sw wav v) traces. ntdec w Decimation factor of the waveform output. For ntdec w=1, traces at every computational time step will be exported. 21 Table 2.5: Station depth specifications. zsw ’dep’ ’fsb’ ’obb’ ’oba’ ’bdi’ (i=0, · · · , 9) Station depth setting Calculate the station location from the given station depth One grid below the free surface (for oceanic areas, the sea surface) One grid below the ocean bottom (seafloor) or ground surface One grid above of the ocean bottom (seafloor) or ground surface Internal velocity discontinuity specified by the velocity model Header name nvhdr ntrace npts Table 2.6: csf headers. Description Format version numbers. Always zero. Number of traces in the file. Number of time samples in the trace. st format Format of the station list file. See Table 2.4. fn stloc Station location filename. wav format Station file format: ’sac’ (usual, recommended) or ’csf’. The csf format Because the SAC format is defined to express the data at one component of one station in a single file, the number of files may become extraordinarily large. In this case, data transfer between computers will become very inefficient. For OpenSWPC version 3.0 or later, users can choose a concatenated SAC format (csf) for the data output by specifying wav format=’csf’. This is a set of SAC binary files connected to a single file, with headers as in Table 2.6. The header consists of three of four-byte floating-point numbers. After the header, SAC-formatted trace records are repeated ntrace times. If this format is selected, the csf files are created at every computation node with a node number in the filename for every component of the traces. 2.6.4 Output Filename Conventions Output data names are determined by the following rules: • Snapshot (odir)/(title).(section).(type).snp • Waveform (odir)/wav/(title).(stnm).(component).sac • Computation time (odir)/wav/(title).tim In the above rules, (section) takes a cross section such as xy or yz. (type) takes v or ps depending on the snapshot data type. (component) takes Vx, Vy, or Vz for the velocity or Ux, Uy, or Uz for the displacement. 22 2.7 2.7.1 Velocity Model Choice of Velocity Model Type Users can choose a uniform (uni), a layered homogeneous medium (lhm), or a NetCDF (grd) file input (grd) velocity type. In addition, a randomly inhomogeneous medium calculated by an external program can be overlaid onto the model. vmodel type Specify the input velocity model. Choose from one of the following. ’uni’ Homogeneous medium with a free surface. The following additional parameters are required: vp0 P-wave velocity [km/s] in the uniform model. vs0 S-wave velocity [km/s] in the uniform model. rho0 Mass density [g/cm3 ] in the uniform model. qp0 QP of the uniform model. qs0 QS of the uniform model. topo0 Topography depth in the uniform model. If this value is greater than zero, seawater is filled from z = 0 to this depth. ’lhm’ Layered Homogeneous Medium. The one-dimensional velocity structure file should be specified as below. fn lhm Medium specification file. Each line specifies the depth of the top of the layer, density, P-wave velocity, S-wave velocity, QP , and QS below that depth. They must be separated by space(s) (see following example). Lines starting with # will be neglected. 1 2 3 4 5 6 7 8 9 10 # depth rho(g/cm ˆ3) vp(km/s) vs(km/s) Qp Qs # ------------------------------------------------------0 2.300 5.50 3.14 600 300 3 2.400 6.00 3.55 600 300 18 2.800 6.70 3.83 600 300 33 3.200 7.80 4.46 600 300 100 3.300 8.00 4.57 600 300 225 3.400 8.40 4.80 600 300 325 3.500 8.60 4.91 600 300 425 3.700 9.30 5.31 600 300 ’grd’ Velocity model input from NetCDF (GMT grd) files. The compilation of OpenSWPC should be performed in accompaniment with the use of the NetCDF library. The following parameters are required. dir grd Directory of the velocity structure (NetCDF) files. 23 fn grdlst List file that specifies the grd files and the associated medium. Each line contains the grd filename (with a single or double quotation mark; recommended), density, P-wave velocity, S-wave velocity, QP , QS , and the layer number integers (0-9) separated by spaces (see following example). Lines starting with # will be neglected. The layer number is used to specify the source or station depth fit to the layer depth. The first NetCDF file will be treated as the ground surface. If the depth of the ground surface is deeper than zero, the depth range from z = 0 to the surface is assumed to be an ocean layer. The grid above the free surface is treated as an air column. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 # grd filename rho vp vs QP QS sw # ------------------------------------------------------’eJIVSM_01_TAB_ .grd ’ 1.80 1.70 0.35 119 70 0 ’eJIVSM_02_BSM_ .grd ’ 1.95 1.80 0.50 170 100 0 ’eJIVSM_03_BSM_ .grd ’ 2.00 2.00 0.60 204 120 0 ’eJIVSM_04_BSM_ .grd ’ 2.05 2.10 0.70 238 140 0 ’eJIVSM_05_BSM_ .grd ’ 2.07 2.20 0.80 272 160 0 ’eJIVSM_06_BSM_ .grd ’ 2.10 2.30 0.90 306 180 0 ’eJIVSM_07_BSM_ .grd ’ 2.15 2.40 1.00 340 200 0 ’eJIVSM_08_BSM_ .grd ’ 2.20 2.70 1.30 442 260 0 ’eJIVSM_09_BSM_ .grd ’ 2.25 3.00 1.50 510 300 0 ’eJIVSM_10_BSM_ .grd ’ 2.30 3.20 1.70 578 340 0 ’eJIVSM_11_BSM_ .grd ’ 2.35 3.50 2.00 680 400 0 ’eJIVSM_12_BSM_ .grd ’ 2.45 4.20 2.40 680 400 0 ’eJIVSM_13_BSM_ .grd ’ 2.60 5.00 2.90 680 400 0 ’eJIVSM_14_BSM_ .grd ’ 2.65 5.50 3.20 680 400 0 ’eJIVSM_15_UPC_ .grd ’ 2.70 5.80 3.40 680 400 0 ’eJIVSM_16_LWC_ .grd ’ 2.80 6.40 3.80 680 400 0 ’eJIVSM_17_CTM_ .grd ’ 3.20 7.50 4.50 850 500 0 ’eJIVSM_18_PH2_ .grd ’ 2.40 5.00 2.90 340 200 1 ’eJIVSM_19_PH3_ .grd ’ 2.90 6.80 4.00 510 300 0 ’eJIVSM_20_PHM_ .grd ’ 3.20 8.00 4.70 850 500 0 ’eJIVSM_21_PA2_ .grd ’ 2.60 5.40 2.80 340 200 2 ’eJIVSM_22_PA3_ .grd ’ 2.80 6.50 3.50 510 300 0 ’eJIVSM_23_PAM_ .grd ’ 3.40 8.10 4.60 850 500 0 node grd MPI node to input the NetCDF data. All NetCDF files are first read by this node, and then, transferred to all nodes via MPI data communication. is ocean In the default (.true.), the depth from z = 0 to the set topography will be treated as an ocean layer. If this parameter is set to .false., the seafloor will be used as a free surface and no seawater will be used. ’user’ A user subroutine defined in src/swpc ∗/m vmodel user.F90 is used for defining velocity model. Recompilation of the code is necessary if this Fortran file is modified. Please refer the comments in the file for input/output of the user subroutine. vcut Cut-off velocity. For the lhm or grd models, a velocity slower than this value will 24 be overwritten by the vcut value. This parameter is used to avoid wavelengths that are too short and violate the wavelength condition (the wavelength is recommended to be longer than 5-10 grids). This substitution will not be performed in the oceanic area. On Treatments of Air and Seawater Layer In OpenSWPC, the air column has a mass density of ρ = 0.001 [g/cm3 ], velocities of VP = VS = 0 [km/s], and intrinsic attenuation parameters of QP = QS = 1010 . In the ocean column, ρ = 1.0 [g/cm3 ], VP = 1.5 [km/s], VS = 0.0 [km/s], and QP = QS = 106 are assumed. The air column is treated as a vacuum with no seismic wave propagation (with zero velocities). However, the mass density must not be zero to avoid division by zero. In the free surface and seafloor, the reduced order of the finite difference is performed according to (Okamoto and Takenaka, 2005; Maeda and Furumura, 2013). These discontinuities are automatically detected as boundaries that change µ and λ from zero to a finite value. 2.7.2 Small-Scale Random Inhomogeneity Users may overlay small-scale velocity inhomogeneities with specified power-law spectra on the background velocity models of ’uni’, ’lhm’, and ’grd’. The small-scale velocity inhomogeneity ξ is defined by external files. From the average velocities VP0 , VS0 , and ρ0 , the fluctuated velocities and density are given as VP = VP0 (1 + ξ) VS = VS0 (1 + ξ) (2.9) ρ = ρ0 (1 + νξ) , where ν = 0.8 is a scaling parameter based on a laboratory experiment (Birch’s law; Sato et al., 2012). Velocity models having this small-scale inhomogeneity are specified by appending rmed to the original velocity models: vmodel type=’uni rmed’, ’lhm rmed’, or ’grd rmed’. For random media generation, the readers are referred to Section 3.4.1. dir rmed A directory name for storing the random media data files. The random media are given as two- or three-dimensional NetCDF files. At each grid location, the velocity fluctuation ξ(I, J, K) is defined. The code automatically reads the corresponding volume from the file; It is not necessary to decompose the NetCDF files into parts for parallel computation. If the computational size (Nx,Ny,Nz) is larger than the random media file size, the media is used repeatedly by applying a circular boundary condition. The simulation codes do not care if the grid sizes of the simulation and the input random media file are identical. Parameters for uni rmed The following parameter is required in addition to the parameters used in vmodel=’uni’. fn rmed0 Name of the random medium file. In this model, the average velocity will be fluctuated based on the input fn rmed0. 25 Parameters for lhm rmed In this model, the small-scale velocity fluctuation is applied to every layer defined by vmodel=’lhm’. It is possible to assign different random velocity models at different layers. The following parameter is substituted in fn lhm: fn lhm rmed List file of the velocity structure. The list file has a similar format to fn lhm; it contains the filenames of the random media files in the rightmost column as in the following example. 1 2 3 4 5 6 7 8 9 10 # depth rho(g/cm ˆ3) vp(km/s) vs(km/s) Qp Qs fn_rmed # ---------------------------------------------------------------------0 2.300 5.50 3.14 600 300 rmedia1 .nc 3 2.400 6.00 3.55 600 300 rmedia1 .nc 18 2.800 6.70 3.83 600 300 rmedia2 .nc 33 3.200 7.80 4.46 600 300 rmedia2 .nc 100 3.300 8.00 4.57 600 300 225 3.400 8.40 4.80 600 300 325 3.500 8.60 4.91 600 300 425 3.700 9.30 5.31 600 300 - In this example, the layers starting from depths of 0 km and 3 km have fluctuations defined in rmedia1.nc, and the layers from 18 km and 33 km are defined in rmedia2.nc. For the layer deeper than 100 km, a dummy filename (-) is given. In this case (i.e., there is no file found), a fluctuation will not be given. The dummy filename is mandatory in this case. Parameters for grd rmed When overlaying the random fluctuations to the layers defined by the model of vmodel=’grd’, it is possible to assign different random media to different layers. The starting depth of the velocity fluctuation can be either the free surface or depths defined by a layer. The filename of the velocity fluctuation is given by the following parameter: fn grdlst rmed A list file that specifies the velocity layer and the random fluctuation files for each layer. The list file has two additional columns at the right: the filename of the random medium and the reference layer number. 1 2 3 4 5 6 7 8 9 10 11 12 13 # grd filename rho vp vs QP QS sw fn_rmed ref # ----------------------------------------------------------------------’eJIVSM_01_TAB_ .grd ’ 1.80 1.70 0.35 119 70 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_02_BSM_ .grd ’ 1.95 1.80 0.50 170 100 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_03_BSM_ .grd ’ 2.00 2.00 0.60 204 120 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_04_BSM_ .grd ’ 2.05 2.10 0.70 238 140 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_05_BSM_ .grd ’ 2.07 2.20 0.80 272 160 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_06_BSM_ .grd ’ 2.10 2.30 0.90 306 180 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_07_BSM_ .grd ’ 2.15 2.40 1.00 340 200 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_08_BSM_ .grd ’ 2.20 2.70 1.30 442 260 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_09_BSM_ .grd ’ 2.25 3.00 1.50 510 300 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_10_BSM_ .grd ’ 2.30 3.20 1.70 578 340 0 ’rmed3d_1 .nc ’ 0 ’eJIVSM_11_BSM_ .grd ’ 2.35 3.50 2.00 680 400 0 ’rmed3d_1 .nc ’ 0 26 14 15 16 17 18 19 20 21 22 23 24 25 ’eJIVSM_12_BSM_ .grd ’ ’eJIVSM_13_BSM_ .grd ’ ’eJIVSM_14_BSM_ .grd ’ ’eJIVSM_15_UPC_ .grd ’ ’eJIVSM_16_LWC_ .grd ’ ’eJIVSM_17_CTM_ .grd ’ ’eJIVSM_18_PH2_ .grd ’ ’eJIVSM_19_PH3_ .grd ’ ’eJIVSM_20_PHM_ .grd ’ ’eJIVSM_21_PA2_ .grd ’ ’eJIVSM_22_PA3_ .grd ’ ’eJIVSM_23_PAM_ .grd ’ 2.45 2.60 2.65 2.70 2.80 3.20 2.40 2.90 3.20 2.60 2.80 3.40 4.20 5.00 5.50 5.80 6.40 7.50 5.00 6.80 8.00 5.40 6.50 8.10 2.40 2.90 3.20 3.40 3.80 4.50 2.90 4.00 4.70 2.80 3.50 4.60 680 680 680 680 680 850 340 510 850 340 510 850 400 400 400 400 400 500 200 300 500 200 300 500 0 0 0 0 0 0 1 0 0 2 0 0 ’rmed3d_1 .nc ’ ’rmed3d_1 .nc ’ ’rmed3d_1 .nc ’ ’rmed3d_1 .nc ’ ’rmed3d_3 .nc ’ ’rmed3d_3 .nc ’ ’rmed3d_2 .nc ’ ’rmed3d_2 .nc ’ ’rmed3d_3 .nc ’ ’rmed3d_2 .nc ’ ’rmed3d_2 .nc ’ ’rmed3d_3 .nc ’ 0 0 0 0 0 0 18 18 18 21 21 21 The reference layer number defines the reference depth plane of the random media. If this number is zero, the depth grid number of the computational model is directly used to assign the random media. This is exactly the same as the behavior of the uni rmed or lhm rmed models. If the nonzero value of the reference layer number NR is specified, the depth of the NR layer is treated as the base plane. The depth grid of the random medium is measured from this depth. Introducing this reference plane, the inclined random media according to the velocity discontinuity (such as the plate boundary) can be specified. In the above example, the 18th and 21st layers are treated as the references of 18–20th and 21–23th layers, respectively. Truncation of Velocity Fluctuations If the magnitude of the velocity fluctuation becomes too large, there can be a spot with nonphysical velocity, such as negative velocity or a velocity too large for the Earth medium. The simulation may be unstable under the following conditions: 1. The fluctuated velocity V = (1 + ξ)V0 exceeds the stability condition for cases with ξ > 0. 2. The velocity has unrealistic negative values for cases with ξ < −1.0. 3. The mass density has negative values for cases with ξ < −1.25. To avoid such situations, OpenSWPC automatically limits the range of the fluctuated velocity to vcut≤ v ≤ 0.95 × vmax , where vcut is an input parameter and vmax is the maximum possible velocity derived from the stability condition. In addition, the following parameter controls the minimum density. rhomin Minimum mass density in g/cm3 . (1.0 g/cm3 by default.) 2.8 2.8.1 Earthquake Source Specification Moment Rate Function This section describes the moment rate functions, Ṁ(t), that can be used in OpenSWPC by choosing the parameter stftype. In the following, all moment rate functions have a duration (or characteristic time) TR and are normalized so that the total moment is 1. 27 1 TR (0 ≤ t ≤ TR ) Box-car function (boxcar) ṁR (t) = Triangle function (triangle) 4t/TR2 T ṁR (t) = −4(t − TR )/TR2 Herrmann function (herrmann) 16t2 /TR3 ṁH (t) = −2(8t2 − 8tTR + TR2 )/TR3 16 (t − TR )2 /T3 R Cosine function (cosine) [ ( )] 1 2πt ṁC (t) = 1 − cos TR TR ( 3π πt sin3 4TR TR ) Küpper Wavelet (kupper) ṁK (t) = t − exp type (texp) ] [ (2π) t 2πt ṁE (t) = exp − TR TR2 2 (2.10) (0 ≤ t ≤ TR /2) (TR /2 < t ≤ TR ) (2.11) (0 ≤ t ≤ TR /4) (TR /4 < t ≤ 3TR /4) (3TR /4 < t ≤ TR ) (2.12) (0 ≤ t ≤ TR ) (2.13) (0 ≤ t ≤ TR ) (2.14) (t ≥ 0) (2.15) Figure 2.7 shows each moment rate function and its Fourier spectrum. The moment rate functions have a roll off of f −1 – f −4 at frequencies of f ≫ 1/TR . To avoid numerical dispersion, the source spectrum should be sufficiently small at the highest target frequency. As this maximum frequency, we adopt fmax = 2/TR for all types of source time functions (the red dotted line in Figure 2.7). If the parameter is appropriately set so that numerical dispersion does not occur at frequencies below fmax , the result should not be contaminated by numerical dispersion. In addition, the uppermost frequency, where the spectrum response of the source time function becomes flat in the frequency domain, is approximately f ≤ 1/(2TR ) (the blue dotted line in Figure 2.7). 2.8.2 Moment Tensor Source The source mechanisms of the faulting are given by a six-component moment tensor or by three parameters of a double couple source (strike, dip, rake). The source locations can be given either by their computational or geographical coordinates. Therefore, there are eight possible formats to describe the source (see Table 2.7). In the program, sources are given as a stress-drop source by using the moment rate function. The moment rate function is chosen from the given six functions (Figure 2.7). They require parameters in the source list file for their starting time T0 , duration TR , and total moment M0 . OpenSWPC can accept multiple point sources as multiple lines in the source list file. There is no fixed limit to the number of sources (in practice, this is determined by the memory size). By gradually changing the starting time and source location, a finite fault rupture can be mimicked. In the source list file, lines starting with # will be ignored. By setting sdep fit, the source depth can be changed so that it fits the medium’s velocity boundary. 28 (a) Boxcar (b) Triangle 10 2 0 100 2 10−1 10 10−2 10−3 1 10 10−1 −2 10−3 1 −4 10−4 10−5 10−5 0 0 10−6 0.0 0.5 1.0 10−6 0.1 1 10 0.0 (c) Herrmann 0.5 1.0 0.1 1 10 0.1 1 10 0.1 1 10 (d) Cosine 10 2 0 100 2 10−1 10 10−2 10−3 1 10 10−1 −2 10−3 1 −4 10−4 10−5 10−5 0 0 10−6 0.0 0.5 1.0 10−6 0.1 1 10 0.0 (e) Kupper 0.5 1.0 (f) t−exp 10 2 0 100 2 10−1 10 10−2 10−3 1 10−1 −2 10−3 1 10−4 10 10−4 −5 10−5 0 0 10−6 0.0 0.5 t/TR 1.0 10−6 0.1 1 10 0.0 f TR 0.5 t/TR 1.0 f TR Figure 2.7: Moment rate functions Ṁ(t) (left) and their Fourier spectra (right). 29 Table 2.7: Format of the source list file. Type ’xym0ij’ ’xym0dc’ ’llm0ij’ ’llm0dc’ ’xymwij’ ’xymwdc’ ’llmwij’ ’llmwdc’ Format x x lon lon x x lon lon y y lat lat y y lat lat z z z z z z z z T0 T0 T0 T0 T0 T0 T0 T0 TR TR TR TR TR TR TR TR M0 M0 M0 M0 MW MW MW MW mxx strike mxx strike mxx strike mxx strike m yy dip m yy dip m yy dip m yy dip mzz rake mzz rake mzz rake mzz rake m yz mxz mxy m yz mxz mxy m yz mxz mxy m yz mxz mxy In this case, the depth in the source list file will be ignored. The layer number should be specified in the fn grd or fn grd rmed list files. stf format Format of the source list file. Choose from ’xym0ij’ or ’llmwdc’ for example. See Table 2.7 for the complete list. stftype Choice of the source time function. Select from ’boxcar’, ’triangle’, ’herrmann’, ’kupper’, ’cosine’, and ’texp’. See Figure 2.7 for these functions. fn stf Filename of the source list. sdep fit Flag to fit the source depth to the velocity discontinuity. ’asis’: do not fit (default). ’bd{i}’(i=1,2,· · · 9): fits to the i-th boundary specified in the rightmost column of fn grdlst. 2.8.3 Body Force Mode A body force source can be used instead of a moment tensor source. In this mode, the three-component force vector ( fx , f y , fz ) should be specified. The force vector is assumed to have a bell-shaped source time function, as in the case of the moment tensor source. Although there is no restriction on the number of body force elements, it is not possible to use both a moment tensor and a body force at the same time. bf mode Flag for the body force mode. If this is .true., the following parameters are used for the body force and the moment tensor source is ignored. stf format Format of the source file. See Table 2.8. stftype Choice of the source time function. Same as the case with a moment tensor source. 30 Table 2.8: File formats of the body force files. Type ’xy’ ’ll’ Format x lon y lat z z T0 T0 TR TR fx fx fy fy fz fz fn stf Filename of the source list file. The format is described in Table 2.8. sdep fit Flag to fit the source depth to a specified velocity discontinuity. Same as the case with a moment tensor source. 2.8.4 Plane Wave Mode A plane wave incident from the bottom can be used as an input source instead of the moment tensor or body force sources. In OpenSWPC, plane wave incidence is achieved by setting the velocity vector and stress tensor components based on the analytic solution of a plane wave propagating upward as the initial condition. The specification of the initial conditions includes the depth of the initial plane wave (pw ztop) and its characteristic length (pw zlen; corresponding to the wavelength), the strike and dip angle of the plane wave (pw strike, pw dip), and the polarization direction (rake angle) in the case of an S-wave (pw rake). See Figure 2.8 for the geometry. The definitions of the strike, dip, and rake parameters follow those of the earthquake source fault geometry of Aki and Richards (2002). For three-dimensional space, pw strike=0 results in the plane dip toward the y-direction (east for phi=0). A rake angle of pw rake=0◦ or pw rake=180◦ will result in pure SH waves whose polarization is parallel to the free surface. The initial plane wave occupies a depth range of pw zlen (km) starting at depths of z =pw ztop at the center of the horizontal coordinate. The depth dependence of the wave amplitude is determined by the source time functions used in the moment rate function as a function of space (Figure 2.8). Via the definition of the source time function, the integration of the initial plane wave along the propagation direction will be normalized to 1. pw mode Flag to use the plane-wave mode. If it is .true., all point-source locations (body force or moment tensor source) will be ignored. pw ztop z-value of the top of the initial plane wave at x = y = 0. pw zlen Characteristic spatial scale of the initial plane wave. pw ps Plane wave type. Choose from ’p’ or ’s’ pw strike Strike direction of initial plane wave in degrees measured from the x-axis. 31 x x, y ke tri s pw_ pw_dip 0 F(z) pw_ztop la riz at io n pw_dip pw_rake pw_z len S- wa ve po pw_ztop y z z Figure 2.8: Geometry of the plane wave specification. (Left) The specification of the uppermost plane and the polarization direction. (Right) The depth cross section of the initial plane wave. pw dip Dip angle of the initial plane wave in degrees. The initial plane wave propagates vertically if this angle is zero. pw rake Polarization direction of initial plane S-wave in degrees measured from the horizontal plane. stftype Source time function type. Same as the cases with the moment tensor or body force sources. The use of the PML absorbing boundary condition (abc type=’pml’; see Section 2.9) is strongly recommended for the case of plane wave incidence. The simple Cerjan’s (abc type=’cerjan’) condition always causes significant contamination by artificial reflections (Figure 2.9). Even when using the PML boundary, the tilted plane wave incidence (with nonzero pw dip angle) causes some amount of artificial reflections. It is highly recommended that the boundary effect be confirmed with snapshot visualization when using this plane wave mode. 2.9 Absorbing Boundary Conditions Users can choose an absorbing boundary condition from the auxiliary differential equation, the complex frequency-shifted perfectly matched layer (ADE CFS-PML Zhang and Shen, 2010), and Cerjan’s sponge condition (Cerjan et al., 1985). 32 Figure 2.9: Snapshots of the absolute values of divergence (red) and rotation (green) for the case of vertical plane S-wave incidence with (a) Cerjan’s condition and (b) PML boundary conditions. The entire computational domain is separated into interior and exterior regions by the thickness of the absorber na, as shown in Figure 2.10. Because this program assumes the existence of a free surface and ignores acoustic waves in the air column, the waves in the top boundary will not be absorbed. At a given horizontal grid location (I,J), the depth grid deeper than kbeg a will be used as the attenuator. For computational efficiency in the PML boundary condition, OpenSWPC does not solve the viscoelastic constitutive equation in the absorber. Note that, in the case of a medium having very small Q values, this may lead to a velocity gap between the interior and exterior regions due to physical dispersion. For Cerjan’s absorbing condition, the parameters suggested by Cerjan et al. (1985) are embedded in the source code. However, these parameters are scaled according to the width of the absorber na. The PML absorber is usually far superior to Cerjan’s sponge in its efficiency in avoiding artificial reflection from the boundaries. However, PML occasionally results in numerical instabilities, particularly for a medium with a strong velocity contrast and after several time steps. In such cases, Cerjan’s sponge always gives a very stable result. abc type Type of the absorbing boundary condition. Choose from ’pml’ or ’cerjan’. na Thickness of the absorbing layer in numbers of grids. Usually, 10-20 grids are chosen. stabilize pml The low velocity layer is removed if this flag is .true., to stabilize PML. 33 Vaccum (no seismic wave) free surface x,y Interior Region z (wave propagation) kbeg_a Exteria Region (Absorber) na Figure 2.10: Schematic of the definition of the absorber region. The red dotted line indicates the location of kbeg a(I,J). 2.10 Checkpointing and Restarting Some large-scale computers limit the computational time of a single job. To achieve longduration computation, OpenSWPC can export all memory contents to files at specific times (checkpointing), and then continue the simulation as another job (restarting). If this function is turned on, OpenSWPC will terminate the computation after an elapsed time of ckp time (in seconds) and will export all memory images. For the next job, OpenSWPC first tries to find the directory cdir to locate the checkpointing file. If there are checkpointing files, OpenSWPC reads them to continue the simulation. Otherwise, OpenSWPC starts the simulation from scratch. After finishing the computation of all time steps, OpenSWPC removes most of the contents of the checkpointing files. However, it does not delete the checkpointing files. This is to avoid unexpectedly starting the computation from the beginning and overwriting the output files. This function is only available for the three-dimensional simulation code (swpc 3d.x). is ckp The flag to use checkpointing/restarting. cdir Output directory name of the checkpointing file. At restart, the checkpointing files are assumed to be in this directory. ckp time Checkpointing time in seconds. ckp interval Investigate if the computation time exceeds ckp time periodically at this interval. Setting this interval step as too small may affect the performance of the computation. 34 ckp seq Sequential output mode. If this flag is .true., the I/O of the checkpointing files is sequentially performed from the zero-th MPI node. If the file system is shared by several computational nodes, this flag effectively improves the I/O performance. 2.11 Reciprocity Mode This mode excites the seismic wave at a specified station location and exports the velocity and/or strain velocity of multiple virtual source locations. Based on the reciprocity theorem, this result corresponds to the body force and/or moment tensor response from virtual source locations observed at specified stations. If the time duration of the source time function is sufficiently short, they can be treated as Green’s functions. If we denote the Green’s tensor, from the virtual source ξ to the receiver r, as Gij (r, t; ξ), this mode simulates the convolution of the spatial derivatives of Green’s tensor with the source time function s(t) as ∂Gix (ξ, t; r) ∂Gix (r, t; ξ) ∗ s(t) = ∗ s(t) ∂ξx ∂ξx ∂Giy (r, t; ξ) ∂Giy (ξ, t; r) (r, t; ξ) ≡ GM2 ∗ s(t) = ∗ s(t) i ∂ξ y ∂ξ y ∂Giz (r, t; ξ) ∂Giz (ξ, t; r) (r, t; ξ) ≡ GM3 ∗ s(t) = ∗ s(t) i ∂ξz ∂ξz ( ) ( ) (2.16) ∂Giy (r, t; ξ) ∂Giz (r, t; ξ) ∂Giy (ξ, t; r) ∂Giz (ξ, t; r) M4 + + Gi (r, t; ξ) ≡ ∗ s(t) = ∗ s(t) ∂ξz ∂ξ y ∂ξz ∂ξ y ) ( ) ( ∂Gix (ξ, t; r) ∂Giz (ξ, t; r) ∂Gix (r, t; ξ) ∂Giz (r, t; ξ) M5 + ∗ s(t) = + ∗ s(t) Gi (r, t; ξ) ≡ ∂ξz ∂ξx ∂ξz ∂ξx ( ) ( ) ∂Gix (r, t; ξ) ∂Giy (ξ, t; r) ∂Gix (r, t; ξ) ∂Giy (ξ, t; r) M6 Gi (r, t; ξ) ≡ + ∗ s(t) = + ∗ s(t), ∂ξ y ∂ξx ∂ξ y ∂ξx (r, t; ξ) ≡ GM1 i which corresponds to the moment tensor response. Optionally, the body-force response GB1 i (r, t; ξ) ≡ Gix (r, t; ξ) ∗ s(t) = Gix (ξ, t; r) ∗ s(t) GB2 i (r, t; ξ) ≡ Giy (r, t; ξ) ∗ s(t) = Giy (ξ, t; r) ∗ s(t) GB3 i (2.17) (r, t; ξ) ≡ Giz (r, t; ξ) ∗ s(t) = Giz (ξ, t; r) ∗ s(t) can be calculated. To use this mode, the users should specify the station name green stnm of the receiver. This station name should be contained in the station list file. OpenSWPC radiates the seismic wave by an excitation force with a direction specified by the green cmp parameter and a source time function of the rise time, green trise. To obtain the full response of all components, three independent simulations with green cmp=’x’, ’y’, and ’z’ are necessary. The virtual source location should be given in the Cartesian or geographical coordinates and depth (the format is described in Table 2.9) with unique integer ID numbers (gid). Multiple virtual source locations can be specified in the simulation. The gids do not need to be sequential. The output file is stored in the directory (odir)/green/(gid) in the SAC format with the name convention (title) (green cmp) mij .sac (for the moment tensor response) or (title) (green cmp) fi .sac (for the body force response). 35 Table 2.9: Virtual source location format for the reciprocity mode. Type ’xyz’ ’llz’ Format x lon y lat z z gid gid The amplitudes of the output files are multiplied by 109 to compare the SAC-formatted files in nm or nm/s units. The vertical component of the output file is changed to be positive upward. However, the derivative with respect to depth is performed according to the original definition of positive downward. green mode Flags to turn the reciprocity mode on. If this is .true., the other earthquake source parameters will be ignored. green stnm Name of the virtual station. This name must be included in the station list. green cmp Component at the virtual receiver. Choose from ’x’, ’y’, or ’z’. green trise Rise time of the source time function convolved with the simulated Green’s function. green bforce If .true., calculate the body force response as well as the moment tensor response. The default setting is .false.. green fmt Format specification of the virtual source location. Choose from ’xyz’ (Cartesian coordinate; default) or ’llz’ (longitude, latitude, and depth). green maxdist The reciprocity wave will only be calculated if the horizontal distance is shorter than this parameter. Specify in units of km. fn glst Name of the virtual source location file. stftype Source time function type. Same as the case with the moment tensor source. ntdec w Temporal decimation factor of the output waveforms. Same as the case with the normal waveform output. 36 x xend cod es N PSV /SH phi pro for 2D jec pro n cro ss sec tio n Station x-z Epicenter xbeg tio jec tio n ybeg 0 yend y Figure 2.11: Cross section for the calculation in the 2D codes on the horizontal (x–y) plane. All stations and epicenters are projected onto the x–z cross section. 2.12 About Two-Dimensional Codes OpenSWPC contains P-SV (swpc psv) and SH (swpc sh) codes, which work with the same parameter file. In these 2D codes, the simulation will be performed along the x − −z cross section of y = 0. The parameters related to the y-direction will be omitted. The MPI partition will, therefore, be 1D, only in the x-direction. Note that all stations and sources outside the cross section will be projected onto the cross section, as schematically shown in Figure 2.11. For plane wave incidence, pw strike and pw rake will be fixed according to the type of code. Only the dip angle (pw dip) can be changed. 2.13 Other Parameters stopwatch mode Measure the computation times at major subroutines and export the accumulated times to (odir)/(title).tim. This function is used for benchmarking and performance tuning. 37 benchmark mode If this flag is .true., the fixed homogeneous medium and single-point moment tensor source will be selected irrespective of the parameter specification. This is used for validation and performance measurements. ipad, jpad, kpad Expand the Fortran array sizes along the x-, y-, and z-directions. In some computer architectures, the computation speed is very sensitive to the array size. In such cases, slightly changing the array size using these parameters may improve the performance. The expanded array will not be used for the simulation. Therefore, the simulation result is not affected by changing this option. 38 Chapter 3 Related Tools 3.1 3.1.1 Snapshot data handling read snp.x Snapshot files in both NetCDF and the originally defined binary format can be extracted or visualized by the program read snp.x. 1 2 read_snp .x -i snapfile [-h] [-ppm|-bmp] [-pall] [-mul var | -mul1 var -mul2 var ...] [-bin] [-asc] [-skip n] -h Print the header information defined in the snapshot, as in the following example. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 > ../ bin/ read_snp .x -i swpc_3d .xz.ps.snp -h [ binary type] : [code type] : [ header version ]: [title] : [date generated ]: [ coordinate ] [data type] [ns1] [ns2] [beg1] [beg2] [ds1] [ds2] [dt] [na1] [na2] [nmed] [nsnp] [clon] [clat] : : : : : : : : : : : : : : : STREAMIO SWPC_3D 3 swpc_3d 1408015126 2014 -08 -14 T11 -18 -46 xz ps 256 256 -63.87500 -9.87500 0.25000 0.25000 0.05000 20 20 3 2 143.50000 42.00000 39 -ppm/-bmp Visualize and export the image files in ppm or bmp format. The ppm or bmp directory will be automatically created in the current directory and image files with sequential numbers will be stored there. If the snapshot file is displacement or velocity, the absolute values of the vertical and horizontal amplitudes will be colored red and green, respectively. For the PS file, the absolute values of the divergence and rotation vector will be colored similarly. If the absolute value option is specified, the black-redyellow-white color palette (similar to the “hot” color palette in GMT) will be adopted. For cross sections along the surface (ob, fs), the topography color map will be overlaid. For other cross sections, the velocity structure in the section will be overlaid. -pall Visualize including the absorbing boundary region. This option works only if it is used with -ppm/-bmp. By default, the absorbing boundary region will be clipped. -mul var|-mul1 var -mul2 var ... Multiply var by the amplitude for visualization. Adjust the visualized color by changing this value. Optionally, by specifying -mul1 or -mul2, for example, one may change the weight of the amplitude by component. -abs Visualize the absolute value of the vector. This only works with the velocity or displacement snapshots. -bin|asc Export the snapshot data to binary (-bin) or ascii (-asc) files. The data file will be created in the automatically created bin or asc directories. The binary formatted data can be directly used in GMT with the xyz2grd module by appending the -bis option. -skip n Skip the first n snapshots for visualization or data exports. 3.1.2 diff snp.x This program takes the difference between two snapshots and exports it to another snapshot file. 1 > diff_snp .x snap1 snap2 diffile The output file format (NetCDF or binary) depends on the input file format. 3.2 3.2.1 Supporting Parameter Settings fdmcond.x The grid width in space and time in the finite difference method is controlled by the stability condition. The wavelength condition will affect the allowed maximum frequency radiated from the source. The tool fdmcond.x can help determine these parameters to satisfy the conditions. After the user specifies several parameters, such as the grid width, maximum frequency (fmax), rise time (Tr), and minimum and maximum velocities in the medium (vmin, vmax), the program can suggest the other parameters. 40 Example 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 > ./ fdmcond .x ---------------------------------------------------------------------FDM CONDITION ---------------------------------------------------------------------- Model Dimension ? --> 3 2) 2D 3) 3D Source Type? --> 3 1) Triangle 2) Herrmann 3) Kupper Parameter 1) dh 2) dh 3) dh 4) dh 5) dh 6) dh 7) fmax 8) Tr 9) vmin Combination ? --> 5 ( space grid), fmax ( space grid), Tr ( space grid), fmax ( space grid), Tr ( space grid), vmin ( space grid), vmin (max freq .) , vmax (rise time) , vmax (min vel .) , vmax (max freq .), (rise time), (max freq .), (rise time), (min vel .), (min vel .), (max vel .), (max vel .), (max vel .), vmax (max vel .) vmax (max vel .) dt (time grid) dt (time grid) vmax (max vel .) dt (time grid) dt (time grid) dt (time grid) dt (time grid) Assumed Parameters : dx = 0.25 dy = 0.25 dz = 0.25 vmin = 0.3 vmax = 8.0 Derived dt fmax Tr 3.2.2 Parameters : <= 0.01546 <= 0.17143 >= 13.41667 mapregion.x The geographical region of the simulation will be automatically determined by the parameters clon, clat, phi, xbeg, ybeg, nx, ny, dx, and dy. The mapregion.x program reads the parameter file and exports the outer edge of the region in longitude and latitude. 1 > mapregion .x -i input.inf -o region .dat 41 If the option -o is omitted, the result will be printed to the standard output on the screen. This program will also estimate the total memory usage in the standard error output. 3.2.3 mapregion.gmt4, mapregion.gmt5 These scripts use mapregion.x to visualize the region by using GMT4 or GMT5. By default, these scripts plot only the region around the Japanese Islands. 3.3 3.3.1 Velocity Structure qmodel tau.x Calculate the frequency dependence of Q−1 and the body wave dispersion from the input parameter file. 1 > qmodel_tau .x -nm [nm] -i [ prm_file ] -f0 [ min_freq ] -f1 [ max_freq ] -nf [ngrid] This discretizes the frequency range from min freq to max freq into ngrid and exports Q−1 ( f ) and physical dispersion. The latter is normalized to 1 at the reference frequency. The parameters related to the viscoelastic body are read from the input parameter file; however, the number of bodies nm should be specified separately because it is hard-coded into the program. 3.3.2 grdsnp.x From the input parameter file, calculate and print the discontinuity of the input NetCDF file in Cartesian coordinates for the simulation (x, y, depth) in the standard output. This program is used to confirm the coordinate transformation and the detailed digital model, and to visualize the model in the computational domain. 1 > grdsnp .x -i [ prm_file ] -g [ grd_file ] 3.4 3.4.1 Generation of Random Media gen rmed3d.x Generate a three-dimensional random medium file. 1 gen_rmed3d .x [-o outfile ] [-nx nx] [-ny ny] [-nz nz] [-epsil epsilon ] [-kappa kappa] [-dx dx] [-dy dy] [-dz dz] [-ax ax] [-ay ay] [-az az] [-ptype ptype] {-seed seed_number } -o outfile Name of the output file. -nx nx -ny ny -nz nz Number of grids in the x-, y-, and z-directions. They must be a power of 2. epsil epsilon Root mean square (RMS) of the velocity fluctuation ε. 42 Figure 3.1: Example of the visualization of a 3D random medium using ParaView. -ax ax -ay ay -az az Characteristic scales in the x-, y-, and z-directions in units of km. -dx dx -dy dy -dz dz Grid width in the x-, y-, and z-directions. They should be identical to the simulation parameters. -ptype ptype Choice of power spectrum density functions (PSDFs) of the random media model in wavenumber space: 1 for Gaussian, 2 for Exponential, and 3 for von Kármán. -kappa kappa The parameter κ for the von Kármán-type PSDF. -seed seed number Specify the seed number of the random variable generation (optional). If this option is not specified, the seed number is automatically generated based on the execution date and time. The random media file will be stored in the NetCDF format. Various software, such as ParaView1 (Figure 3.1) and Panoply2 , can be used for visualization. 1 http://www.paraview.org 2 http://www.giss.nasa.gov/tools/panoply 43 3.4.2 gen rmed2d.x Generate a 2D random media file. Its usage is same as that of gen rmed3d.x, with parameters related to the y-direction omitted. 3.5 3.5.1 Miscellaneous Tools timvis scripts Four scripts, timvis.gmt4, /timvis.gmt5, timvis abs.gmt4, and /timvis abs.gmt5, are used to visualize the elapsed time of the computation obtained with the input parameter stopwatch mode = .true. by using GMT versions 4 and 5. 3.5.2 Geographic Coordinate Conversion The Fortran programs ll2xy.x and xy2ll.x can project and inversely project the geographic and Cartesian coordinates with the same algorithm as OpenSWPC. These tools are provided for OpenSWPC version 3.0 or later. 44 Chapter 4 Additional Materials 4.1 Hints for Parameter Settings The 3D simulation is bounded by the total memory size. The code requires mMP = 116 + 24NM = 188 (NM = 3) bytes (4.1) of memory for the case of mixed precision (MP=DP) with a GNZ viscoelastic body of NM=3. Note that this is a coarse estimate excluding the effect of an absorbing boundary. The computation time can be roughly estimated by the parameter nG , which is defined as the number of spatial and/or temporal grids that one CPU can process in a second. This value depends on the CPU, as shown in Table 4.1. The total computation time can be estimated by tcomp = nx × ny × nz × nt [s], nG × ncore (4.2) where ncore is the number of CPU cores used in the computation. If the estimated time exceeds that provided by the computer system, it is recommended to make the model size smaller and/or to use checkpointing/restarting. 4.2 Hints for Modifying the Code 4.2.1 Defining Your own Velocity Model The velocity structure is defined by the subroutine vmodel ∗, called by the module m medium.F90. These subroutines commonly have the input/output parameters defined in Table 4.2. By creating a Fortran subroutine that returns the medium parameters rho, lam, mu, qp and Table 4.1: Architecture Name Mac Pro 2010 EIC (ERI, UTokyko) The Earth Simulator (3rd gen.) Performance parameter nG . CPU Xeon X5670 2.93GHz Xeon E5-2680 v3 2.5 GHz NEC SX-ACE 45 #core/CPU 6 12 4 nG 6.7 × 106 7.0 × 106 57 × 106 Table 4.2: Input/output specification of subroutines for velocity models. Variable name In/Out Type Description io prm in int I/O number of the input parameter file i0,i1 in int Start/end indices of arrays in x-direction j0,j1 in int Start/end indices of arrays in y-direction k0,k1 in int Start/end indices of arrays in z-direction xc(i0:i1) in real x grid locations yc(i0:i1) in real y grid locations zc(i0:i1) in real z grid locations vcut in real Cut-off velocity rho(k0:k1,i0:i1,j0:j1) out real Mass density [g/cm3 ] lam(k0:k1,i0:i1,j0:j1) out real Lame coefficient λ [g/cm3 ] mu(k0:k1,i0:i1,j0:j1) out real Lame coefficient µ [g/cm3 ] qp(k0:k1,i0:i1,j0:j1) out real QP qs(k0:k1,i0:i1,j0:j1) out real QS bddep(i0:i1,j0:j1,0:NBD) out real Discontinuity boundary depths [km] qs at locations given in the input of the subroutines xc, yc, and zc, it is easy to add a new velocity model. The topography and bathymetry are automatically investigated in the m medium module after calling the vmodel ∗ routine. To make this investigation work properly, the medium parameter mu must be zero in the air and ocean columns and lam must be zero in the air column. The variables bddep(:,:,0) are assumed to be the topography, and are used for the snapshot output. The other values of bddep(:,:,1:NBD) are used to fit the source and/or station location to the discontinuity depths. Providing dummy values of these functions is not necessary. 4.2.2 Defining Your own Source Time Function The source time function is called by the source momentrate Fortran function in m source.F90 based on the choice of stftype. The definitions of the source time functions are given in share/m fdtool.F90. It is easy to add a new source time function here and to add the call to the new function in the m source module. All of the pre-defined source time functions take two time parameters, tbeg and trise. In the source code, they are stored in the array variable srcprm(:). If the new source time function requires more than three parameters, the user can expand the array srcprm(:) to store them. 4.2.3 Appending New Control Parameters In many Fortran modules, the first set-up is performed by subroutines called (modulename) setup during the first computation. Some of the setup modules read parameters from the input parameter file. These parameters are read by the subroutine readini, which is defined in shared/m readini.F90. 46 Bibliography Aki, K., and P. G. Richards (2002), Quantitative Seismology: Theory and Methods, 2nd eidtion ed., University Science Books. Blanch, J. O., J. O. Robertsson, and W. W. Symes (1994), Modeling of a constant Q: methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique, Geophysics, 60, 176–184, doi:10.1111/j.1365-246X.2004.02300.x. Cerjan, C., D. Kosloff, R. Kosloff, and M. Reshef (1985), A nonreflecting boundary condition for discrete acoustic and elastic wave equations, Geophysics, 50(4), 705–708. Koketsu, K., H. Miyake, and H. Suzuki (2012), Japan Integrated Velocity Structure Model Version 1, Proceedings of the 15th World Conference on Earthquake Engineering, p. Paper No.1773. Maeda, T., and T. Furumura (2013), FDM Simulation of Seismic Waves, Ocean Acoustic Waves, and Tsunamis Based on Tsunami-Coupled Equations of Motion, PAGEOPH, 170(12), 109–127, doi:10.1007/s00024-011-0430-z. Maeda, T., T. Furumura, S. Noguchi, S. Takemura, S. Sakai, M. Shinohara, K. Iwai, and S.-J. Lee (2013), Seismic- and tsunami-wave propagation of the 2011 Off the Pacific Coast of Tohoku Earthquake as inferred from the tsunami-coupled finite-difference simulation, Bulletin of the Seismological Society of America, 103(2B), 1456–1472. Okamoto, T., and H. Takenaka (2005), Fluid-solid boundary implementation in the velocitystress finite-difference method, Zisin, 57, 355–364. Robertsson, J. O., J. O. Blanch, and W. W. Symes (1994), Viscoelastic finite-difference modeling, Geophysics, 59(9), 1444–1456, doi:10.1190/1.1443701. Sato, H., M. C. Fehler, and T. Maeda (2012), Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Second Edition, Springer Berlin Heidelberg, Berlin, Heidelberg, doi:10.1007/978-3-642-23029-5. Zhang, W., and Y. Shen (2010), Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling, Geophysics, 75(4), T141– T154, doi:10.1190/1.3463431. 47 Acknowledgments This project was supported by the Collaborative Research Program of the Earthquake Research Institute at the University of Tokyo (2015-B-01), the Core-to-Core Collaborative Research Program of the Earthquake Research Institute at the University of Tokyo, and the Disaster Prevention Research Institute at Kyoto University (2016-K-06). The code was developed through research collaborations with Takashi Furumura, Shunsuke Takemura, Masaru Todoriki, Futoshi Mori, Nana Yoshimitsu, Hiroyuki Kumagai, Hanae Morioka, Aitaro Kato, Issei Doi, Nozomi Kanaya, and Takehi Isse. The author would like to thank Enago for the English language review. 48 Revision History 2015-06-04 First closed version for the ERI/UT joint usage program. 2015-06-10 Revision for the new Earth Simulator. 2015-06-29 Added random media. 2015-07-14 MPI/OpenMP hybrid parallel simulation mode. 2015-12-04 Text revision. 2016-01-14 Body force and reciprocity modes. 2016-02-03 Output in NetCDF format. 2016-05-05 (v1.0) Official open-source release, 2016-06-19 (v2.0) Hybrid parallel simulation for 2D codes. 2016-08-21 (v3.0) Improved reciprocity mode, geographic projection tools, and csf waveform format. 2017-09-21 (v4.0) Minor bugfixes, new binary output for waveform, updated references. 49 Copyright & License This software is provided under the MIT license. Copyright (c) 2013-2017 Takuto Maeda THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF, OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. The author requests that the user cite (at least one of) the following papers in any publications that result from the use of this software, although this is not an obligation. Accompanying Paper • Maeda, T., S. Takemura, and T. Maeda (2017), OpenSWPC: An open-source integrated parallel simulation code for modeling seismic wave propagation in 3D heterogeneous viscoelastic media, Earth Planets Space, 69, 102, doi:10.1186/s40623-017-0687-2. Related Papers • Furumura, T. and L. Chen (2004), Large scale parallel simulation and visualization of 3D seismic wavefield using the Earth Simulator, Computer Modeling of Engineering and Sciences, 6(2), 153-168. • Furumura, T. and L. Chen (2005), Parallel simulation of strong ground motions during recent and historical damaging earthquakes in Tokyo, Japan, Parallel Computing, 31, 149-165. • Furumura, T. Hayakawa, M. Nakamura, K. Koketsu, and T. Baba (2008), Development of long-period ground motions from the Nankai Trough, Japan, earthquakes: Observations and computer simulation of the 1944 Tonankai (Mw8.1) and the 2004 SE Off-Kii Peninsula (Mw7) Earthquakes, Pure Appl. Geophys., 165, 585-607. • Furumura, T. and T. Saito (2009), An integrated simulation of ground motion and tsunami for the 1944 Tonankai earthquake using high-performance super computers, J. Disast. Res., 4(2), 118-126. 50 • Maeda, T., and T. Furumura (2013), FDM simulation of seismic waves, ocean acoustic waves, and tsunamis based on tsunami-coupled equations of motion, Pure Appl. Geophys., 170(1-2), 109-127, doi:10.1007/s00024-011-0430-z. • Noguchi, S., T. Maeda, and T. Furumura (2013), FDM simulation of an anomalous later phase from the Japan Trench subduction zone earthquakes, Pure Appl. Geophys., 170(1-2), 95-108, doi:10.1007/s00024-011-0412-1. • Maeda, T., T. Furumura, S. Noguchi, S. Takemura, S. Sakai, M. Shinohara, K. Iwai, S. J. Lee (2013), Seismic and tsunami wave propagation of the 2011 Off the Pacific Coast of Tohoku Earthquake as inferred from the tsunami-coupled finite difference simulation, Bull. Seism. Soc. Am., 103(2b), 1456-1472, doi:10.1785/0120120118. • Maeda, T., T. Furumura, and K. Obara (2014), Scattering of teleseismic P-waves by the Japan Trench: A significant effect of reverberation in the seawater column, Earth Planet. Sci. Lett., 397(1), 101-110, doi:10.1016/j.epsl.2014.04.037. • Noguchi, S., T. Maeda, and T. Furumura (2016), Ocean-influenced Rayleigh waves from outer-rise earthquakes and their effects on durations of long-period ground motion, Geophys. J. Int., 205(2), 1099-1107, doi:10.1093/gji/ggw074. • Takemura, S., T. Maeda, T. Furumura, and K. Obara (2016), Constraining the source location of the 30 May 2015 (Mw 7.9) Bonin deep-focus earthquake using seismogram envelopes of high-frequency P waveforms: occurrence of deep-focus earthquake at the bottom of a subducting slab, Geophys. Res. Lett., 43, 4297-4302, doi:10.1002/2016GL068437. • Yoshimitsu, N., T. Furumura, and T. Maeda (2016), Geometric effect on a laboratoryscale wavefield inferred from a three-dimensional numerical simulation, J. Appl. Geophys., 132, 184-192, doi:10.1016/j.jappgeo.2016.07.002. (A computational study on seismic wave propagation in the scale of ∼ 10 cm) • Maeda, T., K. Nishida, R. Takagi, and K. Obara (2016), Reconstruction of a 2D seismic wavefield by seismic gradiometry, Prog. Earth Planet. Sci., 3, 31. doi:10.1186/s40645016-0107-4. (A proposal of a seismic wavefield analysis, with verification via a numerical simulation) • Todoriki, M., T. Furumura, and T. Maeda (2017), Effects of seawater on elongated duration of ground motion as well as variation in its amplitude for offshore earthquakes, Geophys. J. Int., 208, 226-233, doi:10.1093/gji/ggw388. • Toya, M., A. Kato, T. Maeda, K. Obara, T. Takeda, and K. Yamaoka (2017), Down-dip variations in a subducting low-velocity zone linked to episodic tremor and slip: a new constraint from ScSp waves, Scientific Reports, 7, 2868, doi:10.1038/s41598-017-03048-6. • Morioka, H., H. Kumagai, and T. Maeda (2017), Theoretical basis of the amplitude source location method for volcano-seismic signals, J. Geophys. Res., 122, 6538-6551. doi:10.1002/2017JB013997. 51
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 52 Creator : TeX output 2017.09.21:1029 Producer : dvipdfmx (20170318) Create Date : 2017:09:21 10:29:57+09:00EXIF Metadata provided by EXIF.tools