Numerical Reduced Order Modeling for Wave Equations in Heterogeneous - PowerPoint PPT Presentation

Numerical Reduced Order Modeling for Wave Equations in Heterogeneous Media Tom Hagstrom Southern Methodist University Support: NSF and DOE .

Maxwell’s equations for a dispersive medium: ∂ 1 ∂t ( E + K e ∗ E ) = ∇ × H ǫ 0 ∂ − 1 ∂t ( H + K m ∗ H ) = ∇ × E µ 0 Typically describe the convolution kernels via their Laplace transform - for isotropic problems they will be scalar but in anisotropic cases they will be matrices. In all the examples we consider the convolutions are “lower order” terms in that: ˆ lim K e,m = 0 . s →∞ We believe that this condition combined with the condition that ˆ K is analytic in ℜ s > 0 implies the well-posedness of the Cauchy problem and finite propagation speed. (Proof? Hille-Yosida + Paley-Wiener?) These are multiscale models in that the convolutional terms model the interaction of the complex medium with the waves.

Dispersion relation s = − iω , isotropic case: ω 2 (1 + ˆ K e ( − iω ))(1 + ˆ K m ( − iω )) = | k | 2 Group velocity - convenient formula from, “On the analysis of perfectly matched layers for a class of dispersive media and applications to negative index metamaterials”, B´ ecache, Joly and Vinoles, Math. Comp. 2018. � (1 + ˆ K e )(1 + ˆ 2 K ) �� k V g = c 0 | k | . � � � � (1 + ˆ ω (1 + ˆ (1 + ˆ ω (1 + ˆ K e ) d K m ) d K m ) K e ) dω dω One of many interesting phenomena in metamaterials is the possibility of reverse modes: V g · V p < 0, where V p is the phase velocity - here is an illustration with a simple Drude model and TE modes: K e = ω 2 K m = ω 2 ˆ ˆ e m s 2 , s 2 , We force it at 3 different frequencies - one with phase and group velocity aligned, one in the band gap, and one with reverse modes.

As in the Drude model, there are many instances where the convolution can be approximated to the desired accuracy by a sum of localizable convolutions - that is we can take ˆ K e,m to be rational functions of s . More directly this will be true if ˆ K e,m can be approximated uniformly by rational functions on an appropriate inversion contour: n e α e,j s + β e j ˆ � K e ≈ , ( s + γ e,j ) 2 + η 2 e,j j =1 n m α m,j s + β m j ˆ � K m ≈ . ( s + γ m,j ) 2 + η 2 m,j j =1 In the time domain this corresponds to approximating K e,m by sums of exponential/trigonometric functions. Such approximations were used, for example, to evaluate exact radiation boundary conditions by Alpert, Greengard and H, (SINUM 2000). A general algorithm for their computation is given in Xu and Jiang, (J. Sci. Comput. 2013). The idea is to use a linearized rational least squares algorithm and slowly add poles. More generally one can construct such approximations after removing a small nonsmooth part of the convolution near t = 0 - corresponds to nonsmooth behavior of ˆ K for large s - we will mention this later on when considering dispersion relations involving fractional derivatives. Another very general approach is based on piecewise (in t ) exponential/trigonometric approximations using quadrature on appropriate contours; see Lubich and Sch¨ adle (SISC 2002); see also Trefethen, Weideman and Schmeizer (BIT 2006). Interpolatory methods (data based - now called machine learning!) - e.g. the AAA algorithm (Nakatsukasa, S´ ete, Trefethen 2018), rational Krylov (G¨ uttel - rktoolbox), AAK variant ( Damle, Beylkin, Haut and Monzon 2018)

Note that if α = γ = 0 one obtains the generalized Lorentz models considered by B´ ecache et al. In any case we can introduce auxiliary variables so that the model becomes a symmetrizable hyperbolic system with a zero order perturbation: n e ∂E 1 � = ∇ × H − α e,j E − γ e,j P j + Q j ∂t ǫ 0 j =1 n m ∂H − 1 � = ∇ × E − α m,j H − γ m,j R j + V j ∂t µ 0 j =1 ∂P j = α e,j E − γ e,j P j + Q j ∂t ∂Q j η e,j E − β 2 = e,j P j − γ e,j Q j ∂t ∂R j = α m,j H − γ m,j R j + V j ∂t ∂P j η m,j H − β 2 = m,j R j − γ m,j V j ∂t

Succinctly with scaling the system can be put in the form: ∂W ∂W ∂W ∂W ∂t = A 1 ∂x + A 2 ∂y + A 3 ∂z + SW, j and in some cases with scaling S + S T ≤ 0. (For the Lorentz case S + S T = 0 - see B´ where A j = A T ecache et al). (Stability analysis for S with nonpositive eigenvalues but not symmetrizable is more involved - symmetrizers in frequency space.) We conclude that the system: • Leads to a well-posed Cauchy problem. • Has the same domain of dependence as the nondispersive Maxwell system. • Is well posed under the same boundary conditions as the nondispersive Maxwell system.

