hawen: time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization

Many applications such as seismic and medical imaging, material sciences, or helioseismology and planetary science, aim to reconstruct properties of a non directly accessible or non-visible interior. For this purpose, they rely on waves whose propagation through a medium interrelates with the physical properties (density, sound speed, etc.) of this medium. The methodology for imaging with waves comprises of two main stages illustrated in Figure 1. In the data acquisition stage (Figure 1a), the medium response to probing waves is recorded (e.g., seismic waves from Earthquakes recorded by ground network). In the second stage, we rely on a reconstruction procedure which iteratively updates an initial model of physical parameters, so that numerical simulations approach the measurements (Figure 1b). This procedure is employed, for instance, for seismic (reconstruction of subsurface layers) and medical (disease diagnostic) imaging, see the introduction of Faucher & Scherzer (2020) and the references therein.


Summary
Many applications such as seismic and medical imaging, material sciences, or helioseismology and planetary science, aim to reconstruct properties of a non directly accessible or non-visible interior. For this purpose, they rely on waves whose propagation through a medium interrelates with the physical properties (density, sound speed, etc.) of this medium. The methodology for imaging with waves comprises of two main stages illustrated in Figure 1. In the data acquisition stage (Figure 1a), the medium response to probing waves is recorded (e.g., seismic waves from Earthquakes recorded by ground network). In the second stage, we rely on a reconstruction procedure which iteratively updates an initial model of physical parameters, so that numerical simulations approach the measurements (Figure 1b). This procedure is employed, for instance, for seismic (reconstruction of subsurface layers) and medical (disease diagnostic) imaging, see the introduction of  and the references therein. b) The reconstruction algorithm starts from initial properties and compares simulations of wave propagation with the measurements, then iteratively updates those properties. The green block corresponds to the forward problem (modeling the propagation of waves) and the orange ones to the iterative correction. hawen solves both the forward and inverse problems associated with timeharmonic waves.

Statement of need
hawen is designed to address large-scale inverse wave problems, and includes 1. simulating time-harmonic waves in heterogeneous media (e.g., visco-acoustic and viscoelastic propagation, modal equations in helioseismology); 2. performing iterative minimization to solve the inverse problem, cf. Figure 1.
It combines MPI and OpenMP parallelism, it is deployed on supercomputers and has been used in studies on seismic imaging (Faucher, Alessandrini, et al., 2020;, as well as in helioseismology (Barucq et al., 2019Barucq, Faucher, Fournier, et al., 2020a, 2020b. The software works with an input parameter file: a text file listing the problem configuration (dimension, choice of viscous model, etc.), such that the computational details are transparent for the users. hawen has its dedicated User's guide (Faucher, 2020) that details its utilization and all the available functionalities, it is available on the software dedicated website, which contains illustrations, installation guide and tutorials.
One specificity of hawen is to implement the hybridizable discontinuous Galerkin (HDG) method (Arnold et al., 2002;Cockburn et al., 2009) for the discretization of the wave equations. It helps reduce the cost of the computations by providing smaller linear systems compared to, e.g., finite elements, depending on the order of approximation Kirby et al., 2012). For seismic applications, current software mostly relies on the spectral element method (Komatitsch & Tromp, 1999;Komatitsch & Vilotte, 1998), such as specfem, or on finite differences discretization.
The software handles unstructured domains (allowing input meshes of different formats) and thus accounts for complex geometry. In addition, hawen relies on the massively parallel sparse direct solver MUMPS (Amestoy et al., 2001(Amestoy et al., , 2006 for the matrix factorization, hence handles problems with many right-hand sides. Despite the technicality of the underlying methods, the purpose of hawen is to provide a unified and evolutive framework to address large-scale visco-elastic and visco-acoustic inverse problems, with a user-friendly interface using the input parameter files.

Modeling the propagation of waves
The first feature of hawen is to simulate the propagation of time-harmonic waves in different types of medium: this is the forward problem in Figure 1. The propagation of waves is characterized by a Partial Differential Equation which depends on the type of the waves (e.g., mechanical or electromagnetic), and on the type of medium considered (e.g., fluid or solid) (Carcione, 2007;Faucher, 2017;Slawinski, 2010). hawen handles visco-acoustic and viscoelastic propagation, as well as waves propagation for helioseismology, these are described in the User's documentation (Faucher, 2020). Our implementation is based upon the HDG method and the computational steps are illustrated in Figure 2.

Parallelism:
the mesh cells are distributed among the mpi processors.

Discretization:
Create the global matrix associated with HDG discretization.

Solver:
Solve the linear system with Mumps.

Save:
Save the solution.

Illustrations:
Disk-domain of radius 1m, decomposed in 3040 triangle cells. The wave speed and density are constant c = 0.1m s −1 , ρ = 1000kg m −2 . The point-source is positioned in the center.
Parallelism: one mpi processor handles a part of the domain, non-overlapping with the other ones.
Real part of the pressure field p solution to the acoustic problem (1) at 0.4Hz.

Figure 2:
Illustration of the computational steps for the numerical resolution of the forward problem. For simplicity, we illustrate with an homogeneous medium.

Inverse problem via iterative minimization
In the forward problem, the objective is to simulate the propagation of waves, given the medium physical properties. Conversely, the objective of the inverse problem is to recover the physical properties, given some observations of waves, see Figure 1. To address the inverse problem, hawen solves the minimization problem: . (1) The misfit function J is defined to evaluate a distance function comparing the observed data d with simulations of wave propagation restricted to the position of the measurements: F(m), where m represents the physical properties. The iterative minimization scheme successively updates the physical properties used for the simulations (Figure 1b), and follows a Newtontype method (Nocedal & Wright, 2006). In the context of seismic imaging, it is the Full Waveform Inversion method, cf. Virieux & Operto (2009).
hawen offers several options to conduct the iterative minimization, such as the choice of misfit function and method to conduct the minimization. These are further listed and detailed in the software documentation (Faucher, 2020).