key: cord-0045891-wd8ydk54
authors: Wolf, Sebastian; Gabriel, Alice-Agnes; Bader, Michael
title: Optimization and Local Time Stepping of an ADER-DG Scheme for Fully Anisotropic Wave Propagation in Complex Geometries
date: 2020-05-22
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50420-5_3
sha: 0d20b3a807e29a3f681b8f54775ecd50d72ee422
doc_id: 45891
cord_uid: wd8ydk54

We present an extension of the earthquake simulation software SeisSol to support seismic wave propagation in fully triclinic anisotropic materials. To our best knowledge, SeisSol is one of the few open-source codes that offer this feature for simulations at petascale performance and beyond. We employ a Discontinuous Galerkin (DG) method with arbitrary high-order derivative (ADER) time stepping. Here, we present a novel implementation of fully physical anisotropy with a two-sided Godunov flux and local time stepping. We validate our implementation on various benchmarks and present convergence analysis with respect to analytic solutions. An application example of seismic waves scattering around the Zugspitze in the Bavarian Alps demonstrates the capabilities of our implementation to solve geophysics problems fast.

To successfully model earthquakes and perform seismic simulations, accurate models for the source dynamics and the propagation of seismic waves are needed. For seismic wave propagation, acoustic, isotropic and anisotropic elastic, attenuating and poroelastic materials are the most important rheologies [21] . Seismic anisotropy affects speed and scattering of seismic waves depending on the direction of propagation and can be found on all scales in the solid Earth. Anisotropy can stem from finely layered or cracked materials, the internal crystal structure of minerals or the alignment of ice crystals in glaciers. Anisotropic material behavior is observed in fault zones [16, 17] and accounted for on global scale in refinements of the Preliminary Reference Earth Model [8] . Effective anisotropy on the scales of seismic wavelengths can be modeled by assuming homogeneous materials with directional dependent properties.

Anisotropy is one of the key seismic properties next to velocity and intrinsic attenuation. Locally, at the scale of earthquake fault zones, large variations in anisotropy reflect the strong material contrasts, extreme strains, and high dissipation of damaged rock. At the global scale, variations in anisotropy provide snapshots of our planet's interior that inform our understanding of plate tectonics. Imaging of anisotropy is also crucial in industry contexts such as exploration or geothermal reservoir development and maintenance. All these applications require efficient forward solvers, ideally accounting for physical anisotropy together with the geometrical complexity of the geological subsurface. Highorder accuracy is crucial to resolve small variations of anisotropy, which are often within only a few percent variation of isotropic material, depending on tectonic context.

Anisotropic material behavior has been successfully included in Finite Difference schemes [11, 22] , pseudo-spectral methods [5] , Spectral Element codes [13] and Discontinuous Galerkin (DG) schemes for the velocity-stress formulation [20] and for the velocity-strain formulation [26] . Only few open-source codes exist which are able to simulate seismic wave propagation in anisotropic materials and which are also tailored to run efficiently on supercomputers. The DG ansatz on unstructured grids allows us to include full physical anisotropy as we do not encounter geometrical restrictions. The DG software SeisSol has undergone end-to-end performance optimization over the last years [3, 10, 25] . However, anisotropic effects have been neglected thus far.

In this paper, we present a novel implementation of fully anisotropic wave propagation that exploits SeisSol's high-performance implementation of elementlocal matrix operations and supports local time stepping. We first lay out the physical model and state the governing equations (Sect. 2). In Sect. 3 these equations are discretized using the DG method combined with arbitrary high-order derivative time stepping (ADER-DG). Our main numerics contribution is to introduce a two-sided numerical flux of the Godunov type in conjunction with a free-surface boundary condition based on solving an inverse Riemann problem. We highlight details of how we implemented theses features into the existing code base, and extended it to make use of local time stepping. Here, the key novelty is a general approach to integrate a numerical eigenvalue solver in SeisSol that replaces analytically derived formulas in the respective precomputation steps. In Sect. 5 we verify our implementation against various analytical solutions and community benchmark problems. We also present an updated reference solution for the AHSP1 benchmark [18] , since our implementation revealed physical inconsistencies in the previous community reference solution. To demonstrate the capability of our code to solve real-world geophysical problems we model anisotropically scattering seismic waves radiating from a point source under the strong topography contrasts of Mount Zugspitze in the Bavarian Alps.

