1. technology and computing
2. consumer electronics
3. tv and video equipment
4. projectors

5/3/2015
Instructor
Dr. Raymond Rumpf
(915) 747‐6958
[email protected]
EE 4395/5390 – Special Topics
Computational Electromagnetics
Lecture #8
Perfectly Matched Layer
 These notes may contain copyrighted material obtained under fair use rules. Distribution of these materials is strictly prohibited 
Lecture 8
Slide 1
Outline
•
•
•
•
•
•
Background Information
The Uniaxial Perfectly Matched Layer (UPML)
Stretched Coordinate PML (SC‐PML)
Calculating the PML Parameters
Incorporating a UPML into Maxwell’s Equations
Incorporating a SC‐PML into Maxwell’s Equations
Lecture 8
Slide 2
1
5/3/2015
Background Information
Lecture 8
Slide 3
Tensors
Tensors are a generalization of a scaling factor where the direction of a vector can be altered in addition to its magnitude.
Scalar Relation 

V

V

aV

a
V
 
 Tensor Relation
a
  xx
 a V  a yx

 azx
Lecture 8
axy
a yy
azy
axz  Vx 

a yz  Vy 
azz  Vz 
Slide 4
2
5/3/2015
Reflectance from a Surface with Loss
Complex Refractive Index
N  n  j
n  ordinary refractive index
  extinction coefficient
Reflectance from a Lossy Surface
1 n   2

R
2
1  n    2
2
air
** Loss contributes to reflections
N  n  j
Lecture 8
Slide 5
Reflection, Transmission and Refraction at an Interface
Angles
inc   ref  1
n1 sin 1  n2 sin  2 Snell’s Law
TE Polarization
2 cos 1  1 cos  2
 2 cos 1  1 cos  2
2 2 cos 1

 2 cos 1  1 cos  2
rTE 
tTE
TM Polarization
2 cos  2  1 cos 1
1 cos 1  2 cos  2
2 2 cos 1

1 cos 1  2 cos  2
rTM 
tTM
Lecture 8
Slide 6
3
5/3/2015
Maxwell’s Equations in Anisotropic Media
Maxwell’s curl equations in anisotropic media are:


  H  j 0  r  E


  E   j0  r  H
These can also be written in a matrix form that makes the tensor aspect of  and  more obvious.
 0  z y   H x 
 xx  xy  xz   Ex 
 





0  x   H y   j 0  yx  yy  yz   E y 
 z

 
 zx  zy  zz   Ez 
0   H z 


 y x
 0  z y   Ex 
  xx  xy  xz   H x 
 





0  x   E y    j0   yx  yy  yz   H y 
 z

 
  zx  zy  zz   H z 
0   Ez 


 y x
Lecture 8
Slide 7
Types of Anisotropic Media
There are three basic types of anisotropic media
 iso
 0

 0
 o
0

 0
 a
0

 0
Lecture 8
0
 iso
0
0
o
0
0
b
0
0 
0 
 iso 
0
0 
 e 
0
0 
 c 
isotropic
uniaxial
Note: terms only arise in the off‐
diagonal positions when the tensor is rotated relative to the coordinate system.
biaxial
Slide 8
4
5/3/2015
Maxwell’s Equations in Doubly‐Diagonally Anisotropic Media
Maxwell’s equations for diagonally anisotropic media can be written as
 z
 0
 
 z
 
 y
0

x
 Hx 
 xx

   H y   j 0  0
0   H z 
 0

y

x
0
 yy
0
0   Ex 
0   E y 
 zz   Ez 
 0
 
 z
 
 y
 z
0

x
  Ex 
  xx

 x   E y    j0  0
0   Ez 
 0

y
0
 yy
0
0  Hx 
0   H y 
 zz   H z 
We can generalize further by incorporating loss.
 0  z y   H x 
 x   xE j
  Ex 
0
0
 


 

E
 
 y   y j
0  x   H y   j 0 
0
0
 z
  Ey 
E

 


 z   z j   Ez 
0   H z 
0
0

 y x
 0  z y   Ex 
 Hx 
  x   xH j
0
0
 

 


H
 
0  x   E y    j0 
0
0
 y   y j
 z
 H y 
