# A quasi-Trefftz method for time-harmonic electromagnetism based on a Flux Reconstruction local solver

Sébastien Pernet, Matthias Rivet, Sébastien Tordeux

► To cite this version:

Sébastien Pernet, Matthias Rivet, Sébastien Tordeux. A quasi-Trefftz method for time-harmonic electromagnetism based on a Flux Reconstruction local solver. NUMELEC 2024, Jul 2024, Toulouse, France. hal-04843838

HAL Id: hal-04843838

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

Submitted on 17 Dec 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A quasi-Trefftz method for time-harmonic electromagnetism based on a Flux Reconstruction local solver

Sébastien Pernet¹, Matthias Rivet¹,², Sébastien Tordeux²

¹ DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France

² EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142

Email: sebastien.pernet@onera.fr, matthias.rivet@onera.fr, sebastien.tordeux@inria.fr

**Résumé**—The simulation of time-harmonic electromagnetic waves in wide domains (in terms of wavelengths) with classic methods, as Discontinuous Galerkin, induces large linear systems: their direct resolution becomes memory-prohibitive while the matrices are poorly adapted to iterative methods. On the contrary, Trefftz methods embody a framework which is able to deal with such configurations. Yet, the classic choice of plane waves as basis functions limits their efficiency and usefulness in heterogeneous or challenging domains. Thus, a quasi-Trefftz approach can be considered by using a Flux Reconstruction method to design an adaptable and robust Trefftz space.

## I. MODEL PROBLEM

For a linear isotropic material which occupies a domain $\Omega \subset \mathbb{R}^3$, the time-harmonic Maxwell's equations in absence of charges correspond to

$$i\kappa \begin{pmatrix} \mathbf{E} \\ \mathbf{H} \end{pmatrix} + \begin{pmatrix} -\nabla \times \mathbf{H} \\ \nabla \times \mathbf{E} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \quad \text{in } \Omega, \quad (1)$$

where $\kappa > 0$ stands for the wavenumber. This can be rewritten in a hyperbolic setting thanks to the flux operators $\mathbf{F}^j \in \mathbb{R}^{6 \times 6}$, for $1 \le j \le 3$, as

$$i\kappa\mathbb{E} + \sum_{j=1}^{3} \frac{\partial\phi^j}{\partial x_j} = \mathbf{0} \quad \text{in } \Omega \text{ with } \phi^j = \mathbf{F}^j\mathbb{E}, \quad (2)$$

where $\mathbb{E} = (\mathbf{E}, \mathbf{H})^t$ and $\phi^j$ denote the electromagnetic field and the fluxes respectively. This equation is completed by an impedance Boundary Condition (BC)

$$\gamma_t[\mathbf{n}_{\partial\Omega}]\mathbf{E} + Z_{\partial\Omega}\gamma_{\times}[\mathbf{n}_{\partial\Omega}]\mathbf{H} = \mathbf{g} \quad \text{on } \partial\Omega, \quad (3)$$

defined on the boundary of $\Omega$, denoted as $\partial\Omega$, for an impedance $Z_{\partial\Omega}: \partial\Omega \to \mathbb{C}^3$ and a tangential field $g \in \mathbf{L}_t^2(\partial\Omega)$. Finally, $\gamma_t[\mathbf{n}_{\partial\Omega}] = (\mathbf{n}_{\partial\Omega} \times \cdot) \times \mathbf{n}_{\partial\Omega}$ and $\gamma_{\times}[\mathbf{n}_{\partial\Omega}] = \mathbf{n}_{\partial\Omega} \times \cdot$ stand for the tangential component and trace respectively, with $\mathbf{n}_{\partial\Omega}$ the outgoing normal from $\Omega$.

We introduce a mesh $\mathcal{T}_h$ of the domain $\Omega$. For any cell $T \in \mathcal{T}_h$, we denote the restriction of the electromagnetic field as $\mathbb{E}^T = (\mathbb{E}^T, \mathbf{H}^T) := \mathbb{E}|_T$ whose tangential electric component and magnetic trace are defined as $\gamma_t \mathbb{E}^T = \gamma_t[\mathbf{n}_T] \mathbb{E}^T$ and $\gamma_{\times}^T \mathbf{H}^T = \gamma_{\times}[\mathbf{n}_T] \mathbf{H}^T$ respectively, where $\mathbf{n}_T$ stands for the outgoing normal of $T$.

## II. QUASI-TREFFTZ METHOD

The Trefftz method relies on a reciprocity formula which encompasses the energy conservation. Technically, it can be described as a DG approach for which the Galerkin space

$\mathbb{X} := \prod_{T \in \mathcal{T}_h} \mathbb{X}_T$ is made of local solutions (i.e. in each mesh cell $T \in \mathcal{T}_h$) of the Maxwell's equations. It then writes [6]: Find $\mathbb{E} = (\mathbf{E}, \mathbf{H}) \in \mathbb{X}$ such that $\forall \mathbb{E}' = (\mathbf{E}', \mathbf{H}') \in \mathbb{X}$,

