# Optimised Correction Polynomial Functions for the Flux Reconstruction Method in Time-Harmonic Electromagnetism

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

► To cite this version:

Matthias Rivet, Sébastien Pernet, Sébastien Tordeux. Optimised Correction Polynomial Functions for the Flux Reconstruction Method in Time-Harmonic Electromagnetism. 2024. hal-04517554

HAL Id: hal-04517554

https://inria.hal.science/hal-04517554

Preprint submitted on 3 Apr 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
# OPTIMISED CORRECTION POLYNOMIAL FUNCTIONS FOR THE FLUX RECONSTRUCTION METHOD IN TIME-HARMONIC ELECTROMAGNETISM

MATTHIAS RIVET<fnref n="1" />, SÉBASTIEN PERNET<fnref n="2" /> AND SÉBASTIEN TORDEUX<fnref n="3" />

**Abstract.** The Flux Reconstruction (FR) method is classically used in the Computational Fluid Dynamics field. However, its use for the simulation of electromagnetic wave propagation is not as developed yet. Following on from the development of a priori error estimates for the 1D wave equations, we introduce optimisation problems to allow an adaptation of the FR correction polynomial functions to the discretisation parameters. Showing notable accuracy gains in 1D, especially in the preasymptotic regime, we generalise this procedure to the 3D Maxwell's equations, leading to similar interesting possibilities to reduce the computational cost for a given accuracy.

**2020 Mathematics Subject Classification.** 35L02, 35Q61, 65N35, 65N15.

03/04/2024.

## 1. INTRODUCTION

This article is interested in the 3D time-harmonic Maxwell's equations equipped with impedance Boundary Conditions (BCs). For an orthogonal parallelepiped domain $\Omega \subset \mathbb{R}^3$ of boundary $\partial\Omega$, their hyperbolic form stands as:

**Problem 1.1** (Time-harmonic Maxwell problem). Find the electromagnetic field $\mathbb{E} := (\mathbf{e}, \mathbf{h}) \in [\mathbb{L}^2(\Omega)]^6$ s.t.