Given these facts stability proofs for many numerical methods are directly generalizable to these models. For example Yang, Huang and Li (J. Sci. Comput. 2016) develop methods based on edge elements. For example we can use standard central or upwind DG schemes, but also other energy stable methods such as dissipative or conservative Hermite methods. For a DG scheme we use the fact that the fluxes are unchanged, so we have: Φ T ∂W ∂W ∂W ∂W � � ∂t − Φ T A 1 ∂x − Φ T A 2 ∂y − Φ T A 3 ∂z − Φ T SW − Φ T ( n · A )( W − W ∗ ) = 0 . Ω j ∂ Ω j Choosing the fluxes exactly as for the nondispersive Maxwell system, as they are the only terms that appear in n · A , and choosing the test function Φ = W we have d � � | W | 2 ≤ � � W T ( S + S T ) W, dt Ω j Ω j j j immediately yielding stability. If we further assume S + S T ≤ 0 we conclude that the L 2 norm is uniformly bounded in time.

There are dispersive models which can’t be accurately represented by rational functions due to irregularity of the convolution kernel for short times (or its transform for large s ). An example is given by Havrilak-Negami dielectric model which involves fractional derivatives: for 0 < α ≤ 1, 0 < β ≤ 1 η ˆ K e = (1 + ( sτ ) α ) β . Incorporation of this model into the Yee scheme is achieved in Causley, Petropoulos and Jiang (JCP 2011). For small t they directly treat the convolution over the current time step using an expansion of K e : � t � α ( k + β ) − 1 K e ( t ) ≈ η � c k ( α, β ) . τ τ k They then build a modification of the Yee scheme to account for the local behavior. For the history, that is the convolution past ∆ t , they use a rational approximation as above computed via quadrature applied to integral representations of the kernel. Recently Baffet and Hesthaven have used the arguments in Alpert, Greengard, H. to directly construct approximations to time fractional derivative operators ∆ t away from the current time and manage to construct a scheme where the pole locations are independent of ∆ t to allow high order adaptive time stepping based on Adams methods. (Math. Comp., J. Sci. Comput. 2017).

All of this assumes that K e,m is known. In many cases this may not be true - as shown, e.g., in work by Bouchitt´ e a first approximation can be obtained using homogenization theory, but formally at least these formulas require a very large scale separation. This leads to the question of whether a numerical approach could be useful - numerical multiscale methods. (Review by Abdulle and Henning 2018) Much work leverages ideas form elliptic problems - doesn’t really take advantage of the fundamental feature of hyperbolic pdes which is local domain-of-dependence. The most famous method is the so-called Heterogeneous Multiscale Method (HMM). (e.g. E and Engquist, Multiscale Methods in Science and Engineering, 2009). The basic idea of the HMM framework is to: • Assume a macroscopic model exists (in our case the dispersive Maxwell system) • Introduce a macro numerical method appropriate to the model • Solve the full microscopic model in microscopic neighborhoods of quadrature/flux nodes of the macroscopic dis- cretization and average over space-time to provide the missing quadrature/flux data. An issue in applying HMM to wave problems is that the local solutions oscillate in time rather than relax. There are some recent papers on the scalar wave equation.

Arjmond and Runborg, Mult. Mod. Sim. 2014 Here they study the ill-posed macro model: ∂ 2 u ∂t 2 = α∂ 2 u ∂x 2 + ǫ 2 β ∂ 4 u ∂x 4 They use the wave equation as the microscopic model with specially prepared initial data and averaging operators and prove that the flux is in their discretization is correct with a small error. Abdulle, Grote and Stohrer, Comptes Rendu Math. 2013 Here they study the well-posed macro model: ∂ 2 u ∂t 2 = α∂ 2 u ∂ 4 u ∂x 2 + ǫ 2 β α ∂x 2 ∂t 2 They solve the Neumann problem in the microscopic level and use the solution both in the quadrature for the mass matrix and the stiffness matrix: ∂ 2 ∂t 2 (( u H , v H ) + ( u H , v H ) M ) + B H ( u H , v H ) = ( f, v H ) w j � � � a ǫ ( x ) ∇ v h · ∇ w h = B H ( v H , w H ) | K δ | K δ T H j w j � � � ( v h − v H,lin ) · ( w h − w H,lin ) = ( v H , w H ) M | K δ | K δ T H j � a ǫ ( x ) ∇ ( v h − v H,lin ) · ∇ z h = 0 . K δ Question : How can we apply this method to Maxwell’s equations where the macroscopic model is the dispersive system we’ve been talking about?