$$\sum_{T \in \mathcal{T}_h} \int_{\partial T} \widehat{\gamma_{\times}^T \mathbf{H}^T} \cdot \overline{\gamma_t \mathbb{E}^T} + \widehat{\gamma_t \mathbb{E}^T} \cdot \overline{\gamma_{\times}^T \mathbf{H}^T} = \ell(\mathbb{E}') \quad (4)$$

where $\widehat{\cdot}$ denotes a numerical trace computation (thanks to an upwinding mechanism or a Riemann solver for example) and $\ell$ stands for a linear form which includes the BCs.

In particular, this approach limits the amount of degrees of freedom, by being only defined on the faces, and is intrinsically adaptable to an iterative resolution through an inherent contraction property [1]. Moreover, the use of physical basis functions, such as plane waves, tends to limit the numerical pollution effects [5].

Yet, the classic plane wave choice is used to leading to numerical dependence phenomena which impact the conditioning properties of the matrix to inverse [2], and prevents adaptability of the basis to local mesh properties (as heterogeneous characteristics of the material or expected singularities which will be hardly reproduced in this case). To face such limitations, we consider a quasi-Trefftz approach [3] : the idea is then to parameterise the space of local solutions by a polynomial approximation of an impedance trace space.

## III. FLUX RECONSTRUCTION LOCAL SOLVER

The choice of the local solver is then of special importance, and we chose to adapt the Flux Reconstruction method (FR) [4], which has already established itself for hyperbolic problems in the Computational Fluid Dynamics community, to this time-harmonic framework for a Cartesian mesh $\mathcal{T}_h^{FR}$. It relies on the strong form of the hyperbolic problem (2) by approximating the solution $\mathbb{E}$ and the fluxes $\phi^j$ by piecewise polynomial functions $\mathbb{E}_h$ and $\tilde{\phi}_h^j$ of degrees $k := (k_1, k_2, k_3)$ and $k^{j,1} = (k_{j'} + \delta_{j,j'})_{1 \le j' \le 3}$. For any mesh cell $T \in \mathcal{T}_h^{FR}$, this leads to

$$i\kappa\mathbb{E}_h + \sum_{j=1}^{3} \frac{\partial\tilde{\phi}_h^j}{\partial x_j} = \mathbf{0} \quad \text{in } T. \quad (5)$$

To take BCs and flux continuity properties into account, the flux variables $\tilde{\phi}_h^j$ result from a polynomial correction of $\phi_h = \mathbf{F}^j \mathbb{E}_h$ which ensures it fits with a given polynomial numerical trace of degree $(k_{j'})_{j' \ne j}$, denoted $\gamma_F \phi_h^j$, on any face $F$ orthogonal to $e_j$.
For example, for any $\mathbf{x} = (x_1, x_2, x_3)$ of a cell $T$ whose 'left' and 'right' boundaries (in the $x_1$ direction) are denoted as $F^+$ and $F^-$ respectively, such a correction takes the form of

$$ \tilde{\phi}_{\mathbf{h}}^{1} = \phi_{\mathbf{h}}^{1} + \delta_{\mathbf{h},F^{+}}^{1,\rightarrow} + \delta_{\mathbf{h},F^{-}}^{1,\leftarrow}, \quad (6) $$

with

$$ \delta_{\mathbf{h},F^{+}}^{1,\rightarrow}(\mathbf{x}) = (\gamma_{F^{+}}\tilde{\phi}_{\mathbf{h}}^{1} - \phi_{\mathbf{h}}^{1}|_{F^{+}})(x_2, x_3) P_{T}^{1,\rightarrow}(x_1), \quad (7) $$

$$ \delta_{\mathbf{h},F^{-}}^{1,\leftarrow}(\mathbf{x}) = (\gamma_{F^{-}}\tilde{\phi}_{\mathbf{h}}^{1} - \phi_{\mathbf{h}}^{1}|_{F^{-}})(x_2, x_3) P_{T}^{1,\leftarrow}(x_1), \quad (8) $$

where $P_T^{1,\rightarrow}$ and $P_T^{1,\leftarrow}$ are flux correction polynomial functions of degree $k_1 + 1$.

This method allows to retrieve classic methods such as nodal Discontinuous Galerkin (DG) and Spectral Difference (SD) in usually less expensive ways, for specific flux correction polynomial functions [4]. Moreover, this framework intrinsically ensures that the polynomial numerical solution strongly solves the hyperbolic system (5).

## IV. NUMERICAL ANALYSIS OF THE FR LOCAL SOLVER

The FR method relies on a choice of correction polynomial functions $P^\rightarrow$ and $P^\leftarrow$, and one may classically consider Radau polynomials or Lagrange polynomials associated to Gauss points (which respectively allow to retrieve a nodal DG and a SD formulation [4]), that we denote as FR_Radau and SD_IG respectively.

Yet, we derived an a priori error estimate of the FR method [7], with explicit dependence with respect to the discretisation parameters, for the 1D Maxwell's equations with incoming BCs (i.e. $Z_{\partial\Omega} = 1$ for $g_1, g_2 \in \mathbb{C}$ on the left and right boundaries respectively) for a uniform mesh of step $h$ of the domain $\Omega = [0, L]$, and symmetric correction polynomials :

$$ \|\mathbb{E} - \mathbb{E}_{\mathbf{h}}\|_{L^2(\Omega)} \lesssim C^{\rightarrow}(\kappa, L, h, P^{\rightarrow}) (|g_1|^2 + |g_2|^2), \quad (9) $$

where $A \lesssim B$ means that there exists a constant $C > 0$, independent of $\kappa, L, h$ and $P^\rightarrow$ such that $A \le BC$. In particular, we obtained the following quasi-optimal behaviour in the asymptotic regime : $C^{\rightarrow}(\kappa, L, h, P^\rightarrow) = O(h^{2k+2})$.

This allowed, for a given discretisation, to look for the correction polynomial function $P^\rightarrow$ which minimises this bounding. We refer to this choice as FR_opti.

The influence of such choices on the h-convergence of the FR method are illustrated in Figure 1, for a domain $\Omega = [0, 3]^3$ with uniform mesh and polynomial degree $k, \kappa = 2\pi$, random positive values of $Z_{\partial\Omega}$ (taken as constant on each face) and an exact solution made of a sum of 5 random planewaves. This highlights the good properties of a high-order accurate method and the optimised correction polynomials show interesting improvements in comparison with the classic choices.

## V. QUASI-TREFFTZ NUMERICAL VALIDATION

In this presentation, we will finally focus on the use of the FR method as a local solver in a quasi-Trefftz approach for the 3D Maxwell's equations. For example, let's consider a domain $\Omega = [0, 1]^3$ uniformly meshed with 5 cubic macrocells in each direction, that we denote as $\mathcal{T}_h$. For a face $F$ of

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Illustration of the h-convergence of the FR method according to the number of degrees of freedom per wavelength $N_{ddl/\lambda}$: evolution of the relative $L^2$ error w.r.t. the norm of the exact solution (Fig. 1a) and FR_Radau error (Fig. 1b), for SD_IG (—*—), FR_Radau (—*—) and FR_opti (—*—), with polynomial degrees $k = 2$ (—*—) and $k = 5$ (—*—). The quasi-optimal slopes are denoted as (---).</figcaption></figure>

a given $T \in \mathcal{T}_h$, we introduce a set of 9 linearly independent polynomial face functions ($f_j^F$)$_{1 \le j \le 9}$: $F \to \mathbb{R}$ and the induced BCs

$$ g_{t,j}^{F} : \partial T \to \mathbb{C}^{3} \text{ with } g_{t,j}^{F} = \begin{cases} f_{j}^{F} e_{t}^{F} & \text{on } F, \\ \mathbf{0} & \text{on } \partial T \setminus F, \end{cases} \quad (10) $$

where $t \in \{1, 2\}$, $e_1^F$ and $e_2^F$ denote the two canonical unitary vectors tangential to $F$. The approximate solutions associated to these 108 BCs are computed in $T$ thanks to a local FR solver and used as basis functions for the Trefftz method. In the presentation, we will focus on the dependencies of this quasi-Trefftz approach in terms of the FR solver and boundary functions. Comparisons with the Trefftz method on matrix conditioning and iterative resolution will be held on challenging configurations to highlight the benefits of the quasi-Trefftz approach.

## RÉFÉRENCES

[1] O. Cessenat, B. Després. "Application of an Ultra Weak Variational Formulation of Elliptic PDEs to the Two-Dimensional Helmholtz Problem", SIAM Journal on Numerical Analysis, vol. 35, pp. 255-299, 1998.

[2] S. Congreve, J. Gedicke, I. Perugia. "Numerical Investigation of the Conditioning for Plane Wave Discontinuous Galerkin Methods", dans *Numerical Mathematics and Advanced Applications ENUMATH 2017*, pp. 493-500, 2019, Springer, Cham.

[3] H. Fure, S. Pernet, M. Sirdey, S. Tordeux, "A discontinuous Galerkin Trefftz type method for solving the two dimensional Maxwell equations", SN Partial Differential Equations and Applications, vol. 1, pp. 1-25, 2020.

[4] H. T. Huynh, "A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods", dans *18th AIAA Computational fluid dynamics conference*, p. 4079, 2007.

[5] F. Ihlenburg, I. Babuška. "Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM", *Computers & Mathematics with Applications*, vol. 30, pp. 9-37, 1995.

[6] S. Pernet, M. Sirdey, S. Tordeux, "Ultra-weak variational formulation for heterogeneous Maxwell problem in the context of high performance computing", *ESAIM: Proceedings and Surveys*, vol. 75, pp. 96-121, 2023.

[7] M. Rivet, S. Pernet, S. Tordeux. "Flux reconstruction method for time-harmonic linear propagation problems: 1D a priori error analysis", HAL, 2023.