$$ \left\{ \begin{array}{ll} i\kappa \mathbf{M} \mathbb{E} + \sum_{j=1}^{3} \frac{\partial \phi^j}{\partial x_j} = \mathbf{0} & \text{with } \phi^j = \mathbf{F}^j \mathbb{E} \quad \text{in } \Omega, \\ \\ \gamma_t[\mathbf{n}_{\partial\Omega}]\mathbf{e} + Z_{\partial\Omega} \gamma_{\times}[\mathbf{n}_{\partial\Omega}]\mathbf{h} = \mathbf{g} & \text{on } \partial\Omega, \end{array} \right. \quad (1) $$

by denoting the fluxes in each direction $j \in [1, 3]$ as $\phi^j$, in addition to

$$ \mathbf{M} = \begin{bmatrix} \varepsilon_r \mathbf{I}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mu_r \mathbf{I}_3 \end{bmatrix} \quad \text{and} \quad \mathbf{F}^j = \begin{bmatrix} \mathbf{0}_3 & -\mathbf{\Gamma}_\times[\mathbf{e}_j] \\ \mathbf{\Gamma}_\times[\mathbf{e}_j] & \mathbf{0}_3 \end{bmatrix} \quad \text{with } \forall \mathbf{n} = (n_1, n_2, n_3) \in \mathbb{R}^3, \mathbf{\Gamma}_\times[\mathbf{n}] = \begin{bmatrix} 0 & -n_3 & n_2 \\ n_3 & 0 & -n_1 \\ -n_2 & n_1 & 0 \end{bmatrix}, \quad (2) $$

<footnote>Keywords and phrases: Flux Reconstruction method, Discontinuous Galerkin, Spectral Differences method, correction polynomials, Maxwell's equations</footnote>

<footnote n="1">¹ DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France & EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142</footnote>

<footnote n="2">² DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France</footnote>

<footnote n="3">³ EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142</footnote>
OPTIMISED CORRECTION POLYNOMIAL FUNCTIONS FOR THE FLUX
RECONSTRUCTION METHOD IN TIME-HARMONIC ELECTROMAGNETISM

MATTHIAS RIVET<fnref n="1" />, SÉBASTIEN PERNET<fnref n="2" /> AND SÉBASTIEN TORDEUX<fnref n="3" />

**Abstract.** The Flux Reconstruction (FR) method is classically used in the Computational Fluid Dynamics field. However, its use for the simulation of electromagnetic wave propagation is not as developed yet. Following on from the development of a priori error estimates for the 1D wave equations, we introduce optimisation problems to allow an adaptation of the FR correction polynomial functions to the discretisation parameters. Showing notable accuracy gains in 1D, especially in the preasymptotic regime, we generalise this procedure to the 3D Maxwell's equations, leading to similar interesting possibilities to reduce the computational cost for a given accuracy.

**2020 Mathematics Subject Classification.** 35L02, 35Q61, 65N35, 65N15.

03/04/2024.

1. INTRODUCTION

This article is interested in the 3D time-harmonic Maxwell's equations equipped with impedance Boundary Conditions (BCs). For an orthogonal parallelepiped domain $\Omega \subset \mathbb{R}^3$ of boundary $\partial\Omega$, their hyperbolic form stands as:

**Problem 1.1** (Time-harmonic Maxwell problem). Find the electromagnetic field $\mathbb{E} := (\mathbf{e}, \mathbf{h}) \in [L^2(\Omega)]^6$ s.t.

$$
\left\{
\begin{array}{ll}
i\kappa \mathbf{M} \mathbb{E} + \sum_{j=1}^{3} \frac{\partial \phi^j}{\partial x_j} = \mathbf{0} & \text{with } \phi^j = \mathbf{F}^j \mathbb{E} \quad \text{in } \Omega, \\
\gamma_t[\mathbf{n}_{\partial\Omega}]\mathbf{e} + Z_{\partial\Omega} \gamma_{\times}[\mathbf{n}_{\partial\Omega}]\mathbf{h} = \mathbf{g} & \text{on } \partial\Omega,
\end{array}
\right.
\tag{1}
$$

by denoting the fluxes in each direction $j \in [1, 3]$ as $\phi^j$, in addition to

$$
\mathbf{M} = \begin{bmatrix} \varepsilon_r \mathbf{I}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mu_r \mathbf{I}_3 \end{bmatrix} \quad \text{and} \quad \mathbf{F}^j = \begin{bmatrix} \mathbf{0}_3 & -\mathbf{\Gamma}_\times[\mathbf{e}_j] \\ \mathbf{\Gamma}_\times[\mathbf{e}_j] & \mathbf{0}_3 \end{bmatrix} \quad \text{with } \forall \mathbf{n} = (n_1, n_2, n_3) \in \mathbb{R}^3, \mathbf{\Gamma}_\times[\mathbf{n}] = \begin{bmatrix} 0 & -n_3 & n_2 \\ n_3 & 0 & -n_1 \\ -n_2 & n_1 & 0 \end{bmatrix}, \tag{2}
$$

<footnote>Keywords and phrases: Flux Reconstruction method, Discontinuous Galerkin, Spectral Differences method, correction polyno-
mials, Maxwell’s equations</footnote>

<footnote n="1">¹ DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France & EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142</footnote>

<footnote n="2">² DTIS, ONERA, Université de Toulouse, 31000, Toulouse, France</footnote>

<footnote n="3">³ EPI Makutu, Inria, Université de Pau et des Pays de l'Adour, TotalEnergies, CNRS UMR 5142</footnote>
where $\varepsilon_r$ and $\mu_r$ are the relative permittivity and permeability, $Z_{\partial\Omega} \in L^\infty(\partial\Omega)$ denotes the boundary impedance that we suppose constant on each face, $\mathbf{g}$ is a tangential vector field of $[L^2(\partial\Omega)]^3$, $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$ is the canonical basis of $\mathbb{R}^3$ and $\mathbf{I}_3$ and $\mathbf{0}_3$ stand for the identity and null matrices of $\mathbb{C}^{3\times3}$.
The tangential component and trace of $\mathbf{w}$ are denoted on $\partial\Omega$ as

$$ \gamma_t[\mathbf{n}_{\partial\Omega}]\mathbf{w} = \mathbf{w} - (\mathbf{w} \cdot \mathbf{n}_{\partial\Omega}) \mathbf{n}_{\partial\Omega} = (\mathbf{n}_{\partial\Omega} \times \mathbf{w}) \times \mathbf{n}_{\partial\Omega} \quad \text{and} \quad \gamma_{\times}[\mathbf{n}_{\partial\Omega}]\mathbf{w} = \mathbf{n}_{\partial\Omega} \times \mathbf{w}, \qquad (3) $$

where $\mathbf{n}_{\partial\Omega}$ stands for the unitary outgoing normal from the domain $\Omega$.

**Remark 1.2.** The question of the well-posedness of Problem 1.1 has been studied in the reference book [4]. Then, a large variety of numerical methods were developed to deal with the simulation of electromagnetic waves propagation, as the Finite Difference Method [10], the Finite Element Method [1] and the Discontinuous Galerkin [2].

Yet, instead of focusing on one of these classical methods, we chose to adapt the Flux Reconstruction method (FR) [3] to this time-harmonic framework for a Cartesian mesh $\mathcal{T}_h$ of $\Omega$. Indeed, the FR method was introduced and is widely used in the Computational Fluid Dynamics community, but is not as common for wave propagation simulation.

It relies on the strong form of the hyperbolic problem (1) by approximating the solution $\mathbb{E}$ and the fluxes $\phi^j$ by piecewise polynomial functions $\mathbb{E}_h$ and $\tilde{\phi}_h^j$ such that for all $T \in \mathcal{T}_h$, $\mathbb{E}_h|_T \in [\mathbb{Q}_k]^6$ and $\tilde{\phi}_h^j|_T \in [\mathbb{Q}_k + x_j \mathbb{Q}_k]^6$, where $\mathbf{k} := (k_1, k_2, k_3)$ and $\mathbb{Q}_k$ is the polynomial space of degree at most $k_j$ in direction $j$, which verify

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

To close the system (4), it remains to define $\tilde{\phi}_h^j$ which is obtained from a polynomial correction of $\phi_h^j = \mathbf{F}^j \mathbb{E}_h$. This correction will be chosen to ensure the respect of the BCs and the flux continuity properties.

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 (corresponding to $x_1 = x_1^{T,+}$ and $x_1 = x_1^{T,-}$ respectively), such a correction takes the form

$$ \tilde{\phi}_{\mathbf{h}}^{1} = \phi_{\mathbf{h}}^{1} + \delta_{\mathbf{h},F^{+}}^{1,\to} + \delta_{\mathbf{h},F^{-}}^{1,\leftarrow} \quad \text{in } T \quad \text{with} \quad \delta_{\mathbf{h},F^{\pm}}^{1,\leftarrow}(\mathbf{x}) = (\hat{\phi}_{\mathbf{h},F^{\pm}}^{1} - \phi_{\mathbf{h}}^{1}|_{F^{\pm}})(x_2, x_3) P_{T}^{1,\leftarrow}(x_1), \qquad (5) $$

where $\hat{\phi}_{\mathbf{h},F^{\pm}}^1$ stands for a Riemann flux on $F^{\pm}$ (see [6, Sections 1-2] for a more general introduction and numerical trace definitions) and $P_T^{1,\leftarrow} = P^{1,\to} \circ \tau_T$, with the reference flux correction polynomial functions $P^{1,\leftarrow}$ lying in the admissible space

$$ \mathbb{P}_{adm}^{k_1} = \{ (P \to, P \leftarrow) \in [\mathbb{P}_{k_1+1}(0,1)]^2, P \to(0) = 1, P \to(1) = 0, P \leftarrow(0) = 0 \text{ and } P \leftarrow(1) = 1 \}, \qquad (6) $$

where $\mathbb{P}_{k_1+1}(0,1)$ stands for the space of polynomial functions of degree at most $k_1+1$, and

$$ \forall x_1 \in [x_1^{T,+}, x_1^{T,-}], \quad \tau_T(x_1) = [x_1 - x_1^{T,+}] [x_1^{T,-} - x_1^{T,+}]^{-1}. \qquad (7) $$

In particular, the FR framework allows to retrieve classic methods such as nodal Discontinuous Galerkin (DG) and Spectral Difference (SD) in less expensive ways, for specific flux correction polynomial functions [3]. In this article, we investigate the question of the optimisation of the correction polynomial functions, aiming at improving the properties of the FR method.
where $\varepsilon_r$ and $\mu_r$ are the relative permittivity and permeability, $Z_{\partial\Omega} \in L^\infty(\partial\Omega)$ denotes the boundary impedance that we suppose constant on each face, $\mathbf{g}$ is a tangential vector field of $[L^2(\partial\Omega)]^3$, $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$ is the canonical basis of $\mathbb{R}^3$ and $\mathbf{I}_3$ and $\mathbf{0}_3$ stand for the identity and null matrices of $\mathbb{C}^{3\times3}$.
The tangential component and trace of $\mathbf{w}$ are denoted on $\partial\Omega$ as

$$ \gamma_t[\mathbf{n}_{\partial\Omega}]\mathbf{w} = \mathbf{w} - (\mathbf{w} \cdot \mathbf{n}_{\partial\Omega}) \mathbf{n}_{\partial\Omega} = (\mathbf{n}_{\partial\Omega} \times \mathbf{w}) \times \mathbf{n}_{\partial\Omega} \quad \text{and} \quad \gamma_{\times}[\mathbf{n}_{\partial\Omega}]\mathbf{w} = \mathbf{n}_{\partial\Omega} \times \mathbf{w}, \qquad (3) $$

where $\mathbf{n}_{\partial\Omega}$ stands for the unitary outgoing normal from the domain $\Omega$.

**Remark 1.2.** The question of the well-posedness of Problem 1.1 has been studied in the reference book [4]. Then, a large variety of numerical methods were developed to deal with the simulation of electromagnetic waves propagation, as the Finite Difference Method [10], the Finite Element Method [1] and the Discontinuous Galerkin [2].

Yet, instead of focusing on one of these classical methods, we chose to adapt the Flux Reconstruction method (FR) [3] to this time-harmonic framework for a Cartesian mesh $\mathcal{T}_h$ of $\Omega$. Indeed, the FR method was introduced and is widely used in the Computational Fluid Dynamics community, but is not as common for wave propagation simulation.

It relies on the strong form of the hyperbolic problem (1) by approximating the solution $\mathbb{E}$ and the fluxes $\phi^j$ by piecewise polynomial functions $\mathbb{E}_h$ and $\tilde{\phi}_h^j$ such that for all $T \in \mathcal{T}_h$, $\mathbb{E}_h|_T \in [\mathbb{Q}_k]^6$ and $\tilde{\phi}_h^j|_T \in [\mathbb{Q}_k + x_j \mathbb{Q}_k]^6$, where $\mathbf{k} := (k_1, k_2, k_3)$ and $\mathbb{Q}_k$ is the polynomial space of degree at most $k_j$ in direction $j$, which verify

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

To close the system (4), it remains to define $\tilde{\phi}_h^j$ which is obtained from a polynomial correction of $\phi_h^j = \mathbf{F}^j \mathbb{E}_h$. This correction will be chosen to ensure the respect of the BCs and the flux continuity properties.

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 (corresponding to $x_1 = x_1^{T,+}$ and $x_1 = x_1^{T,-}$ respectively), such a correction takes the form

$$ \tilde{\phi}_h^1 = \phi_h^1 + \delta_{h,F^+}^{1,\to} + \delta_{h,F^-}^{1,\leftarrow} \quad \text{in } T \quad \text{with} \quad \delta_{h,F^\pm}^{1,\leftarrow}(\mathbf{x}) = (\hat{\phi}_{h,F^\pm}^1 - \phi_h^1|_{F^\pm}) (x_2, x_3) P_T^{1,\leftarrow}(x_1), \qquad (5) $$

where $\hat{\phi}_{h,F^\pm}^1$ stands for a Riemann flux on $F^\pm$ (see [6, Sections 1-2] for a more general introduction and numerical trace definitions) and $P_T^{1,\leftarrow} = P^{1,\to} \circ \tau_T$, with the reference flux correction polynomial functions $P^{1,\to}$ lying in the admissible space

$$ \mathbb{P}_{adm}^{k_1} = \{ (P \to, P \leftarrow) \in [\mathbb{P}_{k_1+1}(0,1)]^2, P \to(0) = 1, P \to(1) = 0, P \leftarrow(0) = 0 \text{ and } P \leftarrow(1) = 1 \}, \qquad (6) $$

where $\mathbb{P}_{k_1+1}(0,1)$ stands for the space of polynomial functions of degree at most $k_1+1$, and

$$ \forall x_1 \in [x_1^{T,+}, x_1^{T,-}], \quad \tau_T(x_1) = [x_1 - x_1^{T,+}] [x_1^{T,-} - x_1^{T,+}]^{-1}. \qquad (7) $$

In particular, the FR framework allows to retrieve classic methods such as nodal Discontinuous Galerkin (DG) and Spectral Difference (SD) in less expensive ways, for specific flux correction polynomial functions [3]. In this article, we investigate the question of the optimisation of the correction polynomial functions, aiming at improving the properties of the FR method.
## 2. A NEW CLASS OF OPTIMISED CORRECTION POLYNOMIAL FUNCTIONS FOR 1D MAXWELL'S EQUATIONS

The Cartesian structure of the discretisation allows to define efficient 3D methods thanks to optimisations that take place in a 1D framework, on which we will focus first. The efficiency of the associated FR methods are then compared to a nodal DG method for general BCs.

### 2.1. Construction strategy based on a specific 1D Maxwell problem

In this part, we consider the homogeneous time-harmonic 1D Maxwell's equations (the 1D wave equations) in a homogeneous material with incoming BCs, which can be written as:

**Problem 2.1** (Time-harmonic 1D wave problem). Let $\Omega := [0, L] \subset \mathbb{R}$ be an interval and $\kappa > 0$ be the wavenumber. Find $\mathbf{y} = (u, v)^T \in [\mathbb{H}^1(\Omega)]^2$ such that

$$i\kappa \mathbf{y} + \frac{d\phi}{dx} = \mathbf{0} \quad \text{with} \quad \phi = \mathbf{F}\mathbf{y} \quad \text{in } \Omega, \quad u(0) - Z_1 v(0) = g_1 \quad \text{and} \quad u(L) + Z_2 v(L) = g_2, \quad (8)$$

with the linear flux operator $\mathbf{F} = \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix}$, $Z_1 = Z_2 = 1$ (modeling incoming BCs) and $(g_1, g_2) \in \mathbb{C}^2$.

We consider the associated FR formulation, equipped with numerical traces detailed in [6, section 3.1] for a uniform mesh of cell size $h$, whose solution is denoted as $\mathbf{y}_h$. In particular, in what follows, we will only consider correction polynomial couples that respect the 'symmetry property' [3,9]

$$\forall x \in [0, 1], P^{\leftarrow}(x) = P^{\rightarrow}(1-x). \quad (9)$$

In [6], we derived a quasi-optimal asymptotic a priori estimate of the error $\varepsilon_h = \mathbf{y} - \mathbf{y}_h$ in the $L^2$-norm, which can be synthesised as:

$$\|\varepsilon_h\|_{L^2(\Omega)} \lesssim \kappa^{-1} (\kappa L) (\kappa h)^{2k+2} C_{asymp}(\kappa, L, P^{\rightarrow}) (|g_1|^2 + |g_2|^2), \quad (10)$$

where the dependencies of $C_{asymp}$ in terms of its parameters are explicitly known, and the notation $A \lesssim B$ means that there exists a constant $C > 0$, independent of $\kappa, L, h$ and $P^{\rightarrow}$, such that we have $A \le B C$. Our optimisation strategy will then involve the choice of the polynomial function that minimizes $C_{asymp}$, i.e.

$$P_{asymp}^{\rightarrow} := \arg \min_{\tilde{P} \in \mathbb{P}_{adm}^k} C_{asymp}(\kappa, L, \tilde{P}). \quad (11)$$

Yet, this bounding is sharp only for very small values of $h$: for coarse meshes, it may be irrelevant to define optimised correction polynomial functions. Consequently, we propose to resort to a preasymptotic estimate, obtained thanks to the integration of the error expression (132), in addition to (103) and (124) of [6]:

**Proposition 2.2.** We have the mesh-dependent error estimate

$$\|\varepsilon_h\|_{L^2(\Omega)} \lesssim C_{refin}(\kappa, h, L, P^{\rightarrow}) (|g_1|^2 + |g_2|^2), \quad (12)$$

where

$$C_{refin}(\kappa, h, L, P^{\rightarrow}) = h \sum_{n=1}^{N} \left[ \begin{array}{l} \left| \frac{\alpha_{\kappa h, k}}{T_{k+1}^{\rightarrow}} \right|^2 (\kappa h)^{2k+2} (1 + \gamma_{\kappa h, k}^{\rightarrow})^{2(n-1)} \int_0^1 |J_{\kappa h}^{\rightarrow}(x)|^2 dx + \left( (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} - 1 \right)^2 \\ + 2 \left| \frac{\alpha_{\kappa h, k}}{T_{k+1}^{\rightarrow}} \right| (\kappa h)^{k+1} (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} \left( (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} - 1 \right) \int_0^1 |J_{\kappa h}^{\rightarrow}(x)| dx \end{array} \right], \quad (13)$$
## 2. A NEW CLASS OF OPTIMISED CORRECTION POLYNOMIAL FUNCTIONS FOR 1D MAXWELL'S EQUATIONS

The Cartesian structure of the discretisation allows to define efficient 3D methods thanks to optimisations that take place in a 1D framework, on which we will focus first. The efficiency of the associated FR methods are then compared to a nodal DG method for general BCs.

### 2.1. Construction strategy based on a specific 1D Maxwell problem

In this part, we consider the homogeneous time-harmonic 1D Maxwell's equations (the 1D wave equations) in a homogeneous material with incoming BCs, which can be written as:

**Problem 2.1** (Time-harmonic 1D wave problem). Let $\Omega := [0, L] \subset \mathbb{R}$ be an interval and $\kappa > 0$ be the wavenumber. Find $\mathbf{y} = (u, v)^T \in [\mathbb{H}^1(\Omega)]^2$ such that

$$i\kappa \mathbf{y} + \frac{d\phi}{dx} = \mathbf{0} \quad \text{with} \quad \phi = \mathbf{F}\mathbf{y} \quad \text{in } \Omega, \quad u(0) - Z_1 v(0) = g_1 \quad \text{and} \quad u(L) + Z_2 v(L) = g_2, \quad (8)$$

with the linear flux operator $\mathbf{F} = \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix}$, $Z_1 = Z_2 = 1$ (modeling incoming BCs) and $(g_1, g_2) \in \mathbb{C}^2$.

