Transcript MA375 - Rice U - Computational and Applied Mathematics

MA557/MA578/CS557
Lecture 21
Spring 2003
Prof. Tim Warburton
[email protected]
1
Recap
•
So far we have derived the advection-(diffusion) equation in a onedimensional domain.
•
We created a numerical scheme based on piecewise polynomial
approximations.
•
Boundary conditions between sub-intervals are imposed by the use of
fluxes.
•
We proved stability and consistency for these schemes.
•
But…. we do not live in a one-dimensional world
•
So – on the march to our three-dimensional world we are going to
investigate PDE’s with two spatial-dimensions.
2
Two-Dimensional Polynomial Space
• For the following we will use a two-dimensional space
which is spanned by the polynomials of total order not
greater than p:
P p    x a y b
0a b p
• Later we will choose an orthonormal basis to discretize
this space.
3
Two-Dimensional Advection
• Recall in one-dimensions we chose an arbitrary section of a
pipe.
• We monitored the flux of fluid into and out of the ends of the
pipe.
• The conservation equation in 1D we derived was:
b
d
C  u  b, t  C  b, t   u  a, t  C  a, t 

dt a
• Now we start with a 2D domain and consider any simply
connected sub-region:
4
Our Domain 

5
Our Sub-region



6
Outward Normal n

n

7
Conservation
• Suppose the fluid has a mean velocity
• This time
u is
u everywhere.
a two-vector.
n=outward pointing normal

u
• The flux of fluid through the boundary is the integral of the
concentration C being advected normal to the surface of w:
flux 
   u  n  CdS

8
Conservation Law
• Now we know how much of the concentrate is being
advected through the boundary of our sub-region.
• The conservation law is now:
d
CdV     u  n  CdS

dt 

• We apply the divergence theorem (beware this relies on
smoothness arguments):
d
CdV     u  n  CdS

dt 

     uC  dV

9
In Component Form
• Assuming ubar is constant this becomes:
u 
u 
v 
  
 x 
 
  
 x 
d
CdV     u   CdV

dt 

 
 
    u  v  CdV
y 
  x
10
Finalize
• For all sub-regions:
 
d
 
CdV    u  v  CdV

dt 
x
y 

 C
C
C 
 
u
v
 dV  0
t
x
y 

• Assuming enough continuity we obtain:
C
C
C
u
v
0
t
x
y
11
Straight To DG
• The pde is:
C
C
C
u
v
0
t
x
y
• We substitute the DG scheme (using Lax-Friedrichs fluxes):
Find C  P     0, T  such that   P  
p
p




 
 C
C
C 
  , t  u x  v y     ,  2

  
where   u  n ,


C   C

C


 C    0





12
+/- Notation
• For the purpose of the numerical scheme notation:
– The + notation implies the limit of C from outside the sub-region
– The - notation implies the limit of C from inside the sub-region
C+
n
C-


C

C

C
 
13
Sample Domain Partition
Suppose we have a domain:
Which we partition into 5 triangles:
5
3
2
4
1
We note the unique edges, and
arbitrarily label their edge limits
with +/-
3- 3+
44+
55+
2+
21+
114
Stability
Set   C







 C
C
C 


C
,

u

v

C
,
 t


x
y  
2







 C  C    0





  u n  
1 d
C 
 C , C    C ,
2 dt
2




    
 C ,


2






 C  C    0





  
1 d

 C , C    C , C 
2 dt
2




    
 C ,


2




 C   C     0




15

j
Suppose We Partition  

  
   
1d

 C , C  j   C , C    C , 

2
2
 j 2 dt

 j 



j




 C  C    0



 j


    
   
1d


C
,
C

C
,
C

C


j 2 dt
j  2   ,  2
j

 j 





 C  C    0



 j

1d
j 2 dt  C , C  j




  
  



C
,
C

C
,
C
 
 

2
2



e 
e

 











 C   C      0
    C  , 

 
2 
unique  



e
edges




 
       

 C  C   
  C , 

 
2 
 


e 
16
Continuing
1d
j 2 dt  C , C  j




  
  



C
,
C

C
,
C
 
 

2
2



e 
e

 


 
       
 C  C      0
   C , 

 
2 
unique  



e
edges e




 
        

 C  C   
  C , 


 
2
 


e 

1d
j 2 dt  C , C  j

    
 
 C  , 
C    C ,

 2   
 2
  

 e 

unique 



edges e
 C ,    C 
e

 


  
C  
  
 e 



Where we have gathered like terms
(note we ignored the boundary terms)
17
Finally -- Stability
1d
j 2 dt  C , C  j

    
    
 C  , 
C    C ,
C  

 2   
 2   
   

 e 

 e 
unique 




edges e
 C ,    C 

e


   





 C  C  , C  C  
   

unique   2 



e
edges e
 



d
C
dt
2

2
L


j

j 




