Yanlai Chen
email: yanlai.chen@umassd.edu

Jan S Hesthaven
email: jan.hesthaven@brown.edu

Yvon Maday
email: maday@ann.jussieu.fr

Jerónimo Rodríguez
email: jeronimo.rodriguez@usc.es

Xueyu Zhu
email: xueyuzhu@brown.edu

Certified Reduced Basis Method for Electromagnetic Scattering and Radar Cross Section Estimation

Keywords: Reduced Basis Method, Electromagnetic Scattering, Radar Cross Section, Empirical Interpolation Method

We study nontrivial applications of the reduced basis method (RBM) for electromagnetic applications with an emphasis on scattering and the estimation of radar cross section (RCS). The method and several extensions are explained with two examples with different characteristics. Parameters that are allowed to vary within the model include frequency, incident angle and measurement angle as well as the geometry of the scatterers. With appropriate applications of the empirical interpolation method (EIM), transformation of the domain, configuration of perfectly matched layer, exponential convergence of the reduced basis solution over the entire parameter domain is achieved. Moreover, we demonstrate that this approach allows for the effective capture of the critical behavior, in this case through shapes that minimize scattering. This further highlights the robustness and quality of the greedy approximation and the reduced basis method approach.

Introduction

In many applications such as computational optimization, control and design, one is required to rapidly yet accurately predict some quantities of interest under the variation of a set of parameters ν ∈ D ⊂ R p . Here, s e (ν) := ℓ e (u e (ν); ν) ∈ C is an output of interest with ℓ e being a functional and u e is the solution of a parametrized partial differential equation (PDE). Let us denote the PDE in its weak form a e (u e (ν), v; ν) = f e (v; ν), ∀ v ∈ X e , [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] where a e and f e are bilinear and linear forms, respectively, associated with the PDE. We shall denote X e as the space of the exact solution u e . In practice, one could use a finite element (FE) discretization to approximate the solution u N (ν) ≃ u e (ν). Hence, given ν ∈ D ⊂ R P , find u N

(ν) ∈ X N satisfying a N (u N (ν), v; ν) = f N (v; ν), ∀v ∈ X N . ( 2
)
An approximate value of the output of interest can be thus computed by s N (ν) := ℓ N (u N (ν); ν) ∈ C. Here X N is the finite element space approximating X e with dim(X N )≡ N , a N (•, •; •) and m N (•; •), m ∈ {f, l} are computable approximations of a(•, •; •) and m(•; •), m ∈ {f, l}, respectively. We assume u N (ν) provides a reference solution, referred to as the truth approximation, that is accurate enough for all ν ∈ D. To ensure this, one must usually choose a very large N . As a result, the cost to solve for the truth approximation is likely very high, especially when the solution is needed for many instances of ν, in which case the cost becomes prohibitive.

The RBM, introduced in [START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF][START_REF] Fink | On the error behavior of the reduced basis technique for nonlinear finite element approximations[END_REF], offers a particularly well suited solution to the challenges of this "many-query" scenario, aiming to improve the simulation efficiency and reduce the overall computational cost. A fundamental observation and assumption utilized by RBM is that the parameter dependent solution u e (ν) is very often not simply an arbitrary member of the infinite-dimensional space X e associated with the PDE, but rather it evolves on a lower-dimensional manifold induced by the parametric dependence. Therefore, one can expect that as ν (∈ D ⊂ R q ) varies, the set of all solutions u e (ν) can be well approximated by a finite and low dimensional vector space. In particular, the RBM method proposes to approximate the solution of the PDE for an arbitrary value of the parameter ν ∈ D as a linear combination of solutions of the same PDE for other adequately chosen parameters ν i , i ∈ {1, . . . , N }, i.e., u N (ν) = ∑ N i=1 c i (ν)u N (ν i ). The truth approximation problem is then replaced by the following problem: find u N (ν) ∈ X N such that

a N (u N (ν), v; ν) = f N (v; ν), ∀v ∈ X N , (3)
where

X N = span { u N (ν i ), i ∈ {1, . . . , N } } .
As can be expected, the choice of the initial set of parameters used to compute the basis functions appropriately is crucial for the method to succeed. This is guided by the combination of a rigorous a posteriori error estimators, also used to certify the quality of the approximations, and a judicious greedy algorithm. We refer to [START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF][START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF] and also the review paper [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF] and the extensive reference therein for detail. Theoretically, exponential a priori convergence result of the reduced basis approximation is confirmed for a one dimensional parametric problem [START_REF] Maday | A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations[END_REF]. More recently, exponential convergence of the greedy algorithm for continuous and coercive problems with parameters in any dimension has been established in [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF], and improved in [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF].

The RBM becomes particularly valuable when the forms involved in the PDE satisfy the so-called affine assumption, namely,

a N (u, v; ν) = Qa ∑ q=1 Θ a q (ν)a N q (u, v), m N (v; ν) = Qm ∑ q=1 Θ m q (ν)m N q (v), (4)
where m ∈ {f, l}. In this case the computational strategy centers around splitting into offline and online parts. In the offline part (computations that are performed once), the reduced basis space and various elements needed for the error estimation is precomputed, at a cost scaling with N . In the online part, where the reduced basis approximation is used to approximate the solution for any new parameter, the computational cost depends only on N (the dimension of the reduced basis space) and Q m , m ∈ {a, f, l} (the so-called complexity of the forms) but not on N (the dimension of the FE space where the truth approximation is computed). Since N ≪ N this introduces a potential for a substantial computational saving.