We consider the associated FR formulation, equipped with numerical traces detailed in [6, section 3.1] for a uniform mesh of cell size $h$, whose solution is denoted as $\mathbf{y}_h$. In particular, in what follows, we will only consider correction polynomial couples that respect the 'symmetry property' [3,9]

$$\forall x \in [0, 1], P^{\leftarrow}(x) = P^{\rightarrow}(1-x). \quad (9)$$

In [6], we derived a quasi-optimal asymptotic a priori estimate of the error $\varepsilon_h = \mathbf{y} - \mathbf{y}_h$ in the $L^2$-norm, which can be synthesised as:

$$\|\varepsilon_h\|_{L^2(\Omega)} \lesssim \kappa^{-1} (\kappa L) (\kappa h)^{2k+2} C_{asymp}(\kappa, L, P^{\rightarrow}) (|g_1|^2 + |g_2|^2), \quad (10)$$

where the dependencies of $C_{asymp}$ in terms of its parameters are explicitly known, and the notation $A \lesssim B$ means that there exists a constant $C > 0$, independent of $\kappa, L, h$ and $P^{\rightarrow}$, such that we have $A \le B C$. Our optimisation strategy will then involve the choice of the polynomial function that minimizes $C_{asymp}$, i.e.

$$P_{asymp}^{\rightarrow} := \arg \min_{\tilde{P} \in \mathbb{P}_{adm}^k} C_{asymp}(\kappa, L, \tilde{P}). \quad (11)$$

