# Numerical methods

Some methodologies that we will developed, extended, verified and applied are briefly outlined:

**Finite differences (FD):**
FD methods are based on discretizations of the governing partial
differential equations on (quasi-)regular grids. Space derivatives and
temporal evolution are carried out through FD approximations of the
partial derivatives. This implies that the spatial operators are of
local nature: only information from neighbouring points is needed to
extrapolate the system in time. This has important consequences for
parallel computing. When implementing FD algorithms on parallel
computers the communication distances are small which enables efficient
algorithms on modern supercomputers. FD algorithms have the advantage
of being conceptually simple. This implies that they are quickly
adapted to specific problems which is one of the reasons why they are
so tremendously popular. They are also the least accurate, unless
clever methods to reduce numerical dispersion are employed. While
spectral elements seem to be an attractive alternative it is not clear
whether modified, optimised FD schemes may not be competitive with
spectral elements. It is one of the goals of this project to answer
this question. FD techniques will continue to play an important role in
earthquake scenario simulations, axi-symmetric global seismic wave
propagation and in reservoir seismology. They will not be applicable to
whole Earth wave propagation. Important extensions of existing FD
methods will include questions on how to model faults with complex
geometries, models for topography and nonlinear elasticity. Recent
modifications may facilitate the implementation of faults and
topography. However, careful comparative studies are necessary to
evaluate their efficiency. Extensions of the FD approach are scheduled
to incorporate large variations of seismic velocities. This may require
the use of a multi-domain approach and space-dependent time steppin.

**Pseudo-spectral methods (PS):**
The accuracy of FD based approximations of partial derivatives is of
low order. An alternative approach is the use of the properties of the
Fourier transform to calculate the derivatives (or interpolations) with
spectral accuracy (i.e. accuracy to machine precision). This has
certain advantages as one needs fewer points per wavelength to achieve
similar accuracy as with the FD approach. In 3D this considerably
reduces the required memory compared to FD. However, there are
drawbacks. The use of spectral techniques leads to the well known
Gibb's phenomenon. When simulating media with strong discontinuities
unphysical high-frequency oscillations appear in the solutions.
Therefore, this approach is suited only for relatively smooth models
(for which one could argue that perturbation methods are more
efficient). However, the potential in the pseudo-spectral approach is
in the combination with FD methods. It is possible to combine PS (e.g.
for horizontal derivatives and laterally smooth velocity variations)
with FD (e.g. for the vertical derivatives which have to account for
material discontinuities). Such an approach has a further advantage in
connection with the implementation on parallel computer architecture.
The *domain decomposition* necessary for parallelization can be
carried out along the FD axis. Therefore, combinations of the PS
approach with other methods will be considered and evaluated in the
course of this project.

**Finite (spectral) elements (F(S)E):**
The FE method originates from considerations of static elasticity in
engineering and is today maybe the most established numerical method in
static and dynamic elasticity. Nevertheless, it has found less
attention in the area of seismology probably due to the fact that its
implementation is considerably more involved than the FE approach. In
addition, the original FE approach involves the solution of (in 3D)
gigantic systems of equations which are not easily solved on parallel
computers. However, during the last decade an extension of the FE
method is entering computational physics and has potential to be the
method of choice in many domains: spectral elements. What are spectral
elements? Instead of approximating an arbitrary function inside an
element by linear (quadratic, cubic) functions, the functions are
approximated using spectral basis functions (e.g. Chebyshev, Lagrange)
which implies that they exactly describe a function on the
corresponding collocation points. A major discovery was the fact that
when using Lagrange polynomials the mass matrix is diagonal and the
system of equations is trivial to solve. This implies that the
extrapolation problem becomes explicit which makes an SE algorithm as
easy to parallelize as an FD approach. Spectral elements have been
applied to local problems and to global problems. Even though it is not
fully clear whether the SE approach is advantageous in all domains of
seismology, we envisage that this method will play an important role in
this project. We intend to carry out careful and objective comparative
studies with the other possible approaches.

**Boundary elements (BE):**
For some problem classes in wave propagation studies it is sufficient
to consider only boundaries instead of volumes. Such problems involve
crack propagation, reflection and transmission of seismic waves at
interfaces embedded in homogeneous media, or diffraction from free
surfaces with topography. There are obvious advantages of this approach
compared to the volume grids used in the techniques mentioned above.
The discretization is reduced to two dimensions. The disadvantage is
that scattering due to short-scale heterogeneities and the interaction
with the behaviour on the boundaries (e.g. cracks, rupture) are not
accounted for. Yet, it is an efficient way of simulating the effects of
boundaries and - as with other approaches which are sub-optimal in the
general case - the BE approach may have interesting properties when
combined with other techniques (e.g. with the FD method to simulate
dynamic rupture on curved boundaries, a project currently under way).
Such combinations will be further investigated.

**Parallel computing:**
As the employment of numerical methods is computational expensive one
focus will be the efficient implementation of algorithms on modern
computer architectures (e.g. parallel processing). Most of the involved
partners have own parallel computational facilities (e.g. Linux
Clusters) or access to national supercomputer centers (with substantial
computer time grants). Some of the techniques above are easier some are
more difficult to parallelize. One of the main objectives of the
proposed RTN is to provide (towards the end of the project) an
extensive archive with programs for seismic wave propagation on many
scales. As computers are always slow computers, those algorithms need
to be provided in parallelized form. It can be considered an advantage
or disadvantage that the software industry was not capable of providing
stable automatic tools to parallelize even relatively basic algorithms
(attempts such as High Performance Fortran failed). This, however, had
the consequence that everybody today is using more or less the same
standard: the message passing interface (MPI). This is a library with
which any Fortran or C program can be linked. MPI programs are - with
some limitations - portable and hardware independent. The (young and
experienced) scientist involved in computational seismology in the
coming years will need to be familiar with the basic MPI concepts and
therefore one of the objectives of this RTN is to provide training in
theory and practice in the field of parallel algorithms and their
implementation.

**Other techniques (particle approaches, unstructured grids):**
Several methodologies that do not fall under the categories above are
being applied by research groups involved in this network. An example
is an approach that uses particle based discrete (non-continuum)
mechanics to model wave propagation in fractured, fluid-saturated rock
under ambient crustal stress conditions. This approach is validated
against specifically designed physical laboratory experiments.
Obviously this is a method that will address problems in reservoir wave
propagation. Furthermore there are attempts to develop numerical
methods for wave propagation that work on arbitrary unstructured
(possible time-evolving) grids. If combined with approaches such as
optimal operators there is a chance that such methods are comparable to
finite elements, while being conceptually more simply.