Transcript MA375 - Rice U - Computational and Applied Mathematics

MA557/MA578/CS557
Lecture 26
Spring 2003
Prof. Tim Warburton
[email protected]
1
Note on Homework 8
When you enter h: 0.10 you are actually specifying
Max triangle size = 0.10
So for the order of accuracy in this case you would
take h=sqrt(0.1)
2
Results
• I ran the advection code in a square domain:
• [-4,4]x[-4,4]
• With hash grids like – here h=1
• Solution used:
2
• dt = .125/((umNq)*(umNq+1)*max(Fconst(:)));
10  x t   y 2 
C e
• 5th order Rk scheme used
3
Error in L2 norm
p\h
2
1
0.5
-------------------------------------------------------------------------------------------1 0.93877444976985 0.38815661859066 0.10197215280871
2 0.30215330223980 0.15653021172476 0.02800152435822
3 0.22043684447109 0.07272803132766 0.00626031860408
4 0.11652413181376 0.02876534888079 0.00131839250177
5 0.11289073566956 0.01460227046098 0.00033161249084
6 0.06133875385613 0.00423786600885 0.00003293355651
7 0.04835896860374 0.00184010232130 0.00001364982037
8 0.03584481563628 0.00037174570177 0.00000091542340
-------------------------------------------------------------------------------------------• Note from h=2 to h=1 the convergence is slow (pre-exponential).
• From h=1 to h=0.5 convergence rate is pretty good.
• We estimate the error as: Err = Ch^s
4
P=8,h1=1,h2=0.5
• Set:
err (h1 )  Ch1
s
err (h2 )  Ch2
s
 err (h1 ) 
log 

err
(
h
)
2 

s
 h1 
log  
 h2 
 0.00037174570177 
log 

• Then
0.00000091542340

  8.67
s
 1 
log 

0.5


• So we get approximately 8.5 order accuracy
5
Just Using the h=1,h=0.5 Errors
P approx. order
---------------------------1 1.93
2 2.48
3 3.54
4 4.45
5 5.46
6 7.01
7 7.07
8 8.67
-----------------------------Experimentally (apart from p=6) we find that the rate
seems to be better than expected (we estimated a rate
of sigma-1=p since the solution is regular)
6
Linear Systems of 2D PDEs
• A general linear pde system will look like:
Find q :   [0, t ] 
d
where d  dim(q) such that:
q
q
q
A B
0
t
x
y
7
Lax-Friedrichs Fluxes
• Suppose we know the maximum wave speed of the
q
q
q
hyperbolic system:
A B
0
t
x
y
 are
• As long as the matrix A, B 
codiagonalizable for any linear combination.
n
n
• i.e. C   A   B Has real eigenvalues for any
real alpha,beta.
• And:
max  | C , x   x for some x 

n

,  ,  |  2   2  1  
• then we can use the following DG scheme:
8
Lax-Friedrichs DG Scheme
For Linear Systems
Find q : P p T   [0, t ] 
d
such that   P p T 
 

 q
q
q    Anx  Bn y  
q      , q    0

  , t  A x  B y     , 
2

T  
 T  2
T
where q   q   q 
9
Discrete LF/DG Scheme
• We choose to discretize the Pp polynomial space on
each triangle using M=(p+1)(p+2)/2 Lagrange
interpolatory polynomials.
• i.e.
P p T   hm 1mM
• The scheme now reads:
Find q m  t  , m  1,..., M  such that:
m M
 h , h 
m 1
n
m T

 mM   hm 
dq m mM   hm 
   A  hn ,
q m .....
 q m    B  hn ,

dt
x T  m1  
y T 
m 1  
   Anx  Bn y    I  
  hn , 
 hm  q m   0
 
 
2

 T

10
Tensor Form
Find q mcj  t  , m  1,..., M , c  1,...,dim(q), j=1,...,#tri such that:
m M
 h ,h 
m 1
n
d dim( q )

d 1
d=dim(q )

d 1
d 1
dt
 ....
 m M  h 
 
A cd    hn , m  q mdj    ....
 m1 
x T j
 


 m M  h 
 

B cd    hn , m  q mdj    .....
y T j
 m1 
 


d dim( q )

m T
j
dq mcj
  A cd nx  B cd n y    I   m M




    hn , hm T q mdj    0

  m1
2


cd
11
Simplified Tensor Form
Find q mcj  t  , m  1,..., M , c  1,...,dim(q), j=1,...,#tri such that:
dq ncj
dt
 ....
d  dim( q )
 mM x
 d=dim(q )
 mM y

A
D
q

B
D
q


cd  
nm mdj 
cd  
nm mdj   .....
d 1
d 1
 m1

 m1