H

 


0   Ez 
0
0
 z   z j   H z 

 y x
Lecture 8
Slide 9
Scattering at a Doubly‐Anisotropic Interface
Refraction into a diagonally anisotropic materials is described by
sin 1  bc sin  2
a 0 0
 r    r    0 b 0
 0 0 c 
Reflection from a diagonally anisotropic material is
rTE 
a cos 1  b cos  2
a cos 1  b cos  2
rTM 
 a cos 1  b cos  2
a cos 1  b cos  2
Sacks, Zachary S., et al. "A perfectly matched anisotropic absorber for use as an absorbing
boundary condition." IEEE Trans. Antennas and Propagation, Vol. 43, No. 12, pp. 1460-1463, 1995.
Lecture 8
Slide 10
5
5/3/2015
Notes on a Single Interface
• It is a change in impedance that causes reflections
• Snell’s Law quantifies the angle of transmission
• Angle of transmission and reflection does not depend on polarization.
• The Fresnel equations quantify the amount of reflection and transmission
• Amount of reflection and transmission depends on the polarization
Lecture 8
Slide 11
Uniaxial Perfectly Matched Layer (UPML)
S. Zachary, D. Kingsland, R. Lee, J. Lee, “A Perfectly Matched Anisotropic Absorber for Use as an Absorbing Boundary Condition,” IEEE Trans. on Ant. and Prop., Vol. 43, No. 12, pp 1460–1463, 1995.
Lecture 8
Slide 12
6
5/3/2015
Boundary Condition Problem
If we model a wave hitting some device or object, it will scatter the applied wave into potentially many directions. We do NOT want these scattered waves to reflect from the boundaries of the grid. We also don’t want them to reenter from the other side of the grid (periodic boundaries).
?
?
How do we prevent this?
x
y
Lecture 8
Slide 13
How We Prevent Reflections in Lab
In the lab, we use anechoic foam to absorb outgoing waves. Lecture 8
Slide 14
7
5/3/2015
Absorbing Boundary Conditions
We can introduce loss at the boundaries of the grid!
 0
Absorbing Boundary
x
y
 0
Lecture 8
Slide 15
Oops!!
But if we introduce loss, we also introduce reflections from the lossy regions!!
1  n    2
R
2
1  n    2
2
 0
x
y
 0
Lecture 8
Slide 16
8
5/3/2015
Match the Impedance
We need to introduce loss to absorb outgoing waves, but we also need to match the impedance to the problem space to prevent reflections.
introduce loss here
r   r  j r
adjust this to control impedance
Lecture 8
Slide 17
More Trouble?
By examining the Fresnel equations, we see that we can only prevent reflections from the interface at one frequency, one angle of incident, and polarization.
rTE 
2 cos 1  1 cos  2
cos  2
 0   2  1
cos 1
2 cos 1  1 cos  2
rTM 
2 cos  2  1 cos 1
cos 1
 0   2  1
1 cos 1  2 cos  2
cos  2
Lecture 8
Slide 18
9
5/3/2015
Anisotropy to the Rescue!!
It turns out we can prevent reflections at all angles and for all polarizations if we allow our absorbing material to be doubly‐
diagonally anisotropic.
  and   
x
y
  and   
Lecture 8
Slide 19
Problem Statement for the PML
   ,  
0%
Free Space
0 ,  0
Lecture 8
100%
1
1
2
z
Slide 20
10
5/3/2015
Designing Anisotropy for Zero Reflection (1 of 3)
We need to perfectly match the impedance of the grid to the impedance of the absorbing region.



everywhere
One easy way to ensure impedance is perfectly matched is:
a 0 0
 s    r    r    0 b 0
 0 0 c 
Lecture 8
Slide 21
Designing Anisotropy for Zero Reflection (2 of 3)
If we choose , then the refraction equation reduces to
bc  1
sin 1  bc sin  2  sin  2
 1   2
No refraction!
The reflection coefficients now reduce to
rTE 
a cos 1  b cos  2
a b

a cos 1  b cos  2
a b
rTM 
 a cos 1  b cos  2  a  b

