Time-Domain Electromagnetic Analysis : From FDTD to Wavelets and Beyond

1. C o s t a s D. S a r r i s The Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan, Ann Arbor Advisor : Prof. Linda P. B. Katehi The University of Toronto, June 28, 2002 Time-Domain ElectromagneticAnalysis : From FDTD to Wavelets and Beyond

2. Outline • Introduction • Time Domain numerical schemes for Maxwell’s equations : Research motivation • and state-of-the art. • Wavelet based numerical schemes : How they started, where they are. • The Multi-Resolution Time-Domain technique. • Recent developments in wavelet-based schemes • Formulation and dispersion analysis. • Numerical interface between FDTD and MRTD : Efficient implementation of boundary and • Perfectly Matched Layer conditions. • Implementation of dynamic mesh adaptivity : The case study of a nonlinear optical pulse • propagation. • Advanced Application : Wireless channel modeling • Cosite Interference in VHF transceiver networks. • Mixed electromagnetic-circuit simulations. • Modeling of cosite interference in a forest environment. • Trends and Conclusions

3. Marching in time scheme Finite Difference – Time Domain (FDTD) Method • In FDTD (K.S. Yee, 1966), the computational domain is divided in Yee’s cells and Maxwell’s equations are solved by marching in time. Example : FDTD discretization of • Simple, versatile and robust algorithm. • Second order accurate in space and time. • Inherently parallelizable. • Ten-twenty points per wavelength necessary. • Small time step necessary. • Multi-wavelength domains result to extremely large scale problems.

4. = + Coarse approximation Successive wavelet refinement Multiresolution Decomposition Example: Mesh Refinement in Haar wavelet MRTD Use of N wavelet levels in direction x, brings about a mesh refinement by a factor : Multi-Resolution - Time Domain (MRTD) Method • In MRTD [Krumpholz and Katehi, MTT-T 1996], field components are expanded in a wavelet basis : • Wavelet coefficients are significant near field variations/discontinuities. • MRTD offers a natural framework for the implementation of an adaptive, moving mesh.

5. High Order Methods FDTD + subgridding Wavelet Based Methods Subgrid at details + Interpolations/ extrapolations Higher accuracy / Larger stencil Higher order Wavelet basis Thresholding Modeling of geometric details Dynamic adaptivity Coarse discretization Computational Efficiency Comparison of Novel Time Domain Schemes

6. Haar scaling function Haar mother wavelet Example : Haar Basis • Finite Domain functions – only nearest neighbor coupling of basis functions ( no stencil effects ). • Extension of the scheme to higher orders is relatively simple. • Haar MRTD with scaling functions only is the FDTD scheme. • Low approximation order scheme ( FDTD and Haar MRTD need same number of degrees of freedom).

7. K-th cell p=0 p=1 p=2 p=3 r=2 p=0 p=1 r=1 p=0 r=0 Example : Haar Basis (cont-d) • Haar scaling basis : • Higher order Haar wavelet functions are produced with translations and dilations of the mother wavelet. • Example : Order r wavelet in k-th cell : Haar wavelets of orders 0, 1, 2

8. Example : Battle-Lemarie basis Battle-Lemarie mother wavelet Battle-Lemarie scaling function

9. Use of Wavelets for Mesh Refinement : Promise • Introduction of one wavelet level refines the mesh by a factor of two with respect to a scaling function based numerical scheme. • Example : Pulses + zero order Haar wavelets = FDTD scheme of half its original cell size per direction. • If D is the scaling cell size, the effective cell size for an R-order MRTD is D / r, r = 2 R+1.

10. Use of Wavelets for Mesh Refinement : Results • Battle-Lemarie wavelet based W-MRTD [Krumpholz and Katehi, MTT-T, April 1996] : • Spurious modes observed, attributed to wavelets. • Use of wavelets caused only incremental improvement in accuracy. • Haar wavelet based MRTD [Goverdhanam et al., Fujii and Hoefer] : • For zero order MRTD, scaling and wavelet terms are given by uncoupled equations. • Incremental accuracy improvement due to wavelets : • “..addition of one resolution of wavelets within a Haar based MRTD framework does not improve significantly the numerical accuracy of the underlying coarse grid scheme.” • Grivet-Talocia, IEEE MGWL, October 2000. These observations are in stark contradiction with Multiresolution Analysis principles, as established in approximation theory and signal processing.