Linear elastic materials are characterized by a stress-strain relation in the form

with stress and strain tensors denoted by σ, ∈ R 3×3 . Symmetry considerations reduce the 81 parameters c ijkl to only 21 independent values [2] . Employing Voigt notation we can write the relation in a matrix-vector manner: ⎛ 

This constitutive relation can be combined with the equations of motion of continuum mechanics which we write in the velocity-strain formulation, where the vector Q = σ 11 , σ 22 , σ 33 , σ 13 , σ 23 , σ 13 , u 1 , u 2 , u 3 T defines the quantities of interest. The combined equation reads:

The Jacobian matrices A d , d = 1, 2, 3, can be deduced from the stress-strain relation and have the form

We observe that the second index is constant for each column. To construct the matrices C 2 and C 3 we replace these by 6, 2, 4 and 5, 4, 3 respectively. The blocks R d are the same as for the isotropic case and are detailed in [20] . The material parameters can vary in space. For better readability the space dependence has been dropped in Eq. 

De la Puente et al. [20] presented the numerics of including anisotropic material effects into ADER-DG seismic wave propagation simulations. We here improve the numerical scheme by a two-sided Godunov flux and a free-surface boundary condition, as well as adaptions necessary for local time stepping with anisotropy. A two-sided flux is physically more accurate and allows for coupling between different rheologies. Local time stepping improves performance drastically.

To solve Eq. (3), we follow a DG ansatz [7] . The underlying geometry is approximated by a mesh of tetrahedral elements τ m . For the discretization polynomial ansatz functions Φ l are defined on a reference element τ ref . On each element the numerical solution Q m p is expanded in terms of the basis functions:

Here the function Ξ m : τ m → τ ref is an affine linear coordinate transformation. By Ψ m l we denote the l th basis function transformed to the m th element. On each element Eq. (3) is multiplied by a test function and integration by parts is applied leading to a semi-discrete formulation:

The Jacobians A d pq are element-wise constant. AlsoQ and its time derivative are constant on each cell. This allows us to pull these quantities out of the integrals. Applying a change of variables to the reference element the integrals can be precomputed. Together with an appropriate flux formulation this leads to a quadrature-free numerical scheme.

In DG schemes continuity across element boundaries is only enforced in a weak sense, via the flux term (n d A d pq Q q ) * in Eq. (6) . Hence, a proper numerical flux, which also takes the underlying physics into account, is essential. De la Puente et al. [20] demonstrated anisotropy with one-sided Rusanov flux and discussed an extension to Godunov fluxes for ADER-DG. Two-sided fluxes capture the correct jump conditions of the Rankine-Hugoniot condition on both sides of the interelement boundaries. They have been introduced to SeisSol in [23] for acoustic and (visco)elastic materials. In the following we give an overview over the most important aspects of using two-sided flux formulations and on generalizing the isotropic flux to a two-sided formulation for the anisotropic case.

In Eq. (6), the surface integral over ∂τ m can be dispatched into four integrals over the four triangular faces of each element. We evaluate the flux for each face individually, so we need to transform the quantities Q as well as the stressstrain relation into face-aligned coordinate systems. For anisotropic materials the stress-strain relation is represented by the matrix H, see Eq. (2). We can express the constitutive behavior in the face-aligned coordinate system via the matrixH = N ·H·N T (cf. [4] ), where N is the so-called Bond matrix (cf. [20] ). We define the matrixà to have the same structure as the matrix A 1 but with entries c ij from the matrixH. At each face we have the Jacobians and approximations on the insideà − , Q − and on the outsideà + , Q + .

The Godunov flux approximates the solution at the element boundary by solving a Riemann problem across the element interfaces. First the equations are transformed to a face-aligned coordinate system. The Rankine-Hugoniot condition states that discontinuities travel with wave speeds given via the eigenvalues of the Jacobian. The differences between the quantities are the corresponding eigenvectors. A detailed derivation can be found in [23, 26] .