a cos 1  b cos  2
a b
These are no longer a function of angle!! 
Lecture 8
Slide 22
11
5/3/2015
Designing Anisotropy for Zero Reflection (3 of 3)
ab
If we further choose , the reflection equations reduce to
rTE 
a b
0
a b
rTM 
 a b
0
a b
Reflection is always zero regardless of frequency, angle of incidence, or polarization!! 
Recall the necessary conditions:
bc  1 and a  b
Lecture 8
Slide 23
The PML Parameters (1 of 3)
So far, we have
a 0 0
 s    r    r    0 b 0
 0 0 c 
ab
1
c
Thus, we can write our PML in terms of just one parameter sz.
 sz
 sz    0
 0
0
sz
0
0
0 
sz1 
sz    j 
This form of tensor is why we call this a uniaxial PML.
This is for a wave travelling in the +z direction incident on a z‐axis boundary.
Lecture 8
Slide 24
12
5/3/2015
The PML Parameters (2 of 3)
We potentially want a PML along all the borders.
 sx1
 sx    0

0
0
sx
0
0

0
sx 
sy

 s y    0
0

0
s y1
0
0

0
s y 
 sz
 sz    0
 0
0
0 
sz1 
0
sz
0
These can be combined into a single tensor quantity.
 s y sz

 sx







s
s
s
s
   x  y  z  0


 0

0
sx sz
sy
0

0 


0 


sx s y 
sz 
Lecture 8
Slide 25
The PML Parameters (3 of 3)
z
The 3D PML can be visualized this way…
x
sz  1
 s y sz

 sx

s   0


 0

Lecture 8
0
sx sz
sy
0

0 


0 


sx s y 
sz 
sy  1
y
sx  1
sx  1
sy  1
sz  1
Slide 26
13
5/3/2015
UPML in Cylindrical and Spherical Coordinates
Cylindrical Coordinates
  sz
 s
 

s   0


 0

0

ss
 z 
0

0 


0 

 s 

 sz 
Spherical Coordinates
 r 2 1
 
  r  sr
 s    0
 0


0
sr
0

0

0

sr 


F. L. Teixeira, W. C. Chew, “Systematic Derivation of Anisotropic PML Absorbing Media in Cylindrical and Spherical Coordiantes,” IEEE Microwave and Guided Wave Letters, Vol. 7, No. 11, pp. 371‐373, 1997.
Lecture 8
Slide 27
Stretched Coordinate Perfectly Matched Layer (SC‐PML)
Lecture 8
Slide 28
14
5/3/2015
The Uniaxial PML
Maxwell’s equations with uniaxial PML are:


  H  k0  r  s  E


  E  k0  r  s  H
 s y sz

 sx

s   0


 0


0 


0 


sx s y 
sz 
0
sx sz
sy
0
Lecture 8
Slide 29
Rearrange the Terms
We can bring the PML tensor to the left side of the equations and associate it with the curl operator.




1
1

 s    E  k0   r  H
 s    H  k0  r  E
The curl operator is now
 sz1s y1sx
1
 s     0
 0


0

 s
  sxy s1z z
 sz 1 
  sx s y y
 
Lecture 8
0
1
1
z y x
s s s
0
 ssxy
 
 
1 
s z z
0
sz
sy
 z
 0
 
0   z
1 1  
sz s y s x    y
0
 
1 
s x x
sx
sz
s
0

x
  
 szy


 
0 

y

x
1 
s y y

1 
s x x
0



Slide 30
15
5/3/2015
“Stretched” Coordinates
Our new curl operator is

0

1
 s     ssxy s1z z
 sz 1 
  sx s y y
s
 s xy
 
 
sz
sy
  
s
 szy
0
 
sx
sz
1 
s z z
 
1 
sx x
1 
s y y
 
1 
s x x
0


The factors sx, sy, and sz are effectively “stretching” the coordinates, but they are “stretching” into a complex space.
Lecture 8
Slide 31
Drop the Other Terms
We drop the non‐stretching terms.

0

 s
 s    sxy s1z z

  ssz s1 y
 x y

sx
sy
 
 
 
0
sz
sy
sx
sz
1 
s z z
 
1 
s x x

sy
sz
  
1 
s y y
 