11. Dispersion Analysis • Dispersion analysis consists in the determination of the numerical wavenumber as a function of the numerical frequency : • In practical terms, it answers the question of how coarse a grid can become for an allowed level of phase errors (important for CAD oriented algorithms). • FDTD dispersion analysis [Taflove, 1994] : • Substitute all terms in finite difference equations with plane wave type expressions, • Formulate a linear homogeneous system with respect to the amplitudes Lx, impose condition that the system have a non-trivial solution. p, q , m : space cell indices n : time step index x : x, y, z • Not directly applicable to MRTD, due to the multi-level character of its expansion basis.

12. A Modified Fourier Dispersion Analysis Method for MRTD1 • Dispersion Analysis has been re-formulated from first principles of Fourier analysis. • In the MRTD finite difference equations, replace all scaling coefficients with FDTD type plane wave expressions (one dimension) : • Observe by simple Fourier calculus arguments that : • Formulate a linear system with respect to scaling spectral amplitudes ONLY. Note : The number of unknowns of the system is independent of MRTD order. Three step approach x : y, z [1] : C. D. Sarris, L. P. B. Katehi, “Some Aspects of Dispersion Analysis of MRTD Schemes”, 2001 ACES Conference.

13. Dispersion Analysis : Gridding Effects Haar MRTD Battle-Lemarie MRTD : Electric field scaling node : Electric field scaling function : Magnetic field scaling node : Magnetic field scaling function • Most MRTD studies in the literature so far have been using a standard half scaling cell offset between electric and magnetic field nodes ( referred to from now on as Formulation I ).

14. Dispersion Analysis : Effect of Gridding on Haar MRTD • Dispersion analyses of arbitrary • order MRTD schemes were • performed, based on the modified • Fourier method for MRTD. • One Dimensional Propagation • results are shown. • Notation : Normalized Wavenumber Normalized Frequency • Zero order MRTD dispersion performance coincides with FDTD of cell size equal to the scaling cell size. • In general, this formulation of Haar MRTD schemes, leads to a decrease in resolution by a factor of two in ALL orders.

15. Gridding Effects (cont-d) • A better approach : Vary the offset between electric and magnetic field nodes so that it remain half • the MRTD equivalent cell, taking into account the mesh refinement brought about by the wavelets : Formulation II (proposed in this work) • Equivalent grid points in 0- • Haar MRTD are determined • under the new convention : • s=0.25 (rmax = 0). • The equivalent grid • points are properly offset. Reference : Sarris and Katehi, IEEE MTT-T, Dec. 2001

16. Dispersion Analysis : Effect of Gridding on Haar MRTD (ii) • Dispersion analysis of arbitrary order MRTD schemes under Formulation II has been performed. • Results for one dimensional wave propagation are shown. • Each wavelet level increases the Nyquist limit of the scheme (turning point in the above curves) by a factor of 2 (FDTD Nyquist point is X / p = 1 ) => Consistence with MRA has been enforced.

17. Applying the same gridding principles, the W-MRTD scheme [Krumpholz and Katehi, 1996] is reformulated. Formulation 1 : Half cell E/H offset Dispersion Analysis : Effect of Gridding on W-MRTD Formulation 2 : Quarter cell E/H offset Nyquist limit of the reformulated W-MRTD scheme is now at 2 p.

18. a (x) (z) a=32 cm Results : Resonant Frequencies for a two-dimensional air-filled cavity • TEz modes have been numerically determined for a square air-filled cavity, using FDTD and MRTD according to the two formulations under comparison. • In all MRTD cases 32 by 32 degrees of freedom have been used. • Results from FDTD of 32 by 32 and 16 by 16 cells are used for comparison. Electric field Ey distribution for TE22, TE32 determined by 4 by 4 Haar MRTD (1 by 1 scaling mesh).