Yet, this bounding is sharp only for very small values of $h$: for coarse meshes, it may be irrelevant to define optimised correction polynomial functions. Consequently, we propose to resort to a preasymptotic estimate, obtained thanks to the integration of the error expression (132), in addition to (103) and (124) of [6]:

**Proposition 2.2.** We have the mesh-dependent error estimate

$$\|\varepsilon_h\|_{L^2(\Omega)} \lesssim C_{refin}(\kappa, h, L, P^{\rightarrow}) (|g_1|^2 + |g_2|^2), \quad (12)$$

where

$$C_{refin}(\kappa, h, L, P^{\rightarrow}) = h \sum_{n=1}^{N} \left[ \begin{array}{l} \left| \frac{\alpha_{\kappa h, k}}{T_{k+1}^{\rightarrow}} \right|^2 (\kappa h)^{2k+2} (1 + \gamma_{\kappa h, k}^{\rightarrow})^{2(n-1)} \int_0^1 |J_{\kappa h}^{\rightarrow}(x)|^2 dx + \left( (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} - 1 \right)^2 \\ + 2 \left| \frac{\alpha_{\kappa h, k}}{T_{k+1}^{\rightarrow}} \right| (\kappa h)^{k+1} (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} \left( (1 + \gamma_{\kappa h, k}^{\rightarrow})^{n-1} - 1 \right) \int_0^1 |J_{\kappa h}^{\rightarrow}(x)| dx \end{array} \right], \quad (13)$$
<table><caption>TABLE 1. Common characteristics to the 1D simulation configuration.</caption><thead><tr><th>κ</th><th>L</th><th>Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th></tr></thead><tbody><tr><td>2π</td><td>10</td><td>0.2 - 1.7 i</td><td>2.4 + 1.5i</td><td>2.3 + 0.4i</td><td>-1.2i</td></tr></tbody></table>