1 
s x x
0
 0
  1 
   sz z
  1 
   s y y

 s1z

z
0
1 
s x x


 s1x x 

0 

1 
s y y
Justification
sy
sx
 
1 
s z z
1 
s z z
Inside the z‐PML, sx = sy = 1. This is valid everywhere except at the extreme corners of the grid where the PMLs overlap.
This also implies that the UPML and SC‐PML have nearly identical performance in terms of reflections, sensitivity to angle of incidence, polarization, etc.
Lecture 8
Slide 32
16
5/3/2015
Theoretical Performance
Given the following choice of PML parameters
s 
1 
1 
1 
aˆ x 
aˆ y 
aˆ z
sx x
s y y
sz z
sx  x   1  j
 x  x
 0
  y
sy  y   1  j y
 0
sz  z   1  j
 z  z
 0
 x 

 Lx 
m
 x  x    x ,max  
 y
 y  y    y ,max  
 Ly




 z 

 Lz 
m
m
 z  z    z ,max  
We choose i,max to achieve a target maximum reflectance R at normal incidence according to
 i ,max  
 m  1 ln R
We typically choose
20 Li
3 m 4
 i ,max 
4
0  i
Lecture 8
Slide 33
Calculating the PML Parameters
Lecture 8
Slide 34
17
5/3/2015
The Perfectly Matched Layer (PML)
The perfectly matched layer (PML) is an absorbing boundary condition (ABC) where the impedance is perfectly matched to the problem space. Reflections entering the lossy regions are prevented because impedance is matched. Reflections from the grid boundary are prevented because the outgoing waves are absorbed.
 y  0
PML
 x   y   z  0
Problem
Space
PML
PML
 x  0
 x  0
PML
 y  0
Lecture 8
Slide 35
How to Calculate the PML Parameters
Maxwell’s Eqs. with PML


  E  k0   r  s  H


  H  k0  r  s  E
 s y sz

 sx

 s   0


 0

0
sx sz
sy
0

0 


0 


sx s y 
sz 
NGRID = [Nx Ny];
NPML = [0 0 20 20];
[sx,sy] = calcpml2d(NGRID,NPML);
Computing PML Parameters
sx  x   ax  x  1  j0 x  x  
s y  y   a y  y  1  j0 y  y  
sz  z   az  z  1  j0 z  z  
ax  x   1  amax   x Lx 
p
a y  y   1  amax   y Ly 
az  z   1  amax   z Lz 
0  amax  5
3 p 5
 1
 max
p
p
 x 

 2 Lx 
 y 
 sin 2 
 y  y    max
 2 Ly 


2 z 
 sin 
 z  z    max

 2 Lz 
 sin 2 
 x  x    max
Writing this function will be in homework 
Lecture 8
Slide 36
18
5/3/2015
Visualizing the PML Loss Terms – 2D
For best performance, the loss terms should increase gradually into the PMLs.
 x  x 
x
 y  y 
y
Lecture 8
Slide 37
Note About x/Lx, y/Ly, and z/Lz
The following ratios provide a single quantity that goes from 0 to 1 as you move through a PML region.
x
Lx
and
y
Ly
and
z
Lz
x, y, z  position within PML
Lx , Ly , Lz  size of PML
We can calculate the same ratio using integer indices from our grid.
nx
nx
x

or
NXHI
Lx NXLO
nx = 1, 2, ... , NXLO
nx = 1, 2, ... , NXHI
0
ny
ny
y
or

NYHI
Ly NYLO
ny = 1, 2, ... , NYLO
x
1
Lx
ny = 1, 2, ... , NYHI
0
nz
nz
z
or