19. Resonant Frequencies for a two-dimensional air-filled cavity (cont-d) • Simulations for the schemes under comparison were run at a time step corresponding to 0.9 of their CFL number. • MRTD schemes of various orders up to 4 per direction (scaling cell size = 32 cm) were tested. • Haar MRTD (Form. I) follows a 16x16 FDTD in accuracy. • Haar MRTD (Form. II) consistently follows a 32x32 FDTD.

20. Example : One dimensional cavity structure. W-MRTD : One Dimensional Case Study Domain defined by 3+2 Battle-Lemarie scaling functions. • Theoretical Dispersion Relationship : • A Battle-Lemarie scaling and a wavelet function are used as initial data (injecting normalized wavenumbers from 0 to 2 p), exciting the first 7 modes. • W-MRTD schemes of formulations I, II are employed. • Both schemes are run at half their CFL number.

21. W-MRTD I, II : Field Spectrum • The spectrum of W-MRTD (form. I) contains a mode that appears to be spurious (but..it is not [Sarris and Katehi, IEEE MGWL Feb. 2001] ! ) . • W-MRTD (form. II) appears to resolve all seven modes.

22. Hybridization of MRTD : Motivation • In typical microwave circuit geometries.. Desired (coarse) MRTD mesh • In the modeling of the PML absorber ( hence in all kinds of open problems ). • An efficient solution to such problems is essential for the application of MRTD to practical microwave structures. • Approach of this work : Combine the versatility of FDTD with the adaptivity of MRTD, by developing a hybrid FDTD / MRTD scheme based on a numerical interface . Reference : Sarris and Katehi, Proc. 2001 IEEE IMS, paper submitted to MTT-T, April 2002.

23. Matching of MRTD / FDTD dispersion properties at the interface Haar MRTD FDTD • Evidently, letting : • or : • the dispersion properties of the two schemes match. FDTD cell size = MRTD effective cell size • A reflectionless propagation through an FDTD / Haar MRTD interface can be simulated without interpolations / extrapolations.

24. Two dimensional FDTD / MRTD Interface • Overlapping of the two regions and use of field equivalence principles to de-couple the two • problems (FDTD and MRTD). • Vehicle for data transfer between the two regions is the Fast Wavelet Transform, at the • optimal cost O(N).

25. FDTD / MRTD interface : Validation • 3 by 4 MRTD interfaced with FDTD (MRTD mesh = 1x1). • Effective cell size : 1cm x 1cm. • The first six cavity modes are extracted via FDTD and the interface.

26. FDTD / MRTD interface : Mode Patterns MRTD mesh 1x1 (order 4 by 4) MRTD mesh 4x2 (order 2 by 3) • TE21 and TE22 are extracted via the FDTD / MRTD interface. • Absolutely smooth mode patterns are observed.

27. FDTD / MRTD interface : Metal fin loaded cavity Geometry and Mesh for the FDTD / 2 x 2 MRTD interface • The pure MRTD modeling of this structure is restricted by the size of the scaling function. • When this size exceeds the width/height of the fin, the fin scaling cell unphysically couples the regions above and below the fin. Domain split is then necessary.

28. Metal fin loaded cavity : Results • A Gaussian excitation with its 3 – dB bandwidth equal to half the cut-off frequency of the TE11 mode of the cavity is used. • Absolutely stable performance of the code is observed (electric field sampled in the cavity for time steps 18,000-20,000 is shown).

29. : 2nd order MRTD equivalent grid points : FDTD grid points MRTD / FDTD PML interface • A common problem of wavelet based numerical methods is that they need many grid points to model boundary conditions. • For MRTD, it may seem necessary that at least the degrees of freedom of one MRTD cell are used to model the PML region. • However, since PML itself is terminated into a PEC, the following mesh configuration can be used to terminate the MRTD mesh in a PML region that extends over a fraction of the grid points of a single MRTD cell : • The termination of a second order MRTD mesh in a four grid point PML is assumed. • One MRTD cell contains 8 equivalent grid points. • All grid points that are beyond the PML region are zeroed out and used in the fast wavelet transform of the interface as such.

30. MRTD / FDTD PML interface : Dipole Radiation • An MRTD mesh is terminated in an FDTD based Uniaxial PML absorber ( 8 cells, theoretical reflection coefficient R = exp(-16), order 4 polynomial variation of conductivity) . • The problem of an infinitesimal y-axis directed current element radiation is simulated via the FDTD / MRTD interface. • Results are compared to FDTD. Mesh in FDTD cells

