\n \n \n \n \n \n

\n Synopsis
		sfewefdm < Fwav.rsf ccc=Fccc.rsf den=Fden.rsf sou=Fsou.rsf rec=Frec.rsf wfl=Fwfl.rsf > Fdat.rsf srctype=0 ani=-1 verb=y snap=n free=n dabc=y abcone=n debug=y cfl=y opot=n nqz=sf_n(az) nqx=sf_n(ax) oqz=sf_o(az) oqx=sf_o(ax) nqy=sf_n(ay) oqy=sf_o(ay) nbell=1 jdata=1 nb=100 jsnap=nt fmax=

\nElastic wave equation finite difference modeling in both 2D and 3D, using an explicit time-domain solver.
\n
\n*** Please see the SConstruct in book/tutorial/ewe for a SConstruct that demonstrates how to use
\npredefined functions for using this program. ***
\n
\nThis program is designed to be as generic as possible, and allows you to use files
\nwith arbitrary models, and arbitrary source and receiver geometries. Source types are
\nas generic as possible. Supports arbitrary types of anisotropy as well.
\n
\nThe downside to the generality, is that the program is not as performant as dedicated solvers
\nthat are less flexible. The program is parallelized using OpenMP, so be sure to use a compatible compiler to take
\nadvantage of the performance boost.
\n=========== Rotated Staggered grid ==========================
\nUx,Uz=====================Ux,Uz
\n|| | ||
\n|| ||
\n|| | ||
\n|| tij ||
\n||- - - - - - |- - - - - - -||
\n|| Cij ||
\n|| | ||
\n|| ||
\n|| | ||
\nUx,Uz=====================Ux,Uz
\n=========== OPTIONS =======================================
\nani - The type of anisotropy for this simulation. Valid options:
\nFor 2D:
\nISO/HTI/VTI = 0
\nTTI = 1
\n
\nFor 3D:
\nISO/HTI/VTI = 0
\nTTI = 1
\n
\nVTI, HTI, and Isotropic media are special cases of ISO/HTI/VTI media.
\nTTI media can be represented using TTI media.
\n
\ncfl - Execute the CFL check. If the CFL check fails, then it will cause the program to fail.
\nThe CFL check will check both the stability and accuracy conditions for both p-waves and
\ns-waves. Depending on the type of anisotropy that you specify, the CFL condition will
\nuse a safety factor (that you can override if necessary).
\n
\nNOTE: the CFL condition will return both minimum and maximum
\nconstraints on the grid given your velocity model, desired frequency content, and other
\nparameters. IT IS POSSIBLE TO HAVE NO STABLE, AND ACCURATE SOLUTIONS FOR A GIVEN
\nMODEL WITH GIVEN PARAMETERS. THE CFL CONDITION WILL WARN YOU IF THIS IS THE CASE.
\n
\nYOU MUST SPECIFY fmax Parameter as well!
\n
\n----- STABILITY ------
\nThe stability condition is related to the maximum wave speed and minimum grid sampling
\nas follows:
\n
\ndt < min(dx,dy,dz) / (sqrt(2)*vmax)
\n
\nGiven a time sampling dt, it is possible to determine the minimum dx,dy,dz for stability.
\nvmax is the MAXIMUM velocity for all waves in the model (usually P-wave velocity).
\n
\nFor elastic FD, the P-wave most greatly influences the stability, as it moves fastest
\non the grid.
\n
\nThe stability condition gives us a LOWER bound on the grid sampling for a given dt.
\n
\n------ ACCURACY -------
\nThe accuracy condition is related to the number of gridpoints per wavelength. Thus,
\n
\nsafety*vmin / fmax > N * sqrt(dx^2+dy^2+dz^2)
\n
\nwhere vmin is the minimum wave velocity in the model (usually S-wave), fmax is some
\nrelative measure of the maximum frequency of your wavelet (usually 1.5*peak for Ricker),
\nN is the number of points desired per wavelength (5), and safety is a safety factor that
\nis dependent on the type of anisotropy.
\n
\nFor elastic FD, the S-wave most greatly impacts the accuracy of the solution, as the S-wave
\nis typically much higher frequency and travels at slower wave speeds, meaning shorter
\nwavelengths.
\n
\nThe accuracy condition places an UPPER bound on the grid sampling.
\n
\n---- SAFETY FACTOR -----
\nThe safety factor depends on the type of anisotropy specified, and attempts to place a lower
\nbound on the slowest S-wave velocity (guess):
\n
\nISO/HTI/VTI - (3/4)
\nTTI - (1/2)
\n
\nYou can also override the safety factor using the safety parameter.
\nsafety- Override the safety factor for the CFL condition. This should be a floating point (0-1.0).
\nfmax - An estimate of the highest frequency content in your wavelet (for Ricker use 1.5*peak)
\n
\nfsrf - Use a free surface at the top boundary (z=0).
\nWARNING: The free surface condition appears to introduce numerical artifacts into the simulation.
\nUSE AT YOUR OWN RISK.
\n
\nsnap - Save snapshots of the wavefield every n iterations of the modeling program.
\n
\njsnap - Number of iterations between snapshots of the wavefield.
\n\t i.e. jsnap=10, means save a snapshot every 10 iterations.
\n\t If you had 1000 total iterations, then you would have 100 snapshots total.
\n\t The default, will output no snapshots.
\n
\njdata - Number of time imterations between exporting the data at the receivers.
\n\t i.e. jdata=1, means save a snapshot every iteration, which should be the default.
\n\t This can be used to change the sampling of the data to something different from
\nthe wavelet/wavefield.
\n
\nverb - Print useful information
\ndebug - Print debugging information. This is more detailed than verbose.
\n
\nsrctype - An integer which determines where the source wavelet is injected
\nin the simulation. Valid options are:
\n0 - Acceleration source
\n1 - Displacement source
\n2 - Stress source
\n3 - Tensor source
\nThe default option is 2: Acceleration source.
\nFor Stress, Displacement and Acceleration sources, your wavelet
\nneeds to have only 3 components (z,x,y).
\nFor a Tensor source, you must specify wavelet components for
\nall 3 (2D) or 6 (3D) tensor components in the following order:
\n2D: tzz, txx, tzx
\n3D: tzz, txx, tyy, tyz, tzx, txy
\n
\nHint: To inject an acoustic source, use a stress source,
\nwith equal components on all three components.
\n
\ndabc - Use a sponge layer to attenuate waves off the edge of the grid. Use this in
\ncombination with the nb parameter.
\nabcone- In addition to the sponge layer, using a severe ramp at the very edge of the expanded
\nsponge layer to severely attenuate zero-incidence waves at the boundaries.
\nIt\'s not clear if this condition actually affects most computations.
\nopot - True: output is second spatial derivative of potentials; False: output wavefield.
\nnbell - Size of gaussian used to linearly interpolate curves. A value of 5 seems to work well.
\nnb - Not listed, but is an important parameter. Allows you to control the size of the sponge
\nlayer for the absorbing boundary condition. If you are getting reflections off the sides,
\nwith dabc=y, then make this number larger (int). This pads the grid by this amount on all sides.
\nFor example:
\n|--------------------------|
\n| ramp layer |
\n|r |--------------------| |
\n|a | nb |r |
\n|m | |~~~~~~~~| |a |
\n|p | | MODEL | |m |
\n| | nb | SPACE | nb |p |
\n| | |~~~~~~~~| | |
\n| | nb | |
\n| |--------------------| |
\n| ramp layer |
\n|--------------------------|
\nnqz, nqx, oqz, oqx, nqy, oqy, - Allows you to set the parameters for the axes. Leave as defaults.
\n
\n=============BOUNDARY CONDITIONS ========================
\n
\nThis code enforces a fixed reflecting boundary condition at the
\nedge of the computational domain. The absorbing sponge is used
\nIN ADDITION to this condition.
\n
\n=============FILE DESCRIPTIONS ========================
\n
\nFdat.rsf - An RSF file containing your data in the following format:
\naxis 1 - source location
\naxis 2 - wavefield component (z,x,y) order
\naxis 3 - Time
\n
\nFwav.rsf - An RSF file containing your wavelet information. For elastic modeling, the wavelet needs
\nto have 3 samples on N1 one for each component Z-X-Y (or just Z-X for 2D). The second
\naxis describes the component as a function of time. The sampling interval, origin time,
\nand number of time samples will be used as the defaults for the modeling code.
\n\t i.e. your wavelet needs to have the same length and parameters that you want to model with!
\n\t Ex:
\n\t 1st axis index
\n\t Z component 0 0 0 0 0 0 0 0 0 1 2 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
\n\t X component 1 0 0 0 0 0 0 0 0 1 2 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
\n\t Y component 2 0 0 0 0 0 0 0 0 1 2 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
\n\t\t\t\t 2nd axis
\nNOTE: For tensor sources, you must have an appropriate number of components. See srctype for more information.
\n
\ncccc - An N+1 dimensional RSF file that contains the values for the stiffness coefficients to be used
\nas the model for the modeling code. So, for 2D, this would be a 3D array of values concatenated
\ntogether in the order as described in the anisotropy section. Each coefficient file contains
\nthe value of that coefficient for every point in space.
\nThe axes for this file are: Axis 1: Z; Axis 2: X; Axis 3: Y;
\n
\nThe stiffness tensor coefficients are defined uniformly as follows, where
\n--x---y---z--(y)-----(y) describes how the coefficients depend on space.
\n|C11 C12 C13 C14 C15 C16|
\n| C22 C23 C24 C25 C26|
\n| C33 C34 C35 C36|
\n| C44 C45 C46|
\n| C55 C56|
\n| C66|
\n
\nThe tensor is assumed to be symmetric.
\n
\nOrder of the coefficients in the N+1 dimensional file...
\n(First coefficient is the first 2D array in the 3D array).
\n2D Anisotropy Modes:
\n
\nISO/HTI/VTI: C11, C33, C55, C13
\n"TTI:" C11, C13, C15, C33, C35, C55
\n***TTI basically allows access to all coefs in 2D, but is not really triclinic media
\n------------------------------------------------------------
\n(First coefficient is the first 3D array in the 4D array).
\n3D Anisotropy Modes:
\n
\nISO/HTI/VTI: C11, C22, C33, C44, C55, C66, C12, C13, C23
\nTTI: C11, C12, C13, C14, C15, C16, C22, C23, C24, C25, C26, C33, C34,
\nC35, C36, C44, C45, C46, C55, C56, C66
\n
\n
\nden - An N dimensional RSF file that contains the valuese for the density to be used for the model.
\nFor 2D, this would be a 2D array.
\n
\nsou, rec -The source and receiver RSF files respectively.
\nThe 1st axis contains the locations for the points like so:
\n[x,y,z]
\nThe second axis is a concatenated list of all points in the list.
\nSo, for an array of receivers, it would look like:
\n[x1,y1,z1]
\n[x2,y2,z2]
\n[x3,y3,z3]
\n[x4,y4,z4]
\n
\nwfl - The name of the file to save the wavefield snapshots to. This will be an N+2
\ndimensional file. The file will be organized as follows:
\n1-2(3) axes, spatial coordinates
\n3(4) axis, wavefield components, in the Z,X,(Y) order
\n4(5) axis, time, sequential snapshots
\n***The parentheses indicate what the axes will be for 3D models.
\n
\ndat - The name of the file to save the receiver data to. The data has the format of:
\n\t spatial coordinates, then the data components of the elastic wavefield in the
\n\t same order as the wavefield. Lastly, time.
\n
\n========== USEFUL COMMANDS ============================= \t
\n
\nTo view the wavefield snapshots (2D case):
\nsfwindow < Fwfl.rsf n3=1 f3=0 | sfgrey gainpanel=a pclip=100 | sfpen
\n
\nTo view the data (2D case):
\nsfwindow < Fdat.rsf n3=1 f3=0 | sfgrey gainpanel=a pclip=100 | sfpen
\n
\n========== TROUBLESHOOTING ===============================
\n
\nIf you aren\'t getting output, or your output is full of Nans, make sure
\nthat you have the proper dimensions for your wavelet files, and that
\nyour input files make sense.
\n
\nMake sure your source and receiver points are located inside the
\nmodel space, otherwise you will get all NaNs and the simulation will
\nrun forever.
\n
\n======= TIPS ========
\n
\nIf the simulation seems to slow down as it\'s running, its a pretty
\ngood indication that the simulation has become unstable and is overflowing
\nwith NaNs.
\n
\n
\nModified by Yuting Duan, Colorado School of Mines,2013-10-22.
\n
\n\n