NZHI
Lz NZLO
Lecture 8
nz = 1, 2, ... , NZLO
x
1
Lx
nz = 1, 2, ... , NZHI
Slide 38
19
5/3/2015
Lecture 8
Nx
Nx-NXHI+1
sx(x) = 1
NXLO
1
% ADD XHI PML
for nx = 1 : NXHI
…
sx(Nx-NXHI+nx,:) = ...
end
sx(x) = ax(x)[1+j0x(x)]
% ADD XLO PML
for nx = 1 : NXLO
…
sx(NXLO-nx+1,:) = ...
end
sx(x) = ax(x)[1+j0x(x)]
Visualizing sx in 2D
Slide 39
Visualizing sy in 2D
1
sy(y) = ay(y)[1+j0y(y)]
% ADD YLO PML
for ny = 1 : NYLO
…
sy(:,NYLO-ny+1) = ...
end
NYLO
sy(y) = 1
Ny – NYHI + 1
Ny
Lecture 8
sy(y) = ay(y)[1+j0y(y)]
% ADD YHI PML
for ny = 1 : NYHI
…
sy(:,Ny-NYHI+ny) = …
end
Slide 40
20
5/3/2015
Procedure for Calculating sx and sy on a 2D Grid
1. Initialize sx and sy to all ones.
s x  x, y   s y  x, y   1
2. Fill in x‐axis PML regions using two for loops.
NXLO
NXHI
3. Fill in y‐axis PML regions using two for loops.
NYLO
NYHI
Lecture 8
Slide 41
Example Data for 2D
NGRID
= [7 4];
NPML
= [2 3 1 2];
[sx,sy] = calcpml2d(NGRID,NPML);
sx =
+ 1.5069i
+ 0.2590i
+ 0.1046i
+ 0.5337i
+ 1.5069i
0.0040
0.0014
0.0010
0.0010
0.0011
0.0019
0.0040
+ 1.5069i
+ 0.2590i
+ 0.1046i
+ 0.5337i
+ 1.5069i
0.0040
0.0014
0.0010
0.0010
0.0011
0.0019
0.0040
+ 1.5069i
+ 0.2590i
+ 1.5069i
+ 0.2590i
+ 0.1046i
+ 0.5337i
+ 1.5069i
0.0040
0.0014
0.0010
0.0010
0.0011
0.0019
0.0040
0.0014
0.0014
0.0014
0.0014
0.0014
0.0014
0.0014
+
+
+
+
+
+
+
0.0040
0.0040
0.0040
0.0040
0.0040
0.0040
0.0040
+
+
+
+
+
+
+
+ 0.1046i
+ 0.5337i
+ 1.5069i
1.0e+03 *
0.0040
0.0040
0.0040
0.0040
0.0040
0.0040
0.0040
Lecture 8
p3
 1
 max
1.0e+03 *
0.0040
0.0014
0.0010
0.0010
0.0011
0.0019
0.0040
sy =
amax  3
+
+
+
+
+
+
+
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
0.2590i
0.2590i
0.2590i
0.2590i
0.2590i
0.2590i
0.2590i
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
1.5069i
Slide 42
21
5/3/2015
PML is Not a Boundary Condition
A numerical boundary condition is the rule you follow when an equation references a field from outside the grid.
The PML does not address this issue.
It is simply a way of incorporating loss while preventing reflections so as to absorb outgoing waves.
Sometimes it is called an absorbing boundary condition, but this is still misleading as the PML is not a true boundary condition.
Lecture 8
Slide 43
Incorporating a UPML into Maxwell’s Equations
Lecture 8
Slide 44
22
5/3/2015
Maxwell’s Equations with a UPML
Maxwell’s equations before the PML is added are


  E  k0   r  H


  H  k0  r  E
 xx
 r    yx
 zx

 xy  xz 

 yy  yz 
 zy  zz 
  xx
 r     yx
  zx

 xy
 yy
 zy
 xz 

 yz 
 zz 
The UPML is incorporated as follows.


  E  k0   r  s  H


  H  k0  r  s  E
 s y sz

 sx

 s   0


 0

0
sx sz
sy
0

0 


0 


sx s y 
sz 
Lecture 8
Slide 45
Vector Expansion
Assuming only diagonal tensors
 xx
 r    0
 0
0
 yy
0
0
0 
 zz 
  xx
 r    0
 0
0
 yy
0
0 
0 
 zz 
Maxwell’s equations expand to
Ex Ez
ss

 k0  yy x z H y
z
x
sy
ss
H z H y

 k0 xx y z Ex
y
z
sx
H x H z
ss

 k0 yy x z E y
z
x
sy
E y
H y
ss
Ez E y

 k0  xx y z H x