31. Dipole Radiation : Results 0 by 0 order MRTD (32x32 cells) 2 by 2 order MRTD (8x8 cells) 4 by 4 order MRTD (2x2 cells) • MRTD results are compared to FDTD for three different schemes. • Field sampling points are at A, B, C respectively. • UPML extends over 4, 1, .25 MRTD scaling cells respectively. • Excellent agreement between FDTD and interface results is observed.

32. Dipole Radiation : Results (cont-d) Time Steps = 20, 100, 200, 600 are shown. 3 by 3 order MRTD mesh (2x2 scaling cells).

33. Characterization of a dielectric slab waveguide (4th order MRTD) • Scaling cell size = 8 mm (=1.28 lmin). • 5 wavelet levels to bring the effective cell size down to lmin /25 (range : 0-30 GHz). • The whole slab is included in one cell. 6 grid point PML (0.1875 of a cell )

34. Wavelet Based Adaptivity : Nonlinear Wave Propagation • As an example, the discretization of a system of equations and its boundary conditions describing a pulse compression in an optical fiber filter is presented. • Fiber index of refraction is field intensity – dependent, assuming the form : • Co - sinusoidal term is due to a grating written within the core of the fiber. • Pulse compression comes as a result of negative dispersion that causes the rear of the pulse to travel faster than the front of it. • Numerical modeling via an S-MRTD scheme (Battle-Lemarie scaling functions) was pursued in [Krumpholz and Katehi, MTT-T, 1997]. Larger – than – FDTD execution time was reported (factor of 1.5), despite the application of coarser grids. • References : H. Winful, Appl. Phys. Lett., Mar. 1986 • Sarris and Katehi, 2002 Proc. IEEE IMS

35. Nonlinear Wave Propagation : Results • Results from Haar MRTD (with no adaptivity) are compared to FDTD and excellent agreement is observed. • In terms of execution time, MRTD appears to be slower than FDTD by ~15%. • This slowdown motivates our next step towards formulating adaptive MRTD schemes.

36. Application of Thresholding in MRTD Two Approaches Absolute Thresholding Condition (“hard” thresholding) Relative Thresholding Condition (“soft” thresholding) • times the scaling value of the cell is the (dynamic) threshold under which a wavelet term contribution to the total field is insignificant. • is the threshold under which a wavelet term contribution to the total field is insignificant. • Both conditions are implemented as if – tests in the code. • Hard thresholding is applied in this work (simpler to implement, works well for • wave propagation problems).

37. Adaptive Meshing in MRTD : Front Tracking • Conventional thresholding approaches are based on testing of thresholding conditions throughout the mesh. • This results insignificant operationoverhead. • Instead, one cantrack the wavefrontspropagating in a certain geometry and apply the tests only within those. • Example : 0 – 10 GHz pulse incidence on a dielectric slab of permittivity 2.2 (dots indicate boundaries of non-thresholded coefficients).

38. Initialize all coefficients as “active” Apply update equations for scaling terms Use only active wavelet coefficients Apply update equations for “active” wavelet terms Apply thresholding every N time steps Determine which active wavelet coefficients remain above the hard threshold. Designate as active all nearest neighbors of active wavelet terms (pivot elements). Front tracking step Operation Savings in Adaptive MRTD : Algorithm : active wavelet coefficients : pivot elements Reference : Sarris and Katehi, Proc. 2002 IEEE IMS, long paper to be submitted to MTT-T, Aug. 2002.

39. Execution Time Savings : Linear Wave Propagation • Timing measurements for the 0-10 GHz • pulse propagation were performed for • FDTD and adaptive MRTD. • For thresholding, the previous algorithm is applied at each time step. • As the domain grows large, MRTD • becomes twice as fast as FDTD.

40. Front Tracking : Nonlinear Wave Propagation • Front evolution along the characteristic line z = ct is observed. • Wavelets track the nonlinear evolution of the wavefront, assuming higher values as the pulse gets compressed (implying higher spatial field derivatives around the front).

