# Hybridizable Discontinuous Galerkin method for solving anisotropic wave propagation

Hélène Barucq, Ha Pham, Florian Faucher

► To cite this version:

Hélène Barucq, Ha Pham, Florian Faucher. Hybridizable Discontinuous Galerkin method for solving anisotropic wave propagation. WAVES 2024 - 16th International Conference on Mathematical and Numerical Aspects of Wave Propagation, Jun 2024, Berlin, Germany. pp.361-362. hal-04898997

HAL Id: hal-04898997

https://inria.hal.science/hal-04898997v1

Submitted on 20 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 - Public Domain Dedication 4.0 International License
Hybridizable Discontinuous Galerkin
method for solving anisotropic wave propagation

Hélène Barucq¹,*, Ha Pham¹, Florian Faucher¹

¹Team Makutu, Inria, University of Pau and Pays de l'Adour, TotalEnergies, CNRS, France.

*Email: helene.barucq@inria.fr

Abstract

In this work, we solve the time-harmonic linear anisotropic elastic equation with Hybridized Discontinuous Galerkin (HDG) method, which is particularly effective in exploiting efficiently advanced and parallel computing architecture to perform large-scale simulations. Following the original HDG method devised by Cockburn and co-authors for second-order elliptic equation, e.g., [2], we propose a formulation based on Voigt notation to easily include general anisotropy and facilitate efficient code implementation. We use a numerical flux for anisotropy as an extension of hybridization technique of upwind flux for isotropy, the so-called Godunov flux in [3]. Numerical experiments in [3] demonstrated the efficiency of this stabilization operator in the definition of numerical trace, for varying degrees of anisotropy. Herein we focus on the parameters representation and the interest of using compliance-based formulation to take advantage on high-order approximation.

1 Introduction

We consider the HDG method to solve the linear anisotropic elastic equation in frequency domain. HDG method retains the attractive features of the Discontinuous Galerkin family, such as being locally conservative, flexible with h-p adaptivity strategy and thus unstructured non-conforming meshes, optimal convergence with non-conforming bases, and highly parallelizable. In our implementation, the advantages of HDG method are exploited further by employing a first-order formulation in $(\mathbf{u}, \boldsymbol{\sigma})$ with constitutive law written in terms of the compliance tensor. This allows for the variation of the model parameters within a mesh cell, and as shown in our numerical investigation, leads to better representation of complex media and more numerically stable resolution for these media. Moreover, the discretized problem is compactly described using Voigt representation, which maximizes matrix operation and provides efficient

book-keeping of physical parameters for any generic anisotropic material. In our work [3], we addressed the question of defining a Godunov-like stabilization matrix $\boldsymbol{\tau}_{\mathbf{u}}$ involved in the transmission conditions between elements. Most of HDG methods work with identity-based stabilization, i.e., choosing $\boldsymbol{\tau}_{\mathbf{u}} = \tau I_d$ with $I_d$ the identity matrix and $\tau$ a scalar scaling factor. Other works investigated $\boldsymbol{\tau}_{\mathbf{u}} = \tau \boldsymbol{\Gamma}$, where $\boldsymbol{\Gamma}$ is the Kelvin-Christoffel (KC) matrix for elastostatics. We extended the hybridization technique of upwind flux for isotropy to anisotropy under the assumption of three distinct wave speeds. The construction is based on exact solutions of the Riemann problem and their Rankine-Huguniot jump conditions. We showed that the Godunov stabilization as derived is optimal, while the other two families require correct scaling factors which might be not universal to match the performance of Godunov stabilization. In this work, we focus on the advantage of writing the elastic system in terms of the compliance matrix $\boldsymbol{S}$ rather than the stiffness tensor $\boldsymbol{C}$ to show how it facilitates non-constant description of background materials within a mesh cell.

2 Problem setting

We consider the propagation of time-harmonic waves in linear elasticity, given in terms of the displacement $\mathbf{u}$. Considering a domain $\Omega$ and an interior source $\mathbf{f}$, the elastic wave equations are given by:

$$-\omega^2 \rho \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{f}, \quad (1a)$$

$$\boldsymbol{S}\boldsymbol{\sigma} = \nabla^5 \mathbf{u}. \quad (1b)$$