In this paper, we consider problems of electromagnetic scattering in which a valuable output of interest is the radar cross section (RCS) [START_REF] Knott | Radar Cross Section[END_REF], measuring the reflective strength of a target when illuminated by an electromagnetic source. The design and optimization of this radar signature is naturally of great practical and tactical importance, involving both material and geometric design options. Other parameters to include are incident and observation wave angles and the frequency. In general the linear and bilinear functionals associated with class of problems do not satisfy the affine assumption [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF]. In consequence, an additional approximation based on the EIM will be developed and applied [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. Reduced basis method has been applied to 3D electromagnetic scattering problem. Inspired by the subdomain residual method, the authors of [START_REF] Pomplun | Accelerated a posteriori error estimation for the reduced basis method with application to 3D electromagnetic scattering problems[END_REF] proposed a faster, but non-rigorous, error estimate to deal with affine expansions with very large number of terms.

What remains of this paper is structured as follows. In Section 2, we state the general problem setting before introducing two specific scattering problems that we shall use as benchmarks to illustrate the capabilities of the reduced basis method for complex scattering and radar cross section prediction applications. These examples will involve parametric variations of frequency, incident and observation angles and geometric variations, or combinations of these. In Section 3, we briefly explain the computational strategy and offer some details on the truth approximation solver and the greedy algorithm used to determine the parameters enabling the construction of the reduced basis space [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF][START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF]. The two test problems do not satisfy the affine assumption and we develop an empirical interpolation method (EIM) that facilitates a reduced problem satisfying (4), ensuring an efficient off-line on-line computational strategy. Numerical results illustrating the superior performance of reduced basis method for the two non-trivial examples are presented in Section 4. Finally, we offer some concluding remarks in Section 5.

Presentation of the problem

In this section, we describe the setup and formulation of two non-trivial benchmark cases when the radar cross section is sought as the output of interest. For the first case, the angular frequency and angle of the incident wave are the parameters, while in the more challenging second case the parameter is the shape of the scatterer. As we shall see, this necessitates the nontrivial geometric transformation and the corresponding application of perfectly matched layer (PML).

In both cases, the electromagnetic waves are TM-polarized, i.e., the electric and magnetic fields satisfy E = (0, 0, E z ) and H = (H x , H y , 0) in the Maxwell's equation. As a source, we consider an incident plane wave

E inc (⃗ x) = e iω⃗ x• ⃗ d .
Here, ω is the angular frequency and

⃗ d ≡ ⃗ d(θ i ) := (-cos(θ i ), -sin(θ i ))
is the direction of propagation that depends on θ i , the angle of the incident wave.

The scattered field, generated by the incident wave, satisfies a set of equations (we assume ε = µ = 1 for the sake of simplicity):

iωE z = ∂H y ∂x - ∂H x ∂y in Ω, (5a)
iωH x = - ∂E z ∂y in Ω, ( 5b
)
iωH y = ∂E z ∂x in Ω, ( 5c
)
E z = -E inc on Γ. ( 5d
)
The system is completed by requiring the scattered field to satisfy the so-called Silver-Müller radiation condition at the infinity. Here, Ω is the domain outside of the scatterer and Γ is the boundary of the scatterer. To enable tangible computations, an exterior boundary sufficiently far away from the scatterer, Γ f ar , is enforced to reduce the infinite domain to a finite one and artificial boundary conditions will be applied on the exterior boundary Γ f ar . In this work, either the fields (E z , H x , H y ) are assumed to satisfy the first order absorbing boundary condition (on the first problem) as

(H × n) • ẑ = √ ε µ E z , on Γ f ar , (6)
or we apply PML [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] on an annulus just inside of Γ f ar (on the second problem).

Given the scattered field H and E, we can calculate the radar cross section associated with this scattering process. We define the surface currents

J = n × H, and M = -n × E, F (ω, θ i , θ r ) = √ e i π 4 √ 8πω (-ωµẑ • J -ωẑ × M • r) e iωr•r ′ dC. ( 7
)
Here r, the unit vector for the observation direction, is equal to (cos(θ r ), sin(θ r )) with θ r being the observation angle. r′ is the vector from scatterer to the integration contour which can be any closed contour surrounding the scatterer.

With the linear output functional (7) computed, the RCS (also called bi-static RCS) is defined as

RCS(ω, θ i , θ r ) = 2π |F (ω, θ i , θ r )| 2 |E inc | 2 . ( 8
)
A case of particular interest in practice is the monostatic RCS defined by RCS(ω, θ i , θ r = θ i ) . Often, these magnitudes (quotients of powers) are measured in decibels (dB).

Open Cavity Benchmark

In the first problem to be considered, the scatterer is assumed to be a perfect electric conductor (PEC) open cavity given by (see Figure 1 left)

Γ := [-1, 1] × {-1} ∪ [-1, 1] × {1} ∪ {-1} × [-1, 1].
The unbounded nature of the computational domain is enabled by a first order absorbing boundary condition (often called the Silver-Müller boundary condition) at the boundary Γ f ar = {⃗ x ∈ R 2 such that |⃗ x| = 3}. The output of interest is the radar cross section given by [START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF]. For this application we consider three parameters (ω, θ i , θ r ) ∈ D = [π/2, 5π/2] × [0, 2π] × [0, 2π], namely, the angular frequency, the angle of incidence and the angle of observation.

This is an interesting benchmark as a example of cavity scattering, being an application area of considerable interest, and because it is well known that even minor changes in the scattering of a cavity can lead to dramatic changes in the scattered fields due to a variety of electromagnetic phenomena,e.g., quasi resonant behavior.

Invisible Pacman Benchmark

The second scatterer we consider is a perfectly conducting 2D cylinder with a cut-out wedge. The configuration is illustrated in Figure 1. Here, θ W denotes the angle of the wedge and is a critical parameter for this case. We will also allow the observation angle θ r and angular frequency ω to vary. A curvilinear PML [START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF] is applied around Γ f ar to simulate the unbounded nature of the computational domain.

An interesting phenomenon of this problem is that the scattered fields and thus the RCS change dramatically with only a small change in the wedge angle θ W . In Figure 2, we show the RCS(ω = 10π, θ i = 0, θ r ) in dB when θ W = 0 • , 18.5 • and 21.5 • . We note in particular the 22dB drop in the monostatic scattering (when θ r = 0) with just 3 degrees of change in θ W rendering the Pacman substantially less visible. One of the challenges in this work is to demonstrate that reduced basis methods are capable of capturing such critical behavior across variations in geometry and frequency. Three cases with different wedge angles θ W are plotted. In the three cases θ i = 0 and ω = 10π.

Geometric Transformation

To apply the RBM with computational advantage, the parametric problem needs to be written on a fixed domain. To address this, we transform the problem on a domain depending on the angle θ W to a parametric problem on a fixed domain through a map from a reference domain with a fixed θ r W to a physical domain with a general θ W .

Let (x, y) and (ρ, θ) be rectangular and polar coordinates on the reference domain, and (x, ŷ) and (ρ, θ) be those on the physical domain. The mapping of the domain is achieved through a sequence of mappings with the first and last being the standard transformation between polar and rectangular coordinates:

( x y ) → ( ρ θ ) → ( ρ = ρ θ = f m (θ) ) → ( x y ) . ( 9
)
The function f m (•) in the second mapping is a continuous, piece-wise smooth function that is one-to-one and such that

f m ([-π, -θ r W ]) = [-π, -θ W ], f m ([-θ r W , θ r W ]) = [-θ W , θ W ] and f m ([θ r W , π]) = [θ W , π].
Let us first discuss the systems of equations on the physical domain. We omit the • on all the variables and operators for simplification of notations. Under the assumption of µ = 1, ϵ = 1, the TM-polarized system of time harmonic Maxwell's equations ( 5) can be rewritten using polar coordinate as

       iωE z = 1 ρ ( ∂(ρH θ ) ∂ρ - ∂Hρ ∂θ
) ,

iωH ρ = -1 ρ ∂Ez ∂θ , iωH θ = ∂Ez ∂ρ . ( 10
)
Next, we incorporate the perfectly matched layer to truncate the computational domain. To achieve this, we follow [START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF] for the derivation under curvilinear coordinate as a change-of-variable technique. For this, we define

ρ = { ρ - ∫ ρ a iσ(s) ω ds ρ ≥ a, ρ ρ < a.
Applying the PML layer is then equivalent to replacing all ρ in (10) by ρ. Here, σ is a piecewise quadratic C 1 -function of ρ (constant along the θ-direction).

It is identically zero in the non-PML region and monotonically increasing along the ρ-direction from the PML/non-PML interface to the exterior boundary. To simplify the system, we let

σ = { 1 ρ ∫ ρ a σ(s)ds ρ ≥ a, 0 ρ < a,
and obtain, after some manipulation, the polar coordinate system

       -(iω + σ)E z + σ-σ iω+σ ∂H θ ∂ρ = 1 ρ ( ∂Hρ ∂θ -∂(ρH θ ) ∂ρ ) , (iω + σ)H ρ = -1 ρ ∂Ez ∂θ , (iω + σ)H θ = ∂Ez ∂ρ . ( 11
)
Converting to rectangular coordinate we recover, after some algebraic manipulations,

           -(iω + σ)(iω + σ)E z + A∇ • ( H x H y ) = 0, ∇E z + A ( H x H y ) = 0. (12)
Here,

A = ( (σ -σ) sin θ cos θ -iω -σ sin 2 θ -σ cos 2 θ iω + σ cos 2 θ + σ sin 2 θ (σ -σ) sin θ cos θ ) = (σ -σ) ( sin θ cos θ sin 2 θ -cos 2 θ -sin θ cos θ ) + ( 0 -iω -σ iω + σ 0 )
Finally, we obtain the parametric system on the reference domain to which we apply the reduced basis method

           -(iω + σ)(iω + σ)E z + A J -T ∇ • ( H x H y ) = 0, J -T ∇E z + A ( H x H y ) = 0. ( 13
)
Here, J is the Jacobian matrix of the mapping, ( x y

) → ( x y
) ,

given as

J = ∂( x, y) ∂(x, y) = 1 x 2 + y 2 ( x x + f ′ m y y xy -f ′ m x y x y -f ′ m xy y y + f ′ m x x
) .

The reduced basis method

In the following, we outline the procedure applied to build the reduced basis.

The truth approximation

We first describe the finite element method to solve the parametric systems ( 5) and ( 13) that we want to approximate for a particular parameter choice. We let

u =   E z H x H y   v =   v 1 v 2 v 3  
and, for a given a mesh T h , we use a discontinuous Galerkin method [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF] to solve both equations as in [START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF] but without the elimination of E z since this can not be done when the first order absorbing boundary condition or the PML are present. Here, we omit the details (see e.g. [START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF][START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[END_REF] for the details) of the formulation and simply denote the resulting bilinear form as

a N (•, •; θ W ).
The truth approximation for the scattered field u N (ν) is then the solution of

a N (u N (θ W ), v; θ W ) = -a N (u I , v; θ W ), ∀ v ∈ X N . ( 14
)
Here, we define the following finite element space

X N = {v ∈ ( L 2 (T h ) ) 3 : for all elements K ∈ T h , v|K ∈ (P k (K)) 3 }, ( 15
)
and u I is the incident field (e iω⃗ x• ⃗ d , -d 2 e iω⃗ x• ⃗ d , d 1 e iω⃗ x• ⃗ d ) T .
Remark 1. The volume terms on the right hand side of [START_REF] Maday | A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations[END_REF] involves the incident field and can be removed in practice since it automatically satisfies the three first equations in [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis[END_REF] at the continuous level.

Construction of the reduced basis. The greedy algorithm

The next step in the approximation is to replace the approximation space (15) for the DG formulation ( 14) by a lower dimensional space of the form (see (3))

X N = span { u N (ν i ), i ∈ {1, . . . , N } } . ( 16
)
As indicated previously, the selection of these parameters ν i , i ∈ {1, . . . , N } is crucial to achieve good approximation properties of the reduced space. This selection is performed through a greedy algorithm, presented for completeness in Figure 3 (see [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF][START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF] for the details). The key point of the algorithm is line 3 where the truth approximation, i.e., solution of the problem ( 14) or ( 2), is compared with the reduced basis solution, recovered through (3)). The value of the parameter for which a bound for this error is maximized is taken to be the next parameter values. In practice, the set D is replaced by a discrete training set Ξ train . Moreover, the a posteriori error estimator ∆ N (ν), evaluated at low computational cost , is assumed to be a rigorous and optimal upper-bound of the error [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF][START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF].

1 ≤ ∆ N (ν) ∥u N (ν) -u N (ν)∥ L 2 ≤ C (constant). X 1 = span { u N (ν 1 ) } (with ν 1 arbitrarily chosen) ; 1 for N = 1, . . . , N max do 2 ν N = argmax ν∈D ∆ N (ν); 3 X N = span { u N (ν i ), i ∈ {1, . . . , N } } ;

The empirical interpolation method

The bilinear and linear forms in our benchmark problems both contain functions that are non-affine with respect to the parameter. For an example, the linear form involved in the computation of the output (7) contains terms of the form ψ(x) g(x; ν) dC, [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] where

g(⃗ x, ν) ≡ exp(-iω ⃗ k θr • ⃗ x) = cos(ω(x cos θ r + y sin θ r )) g1(x,ν) -i sin(ω(x cos θ r + y sin θ r )) g2(x,ν)
.

(18) The affine assumption enables the offline-online computational strategy and is hence a key part of computational efficiency of the scheme. When violated, one idea is to approximate the function g(•; •) in ( 17) by a linear combination of functions with separate dependence of the parameter variable and the spatial variable, i.e.,

g(x; ν) ≈ g M (x; ν) = M ∑ m=1 φ M m (ν)q m (x),
using the empirical interpolation method [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. This leads to an approximate linear form clearly satisfying the affine assumption. The empirical interpolation method is based on i) the expansion of the function g(•; •) over a space accounting for the parameter dependence and ii) a robust interpolation procedure. For completeness of the paper we briefly present these two steps [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF].

Construction of the approximation space.

To build an approximation space well adapted to the parameter dependence we consider

W g M := span (g(•; ν g m ), 1 ≤ m ≤ M ) .
The set of parameters S g M = {ν g m , 1 ≤ m ≤ M } is built hierarchically using the greedy algorithm in Figure 4.

Given a set of parameters S g m (and its associated approximation space W g m ), the next parameter ν g m+1 is selected as the one where the difference between g(•, ν) and its best approximation, measured through a suitable norm ∥ • ∥,

g * M (•, µ) := argmin z∈W g M ∥g(•, µ) -z∥ , ( 19
)
is maximized in W g m In practice, the set D on line 4 is replaced by a training set Ξ g train ⊂ D. In addition, we use a L 2 -type norm in such a way that computing [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF] amounts to the solution of a linear system.

S g 0 = ∅; W g 0 = {0}; for M = 1, . . . , M max do ν g M = argmax ν∈D g * M -1 (•, ν) -g(•, ν) S g M = S g M -1 ∪ {ν g M }; W g M = span (g(•, ν), ν ∈ S g M );
end Setting up the interpolation procedure. The interpolation points

I g 0 = ∅; 1 for M = 1, . . . , M max do 2 Compute g M -1 (•; ν g M ) ∈ W g M -1 interpolating 3 g(•; ν g M ) on the points I g M -1 ; Compute r M (•) = g(•; ν g M ) -g M -1 (•; ν g M ); 4 x g M = arg esssup x∈Γ |r M (•)|, I g M = I g M -1 ∪ {x g M }; 5 q M (•) = r M (•)/r M (x g M );
I g M = {x g m , 1 ≤ m ≤ M } ,
that are selected with the algorithm in Figure 5 lead to an efficient and robust interpolation. In addition, this procedure provides a new basis q m (•), m ∈ {1, • • • , M } for the interpolation space such that the computation of the interpolation coefficients φ M m (ν) amounts to solving a triangular linear system. Notice that the family of interpolation points is hierarchical. Moreover, the algorithm is well defined unless r M (•) ̸ = 0. This would, however, imply that the interpolation procedure would be exact for all ν ∈ D when using M -1 terms on the expansion.

In practice, the search of the interpolation point x g M in line 5 is not performed over Γ but on a finite number of points included in Γ (usually the nodes of a mesh). Moreover, the two algorithms in Figures 4 and5 are interwined, see [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. That is, one more basis function q M (•) is added in each step together with the next interpolation point ν g M . Lastly, the best approximation error in ( 19) is replaced by interpolation error as shown by line 3 of Figure 5.

Numerical Results

In this section, we discuss the numerical results for the two benchmark problems, seeking to confirm the high efficiency and good accuracy of the reduced basis method applied to complex scattering problems.

Open Cavity Benchmark

The governing equations are solved using the discontinuous Galerkin method with third order polynomials over a nonconforming mesh locally refined toward the corners and tips of the scatterer. The total number of unknowns is 7520. In Figure 6, we illustrate the fields when ω = 2π and θ i = π/4.

The linear functional associated with the output of interest does not satisfy the affine assumption as discussed in Section 3.3). For this reason we begin by checking the performance of the EIM used to interpolate the function g 1 (•; •) in (18) (the results for g 2 (•; •) are similar) to recover an affine linear form.

Empirical Interpolation Method

The parametric space is in this case two dimensional with ν = (ω, θ r ) since θ i only appears in the right-hand side of the PDE. Thus we need to interpolate g 1 (•; •) on the boundary Γ where the output is computed. The approximation space for this interpolation is obtained using the L 2 (Γ)-norm. Hence, the computation of the best approximation (see [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF]) amounts to solving a linear system with size equal to the number of degrees of freedom on the boundary Γ. The interpolation points are computed using the algorithm in Figure 5, where the maximum is sought (for simplicity) among the nodes used by the approximation and not in the entire Γ.

1-D parameter experiment. To begin, we only consider one parameter, the observation angle θ r ∈ D := [0, 2π], while keeping the angular frequency fixed. The empirical interpolation method amounts in this case to the approximation of a plane wave with angular frequency ω and an arbitrary measurement angle θ r using a finite number of plane waves with the same angular frequency. The results are displayed in Figure 7 for ω = 0.5π and ω = 2.5π. The greedy search performed on the line 4 of the algorithm in Figure 4 is done over the training set Ξ train = {2πk/360} 360 k=0 . The maximum dimension of the discrete space M max is the smallest integer such that the best approximation error -14 uniformly on Ξ train . In the first column of Figure 7 we illustrate the parameters S g Mmax that are selected by the algorithm, hence defining the approximation space. In the second column we plot the corresponding interpolation points I g Mmax that are selected by the algorithm in Figure 5 for both configurations.

ϵ * M (ν) := ∥g(•, ν) -g * M (•, ν)∥, is less than ϵ = 2.5 × 10
To test the quality of the approximation space W g M and the greedy interpolation procedure, we consider a random set of 1000 parameters Ξ rnd uniformly distributed in D. We compute max

ν∈Ξ rnd ∥g(•; ν) -g * M (•; ν)∥ and max ν∈Ξ rnd ∥g(•; ν) -g M (•; ν)∥, ( 20
)
and show these in the third column of Figure 7. Note that the first quantity is the maximum error committed by the best approximation (projection) over Ξ rnd while the second quantity is the maximum error in which case the interpolation is used over Ξ rnd . Both quantities show exponential decay when the number of elements is increased. Moreover, the quality of the projection and the interpolation is comparable, implying that the interpolation points are well selected and the interpolation procedure is robust. As expected, for a higher angular frequency we need more terms in the approximation to achieve a similar accuracy. 2-D parameter experiment. We interpolate the function g 1 (•; •) for the case where both (ω, θ r ) varies on the set D = [0, 3π] × [0, 2π]. In this case the empirical interpolation method approaches a plane wave with arbitrary (ω, θ r ) in the given range by a linear combination of a finite number of plane waves. The results are displayed in Figure 8. On the left we represent the parameters selected by the greedy process on the line 4 of the algorithm in Figure 4, run over the training set Ξ train = {2πk/360} 360 k=0 × {3πk/100} 100 k=0 . Note that the high frequencies (where the function oscillates more) are preferred. In the middle we represent the interpolation points selected by the algorithm in Figure 5.

To test the quality of the approximation space and the interpolation points we evaluate the algorithm over a set of 40000 parameters randomly chosen in the set [0, 2π] × [0, 3π], and compute the quantities in [START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF]. The results are plotted in Figure 8 on the right. We observe that with 15 terms the error begins to decay exponentially. With about 25 terms a 10 -5 accuracy is achieved. We also observe that the larger the range of frequencies is, the higher the dimension of the approximation space is needed. Hence, in applications where large ranges are required, a computational strategy for splitting this range into smaller ones to control overall online cost is interesting and perhaps even essential.

In the next section, when coupling the EIM with the RBM, we use a large number of terms on the expansion (M = 35) in order to ensure a good approximation uniformly across the parameter space.

RBM with EIM

Since the non-affinity only resides in the linear forms for this problem, we only need to apply the empirical interpolation method on the linear forms. Once we have obtained an affine expansion of the linear forms, the next step of the offline computations consists of computing the reduced basis that will be used in the online part. We recall that the output of interest depends on three parameters (ω,

θ i , θ r ) ∈ D = [π/2, 5π/2] × [0, 2π] × [0, 2π].
We start by splitting the parameter domain [START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF] into three sub-domains (corresponding to low frequencies, medium frequencies and high frequencies)

D = ∪ 3 i=1 D i where D 1 = [ π 2 , 3π 2
] ×[0, 2π] 2 , D 2 = [ 3π 2 , 2π ] ×[0, 2π] 2 and D 3 = [ 2π, 5π 2
] ×[0, 2π] 2 .
This allows us to use different discretization spaces for the truth approximation (for small values of ω a coarser mesh could be used) depending on the parameter sub-domain. We apply the greedy algorithm 3 over the train sets

Ξ train 1 = { π 2 + nπ 2 5 } 2 5 n=0 × { 2lπ 2 5 } 2 5 l=0 × { 2mπ 2 5 } 2 5 m=0 ⊂ D 1 , Ξ train 2 = { 3π 2 + nπ 2•2 5 } 2 5 n=0 × { 2lπ 2 6 } 2 6 l=0 × { 2mπ 2 6 } 2 6 m=0 ⊂ D 2 , Ξ train 3 = { 2π + nπ 2•2 6 } 2 6 n=0 × { 2lπ 2 7 } 2 7 l=0 × { 2mπ 2 7 } 2 7 m=0 ⊂ D 3 .
Note that the problem is more sensitive to parameter changes for large values ω. The density of the train set is selected accordingly. In addition, the construction of the reduced basis is independent on θ r , which appears only in the linear form associated to the output. This allows us to run the greedy algorithm over twodimensional sets.

In Figure 9, we show the parameters selected by the greedy algorithm in the three cases (left), and the evolution of the maximum over the train sets Ξ train i , i ∈ {1, 2, 3} of the a posteriori error estimator (right). We observe exponential convergence of the a posteriori error estimator and see, as expected, that the rate of convergence is lower when we consider larger parameter sets or high values for ω.

Let us now test the quality of the reduced basis space in the online stage.

θ i ) is in [π/2, 3π/2] × [0, 2π] (in blue), [3π/2, 2π] × [0, 2π] (in red) and [2π, 5π/2] × [0, 2π
] (in green). On the left we represent the parameters selected by RBM (the larger the markers the earlier they have been chosen). On the right is shown the evolution of the maximum over the train sets Ξ train i , i ∈ {1, 2, 3} of the a posteriori error estimator when the dimension of the reduced basis is increased. 10 and11. We start by computing the monostatic RCS in dB (i.e., we consider

MonoStatic RCS, Figures

θ r = θ i ) with (ω, θ i ) ∈ [π/2, 5π/2] × [0, 2π].
The results are shown in Figure 10 (left). The logarithm in base 10 of the a posteriori error estimator for the linear output F (•, •, •) is also plotted (right). Note that depending on the sub-domain in which the parameter falls (D 1 , D 2 or D 3 ), a different reduced basis space is used. For both plots we perform the computations over the following grid respectively. We present, in Figure 11, the monostatic RCS for ω ∈ {1.5708, 6.2832, 7.2257} computed with two different number of elements on the reduced basis. We can observe that in both cases the approximation is very good even if the output is less reliable with a low number of elements on the basis. Adding a few elements to the reduced basis space dramatically reduces the width of the confidence interval. Note that these plots correspond to three different cuts of the monostatic RCS shown in Figure 10 for the values of ω selected above.

Ξ test 1 = { π 2 + nπ 100 } 100 n=0 × { ( 2lπ 360 , 2lπ 360 ) } 360 l=0 ⊂ D 1 , Ξ test 2 = { 3π 2 + nπ 2•100 } 100 n=0 × { ( 2lπ 360 , 2lπ 360 ) } 360 l=0 ⊂ D 2 , Ξ test 3 = { 2π + nπ 2•100 } 100 n=0 × { ( 2lπ
We plot, in Figure 12, the bistatic RCS in dB for ω = 7.4693 and all possible angles of incidence and measurement and the logarithm in base 10 of the a posteriori error estimator for the output. For both plots we have performed the computations over the following grid that is not a subset of the train set. In Figure 13 we present the monostatic RCS and the bistatic RCS in dB for θ i ∈ {0.0000, 1.2566}. We observe that the convergence of the confidence region provided by the estimator rapidly shrinks when the dimension is increased. Note that these plots correspond to the restriction of the plot 12 to the solid, dashed and dotted lines. In black, the output computed with the reduced basis method. The blue region is the confidence region provided by the a posteriori error estimator. In the first row, we use 150, 150, and 180 bases respectively and the reduced basis output already provides a good approximation of the actual output. However, the confidence region for some angles is rather wide. In the second row, we enrich the reduced basis (from 150 to 180 and from 180 to 220) to certify the output.

Ξ test = {7.4693} × { 2lπ 360 } 360 l=0 × { 2mπ 360 } 360 m=0 ⊂ D 3 ,

Invisible Pacman Benchmark

For the Pacman scattering problem, we consider two parameters: angular frequency ω and the wedge angle θ W , and apply a reduced basis method to the parametric equations ( 13) derived in Section 2.2. To emphasize the dependence on the parameters, the radar cross section (RCS), equations ( 7) and ( 8), is denoted by RCS(ω, θ W , θ i , θ r ). For an example, RCS(3π, 18.5 • , 0, θ r ) in the caption of a figure means that the plot is the RCS versus observation angle θ r with frequency 3π, wedge angle 18.5 • , and incident angle 0.

The scatterer has a fixed diameter (5 wavelengths when ω = 10π) and we have used a curvilinear PML [START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF] of fixed thickness, corresponding to 2.5 wavelengths away from the scatterer, in an annulus of width 1.5 wavelengths when ω = 10π. The PML has a maximum damping coefficient 60.

Truth approximation

We begin by studying some aspects of the truth solver for the problem (13). We fix θ r W = 18.5 • , the reference angle for the wedge, throughout the paper. When θ W = 0, we have the exact solution to compare our result with. Here, we solve the scattering problem for a cylinder with two frequencies: 3π, and 10π and calculate the RCS. We see, from Figure 14, that the finite element solver is highly reliable with relative error on the level of 10 -3 for higher frequency and 10 -4 for lower frequency. To test the efficiency of this PML applied to this problem, we compare the RCS with that obtained by using only a first order absorbing boundary condition, Silver-Müller boundary condition. It is evident from Figure 15 that the application of the curvilinear PML improves the results. Here, the reference is a proven solver for the Pacman scatterer on a square domain with rectilinear PML and no geometric mapping. Next, we study the effect of three possible forms of mapping f m (•) in ( 9). We show here the resulting monostatic RCS in dB from the parametric solver with θ W = 21.5. In Figure 16, "P 2 P 1 P 2 " indicates that the polynomial degrees of the mapping on [-π, -θ r W ], [-θ r W , θ r W ] and [θ r W , π] are 2, 1, and 2 respectively; similarly "P 2 P 2 " degrees on [-π, 0] and [0, π]; "P 1

P 1 P 1 " degrees on [-π, -θ r W ] ∪ [θ r W , π] and [-θ r W , θ r W ]. Naturally, we require f m (-π) = -π, f m (-θ r W ) = -θ W , f m (θ r W ) = θ W , f m (π)
= π for continuity and also demand C 1 -regularity if there are additional degrees of freedom allowing us to do so (first and second cases). We observe that the "P 1 P 1 P 1 " mapping yields the best overall result. The reason for this is found in the observation that this map ensures the most regularity in the physical domain. Hence, we shall use this mapping in the truth solver.

RBM for one parameter without EIM

The interesting phenomenon described in Section 2.2 and exemplified by Figure 2 suggest the study of the dependence of scattering on varying θ W . Therefore, we fix ω = 10π, treat θ W as our parameter in this section and develop a reduced basis method to the Pacman scattering problem. Here, we set the parameter domain for θ W to be [8.5 • , 28.5 • ]. Contrary to the problem treated in Section 4.1 where only the linear forms were non-affine, the bilinear form for this problem is also non-affine thanks to its dependency on θ W through the geometric transformation. Note that only the parameter θ W is responsible for the non-affinity since the bilinear form is affine in ω (see equation [START_REF] Lovgren | The reduced basis element method: Offline-online decomposition in the nonconforming, nonaffine case[END_REF]). EIM will be applied in the next subsection where it is shown that it works well for this problem and, as expected, does not degrade the quality of the RBM solution. However, the results in this subsection are obtained without EIM for simplicity.

The procedure of building the reduced basis space is standard and we omit the details. The accuracy of the reduced basis solution is certified by the residual-based a posteriori error estimate. The inf-sup number will not become zero since this is an open problem. The resulting error estimate is cheap to obtain online. It also guides the selection of the parameters and thus the building of the reduced basis space in the greedy algorithm in 3. See e.g. [START_REF] Chen | Certified reduced basis methods and output bounds for the harmonic maxwell's equations[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics[END_REF] for details. For this example, this algorithm has been run over the train set Ξ train of size 257.

In Figure 17 we illustrate the 29 parameter instances that the reduced basis method picks together with their orders. The first point is picked randomly to start the greedy algorithm. It is also interesting to note that, in order to capture the critical angles, the greedy algorithm does not necessarily select points clustering around those angles. This underscores the strength of the method and that its ability to capture the critical phenomenon does not rely on sampling around the critical points in the parameter domain.

We plot, in Figure 18, the history of convergence of the RB solutions for the worst of the 120 randomly selected parameter values. Exponential convergence is clearly observed. Moreover, the error estimate decreases exponentially with about the same rate. Thus, the effectivity index is roughly constant across the many magnitudes of decrease in the RBM error. With only 20 bases, we obtain an accuracy at the level of 10 -5 . Using the underlying mesh of our numerical experiment, instead of solving systems of dimension above 200, 000 × 200, 000, we only need to solve problems of dimension 60 × 60.

Next, we study the monostatic scattering as a function of the wedge angle, that is 10log 10 (RCS(10π, θ W , 0, 0)). In Figure 19, we plot the truth approximation (in black) and the reduced basis approximations with up to only 11 basis elements. We see that the curve with 11 bases and the truth cuve are "roughly identical" to naked eyes. We also observe the dramatic effect of adding two correct samples on the output.

As described in section 2.2 and shown here by the monostatic scattering curve in Figure 19, there are several critical angles producing minimal monostatic scattering. To visualize this phenomenon, we show, in Figure 20, the electric fields corresponding to θ W = 14.3 • , 18.5 • and 21.5 • degrees. The difference is easily noticeable around the wedge.

Finally, the error estimate provides bounds for the output computed by the RBM. In Figure 21, we plot the RBM output together with the error bounds for The 29 θ W 's the greedy algorithm of the RBM picks to build the RB space. Top: the higher the vertical line, the earlier that point was picked. Bottom: the points scattered on the monostatic scattering curve 10log 10 (RCS(10π, θ W , 0, 0)), the larger the marker, the earlier it is selected. three different numbers of bases. We see that RBM output with just 17 bases captures all essential behavior with a high level of confidence.

RBM for two parameters with EIM

In this section we consider two parameters in the invisible Pacman problem: the angular frequency ω and the wedge angle θ W . As pointed out in Section 4.2.2, both the bilinear form and the linear forms for this problem are non-affine due to their dependency on θ W through the geometric transformation. The non-affinity attributed to θ W is handled by EIM. To do that, every function depending on θ W in a non-affine way (for example all the components of the Jacobian matrix associated to the transformation) are expanded using a onedimensional empirical interpolation method [START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. EIM is applied in a similar situation in [START_REF] Quarteroni | Numerical solution of parametrized Navier-Stokes equations by reduced basis methods[END_REF][START_REF] Lovgren | The reduced basis element method: Offline-online decomposition in the nonconforming, nonaffine case[END_REF]. The number of expansion terms M is taken in {1, 2, 3}. These numbers are quite small but note that it is applied to one-dimensional functions and there are more than 60 such functions in the linear and bilinear forms. This results in Q a and Q m around 200. In this case, a L ∞ (Ω)-norm best approximation is considered. It is important to note that the RBM procedure still needs to be applied to a two dimensional parameter domain. We first test the method for (ω, θ W ) ∈ [3π, 5π] × [8.5 • , 28.5 • ]. In Figure 22 we show the parameter points selected by the RBM, the worst-case convergence of the approximate solution (over the train set) of the RBM coupled with EIM and its estimate. Exponential convergence is observed. Note that the green curve corresponding to the case M = 1 is cut at the starting point of the plateau where the interpolation error begins to dominate. Next, we compute approximate RCS(ω, θ W , 0, 0) is plotted in Figure 26 with 22 bases used. A slice of this picture with ω = 10π shows comparable accuracy with the "RBM with 11 Bases" curve in Figure 19. Although the complexity of the geometric mapping and the system being in a vector form do lead to a large number of terms for the affine decomposition, the resulting RBM coupled with EIM does provide highly accurate approximations that converge exponentially to the truth approximations.

Conclusion

In this paper, we have developed and applied reduced basis method to two non-trivial electromagnetic scattering problems for which the radar signature is an important output of interest. As we have discussed in detail, this development has required the development and application of a variety of techniques such as empirical interpolation, non-trivial geometric mapping, the inclusion of perfectly matched layers into the formulation.

We have demonstrated that in spite of the need to combine these different elements to ensure computational efficiency, exponential convergence of the reduced basis solution toward the truth finite element approximation is obtained across the entire parameter domains. We also demonstrate that even with very narrow critical behavior, the greedy approach captures this behavior very well, even without requiring a dense sampling close to the critical parameter values.

Even with these advances, several challenges still lie ahead to further mature the reduced basis methods for complex scattering applications. These are mainly related to increasing the online efficiency, seeking to demonstrate substantial speedup with as much as two to three orders of magnitude as seen in other applications. The essential components to achieve this found in reducing the complexity of the affine expansion recovered by the empirical interpolations for geometrical parametric variations and the general computational challenges associated with large parameter ranges or high-dimensional parameter domains. We hope to report on such developments in future work.

Figure 1 :

 1 Figure 1: Geometry of the computational domain for the open cavity problem (left) and invisible pacman problem (right).

Figure 2 :

 2 Figure 2: Radar cross sections for the Pacman, 10log 10 (RCS), versus observation angle θr. Three cases with different wedge angles θ W are plotted. In the three cases θ i = 0 and ω = 10π.

4 end 5 Figure 3 :

 453 Figure 3: Algorithm to build the reduced basis space

9

 9

Figure 4 :

 4 Figure 4: Algorithm to build the approximation space

6 end 7 Figure 5 :

 675 Figure 5: Algorithm to compute the interpolation points

Figure 6 :

 6 Figure 6: Solution of the problem on the open cavity (from left to right: Ez, Hx and Hy) for ω = 2π and θ = π/4.

- 1 - 1

 11 (ω(x cos θ + y sin θ)). ω = 0.5 π for cos(ω(x cos θ + y sin θ)). ω = 0(ω(x cos θ + y sin θ)). ω = 2.5 π for cos(ω(x cos θ + y sin θ)). ω = 2Max. interp. error over Ξ rnd Max. proj. error over Ξ rnd

Figure 7 :

 7 Figure 7: Empirical interpolation of the function g 1 (x; (ω, θr)) = cos(ω(x cos(θr) + y sin θr)) with θr ∈ [0, 2π] when ω = 0.5π (top) and ω = 2.50π (bottom) on Γ. In the first (resp. second) column, we represent the set S g Mmax (resp. the interpolation points I g Mmax ). The larger the point is the earlier it has been chosen. Plotted in the third column are max ν∈Ξ rnd ∥g(•; ν)g * M (•; ν)∥ and max ν∈Ξ rnd ∥g(•; ν) -g M (•; ν)∥.

Figure 8 :

 8 Figure 8: Empirical interpolation of the function g 1 (x; (ω, θr)) = cos(ω(x cos(θr) + y sin θr)) with (ω, θr) ∈ [0, 3π] × [0, 2π] on Γ. In the first (resp. second) column we represent the set S g Mmax (resp. the interpolation points I g Mmax ). The larger the point is the earlier it has been chosen. Plotted in the third column are max ν∈Ξ rnd ∥g(•; ν) -g * M (•; ν)∥ and max ν∈Ξ rnd ∥g(•; ν) -g M (•; ν)∥.

  RB solutionMax of the estimator on the training setω ∈ [π/2,3π/2] ω ∈ [3π/2,2π] ω ∈ [2π, 5π/2]

Figure 9 :

 9 Figure 9: Construction of the reduced basis when (ω,θ i ) is in [π/2, 3π/2] × [0, 2π] (in blue), [3π/2, 2π] × [0, 2π] (in red) and [2π, 5π/2] × [0, 2π] (in green). On the left we represent the parameters selected by RBM (the larger the markers the earlier they have been chosen). On the right is shown the evolution of the maximum over the train sets Ξ train

Figure 10 :

 10 Figure 10: Monostatic RCS (i.e. θr = θ i ) in dB (left) for (ω, θ i ) ∈ D 1 ∪ D 2 ∪ D 3 with D 1 = [π/2, 3π/2] × [0, 2π], D 2 = [3π/2, 2π] × [0, 2π] and D 3 = [2π, 5π/2] × [0, 2π]. For each of these three regions we compute a different reduced basis. Logarithm in base 10 of the estimator for the output of interest is on the right. The computations have been done with a reduced basis of dimension equal to 190 for (ω, θ i ) ∈ D 1 , 220 for (ω, θ i ) ∈ D 2 and 250 for (ω, θ i ) ∈ D 3 .

Figure 11 :

 11 Figure 11: From left to right, monostatic RCS in dB for ω ∈ {1.5708, 6.2832, 7.2257}.In black, the output computed with the reduced basis method. The blue region is the confidence region provided by the a posteriori error estimator. In the first row, we use 150, 150, and 180 bases respectively and the reduced basis output already provides a good approximation of the actual output. However, the confidence region for some angles is rather wide. In the second row, we enrich the reduced basis (from 150 to 180 and from 180 to 220) to certify the output.

Figure 12 :

 12 Figure 12: Bistatic RCS in dB (left) for (θ i , θr) ∈ [0, 2π] 2 and ω = 7.4693, and logarithm in base 10 of the estimator for the output of interest (right). The computations are done with a RB dimension equal to 220.

Figure 13 :

 13 Figure 13: From left to right, monostatic RCS, bistatic RCS in dB with θ i = 0.0000, and bistatic RCS with θ i = 1.2566 in dB. ω = 7.4693 in all cases. Plotted in black is the output computed with the reduced basis method using 180 (top) and 220 (bottom) bases. The blue region is the confidence region provided by the a posteriori error estimator. The RB dimension is 180 (top) and 220 (bottom).

3 Figure 14 :

 314 Figure14: RCS for the cylinder 10log 10 (RCS(ω, 0, 0, θr)) with wave number ω = 3 π for the left and 10π for the right. From top to bottom: exact solution, approximation, relative error.

3 Figure 15 :

 315 Figure 15: Left: comparison of the bistatic RCS (θ W = 18.5 • ), 10log 10 (RCS(10π, 18.5, 0, θr)), with application of curvilinear PML and with a first order absorbing boundary condition. Right: relative error of the computation with curvilinear PML.

Figure 16 :

 16 Figure 16: Comparison of three different mappings

Figure 17 :

 17 Figure17: The 29 θ W 's the greedy algorithm of the RBM picks to build the RB space. Top: the higher the vertical line, the earlier that point was picked. Bottom: the points scattered on the monostatic scattering curve 10log 10 (RCS(10π, θ W , 0, 0)), the larger the marker, the earlier it is selected.

Figure 18 :

 18 Figure 18: The worst case convergence history and the corresponding error estimate of the RBM for 120 randomly selected parameter values.

Figure 19 :

 19 Figure 19: Using only 11 bases, we can obtain a very accurate plot of the monostatic scattering (θ i = θr = 0) respect to the wedge angle θ W .

Figure 20 :

 20 Figure 20: Real part (top), imaginary part (middle), and module (bottom) of the electric field with θ W = 14.3 • (left), 18.5 • (middle) and 21.5 • (right).

Figure 21 :

 21 Figure 21: Error bounds of the RB output, RCS(10π, θ W , 0, 0), for three different number of bases.

Figure 22 :

 22 Figure 22: RBM with EIM for ω ∈ [3π, 5π] and θ W ∈ [8.5 • , 28.5 • ]: (a). The parameter values selected by RBM to build the RB basis. The larger the marker, the earlier it has been selected. (b). Convergence of the error with respect to the dimension of the RB space for different number of magic points. (c). Convergence of the error estimate with respect to the dimension of the RB space for different number of magic points.

Figure 23 :

 23 Figure 23: RCS(ω, θ W , 0, 0) in dB computed by the RBM with EIM for ω ∈ [π, 5π] and θ W ∈ [8.5 • , 28.5 • ]. Top: 3D-plot; middle: top view; bottom: point-wise relative error of RCS(ω, θ W , 0, 0).

Figure 25 :

 25 Figure 25: Comparison of the accuracy of RCS(10π, θ W , 0, 0) in dB with truth approximations computed by a second and a fourth order DG method for θ W ∈ [8.5 • , 28.5 • ].

Figure 26 :

 26 Figure 26: RCS(ω, θ W , 0, 0) in dB computed by the RBM (with 22 bases) with EIM for ω ∈ [9.8π, 10π] and θ W ∈ [8.5 • , 28.5 • ].

Acknowledgements

The second (JSH) and the last (XZ) authors acknowledge partial support by OSD/AFOSR FA9550-09-1-0613.

the monostatic RCS, that is RCS(ω, θ W , 0, 0). Plotted in Figure 23 ]. We plot, in the middle, the top view of RCS. The two regions in this parameter domain for the critical configurations that produces very small RCS is clearly identified. At the bottom is the point-wise relative error between this RB approximation and the truth approximation. Notice that the minimum contour level is set to be 0.01 to show that the approximation is very accurate except for small neighborhoods of parametric configurations rendering the Pacman invisible. In these small neighborhoods, the errors are substantial. We thus perform a new set of computations using 40 RB instead of 30. The result is in Figure 24 showing that the error is indeed getting smaller and closely clustering around a few points.

The results in Figures 22 and 23 are obtained by using second order accurate discontinuous Galerkin method for the efficiency of our calculation. This suffices for frequency in [π, 5π] with the resulting truth approximations being accurate enough. However, when the frequency is increased, we need high-order DG scheme to produce more accurate truth approximation. This is clearly shown by Figure 25 where the RB result is converging perfectly to (an inaccurate) DG result if second order DG is used for ω = 10π.