with

$$ \begin{gathered} \forall l \in [0, k+1], \quad T_l \rightarrow (P^\rightarrow)^{(l)}(0), \quad \alpha_{\kappa h, k} \rightarrow (\kappa h)^{-(k+1)} |T_{k+1}| \left| \sum_{l=0}^{k+1} (-i\kappa h)^{-l} T_l \rightarrow \right|^{-1}, \quad (14) \\ \forall x \in [0, 1], \quad J_{\kappa h}(x) = 1 + \int_0^x (P^\rightarrow)^{(1)}(t) e^{i\kappa ht} dt \quad \text{and} \quad \gamma_{\kappa h, k} \rightarrow \alpha_{\kappa h, k} (\kappa h)^{k+1} \frac{|J_{\kappa h}(1)|}{|T_{k+1}|}. \end{gathered} $$

Similarly to (11), we then introduce the associated optimisation problem

$$ P_{refin} \rightarrow := \arg \min_{\tilde{P} \in \mathbb{P}_{adm}^k} C_{refin}(\kappa, h, L, \tilde{P}). \quad (15) $$

We will refer to the FR methods associated to these 'optimised' polynomials as FR_optiAsymp and FR_opti respectively.

**Remark 2.3.** The optimisation problems (11) and (15) are not easy to solve analytically, and different algorithms from the Scipy [8] library were considered for a numerical estimation of the optimum: the 'minimize' function without gradient, with gradient and with the Conjugate Gradient method [5, pp. 120-122], in addition to a differential evolution algorithm [7]. Yet, convergence issues led us to consider the use of all these optimisation algorithms before selecting the result having the smallest gradient norm value.

## 2.2. Numerical results

After introducing the asymptotic and refined optimised correction polynomial functions, we will illustrate their abilities for a general 1D configuration (i.e. for $Z_1, Z_2 \in \mathbb{C}$).

To do so, we take back Problem 2.1 in a domain $\Omega = [0, L]$, with a uniform mesh of $N \in \mathbb{N}^*$ cells of size $h = L N^{-1}$. In all the following experiments, we consider the parameters summarised in Table 1 with uniform polynomial degree $k$ and flux correction polynomial function in all the cells.

We aim at comparing the optimised correction polynomials to the classic choice of Radau polynomials, that we take as a reference: the associated FR_Radau method can be interpreted as a nodal DG method [3].

Table 2 gathers an illustration of the performances of the optimised correction polynomial functions for the 1D wave equations with BCs described in Table 1 for $L = 10$ (in terms of wavelengths). We present, according to the number of degrees of freedom (dofs) per wavelength, $N_{dof/\lambda}$, the L²-norm of the error relatively to the exact solution and the FR_Radau error. Three regimes can be showcased:

1. for $k = 1$ (that we did not illustrate due to its limited interest), both the optimisations show no perceptible difference with the FR_Radau configuration.

2. for $k = 2$, FR_opti shows a quasi-asymptotic behaviour even for small values of $N_{dof/\lambda}$, which allows to present interesting accuracy gains in the preasymptotic regime. Yet, in the asymptotic limit, FR_opti and FR_Radau tend to exhibit similar performances. The interest of FR_optiAsymp is more debatable and is likely to appear only in the asymptotic regime (as it may be expected).

3. The same conclusions stand for $k \ge 3$ putting aside the fact that the asymptotic behaviour of FR_opti remains better than the one of FR_Radau, especially when the polynomial degree $k$ increases.

Thus, the optimisation process allows up to 90% accuracy gains for the L²-norm, especially in the preasymptotic regime. Moreover, high polynomial degrees allow a 15% gain in the asymptotic regime.
<table><caption>TABLE 1. Common characteristics to the 1D simulation configuration.</caption><thead><tr><th>κ</th><th>L</th><th>Z<sub>1</sub></th><th>Z<sub>2</sub></th><th>g<sub>1</sub></th><th>g<sub>2</sub></th></tr></thead><tbody><tr><td>2π</td><td>10</td><td>0.2 - 1.7 i</td><td>2.4 + 1.5i</td><td>2.3 + 0.4i</td><td>-1.2i</td></tr></tbody></table>

with

$$ \begin{gathered} \forall l \in [0, k+1], \quad T_l \rightarrow (P^\rightarrow)^{(l)}(0), \quad \alpha_{\kappa h, k} \rightarrow (\kappa h)^{-(k+1)} |T_{k+1}| \left| \sum_{l=0}^{k+1} (-i\kappa h)^{-l} T_l \rightarrow \right|^{-1}, \quad (14) \\ \forall x \in [0, 1], \quad J_{\kappa h}(x) = 1 + \int_0^x (P^\rightarrow)^{(1)}(t) e^{i\kappa ht} dt \quad \text{and} \quad \gamma_{\kappa h, k} \rightarrow \alpha_{\kappa h, k} (\kappa h)^{k+1} \frac{|J_{\kappa h}(1)|}{|T_{k+1}|}. \end{gathered} $$