41. Nonlinear Wave Propagation : Execution Time Savings Comparison to FDTD for the example of optical pulse compression • Thresholding windowstands for the number of time steps between two subsequent applications of thresholding to wavelet coefficients. • For stability reasons, if s is theCFL number, the maximum window should not exceed 1/s. • Different thresholding windows are compared and shown to yieldsimilar CPU time performance. • Assuming stable performance, sparser thresholding checks imply smaller operation economy.

42. Relative Error in Peak Intensity [%] Nonlinear Wave Propagation : Thresholding Induced Error • Using the unthresholded MRTD code as a reference, the relative error in the peak intensity at the end of the fiber filter is determined. • CPU time is reduced compared to FDTD by a factor of 30%, with absolute errors limited to the order of 0.1 %. • This case study involves complex operations in a domain largely occupied by the pulse.

43. 256 256 Adaptive Meshing and Parallelization : Dipole Radiation 1 • A critical issue for the performance of parallel codes is load balancing along multiple processors. • Adaptive repartitioning of a certain domain can be based on using wavelet coefficients of cells, as • load measures (weights). • Domain repartitioning assigns cells with large wavelet coefficients to more processors and distributes • inactive cells to the rest. • Case study : Hertz dipole radiation. 1 Joint work with P. Czarnul, University of Michigan, Prof. K. Tomko, University of Cincinnati.

44. Adaptive Meshing and Parallelization : Dipole Radiation • Zero order Haar MRTD, absolute threshold for wavelet coefficients = 1.0 Normalized Ey Field Plot Unthresholded wavelet coefficients (mesh refinement by 4) • Parallelization on 64 processors. • Repartitioning of the domain every 500 time steps (based on Zoltan software).

45. Courtesy : humvee.net Mutually coupled monopoles Power Amplifier Mixing stages Antenna T/R switch Vehicle #2 Low Noise Amplifier Cosite Interference in Ad-Hoc Networks : Problem statement Mobiletransmit - receive modules on army vehicles in proximity of each other, may corrupt each other’s performance (desensitization, intermodulation, cross-modulation etc.). Vehicle #1 Modeling Approach: Unified treatment of electromagnetic propagation across the physical channel and the operation of front – end electronics (including amplifier nonlinearities ).

46. Receiver Modeling : Lumped Elements in Time-Domain Schemes Total input current Device current + Displacement current • For the modeling of lumped element loading of the antennas, the state equation approach of [1] is implemented. • This technique allows for the inclusion of active/nonlinear loads as parts of the input stage of the transceiver. [1] : B. Houshmand, T. Itoh, M. Picket-May : “High-Speed Electronic Circuits with Active and Non-linear Components”, ch. 8 in Advances in Computational Electrodynamics : The FDTD Method, A. Taflove, ed.

47. Idev Receiver Modeling : Lumped Elements in FDTD (cont-d) The previous form of Ampere’s law at the device port can be interpreted as Kirchhoff’s current law for an equivalent circuit : Equivalent FDTD cell capacitance • Electromagnetic wave / device interaction is described by Norton/Thevenin type equivalent circuits (Norton equivalent is shown here).

48. Equivalent circuit for an analog VHF receiver MESFET-transistor 2 section Chebychev band-pass filter Single balanced mixer Large signal model (amplifier operation in the saturation region is included). Equivalent circuit for Yee-cell Impedance matching Impedance matching 3 section Chebychev low-pass filter

49. Gain without an interfering signal Gain with an interfering signal with a power level of 10 dBm Gain with an interfering signal with a power level of 13 dBm Frequency of the interfering signal Cositing Effect on Receiver Performance • MESFET-amplifier gain for a 50MHz input signal with a power level of –10dbm, • under cositing conditions is simulated in HP ADS. • Loss in sensitivity up to 25 dB is observed due to interference-driven • saturation of the amplifier.

50. Monopole antennas on a PEC • Antenna length: 1.2m • Antenna impedance: 36 Ohm • Distance between antennas: 1.5m Simple Cosite Interference Scenario : Receiver Sensitivity • 50 Ohm source: • available power: -20dBm • frequency: 50MHz IF-frequency: 1MHz • MESFET receiver: • tuned to 50MHz • input impedance: 49 Ohm • IF power under no interference is –30.7 dBm.