# 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.