Performance analysis of DG and HDG methods for the simulation of seismic wave propagation in harmonic domain

International audienc

Numerical schemes for the simulation of seismic wave propagation in frequency domain

National audienceFull Waveform Inversion (FWI) is an imaging technique which is widely used for Seismic Imaging. It is an iterative procedure solving 2N harmonic wave equations at each iteration of the algorithm if N sources are used. The number $N$ is usually large (about 1000) and the efficiency of the whole simulation algorithm is directly related to the efficiency of the numerical method used to solve the wave equations.Seismic imaging can be performed by solving time-dependent wave equations but there is an advantage in considering frequency domain. It is indeed not necessary to store the solution at each time step of the forward simulation. This is interesting because seismic imaging involves very large problems with a lot of data. Memory must then be used with attention. The main drawback lies then in solving large linear systems, which represents a challenging task when considering realistic 3D elastic media, despite the recent advances on high performance numerical linear algebra solvers. In this context, the goal of our study is to develop new solvers based on reduced-size matrices without hampering the accuracy of the numerical solution.We consider Discontinuous Galerkin (DG) methods formulated on fully unstructured meshes, which are more convenient than finite difference methods on cartesian grids to handle the topography of the subsurface. DG methods and classical Finite Element (FE) methods mainly differ from discrete functions which are only piecewise continuous in the case of DG approximation. DG methods are then more suitable than Continuous Galerkin (CG) methods to deal with hp-adaptivity. This is a great advantage to DG method which is thus fully adapted to calculations in highly heterogeneous media.Nevertheless, the main drawback of classical DG methods is that they are more expensive in terms of number of unknowns than classical CG methods, especially when arbitrarily high order interpolation of the field components is used. In this case DG methods lead to larger sparse linear systems with a higher number of globally coupled degrees of freedom as compared to CG methods with a same given mesh. In this work we consider a hybridizable DG (HDG) method. The principle of HDG method consists in introducing a Lagrange multiplier representing the trace of the numerical solution on each face of the mesh cells. This new variable exists only on the faces of the mesh and the unknowns of the problem depend on it. This allows us to reduce the number of unknowns of the global linear system. Now the size of the matrix to be inverted only depends on the number of the faces of the mesh and on the number of the degrees of freedom of each face. It is worth noting that for the classical DG method it depends on the number of the cells of the mesh and on the number of the degrees of freedom of each cell. The solution to the initial problem is then recovered thanks to independent elementwise calculation.We have compared the performance of the HDG method with the one of nodal DG methods for the 2D elastic waves propagation in harmonic domain

Performance analysis of DG and HDG methods for the simulation of seismic wave propagation in harmonic domain

International audienc

Discontinuous Galerkin methods for solving Helmholtz elastic wave equations for seismic imaging

International audienceOne of the most used seismic imaging methods is the full wave inversion (FWI) method which is an iterative procedure whose algorithm is the following. The algorithm requires to solve 2N wave equations at each iteration. The number of sources, N, is usually large (about 1000) and the efficiency of the inverse solver is thus directly related to the efficiency of the numerical method used to solve the wave equation. Seismic imaging can be performed in the time domain or in the frequency domain regime. We focus here on the second setting. The drawback of time domain is that it requires to store the solution at each time step of the forward simulation. The difficulties related to frequency domain inversion lie in the solution of huge linear systems, which cannot be achieved today when considering realistic 3D elastic media, even with the progress of high-performance computing. In this context, the goal is to develop new forward solvers that reduce the number of degrees of freedom without hampering the accuracy of the numerical solution.We consider here discontinuous Galerkin (DG) methods which are more convenient than finite difference methods to handle the topography of the subsurface. Moreover, they are more adapted than continuous Galerkin (CG) methods to deal with hp-adaptivity. This last characteristics is crucial to adapt the mesh to the different regions of the subsurface which is generally highly heterogeneous. Nevertheless, the main drawback of classical DG methods is that they are expensive because they require a large number of degrees of freedom as compared to CG methods on a given mesh. In this work we consider a new class of DG method, the hybridizable DG (HDG) method. Instead of solving a linear system involving the degrees of freedom of all volumic cells of the mesh, the principle of HDG consists in introducing a Lagrange multiplier representing the trace of the numerical solution on each face of the mesh. Hence, it reduces the number of unknowns of the global linear systems and the volumic solution is recovered thanks to a local computation on each element.We compare the performances of the HDG method with those of classical nodal DG methods, and we present our first results using a HDG method for the first-order form of the elastic wave propagation equations for 2D realistic test-cases

Hybridizable Discontinuous Galerkin method for solving Helmholtz elastic wave equations

International audienceFull Waveform Inversion (FWI) is an imaging technique which is widely used for Seismic Imaging. It is an iterative procedure solving 2N harmonic wave equations at each iteration of the algorithm if N sources are used. The number N is usually large (about 1000) and the efficiency of the whole simulation algorithm is directly related to the efficiency of the numerical method used to solve the wave equations.Seismic imaging can be performed by solving time-dependent wave equations but there is an advantage in considering frequency domain. It is indeed not necessary to store the solution at each time step of the forward simulation. This is interesting because seismic imaging involves very large problems with a lot of data. Memory must then be used with attention. The main drawback lies then in solving large linear systems, which represents a challenging task when considering realistic 3D elastic media, despite the recent advances on high performance numerical linear algebra solvers. In this context, the goal of our study is to develop new solvers based on reduced-size matrices without hampering the accuracy of the numerical solution.We consider Discontinuous Galerkin (DG) methods formulated on fully unstructured meshes, which are more convenient than finite difference methods on cartesian grids to handle the topography of the subsurface. DG methods and classical Finite Element (FE) methods mainly differ from discrete functions which are only piecewise continuous in the case of DG approximation. DG methods are then more suitable than Continuous Galerkin (CG) methods to deal with hp-adaptivity. This is a great advantage to DG method which is thus fully adapted to calculations in highly heterogeneous media.Nevertheless, the main drawback of classical DG methods is that they are more expensive in terms of number of unknowns than classical CG methods, especially when arbitrarily high order interpolation of the field components is used. In this case DG methods lead to larger sparse linear systems with a higher number of globally coupled degrees of freedom as compared to CG methods with a same given mesh. In this work we consider a hybridizable DG (HDG) method. The principle of HDG method consists in introducing a Lagrange multiplier representing the trace of the numerical solution on each face of the mesh cells. This new variable exists only on the faces of the mesh and the unknowns of the problem depend on it. This allows us to reduce the number of unknowns of the global linear system. Now the size of the matrix to be inverted only depends on the number of the faces of the mesh and on the number of the degrees of freedom of each face. It is worth noting that for the classical DG method it depends on the number of the cells of the mesh and on the number of the degrees of freedom of each cell. The solution to the initial problem is then recovered thanks to independent elementwise calculation. Moreover, the parallelization of the HDG formulation does not induce any additional difficulty in comparison with classical DG.We have compared the performance of the HDG method with the one of nodal DG methods for the 2D elastic waves propagation in harmonic domain