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

HAL
open science

# Accelerated iterative DG finite element solvers for large-scale time-harmonic acoustic problems

Axel Modave

► To cite this version:

Axel Modave. Accelerated iterative DG finite element solvers for large-scale time-harmonic acoustic problems. INTER.NOISE 2024 - 53rd International Congress and Exposition on Noise Control Engineering, Aug 2024, Nantes, France. 10.3397/IN_2024_2877. hal-04739754

HAL Id: hal-04739754

https://hal.science/hal-04739754v1

Submitted on 16 Oct 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Accelerated iterative DG finite element solvers for large-scale
time-harmonic acoustic problems*

A. Modave¹

¹POEMS, CNRS, Inria, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France,
axel.modave@ensta-paris.fr

Abstract

Finite element methods are widely used to solve time-harmonic wave propagation prob-
lems, but solving large cases can be extremely difficult even with the computational power
of parallel computers. In this work, the linear system resulting from the finite element
discretization is solved with iterative solution methods, which are efficient in parallel but
can require a large number of iterations. In standard discontinuous Galerkin (DG) meth-
ods, the numerical solution is discontinuous at the interfaces between the elements. In
hybridizable DG methods, additional unknowns are introduced at the interfaces between
the finite elements, and the physical unknowns are eliminated from the global system, re-
sulting in a hybridized system. We have recently proposed a new strategy, called CHDG,
where the additional unknowns correspond to transmission variables, whereas in the stan-
dard approach they are numerical fluxes. This strategy improves the properties of the
hybridized system for faster iterative solution procedures. In this talk, we present and
study a 3D CHDG implementation with nodal finite element basis functions. The result-
ing scheme has properties amenable to efficient parallel computing. Numerical results are
presented to validate the method, and preliminary 3D computational results are proposed.

# 1 Introduction

Numerical schemes based on finite element methods (FEMs) have proven their ability to
handle realistic time-harmonic acoustic problems. Unfortunately, it is very difficult to develop
FEM-based computational tools that are both fast and reliable. It is even more complicated
for high-frequency cases, where the wavelength is small compared to the characteristic size of
the problem, requiring fine meshes.

In this work, we are interested in solving time-harmonic problems with iterative solution
procedures that can run efficiently on parallel computers. However, standard iterative schemes
are generally slow at solving the linear systems resulting from the finite element discretization
of time-harmonic problems. Many iterations are generally required to obtain acceptable
solutions. We have recently proposed a new finite element method, called CHDG [1], where
the linear system has improved properties for efficient iterative solution procedures. We think

<footnote>*Published in *INTER-NOISE and NOISE-CON Congress and Conference Proceedings*, INTER-NOISE24, Nantes, France. Distributed under **Creative Commons CC-BY 4.0** license.</footnote>
that this method is a good candidate to address large-scale realistic problems on large parallel computers.

The CHDG method is based on a standard discontinuous Galerkin (DG) finite element scheme. An additional variable is introduced on the faces of the elements, and the physical variables are eliminated following the hybridization procedure of the hybridizable DG methods, see e.g. [2]. The additional variable corresponds to a numerical flux in standard HDG methods, and it corresponds to a transmission variable in the CHDG method.

## 2 Numerical methods

To present our approach, we consider the general time-harmonic acoustic problem