d  dim( q ) e 3  A n e  B n e   I 

 mM e

cd x
cd y 




S
q



nm  mdj    0



2
d 1
e 1


cd  m1
where:
M nm   hn , hm T
D
x
nm
k M
 hm 
x
1
x
  hn ,
,
D

M
D



nm

nk km
x T

k 1
D
y
nm
k M
 hm 
y
y
  hn ,
, Dnm    M 1  Dkm

nk
y T
k 1

S
e
nm
  hn , hm T , S
e
e
nm

k M
 M 
k 1
1
nk
S ekm
12
Geometric Factors
• We use the chain rule to compute the Dx and Dy
matrices:
r r s s
D  D  D
x
x
r r s s
y
D  D  D
y
y
x
• In the umSCALAR2d scripts the
umDr  D
r
umDs  Ds
r
s
r
s
rx  , sx  , ry  , sy 
x
x
y
y
13
Surface Terms
• In the Matlab code first we extract the nodes on the edges
of the elements:
• fC = umFtoN*C
• The resulting matrix fC has dimension
(umNfaces*(umP+1))xumNel
• This lists all nodes from element 1, edge 1 then element
1, edge 2 then…
• In order to multiply, say fC, by Se we use
umNtoF*(Fscale.*fC) where:
Fscale = sjac./(umNtoF*jac);
14
Example Matlab Code
(Upwind DG Advection Scalar Eqn)
• This is not strictly a global Lax-Friedrichs scheme
since the lambda is chosen locally!!!
15
Example Matlab Code cont
16
Specific Example
• We will try out this scheme on the 2D, transverse
mode, Maxwell’s equations.
17
Maxwell’s Equations (TM mode)
• In the absence of sources, Maxwell’s equations are:
Hx 
t

Ez
0
y
H y 
Ez

0
t
x
   Ez  H x H y


0
t
y
x
Hx 

H y 
x
  permittivity
• Where
y
0
  permeability
H  magnetic field
E  electric field
18
Free Space..
• For simplicity we will assume that the permeability
and permittivity are =1 we then obtain:
H x Ez

t
y
H y Ez

t
x
Ez H x H y


t
y
x
H x H y

x
y
0
0
0
0
19
Divergence Condition
• We next notice that by taking the x derivative of the
first equation and the y derivative of the second
equation we find that:
  H x H y


t  x
y

0

• i.e. the fourth equation (divergence condition) is a
natural result of the equations – assuming that the
initial condition for the magnetic field H is divergence
free.
20
The Maxwell’s Equations We Will Use
H x Ez

0
t
y
H y Ez

0
t
x
Ez H x H y


0
t
y
x
H y
H x
 x, y , t  0  
 x, y , t  0   0
x
y
21
The Maxwell’s Equations As Conservation Rule
H x 

  0    Ez   0
y
x
t
H y 

   Ez    0   0
y
x
t

Ez 
  H y    H x   0
y
t x
H y
H x
 x, y , t  0   0
 x, y , t  0  
y
x
22
TM Maxwell’s Boundary Condition
• We will use the following PEC boundary condition –
which corresponds to the boundary being a perfectly,
electrical conducting material:

x
H H

x
H y  H y
Ez   Ez
23
The Maxwell’s Equations In Matrix Form
q
q
q
A B
0
t
x
y
 Hx 
0 0 0 
0 0 1
q   H y  , A   0 0 1 , B   0 0 0 






E 
 0 1 0 
1 0 0
 z




• We are now assuming that the magnetic field is going
to be divergence free – so we just forget about it for
the moment.
24
Eigenvalues
• Then:
0 0 0 
0 0 1
0
A   0 0 1 , B   0 0 0   C   0





 0 1 0 
1 0 0






0
0

 
 

0 
• C has eigenvalues which satisfy:
 
0  C  I   0



 
  

  
0
 0  C   I       2   2      
   0,   2   2
• So the matrices are co-diagonalizable for all real alpha, beta and
    1 given     1
2
2
25
Project Time
• Due 04/07/03
• Code up a 2D DG solver for an arbitrary hyperbolic
system based on the 2D DG advection code
provided.
• Test it on TM Maxwell’s
– Set A and B as described, set lambdatilde = 1 and
test on some domain of your choosing.
– Use PEC boundary conditions all round.
• Test it on a second set of named equations of your
own choosing (determine these equations by
research)
• Compute convergence rates.. For smooth solutions.
26
Stability
•
If A,B are both symmetric then it is straight-forward
to prove stability
1) Set phi=q
2) Integrate by parts
3) Sum over all elements
4) Simplify
5) Repeat the approaches previously used..
27