unique
edges e
u n
2
C  L2 e  0
2
Where we have gathered like terms
(note we dropped the boundary terms)
18
Now The Details
• Recall in 1D we were restricted to breaking the interval
[a,b] into segments.
• In 2D there are many more choices.
• Two obvious shapes for sub-regions are the triangle and
quadrilateral.
• There is plenty of software available for mesh generation
using triangles.
• Building quadrilateral meshes is slightly tougher for a
complex domain.
19
Mesh Generation
• For simplicity we will use the Triangle mesh generator.
• This has been built into WinUSEMe – review the previous
notes.
• For now assume that the task of creating a mesh has
been done – and – the mesh is held in a file.
• Note that for the DG we do not need to much information
about how the element are connected – it is mostly
sufficient to know which edge connects to which edge.
• We also need to know which faces lie on a domain
boundary.
20
Reference Triangle
• The following will be our basic triangles:
(-1,1)
s
r
(-1,-1)
(1,-1)
• All straight sided triangles are the image of this triangle
under the map:
1
2
3





v
v
v
 x
 r  s  x  1 r  x
 1 s  x 
 y     2   v1    2   v 2    2   v3 

 y  
 y  
 y 
 
21
Reference Triangle
v , v 
3
x
• The following will be our basic triangles:
(-1,1)
s
r
(-1,-1)
3
y
v , v 
1
x
1
y
v , v 
2
x
(1,-1)
• All straight sided triangles are the image of this triangle
under the map:
1
2
3





v
v
v
 x
 r  s  x  1 r  x
 1 s  x 
 y     2   v1    2   v 2    2   v3 

 y  
 y  
 y 
 
22
2
y
Orthonormal Basis for the Triangle
• Fortunately an orthonormal basis for the triangle has been
discovered (and rediscovered multiple times).
• First we need to know some details about the Jacobi
polynomials. These polynomials are parameterized by two reals:
 ,  and their integer orders n,m such that they satisfy the
orthogonality relationship (for integer alpha,beta):
n   ! n   !

2
 1  x   1  x   ,
 ,
1 2   2  Pn  x  Pm  x dx   nm 2n      1 n! n     !
1


• Basic identity:
 ,
Pn
 ,
Pn
1
n   !


n ! !
  x    1
n
Pn ,   x 
23
Rodrigues’ Formula
• The Jacobi polynomials can be generated in a similar way
to the Legendre polynomials:
Pn ,  x  
1
 1 2n n!1  x  1  x 
n



dn


2 n
1  x  1  x  1  x 
n 
dx
• Note the Legendre polynomials are a special case of the
Jacobi polynomials:
Ln  x   Pn0,0  x 
24

Recurrence Relation
P0 ,  x   1
P1 ,  x  
1
2   1      2  x  1 

2
2  n  1 n      1 2n      Pn1,  x  

2n      2 !  ,

2
2
 2n      1      x
 Pn  x 
 2n      1! 

 2  n    n    2n      2  Pn1,  x 
See:
http://www.math.unm.edu/~timwar/MA578S03/MatlabScripts/JACOBI1D.m
25
Orthonormal Basis for the Triangle
• The following basis is due to Koornwinder (later revived by
Proriol, Dubiner, Owens,….)
2(1  r )
a
1
(1  s )
bs
 1  b  2 n1,0
 nm   r , s   Pn0,0  a  
b 
 Pm
 2 
n
• Part of HW7 is to show that this is an orthogonal basis in
the sense that:
1 s
2
2
1 1ij  kl drds   ik 2i  1 jl 2 j  2i  2
26
Homework 7
Q1) Prove that the Koornwinder-Proriol-Dubiner-Owens… basis is
orthogonal on the triangle: T={-1<=r,s;r+s<=0}:
1 s
2
2
1 1ij  kl drds   ik 2i  1 jl 2 j  2i  2
(HINT: under the mapping (r,s)->(a,b) the image of the triangle T is the
square Q={-1<=a,b<=1} and the Jacobian is (1-b)/2 )
Q2) Determine an orthonormal basis
 
nm
0 n  m  p
for the triangle by normalizing the basis functions of the KPDO.. basis
Q3) Construct a p’th order Vandermonde for the equally spaced points on
the triangle (take half of the square of (p+1)x(p+1) points) using the
orthonormal KPDO basis. In Matlab calculate the condition number of
this Vandermonde basis as a function of p.
27
HW7 cont
Q3cont) Use the lower left triangle of points:
28
HW7 cont
Q3cont) The Vandermonde matrix V will be have
(p+1)(p+2)/2 rows and columns.
Vi nm   nm  ri , si 
p  1 p  2 

where 1  i 
2
and 0  n,m;n+m  p
Q4) Extra credit: repeat 3 using half of the tensor product of
the one-dimensional Chebychev nodes.
For Q3 and (Q4) plot the polynomial order on the horizontal
axis and cond(V) on the vertical axis
29
HW7 cont
Q4cont) Use the lower left triangle of points:
30
Next Class
• We will create the derivative and surface matrices used in
the DG scheme.
• We will transform the modal scheme into a nodal scheme.
• We will use a better set of nodes (improve the condition
number of the Vandermonde matrix).
• Time permitting we will prove consistency.
• For more info on Jacobi polynomials:
http://mathworld.wolfram.com/JacobiPolynomial.html
31