We observe that the eigenvectors of the Jacobianà have the form 

The eigenvectors and corresponding eigenvalues resemble three incoming and outgoing waves. To take different material values on the inside and outside into account an eigendecomposition of both Jacobiansà − andà + is performed and the matrix R is constructed taking the first three columns from R − and the last three columns from R + . With indicator matrices I − = diag(1, 1, 1, 0, 0, 0, 0, 0, 0) and I + = diag(0, 0, 0, 0, 0, 0, 1, 1, 1), we can then compute the flux as:

Here T is a matrix that rotates Q from the global coordinate system to the facealigned coordinate system. We take into account that the first six components of Q ∈ R 9 represent a symmetric tensor and the last three components represent a vector. Both parts can be rotated independently, so T combines the rotation of tensorial and vectorial quantities. Analogous to inter-element boundaries, we also impose boundary conditions via a specialized flux. For a free surface boundary we want to impose s = σn = 0. In Eq. (8) the term RI + R −1 T −1 Q − is identified with the state at the inside of the inter-element boundary. We can use this fact to construct a flux which will yield the free surface boundary. To do so we set the traction s to zero and compute the velocity u consistently:

Superscripts b denote values at the boundary. The values s − and u − are the traction and velocity in the face-aligned coordinate system. The matrices R 11 and R 21 slice out the first three columns of R and the rows corresponding to the traction respectively the velocity components. The flux is obtained as F = Tà − Q b . Note that we did not specify the non-traction components of σ. As these lie in the null space of the Jacobianà the flux is not altered by their value.

To integrate Eq. (6) in time SeisSol employs the ADER method [6] . ADER time stepping expands each element solution locally in time as a Taylor series up to a certain order. The time derivatives are replaced by the discretized spatial derivatives following the Cauchy-Kowalewski theorem. To update one element we therefore only need the values of the element itself and its four neighbors, which fosters efficient parallelization. ADER time stepping inherently supports local time stepping [7] : Smaller elements or elements with high wave speeds will be updated more often than large elements or elements with low wave speeds. Local time stepping is thus crucial for performant applications that use strong adaptive mesh refinement or where meshes suffer from badly shaped elements. Setups with a heterogeneous material can also benefit substantially from local time stepping.

Each element has to satisfy the stability criterion Δt m < 1

l m v m for the time step size Δt m , where l m and v m denote the in-sphere diameter and maximum wave speed of element τ m , N is the order of the method.

In anisotropic materials the wave speeds depend on the direction of propagation. This has not been considered in previous work (e.g. [7] ). For a fixed direction d we define the matrix M (d) ij = d k c iklj d l . We calculate the wave speeds in direction d from the eigenvalues λ i of the matrix M (d) via v i = λ i /ρ resulting in a primary and two secondary waves (cf. [4] ). The element-wise maximum wave speed is the maximum of these speeds over all directions d.

SeisSol's ADER-DG discretization is implemented via element-local matrix chain multiplications, which allows for high performance on modern CPUs [10, 25] . All required matrices are precomputed in the initialization phase and optimized kernels are generated for the matrix chain operations [24] . In the following we present the most important choices we made for our implementation of anisotropy. Concerning the matrices, the compute kernels of the isotropic case can be reused, just the assembly of the Jacobians and flux matrices differs.

For each element we store the material values. In comparison to isotropic materials 22 instead of 3 values have to be stored. This overhead is negligible compared to the storage required for the degrees of freedoms. For example, a discretization of order 6 requires 504 degrees of freedom per element.

Concerning the storage of the precomputed matrices, only the Jacobians A and the flux matrices G + and G − change between the isotropic and the anisotropic case. Based on the sparsity pattern of a matrix, the biggest rectangular non-zero block of the matrix is stored. We store the matrices A, G + and G − as full 9 × 9 matrices for the isotropic as well as for the anisotropic case, thus no overhead is produced. The underlying data structures are not changed.