where $\boldsymbol{\sigma}$ is the stress tensor, $\rho$ is the density, $\boldsymbol{S}$ is the compliance tensor which is the inverse of the stiffness tensor $\boldsymbol{S} := \boldsymbol{C}^{-1}$. This formulation allows us to easily treat physical properties that vary within the cell, avoiding having to derive these properties in the variational formulation. By this way, we represent varying properties on large-cell discretization, hence keeping high-order polynomial discretization, which is critical
for the efficiency of the HDG method compared to CG one. Moreover, a flexible representation of the physical properties on the discretized domain is a key-ingredient for inversion.

### 3 Discrete problem in HDG method

For a mesh with $|T|$ cells $K_e$, the HDG discretization is written for local volume unknown $\mathbf{W}^e$ and global numerical trace $\mathbf{\Lambda}$ as, cf. [3],

$$ \left\{ \begin{array}{l} \mathbf{A}^e \mathbf{W}^e + \mathcal{D}^e \mathcal{R}_e \mathbf{\Lambda} = \mathbf{S}^e, \forall e, \\ \displaystyle\sum_{e=1}^{|T|} \mathcal{R}_e^t (\mathcal{B}^e \mathbf{W}^e + \mathcal{L}^e \mathcal{R}_e \mathbf{\Lambda}) = \displaystyle\sum_{e=1}^{|T|} \mathcal{R}_e^t \mathbf{s}^e. \end{array} \right. $$

The problem can be reduced in terms of $\mathbf{\Lambda}$:

$$ \mathbb{K} \mathbf{\Lambda} = \mathcal{S}, \text{ with } \mathcal{S} := \sum_{e=1}^{|T|} \mathcal{R}_e^t (\mathbf{s}^e - \mathcal{B}_e \mathcal{A}_e^{-1} \mathbf{S}^e), $$

$$ \text{and } \mathbb{K} := \sum_{e=1}^{|T_h|} \mathcal{R}_e^t \mathbb{K}^e \mathcal{R}_e, \qquad (3a) $$

$$ \text{where } \mathbb{K}^e := -\mathcal{B}_e (\mathcal{A}^e)^{-1} \mathcal{D}^e + \mathcal{L}_e. \qquad (3b) $$

One solves this global problem in terms of the trace $\mathbf{\Lambda}$; the values of the volume unknowns $\mathbf{W}^e$ are retrieved element-by-element as follows:

$$ \mathbf{W}^e = -(\mathcal{A}^e)^{-1} (\mathcal{D}^e \mathcal{R}_e \mathbf{\Lambda} + \mathbf{S}^e), \forall e. $$

### 4 Numerical experiment

We consider a background model inspired from the standard model S in helioseismology, [1]. It provides us with radially symmetric density and wave speed, employed here as P-wave speed. The domain is a disk of radius 1. The physical parameters are very challenging because the density is exponentially decaying near the boundary, and the wave speed varies sharply. Hence, a piecewise-constant approximation requires mesh refinement near the boundary. Using the discretization with the compliance tensor, we instead represent the physical parameters using a basis of Lagrange functions on each cell. In Figure 1, we compare the density model near surface when represented with piecewise-constant or with Lagrange basis of order 3. The piecewise-constant representation loses the spherical nature of the medium, and the Lagrange polynomial representation preserves the spherical nature of the model. In Figure 2, we see that the piecewise-constant representation leads to noisy wavefields. This illustrates the interest in using the formulation with compliance assisted with a Lagrange polynomial representation.

<figure><img src="image_2.png" /><figcaption>Figure 1: Solar-like density represented on the mesh near surface. The piecewise-constant representation (top) is unable to preserve the radial symmetry, contrary to the piecewise polynomial representation (bottom).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Density-scaled wavefields $\sqrt{\rho} \text{Re}(u_z)$ at 4 mHz using solar-like models with piecewise polynomial representation of the models (left) or piecewise constant (right).</figcaption></figure>

### References

[1] J. Christensen-Dalsgaard, et al., The current state of solar modeling, Science, 1996.

[2] B. Cockburn, J. Gopalakrishnan, and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Num. An., 2009.

[3] H. Pham, F. Faucher and H. Barucq, On the implementation of Hybridizable Discontinuous Galerkin discretization for linear anisotropic elastic wave equation: Voigt-notation and stabilization, Research Report hal-04356602 RR-9533, 2023.