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

HAL
open science

Reduced-order modeling for time-domain
electromagnetics

Stéphane Lanteri, Ting-Zhu Huang, Kun Li, Liang Li, Alan Youssef

► To cite this version:

Stéphane Lanteri, Ting-Zhu Huang, Kun Li, Liang Li, Alan Youssef. Reduced-order modeling for time-domain electromagnetics. Waves 2024 - The 16th International Conference on Mathematical and Numerical Aspects of Wave Propagation, Jun 2024, Berlin, Germany. pp.25-30, 10.17617/3.MBE4AA. hal-04907783

HAL Id: hal-04907783

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

Submitted on 23 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.

Public Domain
Reduced-order modeling for time-domain electromagnetics

Stéphane Lanteri<fnref n="1" />,*, Ting-Zhu Huang<fnref n="2" />, Kun Li<fnref n="3" />, Liang Li<fnref n="2" />, Alan Youssef<fnref n="1" />

**Abstract**

Reduced-order modeling (ROM) allows to construct simplifications of a high-fidelity model of a complex physical problem. The resulting lower fidelity (also referred as surrogate) model can capture the salient features of the source model so that one can quickly study a system's dominant effects using minimal computational resources. In this contribution, we present our works and achievements toward the development of ROM for time-domain electromagnetics.

**Keywords:** time-domain electromagnetics, reduced-order modeling, Proper Orthogonal Decomposition, Discontinuous Galerkin method.

1 Introduction

A variety of physical systems occurring in nature and in engineering are usually described by partial differential equations (PDEs) whose numerical solution requires considerable computational effort, especially in the case of parameterized PDEs, which often need to be solved for a range of parameter values. Such a situation is for instance encountered when studying electromagnetic wave propagation problems in complex geometries for different input parameters including frequency, direction of the incident wave, geometric dimensions, as well as material properties, which are described by the system of time-domain Maxwell's equations. The Discontinuous Galerkin Time-Domain (DGTD) method has progressively emerged as a popular method to solve this system. Although the DGTD method is highly accurate, it forms a high dimensional model (HDM) due to the duplication of the degrees of freedom (DoF) on the boundaries of the mesh elements, and it is often expensive in terms of both CPU time and memory demands. Hence, ROM methods are required, which allows for rapidly, inexpensively, and effectively providing the numerical solution of the problem for a range of parameter values

while maintaining a sufficiently high accuracy. The reduced basis (RB) method [4] is a widely used family of ROM techniques. Here, we report on our contributions on the study of RB-based ROM for time-domain electromagnetics, and more particularly on the applicability of Proper Orthogonal Decomposition (POD) in this context. In the last five years, we have proposed several variants of POD-based ROM, from an intrusive approach to purely data-driven non-intrusive ones. In all cases, the starting point is a Full-Order Model (FOM) relying on a high-order DGTD method for solving the system of time-domain Maxwell equations for heterogeneous media [3].

2 The POD method with snapshots

The POD method with snapshots introduced by Sirovich [15], also known as Karhunen-Loève decomposition, principal component analysis, or singular value decomposition, is probably the most widespread ROM technique at least in terms of applications [2]. The POD method begins with a set of snapshot vectors extracted from high-fidelity numerical simulations or experiments; the POD basis vectors are then obtained by an eigensystem of a correlation matrix generated by the snapshot matrix whose columns are the snapshot vectors. In the calssical POD method, a ROM is created by projecting (Galerkin projection) the semi-discrete numerical scheme used for obtaining the snapshots onto a low-dimensional space spanned by the POD basis. This approach is often referred as an *intrusive* ROM technique since it requires having access to the formulation of the original numerical scheme. The alternative approach that avoids this constraint consists in defining an appropriate decoupling of the construction of the POD basis and the prediction of the solution in a purely data-driven workflow. As a result, one obtains a *non-intrusive* ROM technique with decoupled offline (training) and online (inference) phases.

<footnote n="1">¹Université Côte d'Azur, Inria, CNRS, LJAD, Sophia Antipolis Cedex, France</footnote>

<footnote n="2">²University of Electronic Science and Technology of China Chengdu, China</footnote>

<footnote n="3">³Southwest University of Finance & Economics, Chengdu, China</footnote>

<footnote>*Email: stephane.lanteri@inria.fr</footnote>
Some contributions based on the POD method are presented in [11–14] for solving the time-domain Maxwell’s equations, e.g., in [12], a reduced-order FD method based on the POD method is studied for the time-domain Maxwell’s equations with respect to changes in the excitation frequencies and change of parameters; in [11], the authors propose a reduced-order Finite Difference Time-Domain method based on the POD method for solving the 2d Maxwell’s equations in a lossy medium. However, compared to other physical fields, the development of ROM methods for time-domain electromagnetics has been rather scarce.