Unlike as with isotropic materials the eigenstructure of the Jacobians given in Eq. (7) is hard to express analytically. In the applied scheme the eigendecomposition has to be calculated once for every element. We use the open source software package Eigen3 [9] to obtain the eigenvectors numerically. We chose this approach for three reasons:

(i) Even if an analytic expression for the eigenvalues is available it is lengthy and hence the implementation is error-prone. (ii) From a software engineering point of view the use of the numerical eigenvalue solver replaces a lot of code that was previously needed for each material model individually. We unified the formulation of the Riemann solver for all material models (isotropic, anisotropic and viscoelastic). We use templating to distinguish the assembly of the Jacobians for each model. From then on we can use the same code to precalculate G + and G − . We expect that these software engineering choices make it easy to include additional material models into SeisSol in the future. Also coupling between different material models within the same simulation can be obtained with little overhead. (iii) The question of accuracy and stability of the numerical solver may arise.

But stability of an analytically derived formula is also not guaranteed. Round-off errors and cancellation could drastically influence the accuracy of the derived eigenvectors. With our choice for using a stable numerical solver instead, we circumvent this problem.

Local time stepping is implemented with a clustered scheme to meet the requirements of modern supercomputers [3] . To cluster the elements the required time step for each element has to be known in advance. To obtain the maximum wave speed for one element, we would have to find the maximum wave speed over all directions. This boils down to solving an optimization problem which involves the calculation of eigenvalues of arbitrary matrices. Solving this optimization problem analytically results in lengthy calculations. In practice, the time step is relaxed by a security factor to meet the CFL condition, so the maximum wave speed does not have to be computed exactly. We sample the wave speeds for several directions d and take their maximum as the maximum wave speed v m .

Planar wave analytic descriptions are widely used in wave propagation problems.

Here we present a numerical convergence study to analyze the correct implementation to confirm its expected convergence properties. To this end, we verify our implementation solving the 3-D, anisotropic, seismic wave equations in the form of periodic, sinusoidal waves in a unit-cube as explained in [19] . The computational domain is the unit cube [−1, 1] 3 with periodic boundary conditions. The ansatz for our plane-wave solution is

where ω denotes the frequency and k is the vector of wave numbers. When we combine this with Eq. (3) we see that the initial condition Q 0 has to be a solution of the eigenvalue problem

In the case of linear elasticity there is a zero eigenvalue with multiplicity 3. The other eigenvalues appear pairwise with different signs and correspond to the P wave and two S waves. For isotropic materials the two S waves coincide, whereas for anisotropic media a slow and a fast S wave can be distinguished. For linear PDEs a linear combination of several solutions is a solution again. To take the directional dependence of anisotropic materials into account we superimpose three planar waves with wave number vectors k 1 = π, 0, 0 , k 2 = 0, π, 0 and k 3 = 0, 0, π . For each direction a P wave traveling in the direction of k l and an S wave traveling in the opposite direction has been chosen. When we denote the eigenvectors of the matrix A d k l d with R l and the corresponding eigenvalues with ω l , the analytic solution can be written as

The computational domain is discretized into cubes of edge length h = 

We compare the numerical solution to the analytic solution at time t = 0.02. Figure 1 shows the convergence behavior for the stress component σ 11 in the L 2 -norm. We clearly observe the expected convergence orders. The plots for the L 1 -and L ∞ -norm are comparable. All other quantities also show the expected convergence rates.

The community benchmark LOH1 [18] is designed for isotropic elastic materials. We here use it to validate backwards compatibility, as isotropic elasticity is a special case of anisotropic elasticity. The setup consists of a layered half space On all other surfaces than the free surface (x 3 = 0) we impose absorbing boundary conditions. The mesh is refined around the source and coarsened away from it. The characteristic length is approximately 300 m in the vicinity of the source and grows up to 2000 m towards the boundary. In total the mesh consists of 2.57 million cells. The mesh is large enough that waves do not leave the domain during the computational time and we do not observe problems with artificial reflections. At the same time we keep the computational effort reasonable due to the coarsening towards the absorbing boundaries.