\n \n \n \n \n \n

\n Parameters
		\n \n \n \n   bool abcone=n [y/n] \tuse sharp brake at end of boundary layer \n \n \n\n \n \n \n   int ani=-1 \tAnisotropy type, see comments \n \n \n\n \n \n \n   file ccc= \tauxiliary input file name \n \n \n\n \n \n \n   bool cfl=y [y/n] \tuse CFL check, will cause program to fail if not satisfied \n \n \n\n \n \n \n   bool dabc=y [y/n] \tuse sponge absorbing BC \n \n \n\n \n \n \n   bool debug=y [y/n] \tprint debugging info \n \n \n\n \n \n \n   file den= \tauxiliary input file name \n \n \n\n \n \n \n   float fmax= \t \n \n \n\n \n \n \n   bool free=n [y/n] \tfree surface flag \n \n \n\n \n \n \n   int jdata=1 \t \n \n \n\n \n \n \n   int jsnap=nt \t \n \n \n\n \n \n \n   int nb=100 \tpadding size for absorbing boundary \n \n \n\n \n \n \n   int nbell=1 \tbell size \n \n \n\n \n \n \n   int nqx=sf_n(ax) \t \n \n \n\n \n \n \n   int nqy=sf_n(ay) \t \n \n \n\n \n \n \n   int nqz=sf_n(az) \t \n \n \n\n \n \n \n   bool opot=n [y/n] \toutput potentials \n \n \n\n \n \n \n   float oqx=sf_o(ax) \t \n \n \n\n \n \n \n   float oqy=sf_o(ay) \t \n \n \n\n \n \n \n   float oqz=sf_o(az) \t \n \n \n\n \n \n \n   file rec= \tauxiliary input file name \n \n \n\n \n \n \n   bool snap=n [y/n] \twavefield snapshots flag \n \n \n\n \n \n \n   file sou= \tauxiliary input file name \n \n \n\n \n \n \n   int srctype=0 \tsource type, see comments \n \n \n\n \n \n \n   bool verb=y [y/n] \tverbosity flag \n \n \n\n \n \n \n   file wfl= \tauxiliary output file name \n \n \n

'