### 3 High-fidelity model

We consider the system of time-domain Maxwell’s equations (normalized form) modeling the interaction of an incident wave with a scatterer $\Omega_s$ included in $\Omega \subset \mathbb{R}^3$

$$ \left\{ \begin{array}{l} \mu_r \frac{\partial \mathbf{H}(\mathbf{x}, t)}{\partial t} + \nabla \times \mathbf{E}(\mathbf{x}, t) = 0, \\ \varepsilon_r \frac{\partial \mathbf{E}(\mathbf{x}, t)}{\partial t} - \nabla \times \mathbf{H}(\mathbf{x}, t) = 0, \end{array} \right. \qquad (1) $$

$\forall (\mathbf{x}, t) \in \Omega \times \mathcal{T}$ and where $\mathbf{E} = (E_x, E_y, E_z)^T$ and $\mathbf{H} = (H_x, H_y, H_z)^T$ respectively denote the electric and magnetic fields (the symbol $\mathcal{T}$ denotes the transposition); $\mathcal{T} \subset \mathbb{R}^+$ is the time domain; $\varepsilon_r$ and $\mu_r$ are the relative permittivity and permeability, respectively. In particular, in the following, we may consider $\theta = (\varepsilon_r^{(1)}, \varepsilon_r^{(2)}, \dots, \varepsilon_r^{(d)}) \in \mathcal{P} \subset \mathbb{R}^d$ as problem’s parameters with $\varepsilon_r^{(i)}$ ($i = 1, 2, \dots, d$) being the relative permittivity in the $i$-th domain of a multilayer structure $\Omega_s$, $\mathcal{P}$ being the parameter domain, and $d$ being the number of parameters. Our goal is to solve the system (1) in $\Omega$ with boundary $\partial\Omega = \Gamma_a$, where we impose the first-order Silver-Müller absorbing boundary condition (ABC)

$$ \mathbf{n} \times \mathbf{E} + Z \mathbf{n} \times (\mathbf{n} \times \mathbf{H}) = \mathbf{n} \times \mathbf{E}^{\text{inc}} + Z \mathbf{n} \times (\mathbf{n} \times \mathbf{H}^{\text{inc}}). $$

Here, $\mathbf{n}$ denotes the unit normal vector pointing outward to $\partial\Omega$, $\mathbf{E}^{\text{inc}}$ and $\mathbf{H}^{\text{inc}}$ stand for the incident fields, and $Z = \sqrt{\mu_r/\varepsilon_r}$. System (1) with the boundary condition is completed with initial conditions $\mathbf{E}(\mathbf{x}, 0) = \mathbf{E}_0(\mathbf{x})$, and $\mathbf{H}(\mathbf{x}, 0) = \mathbf{H}_0(\mathbf{x})$ for $\mathbf{x} \in \Omega \subset \mathbb{R}^3$, where $\mathbf{E}_0$ and $\mathbf{H}_0$ denote some given functions. In a Discontinuous Galerkin (DG) framework the computational domain $\Omega$ is divided into elements

$\Omega_h = \bigcup_{i=1}^{N_{\Omega}} K_i$ with $N_{\Omega}$ being the number of elements. Then the DG method seeks an approximate solution $(\mathbf{E}_h, \mathbf{H}_h)$ in the space $\mathbf{V}_h \times \mathbf{V}_h$ that satisfies for all $K_i$ in $\Omega_h$ and $\forall \mathbf{v} \in \mathbf{V}_h$

$$ \left\{ \begin{array}{l} (\mu_r \frac{\partial \mathbf{H}_h}{\partial t}, \mathbf{v})_{K_i} + (\mathbf{E}_h, \mathrm{curl} \mathbf{v})_{K_i} \\ \qquad - \langle \mathbf{E}_h^* \times \mathbf{n}, \mathbf{v} \rangle_{\partial K_i} = 0, \\ \\ (\varepsilon_r \frac{\partial \mathbf{E}_h}{\partial t}, \mathbf{v})_{K_i} - (\mathbf{H}_h, \mathrm{curl} \mathbf{v})_{K_i} \\ \qquad + \langle \mathbf{H}_h^* \times \mathbf{n}, \mathbf{v} \rangle_{\partial K_i} = 0. \end{array} \right. \qquad (2) $$

Here, $\mathbf{E}_h^*$ and $\mathbf{H}_h^*$ are the so-called numerical traces used to communicate information between adjacent elements, and $\mathbf{V}_h = \{\mathbf{v} \in (L^2(\Omega))^3 | \mathbf{v}|_{K_i} \in (\mathbb{P}_{p_i}(K_i))^3, 1 \le i \le N_{\Omega}\}$ where $\mathbb{P}_{p_i}(K_i)$ is the space of nodal polynomials of degree at most $p_i$ inside the element $K_i$. By gathering the electric and magnetic vectors DoF in each element into column vectors, we derive the following global semi-discrete system

$$ \left\{ \begin{array}{l} \mu_r \mathbf{M} \frac{\partial \mathbf{H}_h}{\partial t} = - \mathbf{K} \mathbf{E}_h + \mathbf{S}^i \mathbf{E}_h + \mathbf{S}^e \mathbf{H}_h + \mathbf{B}^e(t), \\ \varepsilon_r \mathbf{M} \frac{\partial \mathbf{E}_h}{\partial t} = \mathbf{K} \mathbf{H}_h - \mathbf{S}^i \mathbf{H}_h - \mathbf{S}^h \mathbf{E}_h - \mathbf{B}^h(t). \end{array} \right. \qquad (3) $$

We refer to [3] for more details on matrices $\mathbf{M}$, $\mathbf{K}$, $\mathbf{S}^i$ and $\mathbf{S}^e$. Moreover, $\mathbf{E}_h^*$ and $\mathbf{H}_h^*$ are computed using a centered scheme and system (3) is time-integrated with a second-order Leap-Frog (LF2) scheme.

### 4 Projected POD-DGTD method

Our first contribution on the design of POD-based ROM for system (1) is introduced in [9]. In this work, the POD basis is created by the eigensystem of the correlation matrix, which is generated by the snapshot matrix whose columns are the snapshot vectors extracted from high-identity DGTD simulations. A POD-DGTD formulation with lower dimension is established by applying a Galerkin projection to the semi-discrete DG scheme (3). We equidistantly extract $\ell$ vectors $\mathbf{u}_h^{n_i}$ ($0 \le n_1 \le \dots \le n_\ell \le N_t, \ell \ll N_t$) from the transient DGTD solutions $\{\mathbf{u}_h^n\}_{n=0}^{N_t}$ ($\mathbf{u}_h^n = \mathbf{E}_h^{(n)}, \mathbf{H}_h^{(n+\frac{1}{2})}$). Further, we formulate two $N \times \ell$ ($N = 3N_{dof}$) snapshot matrices

$$ \mathbf{A}_{\mathbf{u}} = \begin{pmatrix} \mathbf{u}_{h,1}^{n_1} & \mathbf{u}_{h,1}^{n_2} & \cdots & \mathbf{u}_{h,1}^{n_\ell} \\ \mathbf{u}_{h,2}^{n_1} & \mathbf{u}_{h,2}^{n_2} & \cdots & \mathbf{u}_{h,2}^{n_\ell} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{u}_{h,N}^{n_1} & \mathbf{u}_{h,N}^{n_2} & \cdots & \mathbf{u}_{h,N}^{n_\ell} \end{pmatrix}, \quad \mathbf{u} = \mathbf{E}, \mathbf{H}. $$
Let

$$ \mathbf{U}_{\mathbf{u}}^{T} \mathbf{A}_{\mathbf{u}} \mathbf{V}_{\mathbf{u}} = \begin{pmatrix} \sum_{r_{\mathbf{u}} \times r_{\mathbf{u}}}^{\mathbf{u}} & \mathbf{O}_{r_{\mathbf{u}} \times (\ell - r_{\mathbf{u}})} \\ \mathbf{O}_{(N - r_{\mathbf{u}}) \times r_{\mathbf{u}}} & \mathbf{O}_{(N - r_{\mathbf{u}}) \times (\ell - r_{\mathbf{u}})} \end{pmatrix}, $$

where $\mathbf{U}_{\mathbf{u}}$ and $\mathbf{V}_{\mathbf{u}}$ are $N \times N$ and $\ell \times \ell$ unitary matrices, $r_{\mathbf{u}}$ is the rank of $\mathbf{A}_{\mathbf{u}}$, and $\sum_{r_{\mathbf{u}} \times r_{\mathbf{u}}}^{\mathbf{u}} = \text{diag}(\sigma_{\mathbf{u},1}, \sigma_{\mathbf{u},2}, \dots, \sigma_{\mathbf{u},r_{\mathbf{u}}})$ with $\sigma_{\mathbf{u},1} \ge \sigma_{\mathbf{u},2} \ge \dots \ge \sigma_{\mathbf{u},r_{\mathbf{u}}} \ge 0$ being the singular values of $\mathbf{A}_{\mathbf{u}}$. Let $\mathbf{U}_{\mathbf{u}} = (\psi_{\mathbf{u},1}, \psi_{\mathbf{u},2}, \dots, \psi_{\mathbf{u},N})$ and $\mathbf{V}_{\mathbf{u}} = (\phi_{\mathbf{u},1}, \phi_{\mathbf{u},2}, \dots, \phi_{\mathbf{u},\ell})$, the POD basis of dimension $d_{\mathbf{u}}$ ($d_{\mathbf{u}} \le r_{\mathbf{u}}$) is the set $\{\psi_{\mathbf{u},i}\}_{i=1}^{d_{\mathbf{u}}}$ in the ROM context. It is well known that

$$ \mathbf{A}_{\mathbf{u}} \phi_{\mathbf{u},i} = \sigma_{\mathbf{u},i} \psi_{\mathbf{u},i}, \quad \text{and} \quad \mathbf{A}_{\mathbf{u}}^{T} \psi_{\mathbf{u},i} = \sigma_{\mathbf{u},i} \phi_{\mathbf{u},i}, \quad (4) $$

so that $\mathbf{A}_{\mathbf{u}}\mathbf{A}_{\mathbf{u}}^{T}\psi_{\mathbf{u},i} = \sigma_{\mathbf{u},i}^{2}\psi_{\mathbf{u},i}$ and $\mathbf{A}_{\mathbf{u}}^{T}\mathbf{A}_{\mathbf{u}}\phi_{\mathbf{u},i} = \sigma_{\mathbf{u},i}^{2}\phi_{\mathbf{u},i}$ where $i = 1, 2, \dots, r_{\mathbf{u}}$. Let $\mathbf{\Psi}_{\mathbf{u}} = (\psi_{\mathbf{u},1}, \psi_{\mathbf{u},2}, \dots, \psi_{\mathbf{u},d_{\mathbf{u}}})$ ($\mathbf{u} = \mathbf{E}, \mathbf{H}$). Using a Galerkin ansatz, the reduced-basis approximation of the fields takes are given by

$$ \underline{\mathbf{E}}_h \approx \underline{\mathbf{E}}_h^d = \Psi_E \alpha_E(t), \quad \underline{\mathbf{H}}_h \approx \underline{\mathbf{H}}_h^d = \Psi_H \alpha_H(t), \quad (5) $$

where $\alpha_E(t) \in \mathbb{R}^{d_E}$ and $\alpha_H(t) \in \mathbb{R}^{d_H}$ respectively denote a general time-varying solution for $\mathbf{E}$ and $\mathbf{H}$; $\underline{\mathbf{E}}_h^d$ and $\underline{\mathbf{H}}_h^d$ denote the reduced-order solutions. Applying the Galerkin projection to system (3), we obtain the POD-DGTD formulation

$$ \left\{ \begin{array}{l} \mu_r \mathbf{M}^{\mathbf{H}} \alpha'_{\mathbf{H}}(t) = \mathbf{K}^{\mathbf{H},\mathbf{E}} \alpha_{\mathbf{E}}(t) \\ + \mathbf{S}_i^{\mathbf{H},\mathbf{E}} \alpha_{\mathbf{E}}(t) + \mathbf{S}_e^{\mathbf{H},\mathbf{E}} \widehat{\alpha}_{\mathbf{E}}(t) + \mathbf{B}_e^{\mathbf{H}}(t), \\ \varepsilon_r \mathbf{M}^{\mathbf{E}} \alpha'_{\mathbf{E}}(t) = \mathbf{K}^{\mathbf{E},\mathbf{H}} \alpha_{\mathbf{H}}(t) \\ - \mathbf{S}_i^{\mathbf{E},\mathbf{H}} \alpha_{\mathbf{H}}(t) - \mathbf{S}_h^{\mathbf{E},\mathbf{H}} \widehat{\alpha}_{\mathbf{H}}(t) - \mathbf{B}_h^{\mathbf{E}}(t), \end{array} \right. \quad (6) $$

where the symbol ' is a time derivative, and the matrix $\mathbf{M}^{\mathbf{u}}$, $\mathbf{K}^{\mathbf{u},\mathbf{w}}$ and $\mathbf{S}_i^{\mathbf{u},\mathbf{w}}$ ($\mathbf{u}, \mathbf{w} = \mathbf{H}, \mathbf{E}$) are respectively given by $\mathbf{M}^{\mathbf{u}} = \Psi_{\mathbf{u}}^{T}\mathbf{M}\Psi_{\mathbf{u}}$ (of size $d_{\mathbf{u}} \times d_{\mathbf{u}}$), $\mathbf{K}^{\mathbf{u},\mathbf{w}} = \Psi_{\mathbf{u}}^{T}\mathbf{K}\Psi_{\mathbf{w}}$ (of size $d_{\mathbf{u}} \times d_{\mathbf{w}}$) and $\mathbf{S}_i^{\mathbf{u},\mathbf{w}} = \Psi_{\mathbf{u}}^{T}\mathbf{S}_i^{\mathbf{u}}\Psi_{\mathbf{w}}$ (of size $d_{\mathbf{u}} \times d_{\mathbf{w}}$); we refer to [9] for the definition of the other matrices in (6). We then obtain the fully discrete POD-DGTD formulation based on the $\text{LF}_2$ time scheme as follows

$$ \left\{ \begin{array}{l} \varepsilon_r \mathbf{M}^{\mathbf{E}} \frac{\delta \alpha_{\mathbf{E}}}{\Delta t} = \mathbf{K}^{\mathbf{E},\mathbf{H}} \alpha_{\mathbf{H}}^{(n+\frac{1}{2})} \\ - \mathbf{S}_i^{\mathbf{E},\mathbf{H}} \alpha_{\mathbf{H}}^{(n+\frac{1}{2})} - \mathbf{S}_h^{\mathbf{E},\mathbf{H}} \widehat{\alpha}_{\mathbf{H}}^{(n+\frac{1}{2})} - \mathbf{B}_h^{\mathbf{E}}(n\Delta t), \\ \mu_r \mathbf{M}^{\mathbf{H}} \frac{\delta \alpha_{\mathbf{H}}}{\Delta t} = - \mathbf{K}^{\mathbf{H},\mathbf{E}} \alpha_{\mathbf{E}}^{(n+1)} \\ + \mathbf{S}_i^{\mathbf{H},\mathbf{E}} \alpha_{\mathbf{E}}^{(n+1)} + \mathbf{S}_e^{\mathbf{H},\mathbf{E}} \widehat{\alpha}_{\mathbf{E}}^{(n+1)} + \mathbf{B}_e^{\mathbf{H}}((n+\frac{1}{2})\Delta t), \end{array} \right. $$

with $\delta\alpha_{\mathbf{E}} = \alpha_{\mathbf{E}}^{(n+1)} - \alpha_{\mathbf{E}}^{(n)}$, $\delta\alpha_{\mathbf{H}} = \alpha_{\mathbf{H}}^{(n+\frac{3}{2})} - \alpha_{\mathbf{H}}^{(n+\frac{1}{2})}$ (see [9] for the definition of $\widehat{\alpha}_{\mathbf{H}}^{(n+\frac{1}{2})}$ and $\widehat{\alpha}_{\mathbf{E}}^{(n+1)}$). The choice of the dimension $d_{\mathbf{u}}$ of the POD basis is not obvious. We show how to define $d_{\mathbf{u}}$ in this remark. The POD energy or error in the POD basis [2] can be defined by

$$ \vartheta(\psi_{\mathbf{u},1}, \psi_{\mathbf{u},2}, \dots, \psi_{\mathbf{u},d_{\mathbf{u}}}) = \sum_{j=d_{\mathbf{u}}+1}^{r_{\mathbf{u}}} \lambda_{\mathbf{u},j}, \quad \mathbf{u} = \mathbf{E}, \mathbf{H}. $$

If we want the POD energy to be less than some prescribed tolerance $\rho$, i.e., that

$$ \vartheta(\psi_{\mathbf{u},1}, \psi_{\mathbf{u},2}, \dots, \psi_{\mathbf{u},d_{\mathbf{u}}}) \le \rho, $$

then we can choose $d_{\mathbf{u}}$ to be the smallest integer such that $I(d_{\mathbf{u}}) \ge 1 - \rho$ with $I(d_{\mathbf{u}}) = \sum_{j=1}^{d_{\mathbf{u}}} \lambda_{\mathbf{u},j} / \sum_{j=1}^{r_{\mathbf{u}}} \lambda_{\mathbf{u},j}$. To illustrate the capabilities of the POD-DGTD method, we consider a two-dimensional (2d) problem given by the scattering of a modulated Gaussian by an airfoil profile. We set $\varepsilon_r = \mu_r = 1$. The computational domain is artificially truncated by $\Omega = [-1 \text{ m}, 2 \text{ m}] \times [-1 \text{ m}, 1 \text{ m}]$. The boundary of the rectangle is assumed to be the first order Silver-Müller ABC, and the surface of the airfoil is assumed to be the PEC. In the present case, there is no incident field and a local source term is added, which is given by

$$ J_z^s(\vec{x}, t) = A e^{-(\frac{t-4\tau}{\tau})^2} \sin(2\pi f_c(t - 4\tau)) \frac{g(\vec{x})}{\|g(\vec{x})\|}, $$

where $A$ is the amplitude of the signal, $\tau = 1.7 \ 10^3 \text{ps}$, $f_c = 1.2 \text{ GHz}$ and $g(\vec{x})$ is a 2d Gaussian function $g(\vec{x}) = e^{-\beta((x-x_0)^2+(y-y_0)^2)}$ with $(x_0, y_0) = (-0.3,0.0)$ the center of the Gaussian spatial support. Here, the parameter $\beta$ has been chosen such that the source term $J_z^s$ is strongly localized. We set $\beta = 2.5 \ 10^3$. The total simulation time, $T_f$, is set to $T_f = 10 \text{ m}$ (normalized unit). The simulations are performed on a mesh with 8,436 nodes and 16,460 elements, and mesh size $h = 3.993 \ 10^{-2}$. A performance comparison between the DGTD and POD-DGTD methods is given in Table 1. The $L^2$-error is the maximum global error over all the time steps, $T_{\text{cons}}$ denotes the CPU time for the construction of the global system, $T_{\text{sol}}$ denotes the CPU time for the DGTD or POD-DGTD methods, and $T_{\text{POD}}$ denotes the CPU time for the construction of the POD basis. We define Speed-up
<table><caption>Table 1: Comparisons between the DGTD and POD-DGTD methods in terms of CPU time (in seconds) and $L_2$-error for the scattering of a modulated Gaussian by an airfoil profile.</caption><thead><tr><th rowspan="2">P<sub>i</sub> T<sub>cons</sub></th><th colspan="2">DGTD</th><th colspan="4">POD-DGTD</th><th rowspan="2">Speed-up</th><th rowspan="2">L<sup>2</sup>-error</th></tr><tr><th>T<sub>sol</sub></th><th>d<sub>H<sub>x</sub></sub></th><th>d<sub>H<sub>y</sub></sub></th><th>d<sub>E<sub>z</sub></sub></th><th>T<sub>POD</sub></th><th>T<sub>sol</sub></th></tr></thead><tbody><tr><td>1</td><td>115.6</td><td>119.8</td><td>19</td><td>19</td><td>19</td><td>1.4</td><td>31.8</td><td><b>3.6</b></td><td>9.324 10<sup>-2</sup></td></tr><tr><td>2</td><td>295.9</td><td>851.5</td><td>15</td><td>15</td><td>15</td><td>1.9</td><td>73.7</td><td><b>11.2</b></td><td>6.254 10<sup>-2</sup></td></tr></tbody></table>

as Time1/Time2 where Time1 denotes $T_{sol}$ for DGTD solution and Time2 represents the sum of $T_{sol}$ for POD-DGTD solution and $T_{POD}$. The time evolution of the fields $H_x$, $H_y$ and $E_z$ at a selected point of the mesh for the DGTD and POD-DGTD methods with a $P_2$ approximation are shown in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Figure 1: Scattering of a modulated Gaussian by an airfoil profile: distribution of the intensity of the electric field, $|E_z|$, at time $t = 2.11$ m (top) and 3.20 m (bottom).</figcaption></figure>

## 5 Data-driven POD-CSI method

As mentioned previously, one of the major drawbacks of intrusive ROM methods such as the one presented in section 4 is that they require access to the HDM. Non-intrusive ROM techniques are currently developing rapidly with the emergence of various viewpoints [17]. In a non-intrusive ROM, the reduced-order coefficients for a new time (and parameter) value are obtained via an approximate mapping. The high-

fidelity solver is only used to generate the snapshots and the training dataset, and there is a complete decoupling between the online evaluation and the offline training. Several strategies have been proposed such as the regression-based non-intrusive RB methods including the tensor decomposition-based and the Gaussian processes-based regression, and the interpolation-based non-intrusive RB methods comprising the radial basis function (RBF) interpolations, the polynomial interpolations and the cubic spline interpolations (CSI).

In [6] we present a non-intrusive ROM for the solution of parameterized electromagnetic scattering problems. In this approach, in addition to the time parameter, the dielectric permittivity $\epsilon_r$ is also considered as a reduction parameter. As with the POD-DGTD method, snapshot vectors are extracted from high-order DGTD simulations. Because the second dimension of the snapshot matrix is large, we introduce a two-step or nested POD method to extract time- and parameter-independent POD basis functions. SVD is used to compute the principal components of the projection coefficient matrices (also referred to as the reduced coefficient matrices) of full-order solutions onto the RB subspace. A CSI approach is proposed to approximate the dominating time- and parameter-modes of the reduced coefficient matrices without resorting to Galerkin projection. The generation of snapshot vectors, the construction of POD basis functions and the approximation of reduced coefficient matrices based on the CSI method are completed during the offline stage. The reduced-order solutions for new time and parameter values can be rapidly recovered via outputs from the interpolation models in the online stage. In particular, the offline and online stages of the proposed RB method, termed as the POD-CSI method, are completely decoupled, which ensures the computational validity of this POD-CSI method. We illustrate the capabilities of the proposed POD-CSI method in the case of the scattering of plane wave by a multilayer medium where each layer has a cylindrical geometry. The computational domain is artificially truncated by a square with 6.4 m side length on which the first order Silver-Müller ABC is applied. Simulations are performed using a mesh with 2,049 nodes and 4,016 elements of which 865 elements are located inside the mul-
tilayer medium Table 2 summarizes the distribution and range of material parameters considered for this study. The excitation is a

<table><caption>Table 2: Scattering of plane wave by a multilayer heterogeneous medium: the distribution and range of material parameters ($r_i$ is the radius of the $i$th cylinder).</caption><thead><tr><th>Layer <i>i</i></th><th><i>P</i><sup><i>i</i></sup></th><th><i>μ</i><sub><i>r</i></sub><sup>(<i>i</i>)</sup></th><th><i>r</i><sub><i>i</i></sub></th></tr></thead><tbody><tr><td>1</td><td>ε<sub><i>r</i></sub><sup>(1)</sup> ∈ [5.0, 5.6]</td><td>1.0</td><td>0.15</td></tr><tr><td>2</td><td>ε<sub><i>r</i></sub><sup>(2)</sup> ∈ [3.25, 3.75]</td><td>1.0</td><td>0.3</td></tr><tr><td>3</td><td>ε<sub><i>r</i></sub><sup>(3)</sup> ∈ [2.0, 2.5]</td><td>1.0</td><td>0.45</td></tr><tr><td>4</td><td>ε<sub><i>r</i></sub><sup>(4)</sup> ∈ [1.25, 1.75]</td><td>1.0</td><td>0.6</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Scattering of plane wave by a multilayer heterogeneous medium: 1d plots of the real part of DGTD and POD-CSI solutions (DFT) of $H_y$ along $y = 0$ (left), and bistatic RCS (right) based on the DGTD and POD-CSI solutions with the testing parameters $\theta = \{(5.5, 3.7, 2.4, 1.7)\}$.</figcaption></figure>

plane wave such that $E_z^{\text{inc}} = \cos(-kx + \omega t)$ and $H_y^{\text{inc}} = -\cos(-kx + \omega t)$ with $\omega = 2\pi f$ being the angular frequency of the incident wave frequency $f = 300$ MHz, and $k$ being the wave number in vacuum. The physical simulation time is set to 50 periods of the incident wave oscillation for each 4-dimensional parameter $\theta \in P$ with time step $\Delta t = 0.0038$ m using the DGTD method with a $\mathbb{P}_2$ approximation, where $P = P^1 \times P^2 \times P^3 \times P^4 \subset \mathbb{R}^4$ with $\varepsilon_r^{(i)} \in P^i$ ($i = 1, 2, \dots, 4$). We use a grid sampling of tensor product with $N_{\Delta p} = 2$ and 3 uniform points for each parameter to form the training parameter samples $P_h^{\text{tr}}$, respectively, resulting in $N_p = 16$ and 81 points. Each choice of the parameter is sampled for $N_t = 253$ snapshots in time at last period. We use the two-step POD method to extract the POD basis functions from all the snapshots. The performance results obtained by the POD-CSI and DGTD methods with $\mathbb{P}_2$ approximation are summarized in Table 3. We obtain a speed-up of 403, which

show the significantly enhanced efficiency of the POD-CSI method.

<table><caption>Table 3: Scattering of plane wave by a multilayer heterogeneous medium: the comparison between the POD-CSI (offline and online) and DGTD methods in terms of CPU time. The unit of time cost is second.</caption><thead><tr><th colspan="3">Offline</th></tr><tr><th>Snapshots</th><th>Two-step POD</th><th>CSI</th></tr></thead><tbody><tr><td>4.866 × 10<sup>4</sup></td><td>1.699 × 10<sup>2</sup></td><td>1.607 × 10<sup>2</sup></td></tr></tbody><thead><tr><th colspan="2">Online</th></tr><tr><th>POD-CSI</th><th>DGTD</th></tr></thead><tbody><tr><td>1.49</td><td>6.007 × 10<sup>2</sup></td></tr></tbody></table>

## 6 Further works and contributions

In addition to the two contributions discussed in the previous sections, we have investigated several directions to (1) improve the performance of POD-based ROM and, (2) deal with more complex problems. For instance, in [7], we propose a POD-based ROM with an adaptive snapshot selection strategy. Inspired from [1], we derive an error bound for the POD-DGTD method, which is adapted to our particular modeling and discretization settings. A snapshot choosing rule aiming at keeping the error estimate close to a target selection error tolerance is proposed, and an incremental singular value decomposition (ISVD) algorithm is used to update the SVD on-the-fly when a new snapshot is available. For problems relevant to nanaophotonics, in [8], we extend the POD-DGTD method to deal with the system of 3d time-domain Maxwell equations coupled a Drude dispersion model, which is employed to describe the interaction of light with nanometer scale metallic structures. Finally in [10] we extend the performance of the non-intrusive ROM method introduced in [6] by combining the POD method with the Dynamic Mode Decomposition (DMD) method. Here, the DMD method allows to achieve a further dimension reduction of the ROM method, and is used to predict the reduced-order coefficient vectors for future time instants on the previous training parameter set. Moreover, the radial basis function (RBF) is employed for approximating the reduced-order coefficient vectors at a given untrained parameter in the future time instants, leading to the applicability of DMD
method to parameterized problems.

References

[1] D. Amsallem and U. Hetmaniuk. Error estimates for Galerkin reduced-order models of the semi-discrete wave equation. *ESAIM Math. Model. Numer. Anal.*, 48(1):135-163, 2014.

[2] J. Burkardt, M. Gunzburger, and H. Lee. POD and CVT-based reduced-order modeling of Navier-Stokes flows. *Comput. Methods Appl. Mech. Engrg.*, 196(1):337-355, 2006.

[3] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno. Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes. *ESAIM: Math. Model. Numer. Anal.*, 39(6):1149-1176, 2005.

[4] J.S. Hesthaven, G. Rozza, and B. Stamm. *Certified reduced basis methods for parametrized partial differential equations*. SpringerBriefs in Mathematics. Springer, 2016.

[5] J.S. Hesthaven and T. Warburton. *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[6] K. Li, T.-H. Huang, L. Li, and S. Lanteri. Non-intrusive reduced-order modeling of parameterized electromagnetic scattering problems using cubic spline interpolation. *J. Sci. Comput.*, 87(2):1-29, 2021.

[7] K. Li, T.-Z. Huang, L. Li, and S. Lanteri. POD-based model order reduction with an adaptive snapshot selection for a discontinuous Galerkin approximation of the time-domain Maxwell's equations. *J. Comput. Phys.*, 396:106-128, 2019.

[8] K. Li, T.-Z. Huang, L. Li, and S. Lanteri. Simulation of the interaction of light with 3-D metallic nanostructures using a proper orthogonal decomposition-Galerkin reduced-order discontinuous Galerkin time-domain method. *Numer. Methods Partial Differ. Eq.*, 39(2):932-954, 2023.

[9] K. Li, T.-Z. Huang, L. Li, S. Lanteri, L. Xu, and B. Li. A reduced-order discontinuous Galerkin method based on POD for electromagnetic simulation. *IEEE Trans. Antennas Propag.*, 66(1):242-254, 2018.

[10] K. Li, Y. Li, L. Li, and S. Lanteri. Surrogate modeling of time-domain electromagnetic wave propagation via dynamic mode decomposition and radial basis function. *J. Comput. Phys.*, 491(112354), 2023.

[11] Z. Luo and J. Gao. A POD reduced-order finite difference time-domain extrapolating scheme for the 2D Maxwell equations in a lossy medium. *J. Math. Anal. Appl.*, 444(1):433-451, 2016.

[12] R. Mancini and S. Volkwein. An inverse scattering problem for the time-dependent Maxwell equations: nonlinear optimization and model-order reduction. *Numer. Linear Algebra Appl.*, 20(4):689-711, 2013.

[13] V. Pereyra. Wave equation simulation using a compressed modeler. *Amer. J. Comput. Math.*, 3(3):231-241, 2013.

[14] V. Pereyra and B. Kaelin. Fast wave propagation by model order reduction. *Electron. Trans. Numer. Anal.*, 30(7):406-419, 2008.

[15] L. Sirovich. Turbulence and the dynamics of coherent structures. Part 1: coherent structures. *Quart. Appl. Math.*, 45(3):561-571, 1986.

[16] K. Yee. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. *IEEE Trans. Antennas Propag.*, 14(3):302-307, 1966.

[17] J. Yu, C. Yan, and M. Guo. Non-intrusive reduced-order modeling for fluid problems: a brief review. *Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering*, 233(16):5896-5912, 2019.