We ran the simulation with convergence order 6 up to a final time of 5 s and compared our solutions to the reference solution using envelope and pulse misfit [14] . We recorded the velocities at all nine receivers. In Table 1 we present the results for the third, sixth and ninth receiver, which are farthest away from the source (the longer the wave is propagated, the less accurately it is typically resolved due to numerical dissipation and dispersion). With a maximal envelope misfit of 1.32% and a maximal phase misfit of 0.20% we are well within the highest level of accuracy defined in the LOH1 benchmark description. To investigate the computational overhead of incorporating anisotropic effects, we compare the execution times for the setup executed with the isotropic and the anisotropic implementation. We ran each simulation 5 times on 100 nodes of SuperMUC-NG (Intel Skylake Xeon Platinum 8174, 48 cores per node) and averaged the total wall time. The average runtime of the anisotropic implementation was 185.0 s, which is only a little slower than the isotropic implementation, which averaged at 184.6 s. The difference between the runtimes may be explained by the more complicated initialization phase for the anisotropic case. We also point out that the deviation in run times in between different runs was larger than the difference between the averages. The computations achieved an average of 990.3 GFlop/s in the isotropic case and 994.7 GFlop/s in the anisotropic case.

The SISMOWINE test suite [18] proposes the following test case for seismic wave propagation in anisotropic materials: The geometry is a homogeneous full space.

The homogeneous material has a density of ρ = 2700 kg m 3 . The elastic response is characterized by the elastic tensor c ij = 0 expect for 

The source is identical to the LOH1 benchmark. We again refine the mesh around the source with a characteristic edge length of about 300 m and coarsen away from the source. This results in a total of 3.98 million cells. To simulate a full space, absorbing boundary conditions are imposed on all six surfaces. Just as for the LOH1 test case we do not encounter artificial reflections from the boundaries. Figure 2 shows an exemplary comparison between our new implementation and the reference solution of SISMOWINE. One can clearly see that the reference solution does not feature a second shear wave whereas our implementation does. Since shear wave splitting is a well-known physical feature of wave propagation in anisotropic media [2] we assume an error in the proposed reference solution. The correctness of our calculation is confirmed in comparison to an analytical reference: We used the open source tool christoffel [12] to compute the wave speeds for the given material depending on the direction of wave propagation. The shown receiver 6 is located at (7348, 7348, 0), 10392 m away from the source. The calculated arrival time for the P wave is 1.73 s, the slow and the fast S wave arrive after 2.59 s and 3.00 s respectively. We observe that the arrival times align very well with the results calculated by the new SeisSol implementation. In agreement with the original authors and the maintainers of the SISMOWINE project, our here presented solution has been accepted as the new reference solution due to its physical plausibility.

For transversally isotropic materials analytical solutions can be found along the axis of symmetry [5] . As mentioned earlier the representation of the elastic tensor depends on the chosen coordinate system. By tilting the symmetry axis of the transversally isotropic material we can generate an almost densely filled elastic tensor. We take the material characterized by the tensor c ij = 0 except for 

and density ρ = 2590 kg m 3 and tilt it around the x axis about 30 • . We consider the computational domain Ω = [0, 2500] 3 . The source is placed at (1250, 1562.5, 937.5) and a receiver is placed at (1250, 1198.05, 1568.75), which is along the symmetry axis of the tilted material 728.9 m away from the source. The source either acts along the axis of symmetry or orthogonal to the axis. The time history is a Ricker wavelet with dominant frequency f 0 = 16.0 Hz and onset time t 0 = 0.07 s. The whole simulation was run for 0.6 s on a mesh which is refined in a sphere around the source and in a cylinder along the axis of symmetry. In the most refined region the characteristic length is 5 m and grows towards the boundary. In total the mesh consisted of 6.10 million elements.

We compare our solution obtained at the receiver with the analytic solution using envelope and pulse misfit. For the horizontal source we obtained a maximal envelope misfit of 2.09% and a pulse misfit of 0.49%. For the vertical source the misfits were 2.16% and 0.28% respectively. For both source types the numerical solution fits the analytic solution very well. 