Similarly to (11), we then introduce the associated optimisation problem

$$ P_{refin} \rightarrow := \arg \min_{\tilde{P} \in \mathbb{P}_{adm}^k} C_{refin}(\kappa, h, L, \tilde{P}). \quad (15) $$

We will refer to the FR methods associated to these 'optimised' polynomials as FR_optiAsymp and FR_opti respectively.

**Remark 2.3.** The optimisation problems (11) and (15) are not easy to solve analytically, and different algorithms from the Scipy [8] library were considered for a numerical estimation of the optimum: the 'minimize' function without gradient, with gradient and with the Conjugate Gradient method [5, pp. 120-122], in addition to a differential evolution algorithm [7]. Yet, convergence issues led us to consider the use of all these optimisation algorithms before selecting the result having the smallest gradient norm value.

## 2.2. Numerical results

After introducing the asymptotic and refined optimised correction polynomial functions, we will illustrate their abilities for a general 1D configuration (i.e. for $Z_1, Z_2 \in \mathbb{C}$).

To do so, we take back Problem 2.1 in a domain $\Omega = [0, L]$, with a uniform mesh of $N \in \mathbb{N}^*$ cells of size $h = L N^{-1}$. In all the following experiments, we consider the parameters summarised in Table 1 with uniform polynomial degree $k$ and flux correction polynomial function in all the cells.

We aim at comparing the optimised correction polynomials to the classic choice of Radau polynomials, that we take as a reference: the associated FR_Radau method can be interpreted as a nodal DG method [3].

Table 2 gathers an illustration of the performances of the optimised correction polynomial functions for the 1D wave equations with BCs described in Table 1 for $L = 10$ (in terms of wavelengths). We present, according to the number of degrees of freedom (dofs) per wavelength, $N_{dof/\lambda}$, the L²-norm of the error relatively to the exact solution and the FR_Radau error. Three regimes can be showcased:

1. for $k = 1$ (that we did not illustrate due to its limited interest), both the optimisations show no perceptible difference with the FR_Radau configuration.

2. for $k = 2$, FR_opti shows a quasi-asymptotic behaviour even for small values of $N_{dof/\lambda}$, which allows to present interesting accuracy gains in the preasymptotic regime. Yet, in the asymptotic limit, FR_opti and FR_Radau tend to exhibit similar performances. The interest of FR_optiAsymp is more debatable and is likely to appear only in the asymptotic regime (as it may be expected).

3. The same conclusions stand for $k \ge 3$ putting aside the fact that the asymptotic behaviour of FR_opti remains better than the one of FR_Radau, especially when the polynomial degree $k$ increases.

Thus, the optimisation process allows up to 90% accuracy gains for the L²-norm, especially in the preasymptotic regime. Moreover, high polynomial degrees allow a 15% gain in the asymptotic regime.
<figure><img src="image_1.png" /><figcaption>TABLE 2. Performances of the correction polynomial optimisations on the FR method according to the number of dofs per wavelength $N_{dof/\lambda}$ for $L = 10$ and general BCs: L²-error with respect to the exact solution (row 1) and with respect to the FR\_Radau error (row 2), for different polynomial degrees $k \in [2, 5]$ and corrections (FR\_Radau (—★—), FR\_opti (—+—), FR\_optiAsymp (—▲—)).</figcaption></figure>

**Remark 2.4.** We highlight the fact the results are presented for general BCs, even if the optimisation process leans on incoming BCs. Moreover, it remains valid for smaller domains, with $L = 1$ for example. The asymptotic regime appearing quicker then, the remarkable preasymptotic behaviour of FR\_opti is less visible and we focused here on a larger configuration. In particular, we point out that the early quasi-optimal behaviour of FR\_opti appears in spite of a new optimisation for every configuration (contrary to classic choices, where the correction polynomial functions do not evolve with the mesh refinement).

**Remark 2.5.** One may observe a discontinuity in the behaviour of FR\_opti for $N_{dof/\lambda} = 160$ and $k = 5$: this is likely to be due to issues in the optimisation process to determine the correction polynomial function. Yet, we point out that this configuration is far from real use configurations. On the contrary, the preasymptotic gain would allow to obtain relatively precise solutions at a lower cost in more usual cases of interest.

### 3. DIRECTION-WISE OPTIMISATION FOR 3D MAXWELL'S EQUATIONS

As the FR method for the 3D Maxwell's equations relies on the choice of correction polynomial functions in each direction, we consider the direct extension to 3D by optimising according to the parameters of each direction. Due to its more regular performances in 1D, we will restrict to the FR\_opti method.

Then, we consider the same approach with parameters, and then correction polynomial functions, relative to each direction (as expected in (5)), leading to the three optimisation problems:

$$ \forall j \in [1, 3], \quad \min_{P_{j, \rightarrow} \in \mathbb{F}_{\text{adm}}^{k_j}} C_{\text{refin}} \rightarrow (\kappa_j, h_j, L_j, P^{j, \rightarrow}). \quad (16) $$
<figure><img src="image_2.png" /><figcaption>TABLE 2. Performances of the correction polynomial optimisations on the FR method according to the number of dofs per wavelength $N_{dof/\lambda}$ for $L = 10$ and general BCs: L²-error with respect to the exact solution (row 1) and with respect to the FR\_Radau error (row 2), for different polynomial degrees $k \in [2, 5]$ and corrections (FR\_Radau (—*—), FR\_opti (—+—), FR\_optiAsymp (—▲—)).</figcaption></figure>

**Remark 2.4.** We highlight the fact the results are presented for general BCs, even if the optimisation process leans on incoming BCs. Moreover, it remains valid for smaller domains, with $L = 1$ for example. The asymptotic regime appearing quicker then, the remarkable preasymptotic behaviour of FR\_opti is less visible and we focused here on a larger configuration. In particular, we point out that the early quasi-optimal behaviour of FR\_opti appears in spite of a new optimisation for every configuration (contrary to classic choices, where the correction polynomial functions do not evolve with the mesh refinement).