$$ \left\{ \begin{array}{ll} -\iota\kappa p + \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega, \\ -\iota\kappa \mathbf{u} + \nabla p = \mathbf{0} & \text{in } \Omega, \\ p - \mathbf{n} \cdot \mathbf{u} = s & \text{on } \partial\Omega, \end{array} \right. \qquad (2.1) $$

where $\Omega$ is the computational domain, $\partial\Omega$ is the domain boundary, $p(\mathbf{x})$ and $\mathbf{u}(\mathbf{x})$ are the pressure and velocity fields (respectively), $\kappa$ is the (real positive) wavenumber defined on $\partial\Omega$, $\mathbf{n}(\mathbf{x})$ is the outward unit normal, and $s(\mathbf{x})$ is a surface source. Other boundary conditions and source terms can easily be considered, see [1].

### 2.1 Discontinuous Galerkin finite element scheme with upwind fluxes

We consider a mesh of the domain $\Omega$ made of tetrahedral elements. The elements and the faces are denoted $K$ and $F$, respectively.

On each element $K$, the physical fields are approximated by complex polynomial functions, i.e. $q_K \in \mathbb{P}_p(K)$ and $v_K \in [\mathbb{P}_p(K)]^3$, where $p$ is the maximum polynomial degree. The finite element scheme is based on a standard discontinuous Galerkin method described in [3]. With this scheme, the physical fields verify the following local element-wise problem.

$$ \left\{ \begin{array}{ll} \text{Find } p_K \in \mathbb{P}_p(K) \text{ and } \mathbf{u}_K \in [\mathbb{P}_p(K)]^3 \text{ such that} \\ \begin{cases} -\iota\kappa(p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_F \langle \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F, q_K \rangle_F = 0, \\ -\iota\kappa(\mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K + \sum_F \langle \hat{p}_F, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F = 0, \end{cases} & (2.2) \\ \text{for all } q_K \in \mathbb{P}_p(K) \text{ and } \mathbf{v}_K \in [\mathbb{P}_p(K)]^3, \end{array} \right. $$

where $\mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F$ and $\hat{p}_F$ are numerical fluxes. We use the notation $(x, y)_K = \int_K u\bar{y} \, d\mathbf{x}$ and $\langle x, y \rangle_F = \int_F u\bar{y} \, d\mathbf{x}$, and the sum $\sum_F$ is written for each face $F$ of $K$. The global problem is composed of all local problems (2.2).

The numerical fluxes are used to impose the boundary condition and weakly enforce the
continuity of the fields at the interfaces between elements. We use upwind fluxes defined as

$$ \hat{p}_F := \begin{cases} \frac{1}{2}(p_K + p_{K'}) + \frac{1}{2} \mathbf{n}_{K,F} \cdot (\mathbf{u}_K - \mathbf{u}_{K'}) & \text{if } F \not\subset \partial\Omega, \\ (p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K + s)/2 & \text{if } F \subset \partial\Omega, \end{cases} \\ \mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F := \begin{cases} \frac{1}{2} \mathbf{n}_{K,F} \cdot (\mathbf{u}_K + \mathbf{u}_{K'}) + \frac{1}{2} (p_K - p_{K'}) & \text{if } F \not\subset \partial\Omega, \\ (p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K - s)/2 & \text{if } F \subset \partial\Omega, \end{cases} $$

where $\mathbf{n}_{K,F}$ is the unit outward normal to $K$ on face $F$. If $F$ is an interior face (i.e. $F \not\subset \partial\Omega$), then $K'$ is the neighboring element sharing $F$ with $K$.

## 2.2 Hybridization with transmission variables

Instead of directly solving the global problem, we introduce an additional variable on all the faces of the mesh, and we eliminate the physical variables. Then, we solve the resulting *hybridized problem*, where the unknown corresponds to the additional variable. The original physical variables can be recovered in a local post-processing step. The goal is to choose the additional variable in such a way that the hybridized problem has good properties for fast iterative solution methods.

For each face $F$ of each element $K$, we introduce the *outgoing and incoming transmission variables*, denoted $g_{K,F}^{\oplus}$ and $g_{K,F}^{\ominus}$, respectively. They are defined as

$$ g_{K,F}^{\oplus} := p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K, \\ g_{K,F}^{\ominus} := \begin{cases} g_{K',F}^{\oplus} & \text{if } F \not\subset \partial\Omega, \\ s_R & \text{if } F \subset \partial\Omega, \end{cases} $$

where $K'$ is the neighbor of $K$ at face $F$. We have $\hat{p}_F = (g_{K,F}^{\oplus} + g_{K,F}^{\ominus})/2$ and $\mathbf{n}_{K,F} \cdot \hat{\mathbf{u}}_F = (g_{K,F}^{\oplus} - g_{K,F}^{\ominus})/2$. The transmission variables can be interpreted as outgoing/incoming data traveling across the face.

In the CHDG method, the incoming transmission variable (on each face of each element) is defined as an additional unknown of the problem. The global problem is then composed of local problems defined on all the elements, and additional equations defined on all the faces. For each element $K$, the local element-wise problem is

$$ \left| \begin{array}{l} \text{Find } p_K \in \mathbb{P}_p(K) \text{ and } \mathbf{u}_K \in [\mathbb{P}_p(K)]^3 \text{ such that} \\ \left\{ \begin{array}{l} -i\kappa(p_K, q_K)_K - (\mathbf{u}_K, \nabla q_K)_K + \sum_F \frac{1}{2} \langle p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K - g_{K,F}^{\ominus}, q_K \rangle_F = 0, \\ -i\kappa(\mathbf{u}_K, \mathbf{v}_K)_K - (p_K, \nabla \cdot \mathbf{v}_K)_K + \sum_F \frac{1}{2} \langle p_K + \mathbf{n}_{K,F} \cdot \mathbf{u}_K + g_{K,F}^{\oplus}, \mathbf{n}_{K,F} \cdot \mathbf{v}_K \rangle_F = 0, \end{array} \right. \\ \text{for all } q_K \in \mathbb{P}_p(K) \text{ and } \mathbf{v}_K \in [\mathbb{P}_p(K)]^3, \end{array} \right. \qquad (2.3) $$

For each face $F$ of each element $K$, the additional equations corresponds to the problem

$$ \left| \begin{array}{l} \text{Find } g_{K,F}^{\ominus} \in \mathbb{P}_p(F) \text{ such that} \\ \left\{ \begin{array}{l} \langle g_{K,F}^{\ominus}, \xi_{K,F} \rangle_F - \langle p_{K'} + \mathbf{n}_{K',F} \cdot \mathbf{u}_{K'}, \xi_{K,F} \rangle_F = 0, & \text{if } F \not\subset \partial\Omega, \\ \langle g_{K,F}^{\ominus}, \xi_{K,F} \rangle_F = \langle s, \xi_{K,F} \rangle_F, & \text{if } F \subset \partial\Omega, \end{array} \right. \\ \text{for all } \xi_{K,F} \in \mathbb{P}_p(F). \end{array} \right. \qquad (2.4) $$
After eliminating all physical variables, we obtain a hybridized problem where the incoming transmission variable is the only unknown. We introduce the field $g_h^{\ominus}$ with $g_h^{\ominus}|_{K,F} := g_{K,F}^{\ominus}$. The hybridized problem can be written as

$$ \begin{gathered} \text{Find } g_h^{\ominus} \in \bigcup_K \bigcup_F \mathbb{P}_p(F) \text{ such that} \\ \langle g_h^{\ominus}, \xi_h \rangle \partial \tau_h - \langle \Pi S g_h^{\ominus}, \xi_h \rangle \partial \tau_h = \langle b, \xi_h \rangle \partial \tau_h, \tag{2.5} \\ \text{for all } \xi_h \in \bigcup_K \bigcup_F \mathbb{P}_p(F), \end{gathered} $$

with $\langle \cdot, \cdot \rangle \partial \tau_h := \sum_K \sum_F \langle \cdot, \cdot \rangle_F$, an exchange operator $\Pi$, a scattering operator $S$, and a global source term $b$ depending on $s$. The operator $S$ involves the solution of the local problem (2.3) for all the elements. The operator $\Pi$ involves the transmission of data between neighboring elements. We refer to [1] for more details.

## 2.3 Nodal finite element discretization

The numerical fields are represented with nodal basis functions. For each face $F$ of each element $K$, we have

$$ g_{K,F}^{\ominus}(\boldsymbol{x}) = \sum_{n=1}^{N_{fp}} g_{K,F,n}^{\ominus} \ell_{K,F,n}(\boldsymbol{x}), \tag{2.6} $$

where $\{\ell_{K,F,n}\}_n$ are the Lagrange functions associated to finite element nodes defined on $F$, and $\{g_{K,F,n}^{\ominus}\}_n$ are the nodal values of $g_{K,F}^{\ominus}$. See [3] for more details on nodal finite elements.

Using the discrete representation (2.6) and the Lagrange functions as test functions in the hybridized problem (2.5) leads to the linear system

$$ [\mathbf{M} - \Pi \mathbf{S}] \mathbf{g} = \mathbf{b}, \tag{2.7} $$

where $\mathbf{g}$ contains the nodal values of $g_h^{\ominus}$ on the entire mesh, $\mathbf{M}$ contains the mass matrices associated to the faces, and $\Pi$, $S$ and $b$ are discrete versions of $\Pi$, $S$ and $b$. The matrices $\mathbf{M}$ and $S$ are block diagonal, where the size of each block depends on the number of nodes per face. The matrix $\Pi$ is sparse. The application of the matrix $S$ requires the solution of the discrete version of the local element-wise problem (2.3) for all the elements. For these local problems, the physical fields are also discretized with nodal basis functions.

The finite element matrices involved in the local and hybridized problems are computed by using strategies described in [3]. These strategies use basis transformations between the nodal basis functions and orthonormal basis functions on each element.

## 2.4 Iterative procedures

The discrete problem (2.7) can be solved with the block Jacobi Over-Relaxation (JOR) iterative method, i.e.

$$ \mathbf{g}^{\ell+1} = \mathbf{g}^{\ell} + \alpha \mathbf{M}^{-1} (\mathbf{b} - [\mathbf{I} + \Pi \mathbf{S}] \mathbf{g}^{\ell}) $$

with the relaxation parameter $\alpha \in ]0, 1]$. Because $\Pi S$ is a strict contraction, this scheme converges even without relaxation, i.e. $\alpha = 1$.
Krylov iterative methods are widely used to solve Helmholtz problems. In this work, we have implemented the CGNR method, which consists in a conjugate gradient method applied to the normal equation $A^*Ax = A^*b$ for a given system $Ax = b$. In [1], we observe that CGNR is almost as efficient as GMRES for solving the CHDG system with a set of 2D benchmarks. Here, we have tested left and symmetric preconditioning with $M$. It turns out that the second case can be interpreted as rewriting the system in the orthonormal polynomial basis.

# 3 Implementation and preliminary results

## 3.1 Programming and computational aspects

The CHDG method and the iterative solution procedures have been implemented in a dedicated C++ code. The OpenMP library is used for shared memory parallel computing, and mkl is used for solving the local linear systems.

The local problems are constructed and solved in parallel at each iteration. The local matrices are computed by using strategies similar to those described in [3]. From an algorithmic point of view, data exchange between elements is similar to that in parallel DG time-domain solvers, see e.g. [4].

## 3.2 Validation results

In order to validate the code, we consider a free space benchmark (with a non-homogeneous Robin boundary condition and $\kappa = 10\pi$) and a cavity benchmark (with a homogeneous Dirichlet boundary condition, unit volume source term and $\kappa = 5.5\sqrt{2\pi}$). In both cases, the computational domain is $[0, 1]^3$, the characteristic size of an element is $h = 1/6$ and the polynomial degree is $p = 6$. The reference solutions are represented on Figure 1.

Figure 2 shows the decay of the relative error as a function of the number of iterations for several iterative methods and both benchmarks. The error corresponds to the difference between the analytic reference solution and the numerical solution obtained at each iteration. Therefore, this error should converge to the finite element error that would be obtained with a direct solver. The best approximation error, which corresponds to the lowest error achievable with the finite element mesh, is shown with dashed lines.

We observe that the Block Jacobi iterative method (with and without relaxation) does not perform well for the cavity benchmark, and it is efficient for the free-space benchmark. This is in agreement with the 2D results obtained in [1]. In contrast, the CGNR method with symmetric preconditioning is efficient in both cases. This is the most robust combination.

## 3.3 A first computational result

We consider the scattering of the plane wave $p_{\text{inc}}(x) = e^{i\kappa x}$ by a sphere (radius 1) in a spherical computational domain (radius 1.5). Both spheres are centered at the origin. The Dirichlet boundary condition $p(x) = -p_{\text{inc}}(x)$ is prescribed at the boundary of the inner sphere, and an absorbing boundary condition is used at the boundary of the outer sphere. The solution
<figure><img src="image_3.png" /><figcaption>Figure 1: Validation benchmarks. Real part of the reference solutions. Surfaces of isovalues are shown inside the mesh of the boundary.</figcaption></figure>

corresponds to the acoustic scattered field. The wavenumber is $κ = 10π$, the characteristic size of an element is $h = 0.1$, and the polynomial degree is $p = 4$. The mesh generation and the post-processing have been done with **gmsh** [5].

The mesh is composed of 49.174 tetrahedral elements. The global problem corresponds to 6.884.360 degrees of freedom (DOFs), with 1.721.090 DOFs for the pressure field. The hybridized problem has 2.950.440 DOFs.

Snapshots of the numerical solution obtained during the iterations of the CGNR procedure (with symmetric preconditioning) are shown in Figure 3. The quality of the solution obtained after 100 iterations is already good. The runtime per iteration is about 22 seconds on an Intel Xeon CPU Gold 6230 20-core (2.1 GHz) processor. OpenMP was used with 20 threads.

In future work, we plan to extend the code to distributed memory parallel computing with MPI and to tackle more challenging cases.
<figure><img src="image_4.png" /><figcaption>Figure 2: Validation benchmarks. Error on $u(x)$ in $L^2$-norm vs iteration with different iterative schemes. The numerical solution is compared to the analytic reference solution. The dashed lines correspond to the best approximation error.</figcaption></figure>

## Acknowledgements

This work was partially supported by the ANR JCJC project WavesDG (ANR-21-CE46-0010) and by AID via CIEDS (project 2022 – *ElectroMath*). The author would like to thank the IDCS research facility of École Polytechnique for providing access to the Cholesky cluster.

## References

[1] Axel Modave and Théophile Chaumont-Frelet. A hybridizable discontinuous Galerkin method with characteristic variables for Helmholtz problems. *Journal of Computational Physics*, 493:112459, 2023.

[2] Bernardo Cockburn. Static condensation, hybridization, and the devising of the HDG methods. *Building bridges: connections and challenges in modern approaches to numerical partial differential equations*, pages 129–177, 2016.

[3] Jan S Hesthaven and Tim Warburton. *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[4] Axel Modave, Amik St-Cyr, and Tim Warburton. GPU performance analysis of a nodal discontinuous Galerkin method for acoustic and elastic models. *Computers & Geosciences*, 91:64–76, 2016.

[5] Christophe Geuzaine and Jean-François Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. *International journal for numerical methods in engineering*, 79(11):1309–1331, 2009.
<figure><img src="image_5.png" /><figcaption>Figure 3: Scattering benchmark. Numerical solution during iterations using the CGNR procedure. The solution corresponds to the real part of the pressure field. Surfaces of isovalues are shown inside the mesh of the boundary. Only half of the spheres are shown.</figcaption></figure>