Title: Computational Electromagnetics
1Computational Electromagnetics
- Daniel S. Katz
- Director for Cyberinfrastructure Development, CCT
- Adjunct Associate Professor,
- Electrical and Computer Engineering
- Louisiana State University
2Electromagnetics
- Maxwells Equations
- Lots of versions, pick the right set for your
problem and methodology - Wavelength and frequency are inversely related
- An object of size 1 m is one wavelength long at
300 MHz or 2 wavelengths long at 150 MHz - Either frequency or size can be scaled as needed
- Radar Cross Section (RCS) as an example problem
- A plane wave at some incident angle and some
frequency illuminates a target. - Monostatic or Backscatter RCS what energy comes
back? - Bistatic RCS what energy goes off in another
direction? - Toy problem sphere
- Radius r volume 4/3 P r3 surface area 4 P
r2
3CEM methods
- Maxwells Equations-based methods
- Optics-based methods
4 5Outline
- Method of Moments (MoM)
- Finite Difference Time Domain (FDTD)
6Corner Reflector
- Surface currents displayed
- Tone burst has hit reflector
- Burst incident at 45 to left, right sides 10
to bottom - Frequency 10 GHz
- Dimension 5 cm on edge
Photograph of computer monitor, 1989
7Dielectric Lens
- Dielectric lenses can be made in different
materials with different properties - Above quantum well infrared photodetector (QWIP),
can increase QWIP sensitivity by 14x - 20 THz plane wave incident downward
Scanning Electron Microscope (SEM) Image of a
portion of the dielectric lens (credit Dan
Wilson, JPL Microdevices Laboratory)
8Method of Moments
- Maxwells Equations in Integral Form
(Stratton-Chu Equations) - Total field (to which boundary conditions are
applied) is summation of incident field (defined
by incident wave only) and scattered field (the
result of the incident wave interacting with the
scattering object)
9Method of Moments (2)
- Electric Field Integral Equation (EFIE)
- Continuity relationship for current and charge
density
10Method of Moments (3)
- For a perfect electric conductor, tangential
component of E on the surface is zero boundary
condition - And M is zero
- Can rewrite EFIE as operator on
11Method of Moments (4)
- Represent surface current as summation of element
currents and basis functions - Enforce the boundary condition on the surface (at
N unique locations match points) - At each (N) match point, the total field is the
sum of - the incident field, and
- the field caused by currents located on the other
points - This leads to a coupled system of equations by
which every part of the body affects every other
part (Maxwells equations are elliptic) - Sample the region around the match points and
average by using a weight function
12Method of Moments (5)
- Define inner product integral as
- (Integration is over the pth surface area which
surrounds the pth BC match point)
- jn current at location n
- vn forcing function at location n
- znm interaction between locations n and m
13Method of Moments (6)
- Assume incident place wave
- Impedance matrix terms
- Scattered field
14Sequential MoM
- jn current at location n
- vn forcing function at location n
- znm interaction between locations n and m
- Based on geometry, calculate
- LU decompose
- Based on input plane wave, calculate
- Backsolve for
- Use to find far fields, or RCS, using
equivalence principle
(image from http//www.cgal.org/)
15Sequential MoM
- How big is the matrix? (What is N?)
- Depends on the geometry
- Need to have enough points to define the object
- Also want at least 5 points per linear wavelength
- For a sphere of radius 1 wavelength, N 200
- For a sphere of radius 2 wavelengths, N 8000
- N is O(x2), x is linear size of object
- Matrix is size O(N2) or O(x4)
- Filling matrix is O(N2) or O(x4)
- Solving the matrix is O(N3) or O(x6)
- Backsolve for one RHS is O(N2) or O(x4)
- Calculating one RCS is O(N) or O(x2)
16Sequential MoM vs. problem
- One LU decomposition suffices for an object at a
frequency - Multiple incident angles means multiple RHSs
- Multiple scattering angles means using jns
multiple times - Change frequency, redo LU decomposition
- Time saving tricks
- Use symmetry 8 way symmetry (sphere) reduces
matrix size to 1/8th, and solution time to
1/512th - Easy to calculate znm and zmn at the same time,
reduce matrix fill time by almost ½
17Parallel MoM
- Store matrix in block cyclic format for good
performance in solve - Good point
- Can use standard parallel linear algebra
libraries - ScaLAPACK, PLAPACK, etc.
- Hard to fill matrix and RHS
- Easiest to have all the geometry on all
processors - Divide patches by procs, fill rows of matrix that
go with each procs patches, same for RHS - Problems
- Still hard to fill the block-cyclic matrix
- Either lots of communication or lots of
duplicated computation - Probably have to compute znm and zmn
independently
18Direct vs. Iterative Matrix Solution
- Given BAX, where A is a matrix and B and X are
vectors, a direct solution for X is basically
finding A-1B - This is the LU decomposition approach, more or
less - O(N3)
- Iterative schemes are also possible, where the
core of the iteration loop is multiplying A times
a vector y - Each iteration (matrix vector multiplication)
is O(N2) - May need a lot of iterations to converge
- For dense matrices, this usually doesnt make
sense - What if A is sparse?
- Matrix vector multiplication with a sparse matrix
may be low cost, particularly if the number of
non-zero elements is small and only these
elements of the matrix are stored - O(n_non_zeros)
19Fast Multipole Method
- Introduced by Greengard and Rokhlin in mid 1980s
- Group patches that are close to each other
- Leads to errors, but the scale of the errors can
be predefined - Similar to some n-body solution methods
- Use an iterative solver for the matrix
- Dont have to store all elements, can build some
on the fly - Use FMM hierarchically to get solver thats O(n
log n) - Classical MoM is basically dead for large
problems
20Outline
- Method of Moments (MoM)
- Finite Difference Time Domain (FDTD)
21Maxwells Equations in Curl Form
- Maxwells (curl) Equations
- Electric field vector
- Magnetic flux density vector
- Electric flux density vector
- Magnetic field vector
- In linear, isotropic, non-dispersive media
22Maxwells Equations in Curl Form
- Writing out the vector components
(image from wikipedia)
231-D Maxwells Equations
- Assume E only has a z component, and that
everything is constant in y
24Central Differencing
- Consider Taylor series expansion of f(x) expanded
around x0 with an offset of d/2 - Subtract the second from the first
- Notes
- This is second-order accurate
- No sample at x0
251-D FDTD
261-D FDTD details
- Non-rigorously
- Energy should not propagate more than one spatial
step in each temporal step -
- Computer implementation
271-D FDTD Code
- Define media (ca, cb)
- Initialize fields to zero
- Loop over time (n 1 to nmax)
- Loop over space for ez (i0 to imax)
- ezi cai(hyi-hyi-1)
- Loop over space for hy (i1 to imax-1)
- hyi cbi(ezi1-ezi)
281-D FDTD Code - BC
- What about Ez0 and Ezimax?
- We need boundary conditions to ensure that waves
propagate past these points without reflecting - Simple choice, if dt/dxc
- Ezn0 Ezn-1 1
- Mathematic/geometric option in 2d and 3d
- Mur RBC (1981) Mur RBC
- Model absorbing material (virtual range)
- Berenger (1994) Berenger PML
291-D FDTD Code - Inputs
- How to input energy into the system?
- Use a hard source
- ez10 cssin(omegadttimestep)
- Simple, but leads to reflections
- Use a soft source
- Amperes Law
- Apply finite differences
- Separate into normal update and additive source
- ezi cai(hyi-hyi-1)
- ez10 cssin(omegadttimestep)
301-D FDTD Code - Scatterers
- How to find scattered field?
- Use a total field / scattered field formulation
for the main grid - Compute two 1-D grids, one for the incident field
and one for the total/scattered field - Incident grid is homogeneous TF/SF grid has
scatterer geometry - Add/subtract incident field on total
field/scattered field boundaries
311-D FDTD Code Scatterers (2)
- eztotal50 ca50(hytotal50-hytotal49)
- Correct update from difference equation, but
doesnt match grid - hytotal49 hyinc49 hyscat49
- eztotal50 ca50(hytotal50-hyscat49)
(normal update) - eztotal50 - ca50hyinc49 (special update
for TF/SF interface) - Similar changes needed for hy49 update, and ez
and hy at TF/SF interface on right side of grid
32Ghost Cells (2D)
- Parallel Implementation
- Need to update these cells on a given processor,
using second order central differences (one cell
on each side) - In order to update outer cells, need cells one
step further away - These have to be communicated from neighboring
processors
j
i
33Ghost Cells (2D) saving communication?
- What if we communicate more cells every other
time step? - Seems like less communication calls, same amount
of communication total - Better for small messages, where latency matters
more than bandwidth - But there are problems with this idea
j
i
34Load Balancing
- How to divide this domain for 4 procs?
- OpenMP worry about
- Work
- MPI worry about
- Memory
- Work
- Communication
j
i
353-D Decomposition
Standard Domain Decomposition
Required Ghost Cells
One plane of ghost cells must be communicated to
each neighboring processor each time step.
363-D Decomposition (2)
Standard Domain Decomposition
Boundary Decomposition
Data at four faces must be redistributed twice
each time step!!
37Parallel FDTD Modeling Example Periodic
Plasmonic System
wraparound boundary conditions for side domain
walls. PML for top and bottom boundary
3D FDTD domain of unit cell and domain
decomposition
Credit Tae-Woo Lee
38Frequency Domain Results
- FDTD is time domain how to get frequency domain
results for scattering calculations? - Keep track of fields at some locus in the
scattered field domain - Use FFT or similar to transform from time domain
to frequency domain (see first part of notes) - Use equivalence principal to transform to
far-field (RCS) - Can do this for multiple bistatic angles after
one FDTD run - Need to run long enough to reach steady state
- Approx. enough time for fields to move from one
corner of the domain to the opposite corner 4
times - What about multiple frequencies?
- Can use incident pulse with good frequency
content, FFT-1
39FDTD Staircasing
- Want Get
- Can modify update equations near surface for
smoothing
40Engine Inlet
- Aluminum inlet in box coated with radar absorbing
material - 2-D slice of 3-D problem
- Inlet dimensions 5 x 5 x 30 cm
- Scattered electric field magnitude displayed
- 10 GHz plane wave incident from right side
- Slice after 120 cycles of incident wave
Photographs of computer monitor, 1989
41Credits
- MoM material
- Shaeffer, MOM3D Method of Moments Code Manual,
1992 - Wikipedia
- Umashankar, UIC
- FMM material
- Wikipedia
- FDTD material
- Allen Taflove, Northwestern, John Schneider, WSU,
Tae-Woo Lee, LSU