**Remark 2.5.** One may observe a discontinuity in the behaviour of FR\_opti for $N_{dof/\lambda} = 160$ and $k = 5$: this is likely to be due to issues in the optimisation process to determine the correction polynomial function. Yet, we point out that this configuration is far from real use configurations. On the contrary, the preasymptotic gain would allow to obtain relatively precise solutions at a lower cost in more usual cases of interest.

### 3. DIRECTION-WISE OPTIMISATION FOR 3D MAXWELL'S EQUATIONS

As the FR method for the 3D Maxwell's equations relies on the choice of correction polynomial functions in each direction, we consider the direct extension to 3D by optimising according to the parameters of each direction. Due to its more regular performances in 1D, we will restrict to the FR\_opti method.

Then, we consider the same approach with parameters, and then correction polynomial functions, relative to each direction (as expected in (5)), leading to the three optimisation problems:

$$ \forall j \in [1, 3], \quad \min_{P_{j, \rightarrow} \in \mathbb{F}_{\text{adm}}^{k_j}} C_{\text{refin}} \rightarrow (\kappa_j, h_j, L_j, P^{j, \rightarrow}). \quad (16) $$
<figure><img src="image_2.png" /><figcaption>TABLE 3. Performances of the FR methods according to the number of dof per wavelength for $L = 5$: L²-error with respect to the exact solution (row 1) and the FR_Radau error (row 2), for different polynomial degrees $k \in [2, 5]$ and corrections (FR_Radau (—*—), FR_opti1(—-), FR_optiSqrt(—), FR_optiMean (—)).</figcaption></figure>

Yet, while the domain length $L_j$ or mesh step $h_j$ have a clear sense in each direction $j \in [1, 3]$, this is not the case for the wave number $\kappa_j$: only the global value $\kappa$ is defined. Thus, we propose to introduce correction coefficients $\delta_\kappa$ such that the considered wave number for each direction is computed as $\kappa_j = \delta_\kappa \kappa$:

1. $\delta_\kappa^1 = 1$ allows to consider the original value of $\kappa$ in each direction, not taking into account the volumic propagation: it is referred as FR_opti1.

2. $\delta_\kappa^2 = 3^{-1/2}$ is associated to a plane wave whose propagation would be along the diagonal of a cubic volume and could be seen as the most extreme case on average: it is referred as FR_optiSqrt.

3. $\delta_\kappa^3 = 0.5 (\delta_\kappa^1 + \delta_\kappa^2)$ is the mean of the two previous cases: it is referred as FR_optiMean.

The numerical experiments are realised in a domain $\Omega = [0, L]^3$, with a uniform mesh of $N \in \mathbb{N}^*$ cells of size $h = L N^{-1}$ in each direction. In all the following experiments, we consider $\kappa = 2\pi$ and $\varepsilon_r = \mu_r = 1$ with uniform polynomial degree $k$ and flux correction polynomial functions in all the cells and every direction (even if direction-wise choices could have been made). For an exact solution $E_{exa} = (e_{exa}, h_{exa})$ made of a sum of six random plane waves, we consider the associated BCs with random positive values of $Z_{\partial\Omega}$ constant on each face:

$$g = \gamma_t[\mathbf{n}_{\partial\Omega}]e_{exa} + Z_{\partial\Omega} \gamma_{\times}[\mathbf{n}_{\partial\Omega}]h_{exa} \quad \text{on } \partial\Omega. \qquad (17)$$

Then, Table 3 gathers an illustration of the performances of the optimised correction polynomial functions for the 3D Maxwell's equations with general BCs for $L = 5$ (in terms of wavelengths). We present, according to the number of dofs per wavelength, the L²-norm of the error relatively to the exact solution and the FR_Radau error for $k \ge 2$, as the 1D conclusion for $k = 1$ remains valid in 3D.

One can observe that the general observations of the 1D case remain pertinent: the optimised correction allows accuracy gains, especially in the preasymptotic regime. Yet, concerning the wavenumber modification, FR_optiSqrt shows irregular gains, in particular in the asymptotic regime, while FR_opti1 shows limited abilities for coarse meshes at low polynomial order. Finally, the FR_optiMean seems to gather the advantages of these two approaches, showing interesting gains for every mesh refinement and polynomial degree: it appears as a robust and pertinent approach which offers as much as a 40% accuracy gain with respect to FR_Radau.
<figure><img src="image_3.png" /><figcaption>TABLE 3. Performances of the FR methods according to the number of dof per wavelength for $L = 5$: L²-error with respect to the exact solution (row 1) and the FR_Radau error (row 2), for different polynomial degrees $k \in [2, 5]$ and corrections (FR_Radau (—*—), FR_opti1(—-), FR_optiSqrt(—), FR_optiMean (—)).</figcaption></figure>

Yet, while the domain length $L_j$ or mesh step $h_j$ have a clear sense in each direction $j \in [1, 3]$, this is not the case for the wave number $\kappa_j$: only the global value $\kappa$ is defined. Thus, we propose to introduce correction coefficients $\delta_\kappa$ such that the considered wave number for each direction is computed as $\kappa_j = \delta_\kappa \kappa$:

1. $\delta_\kappa^1 = 1$ allows to consider the original value of $\kappa$ in each direction, not taking into account the volumic propagation: it is referred as FR_opti1.

2. $\delta_\kappa^2 = 3^{-1/2}$ is associated to a plane wave whose propagation would be along the diagonal of a cubic volume and could be seen as the most extreme case on average: it is referred as FR_optiSqrt.

3. $\delta_\kappa^3 = 0.5 (\delta_\kappa^1 + \delta_\kappa^2)$ is the mean of the two previous cases: it is referred as FR_optiMean.

The numerical experiments are realised in a domain $\Omega = [0, L]^3$, with a uniform mesh of $N \in \mathbb{N}^*$ cells of size $h = L N^{-1}$ in each direction. In all the following experiments, we consider $\kappa = 2\pi$ and $\varepsilon_r = \mu_r = 1$ with uniform polynomial degree $k$ and flux correction polynomial functions in all the cells and every direction (even if direction-wise choices could have been made). For an exact solution $E_{exa} = (e_{exa}, h_{exa})$ made of a sum of six random plane waves, we consider the associated BCs with random positive values of $Z_{\partial\Omega}$ constant on each face:

$$g = \gamma_t[\mathbf{n}_{\partial\Omega}]e_{exa} + Z_{\partial\Omega} \gamma_{\times}[\mathbf{n}_{\partial\Omega}]h_{exa} \quad \text{on } \partial\Omega. \qquad (17)$$

Then, Table 3 gathers an illustration of the performances of the optimised correction polynomial functions for the 3D Maxwell's equations with general BCs for $L = 5$ (in terms of wavelengths). We present, according to the number of dofs per wavelength, the L²-norm of the error relatively to the exact solution and the FR_Radau error for $k \ge 2$, as the 1D conclusion for $k = 1$ remains valid in 3D.

One can observe that the general observations of the 1D case remain pertinent: the optimised correction allows accuracy gains, especially in the preasymptotic regime. Yet, concerning the wavenumber modification, FR_optiSqrt shows irregular gains, in particular in the asymptotic regime, while FR_opti1 shows limited abilities for coarse meshes at low polynomial order. Finally, the FR_optiMean seems to gather the advantages of these two approaches, showing interesting gains for every mesh refinement and polynomial degree: it appears as a robust and pertinent approach which offers as much as a 40% accuracy gain with respect to FR_Radau.
**Remark 3.1.** We highlight the fact an optimisation is needed for each direction (even if a unique one can be realised in the case of an homogeneous mesh) and every configuration. Yet, this optimisation process has needed negligible computation time with respect to the linear system resolution in all the simulations we conducted (in particular, the differential evolution algorithm, which is the most consuming optimisation technique we considered, has only shown interest for the asymptotic configurations we faced in 1D and its use could be avoided for reasonable 3D cases).

## 4. CONCLUSION

We introduced optimisation processes for the correction polynomial functions on which relies the Flux Reconstruction method for the 1D wave equations. Leaning on a priori error estimates, they allowed noticeable gains in the $L^2$-norm for all polynomial degrees and mesh refinements, in particular in the preasymptotic regime. As the 3D FR approach leans on correction in each direction, such an optimisation was realised direction-wise with different propositions for the wanenumber modification. A reasonable choice of such a modification allowed to keep the interesting 1D gains in this 3D framework, and could reduce the computational cost for this electromagnetic wave propagation problem in general configurations.

## REFERENCES

[1] Gary Cohen and Sébastien Pernet. *Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations*. Scientific Computation. Springer Netherlands, Dordrecht, 2017.

[2] Daniele Antonio Di Pietro and Alexandre Ern. *Mathematical Aspects of Discontinuous Galerkin Methods*, volume 69 of *Mathématiques et Applications*. Springer, Berlin, Heidelberg, 2012.

[3] H.T. Huynh. A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods. *AIAA Paper AIAA 20074079*, pages 1–42, January 2007.

[4] Peter Monk. *Finite element methods for Maxwell's equations*. Oxford University Press, 2003.

[5] Jorge Nocedal and Stephen J Wright. *Numerical optimization*. Springer New York, NY, 2006.

[6] Matthias Rivet, Sébastien Pernet, and Sébastien Tordeux. Flux reconstruction method for time-harmonic linear propagation problems: 1D a priori error analysis. working paper or preprint, December 2023.

[7] Rainer Storn and Kenneth Price. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. *Journal of global optimization*, 11:341–359, 1997.

[8] Pauli Virtanen and et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

[9] Z. J. Wang and Haiyang Gao. A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. *Journal of Computational Physics*, 228(21):8161–8186, November 2009.

[10] Kane Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. *IEEE Transactions on Antennas and Propagation*, 14(3):302–307, May 1966. Conference Name: IEEE Transactions on Antennas and Propagation.
**Remark 3.1.** We highlight the fact an optimisation is needed for each direction (even if a unique one can be realised in the case of an homogeneous mesh) and every configuration. Yet, this optimisation process has needed negligible computation time with respect to the linear system resolution in all the simulations we conducted (in particular, the differential evolution algorithm, which is the most consuming optimisation technique we considered, has only shown interest for the asymptotic configurations we faced in 1D and its use could be avoided for reasonable 3D cases).

## 4. CONCLUSION

We introduced optimisation processes for the correction polynomial functions on which relies the Flux Reconstruction method for the 1D wave equations. Leaning on a priori error estimates, they allowed noticeable gains in the $L^2$-norm for all polynomial degrees and mesh refinements, in particular in the preasymptotic regime. As the 3D FR approach leans on correction in each direction, such an optimisation was realised direction-wise with different propositions for the wanenumber modification. A reasonable choice of such a modification allowed to keep the interesting 1D gains in this 3D framework, and could reduce the computational cost for this electromagnetic wave propagation problem in general configurations.

## REFERENCES

[1] Gary Cohen and Sébastien Pernet. *Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations*. Scientific Computation. Springer Netherlands, Dordrecht, 2017.

[2] Daniele Antonio Di Pietro and Alexandre Ern. *Mathematical Aspects of Discontinuous Galerkin Methods*, volume 69 of *Mathématiques et Applications*. Springer, Berlin, Heidelberg, 2012.

[3] H.T. Huynh. A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods. *AIAA Paper AIAA 20074079*, pages 1–42, January 2007.

[4] Peter Monk. *Finite element methods for Maxwell's equations*. Oxford University Press, 2003.

[5] Jorge Nocedal and Stephen J Wright. *Numerical optimization*. Springer New York, NY, 2006.

[6] Matthias Rivet, Sébastien Pernet, and Sébastien Tordeux. Flux reconstruction method for time-harmonic linear propagation problems: 1D a priori error analysis. working paper or preprint, December 2023.

[7] Rainer Storn and Kenneth Price. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. *Journal of global optimization*, 11:341–359, 1997.

[8] Pauli Virtanen and et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020.

[9] Z. J. Wang and Haiyang Gao. A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. *Journal of Computational Physics*, 228(21):8161–8186, November 2009.

[10] Kane Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. *IEEE Transactions on Antennas and Propagation*, 14(3):302–307, May 1966. Conference Name: IEEE Transactions on Antennas and Propagation.