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=

Elastic wave equation finite difference modeling in both 2D and 3D, using an explicit time-domain solver.

*** Please see the SConstruct in book/tutorial/ewe for a SConstruct that demonstrates how to use
predefined functions for using this program. ***

This program is designed to be as generic as possible, and allows you to use files
with arbitrary models, and arbitrary source and receiver geometries. Source types are
as generic as possible. Supports arbitrary types of anisotropy as well.

The downside to the generality, is that the program is not as performant as dedicated solvers
that are less flexible. The program is parallelized using OpenMP, so be sure to use a compatible compiler to take
advantage of the performance boost.
=========== Rotated Staggered grid ==========================
Ux,Uz=====================Ux,Uz
|| | ||
|| ||
|| | ||
|| tij ||
||- - - - - - |- - - - - - -||
|| Cij ||
|| | ||
|| ||
|| | ||
Ux,Uz=====================Ux,Uz
=========== OPTIONS =======================================
ani - The type of anisotropy for this simulation. Valid options:
For 2D:
ISO/HTI/VTI = 0
TTI = 1

For 3D:
ISO/HTI/VTI = 0
TTI = 1

VTI, HTI, and Isotropic media are special cases of ISO/HTI/VTI media.
TTI media can be represented using TTI media.

cfl - Execute the CFL check. If the CFL check fails, then it will cause the program to fail.
The CFL check will check both the stability and accuracy conditions for both p-waves and
s-waves. Depending on the type of anisotropy that you specify, the CFL condition will
use a safety factor (that you can override if necessary).

NOTE: the CFL condition will return both minimum and maximum
constraints on the grid given your velocity model, desired frequency content, and other
parameters. IT IS POSSIBLE TO HAVE NO STABLE, AND ACCURATE SOLUTIONS FOR A GIVEN
MODEL WITH GIVEN PARAMETERS. THE CFL CONDITION WILL WARN YOU IF THIS IS THE CASE.

YOU MUST SPECIFY fmax Parameter as well!

----- STABILITY ------
The stability condition is related to the maximum wave speed and minimum grid sampling
as follows:

dt < min(dx,dy,dz) / (sqrt(2)*vmax)

Given a time sampling dt, it is possible to determine the minimum dx,dy,dz for stability.
vmax is the MAXIMUM velocity for all waves in the model (usually P-wave velocity).

For elastic FD, the P-wave most greatly influences the stability, as it moves fastest
on the grid.

The stability condition gives us a LOWER bound on the grid sampling for a given dt.

------ ACCURACY -------
The accuracy condition is related to the number of gridpoints per wavelength. Thus,

safety*vmin / fmax > N * sqrt(dx^2+dy^2+dz^2)

where vmin is the minimum wave velocity in the model (usually S-wave), fmax is some
relative measure of the maximum frequency of your wavelet (usually 1.5*peak for Ricker),
N is the number of points desired per wavelength (5), and safety is a safety factor that
is dependent on the type of anisotropy.

For elastic FD, the S-wave most greatly impacts the accuracy of the solution, as the S-wave
is typically much higher frequency and travels at slower wave speeds, meaning shorter
wavelengths.

The accuracy condition places an UPPER bound on the grid sampling.

---- SAFETY FACTOR -----
The safety factor depends on the type of anisotropy specified, and attempts to place a lower
bound on the slowest S-wave velocity (guess):

ISO/HTI/VTI - (3/4)
TTI - (1/2)

You can also override the safety factor using the safety parameter.
safety- Override the safety factor for the CFL condition. This should be a floating point (0-1.0).
fmax - An estimate of the highest frequency content in your wavelet (for Ricker use 1.5*peak)

fsrf - Use a free surface at the top boundary (z=0).
WARNING: The free surface condition appears to introduce numerical artifacts into the simulation.
USE AT YOUR OWN RISK.

snap - Save snapshots of the wavefield every n iterations of the modeling program.

jsnap - Number of iterations between snapshots of the wavefield.
i.e. jsnap=10, means save a snapshot every 10 iterations.
If you had 1000 total iterations, then you would have 100 snapshots total.
The default, will output no snapshots.

jdata - Number of time imterations between exporting the data at the receivers.
i.e. jdata=1, means save a snapshot every iteration, which should be the default.
This can be used to change the sampling of the data to something different from
the wavelet/wavefield.

verb - Print useful information
debug - Print debugging information. This is more detailed than verbose.

srctype - An integer which determines where the source wavelet is injected
in the simulation. Valid options are:
0 - Acceleration source
1 - Displacement source
2 - Stress source
3 - Tensor source
The default option is 2: Acceleration source.
For Stress, Displacement and Acceleration sources, your wavelet
needs to have only 3 components (z,x,y).
For a Tensor source, you must specify wavelet components for
all 3 (2D) or 6 (3D) tensor components in the following order:
2D: tzz, txx, tzx
3D: tzz, txx, tyy, tyz, tzx, txy

Hint: To inject an acoustic source, use a stress source,
with equal components on all three components.

dabc - Use a sponge layer to attenuate waves off the edge of the grid. Use this in
combination with the nb parameter.
abcone- In addition to the sponge layer, using a severe ramp at the very edge of the expanded
sponge layer to severely attenuate zero-incidence waves at the boundaries.
It's not clear if this condition actually affects most computations.
opot - True: output is second spatial derivative of potentials; False: output wavefield.
nbell - Size of gaussian used to linearly interpolate curves. A value of 5 seems to work well.
nb - Not listed, but is an important parameter. Allows you to control the size of the sponge
layer for the absorbing boundary condition. If you are getting reflections off the sides,
with dabc=y, then make this number larger (int). This pads the grid by this amount on all sides.
For example:
|--------------------------|
| ramp layer |
|r |--------------------| |
|a | nb |r |
|m | |~~~~~~~~| |a |
|p | | MODEL | |m |
| | nb | SPACE | nb |p |
| | |~~~~~~~~| | |
| | nb | |
| |--------------------| |
| ramp layer |
|--------------------------|
nqz, nqx, oqz, oqx, nqy, oqy, - Allows you to set the parameters for the axes. Leave as defaults.