y
z
sx
x
Lecture 8

ss
Ex
 k0  zz x y H z
y
sz
x

ss
H x
 k0 zz x y Ez
y
sz
Slide 46
23
5/3/2015
Absorb UPML into  and (3D Grid)
We can absorb the UPML parameters into the material functions.
 xx   xx
s y sz
 xx   xx
sx
sx
ss
 yy   yy x z
sy
ss
 yy   yy x z
sy
 zz   zz
s y sz
sx s y
 zz   zz
sz
sx s y
sz
We can now write Maxwell’s equations as
Ez E y
 H x

 k0  xx
y
z
Ex Ez

 k0  yy H y
z
x
E y Ex

 k0  zz H z
x
y
H z H y

 k0 xx Ex
y
z
H x H z

 k0 yy E y
z
x
H y H x

 k0 zz Ez
x
y
This means we can formulate a code as if there was no PML. All we have to do is modify the materials being modeled near the boundaries.
Lecture 8
Slide 47
Absorb UPML into  and (2D Grid)
Let z be the uniform direction, then d/dz = 0 and sz = 1.
We can still absorb the UPML parameters into the material functions.
 xx   xx
sy
 xx   xx
sx
sy
sx
sx
sy
s
 yy   yy x
sy
 yy   yy
 zz   zz sx s y
 zz   zz sx s y
We can now write Maxwell’s equations as
H Mode
E Mode
H y
H x

 k0 zz Ez
x
y
Ez
 H x
 k0  xx
y
E
 z  k0  yy H y
x
Lecture 8
E y
E
 x  k0  zz H z
x
y
H z
 Ex
 k0 xx
y
H z

 k0 yy E y
x
Slide 48
24
5/3/2015
Incorporating a SC‐PML into Maxwell’s Equations
Lecture 8
Slide 49
Maxwell’s Equations with a SC‐PML
Maxwell’s equations before the PML is added are


  E  k0   r  H


  H  k0  r  E
 xx  xy  xz 


 r   yx  yy  yz 
 zx  zy  zz 


  xx


 r    yx
  zx

 xy
 yy
 zy
 xz 

 yz 
 zz 
The SC‐PML is incorporated as follows.


 s  E   j    H


 s  H  j   E
Lecture 8
 0

 s    s1z z
 1 
  s y y
 s1z

z
0
1 
s x x


 s1x x 

0 

1 
s y y
Slide 50
25
5/3/2015
Vector Expansion
Maxwell’s equations with a SC‐PML expand to
Fully Anisotropic

1 H z 1 H y

 k0   xx Ex   xy E y   xz Ez 
s y y sz z
Diagonally Anisotropic
1 H z 1 H y

 k0 xx Ex
s y y sz z
1 H x 1 H z

 k0   yx Ex   yy E y   yz Ez 
sz z sx x
1 H y 1 H x

 k0   zx Ex   zy E y   zz Ez 
sx x s y y
1 H x 1 H z

 k0 yy E y
sz z sx x
1 H y 1 H x

 k0 zz Ez
sx x s y y
1 Ez 1 E y

 k0  xx H x   xy H y   xz H z
s y y sz z

1 Ez 1 E y

 k0  xx H x
s y y sz z
1 Ex 1 Ez

 k0  yx H x   yy H y   yz H z
sz z sx x

1 Ex 1 Ez

 k0  yy H y
sz z sx x
1 E y 1 Ex

 k0  zx H x   zy H y   zz H z
sx x s y y

1 E y 1 Ex

 k0  zz H z
sx x s y y



Lecture 8
Slide 51
UPML Vs. SC‐PML
Lecture 8
Slide 52
26
5/3/2015
UPML Vs. SC‐PML
Uniaxial PML
Stretched‐Coordinate PML
Benefits
Benefits
• Has a physical interpretation
• Models can be formulated and implemented without considering the PML
• Less computationally intensive
• More efficient implementation in the time‐domain
• Matrices are better conditioned.
Drawbacks
Drawbacks
• Can be more computational intensive to implement
• Resulting matrices are less well conditioned
• Must be accounted for in the formulation and implementation of the numerical method.
• Not intuitive to understand
Lecture 8
Slide 53
27