Accurate numerical simulation of scattering seismic waves by complex geometries are critical for assessing and quantifying seismic hazards. In the context of regional scale seismology wave propagation simulations, many numerical methods are challenged by geometric restrictions or low-order accuracy. We here spatially discretize the surface topography around Mount Zugspitze [1] in an unstructured tetrahedral computational mesh of size 90 km × 90 km up to 70 km depth with a resolution of 600 m at the surface. The mesh contains 1.47 million cells. We chose a discretization of order 6 which results in 740 million degrees of freedom. This means we can resolve frequencies up to 4.2 Hz with an estimated envelope misfit smaller than 1% [15] . A kinematic point source with the same parameters as for the LOH1 test case is placed in the center of the domain at 10 km depth. We visually compare the wave field scattered by topography in an isotropic material with parameters ρ = 2670 kg m 3 , λ = 36.4 GPa, μ = 29.8 GPa with an anisotropic material. The elastic tensor is chosen such that in EWdirection the P wave speed is 6000 m s for both materials. To illustrate the effects of anisotropy the P wave speed in NS-direction of the anisotropic material is 5% lower.

In Fig. 3 snapshots of the vertical velocity field on the free surface are plotted. A circular shape for the isotropic example and an elliptic shape for the anisotropic part illustrate the effects of anisotropic materials on wave propagation under strong velocity contrast. We compare this simulation with and without local time stepping: moving from global to local time stepping drastically reduced the execution time from 5990 s to 210 s. The simulation was run on 50 nodes of SuperMUC-NG. The computational more intense version with global time stepping achieved 1.414 TFlop/s, the version with local time stepping achieved 1.015 TFlop/s. This shows that local time stepping is crucial to obtain fast simulations when the element sizes vary a lot, such as in the case of surface topography.

The earthquake simulation code SeisSol has been successfully extended to take general anisotropic materials into account. A two-sided Godunov flux for anisotropic media has been derived and implemented. Together with the formulation of the free-surface boundary condition as the solution of an inverse Riemann problem it fits well with the other rheological models. Necessary changes to include local time stepping have been described and implemented.

The scheme has been validated against various benchmarks. The expected convergence rates are demonstrated in comparison to analytic solutions. The mismatch between our results on a community benchmark have been discussed with the maintainers and led to an update of the reference solution.

As anisotropy is non-neglectable to describe the Earth's subsurface structure we expect a wide range of applications. Besides the importance of seismic anisotropy for exploration purposes, earthquake fault zones may be characterised by pervasive anisotropy. Earthquake monitoring and forecasting can be built upon this observation.

Quantitative Seismology, 2nd edn

Petascale local time stepping for the ADER-DG finite element method

Wave Fields in Real Media

A spectral scheme for wave propagation simulation in 3-D elastic-anisotropic media

An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -II. The three-dimensional isotropic case

An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -V. Local time stepping and p-adaptivity

Preliminary reference earth model

Eigen v3

Petascale high order dynamic rupture earthquake simulations on heterogeneous supercomputers

Anisotropic wave propagation through finitedifference grids

Solving the christoffel equation: phase and group velocities

Simulation of anisotropic wave propagation based upon a spectral element method

Time-frequency misfit and goodness-offit criteria for quantitative comparison of time signals

Quantitative accuracy analysis of the discontinuous Galerkin method for seismic wave propagation

Observation and modelling of fault-zone fracture seismic anisotropy -I. P, SV and SH travel times

Seismic anisotropy in central north Anatolian fault zone and its implications on crustal deformation

Comparison of numerical methods for seismic wave propagation and source dynamics -the SPICE code validation

Discontinuous Galerkin methods for wave propagation in poroelastic media

An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV

An Introduction to Seismology, Earthquakes, and Earth Structure

3D seismic wavefield modeling in generally anisotropic media with a topographic free surface by the curvilinear grid finite-difference method 3D seismic wavefield modeling in generally anisotropic media

Flexible model extension and optimisation for earthquake simulations at extreme scale

Yet another tensor toolbox for discontinuous Galerkin methods and other applications

Extreme scale multi-physics simulations of the tsunamigenic 2004 Sumatra megathrust earthquake

Full-anisotropic poroelastic wave modeling: a discontinuous Galerkin algorithm with a generalized wave impedance