=============BOUNDARY CONDITIONS ========================

This code enforces a fixed reflecting boundary condition at the
edge of the computational domain. The absorbing sponge is used
IN ADDITION to this condition.

=============FILE DESCRIPTIONS ========================

Fdat.rsf - An RSF file containing your data in the following format:
axis 1 - source location
axis 2 - wavefield component (z,x,y) order
axis 3 - Time

Fwav.rsf - An RSF file containing your wavelet information. For elastic modeling, the wavelet needs
to have 3 samples on N1 one for each component Z-X-Y (or just Z-X for 2D). The second
axis describes the component as a function of time. The sampling interval, origin time,
and number of time samples will be used as the defaults for the modeling code.
i.e. your wavelet needs to have the same length and parameters that you want to model with!
Ex:
1st axis index
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
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
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
2nd axis
NOTE: For tensor sources, you must have an appropriate number of components. See srctype for more information.

cccc - An N+1 dimensional RSF file that contains the values for the stiffness coefficients to be used
as the model for the modeling code. So, for 2D, this would be a 3D array of values concatenated
together in the order as described in the anisotropy section. Each coefficient file contains
the value of that coefficient for every point in space.
The axes for this file are: Axis 1: Z; Axis 2: X; Axis 3: Y;

The stiffness tensor coefficients are defined uniformly as follows, where
--x---y---z--(y)-----(y) describes how the coefficients depend on space.
|C11 C12 C13 C14 C15 C16|
| C22 C23 C24 C25 C26|
| C33 C34 C35 C36|
| C44 C45 C46|
| C55 C56|
| C66|

The tensor is assumed to be symmetric.

Order of the coefficients in the N+1 dimensional file...
(First coefficient is the first 2D array in the 3D array).
2D Anisotropy Modes:

ISO/HTI/VTI: C11, C33, C55, C13
"TTI:" C11, C13, C15, C33, C35, C55
***TTI basically allows access to all coefs in 2D, but is not really triclinic media
------------------------------------------------------------
(First coefficient is the first 3D array in the 4D array).
3D Anisotropy Modes:

ISO/HTI/VTI: C11, C22, C33, C44, C55, C66, C12, C13, C23
TTI: C11, C12, C13, C14, C15, C16, C22, C23, C24, C25, C26, C33, C34,
C35, C36, C44, C45, C46, C55, C56, C66


den - An N dimensional RSF file that contains the valuese for the density to be used for the model.
For 2D, this would be a 2D array.

sou, rec -The source and receiver RSF files respectively.
The 1st axis contains the locations for the points like so:
[x,y,z]
The second axis is a concatenated list of all points in the list.
So, for an array of receivers, it would look like:
[x1,y1,z1]
[x2,y2,z2]
[x3,y3,z3]
[x4,y4,z4]

wfl - The name of the file to save the wavefield snapshots to. This will be an N+2
dimensional file. The file will be organized as follows:
1-2(3) axes, spatial coordinates
3(4) axis, wavefield components, in the Z,X,(Y) order
4(5) axis, time, sequential snapshots
***The parentheses indicate what the axes will be for 3D models.

dat - The name of the file to save the receiver data to. The data has the format of:
spatial coordinates, then the data components of the elastic wavefield in the
same order as the wavefield. Lastly, time.

========== USEFUL COMMANDS =============================

To view the wavefield snapshots (2D case):
sfwindow < Fwfl.rsf n3=1 f3=0 | sfgrey gainpanel=a pclip=100 | sfpen

To view the data (2D case):
sfwindow < Fdat.rsf n3=1 f3=0 | sfgrey gainpanel=a pclip=100 | sfpen

========== TROUBLESHOOTING ===============================

If you aren't getting output, or your output is full of Nans, make sure
that you have the proper dimensions for your wavelet files, and that
your input files make sense.

Make sure your source and receiver points are located inside the
model space, otherwise you will get all NaNs and the simulation will
run forever.

======= TIPS ========

If the simulation seems to slow down as it's running, its a pretty
good indication that the simulation has become unstable and is overflowing
with NaNs.


Modified by Yuting Duan, Colorado School of Mines,2013-10-22.

Parameters
		bool abcone=n [y/n] use sharp brake at end of boundary layer   int ani=-1 Anisotropy type, see comments   file ccc= auxiliary input file name   bool cfl=y [y/n] use CFL check, will cause program to fail if not satisfied   bool dabc=y [y/n] use sponge absorbing BC   bool debug=y [y/n] print debugging info   file den= auxiliary input file name   float fmax=   bool free=n [y/n] free surface flag   int jdata=1   int jsnap=nt   int nb=100 padding size for absorbing boundary   int nbell=1 bell size   int nqx=sf_n(ax)   int nqy=sf_n(ay)   int nqz=sf_n(az)   bool opot=n [y/n] output potentials   float oqx=sf_o(ax)   float oqy=sf_o(ay)   float oqz=sf_o(az)   file rec= auxiliary input file name   bool snap=n [y/n] wavefield snapshots flag   file sou= auxiliary input file name   int srctype=0 source type, see comments   bool verb=y [y/n] verbosity flag   file wfl= auxiliary output file name