Spring 2003

1
MA557/MA578/CS557Lecture 32

• Spring 2003
• Prof. Tim Warburton
• timwar_at_math.unm.edu

2
Heat Equation

• Recall Lecture 17 where we derived the
  advection-diffusion equation in 1D.
• The same type of derivation for the diffusion
  process of a concentrate C yields in 2D

3
Discretization of the Laplacian

• The LDG, semidiscrete, equation satisfies
• Zero Dirichlet boundary conditions (concentration
  pinned to zero at boundary)
• Zero Neumann boundary conditions (insulated
  boundaries)

4
Heat Equation

• Our first approach will be an explicit time
  integration of the equation (using zero Neumann
  bc)

5
Explicit Numerical Heat Equation

• The trouble is that the discrete diffusion matrix
  has a large spectral radius
• We do know that (using the integration by parts
  formula for the DG operator) that the operator is
  non-positive.
• We also need to fit the eigenspectrum in the
  appropriate stability region

6
(No Transcript)
7
We also know that all the eigenvalues are
non-positive real numbers and dtrho(A) must fit
in the stability region..
8
Comment

• We have shown that using an explicit time
  integration scheme may require very small time
  steps
• The root cause of this is that information
  diffuses through the data domain very quickly..
• We will investigate some schemes which solve
  globally coupled systems.

9
Backwards Euler Time Integration

• Suppose we consider the backwards Euler scheme
• Iterator

10
Backwards Euler Time Integration

• Suppose we consider the backwards Euler scheme
• Thus the only conditions for stability are
• However, accuracy depends on dt. The
  discretization of the time derivative has a first
  order truncation error

11
Crank-Nicholson

• We replace the right hand side of the discrete
  scheme with

12
Diagonalize

• If we apply a change of basis so that the
  inv(M)L is diagonal then we decouple the
  iterator into a sequence of independent scalar
  iterators

13
Quick Result

• We can bound the iteration multiplier
• As long as both the eigenvalues and dt are
  non-negative then the multiplier is bounded by
  one and the scheme is stable.
• We have previously shown that L is a non-negative
  operator.

14
Crank-Nicholson Accuracy

• We can briefly compute the temporal accuracy

15
ESDIRK4

• A new implicit RK scheme with a lower triangle
  Butcher tableau has been proposed by Carpenter et
  al. http//fun3d.larc.nasa.gov/papers/carpenter.pd
  f
• The linear system to invert at each sub stage is
  the same.
• The RK scheme comes with an embedded scheme to
  provide sth order and (s-1)th order time
  approximations

16
4th Order ESDIRK4 Implicit RK Schemeinitializer
5 stages

• Iterator (Ctilde are intermediate values of C,
  Chat6 is the embedded 3rd order approximation of
  Cn1)

17
4th Order ESDIRK4 Implicit RK Scheme

• Iterator (Ctilde are intermediate values of C,
  Chat5 is the embedded 3rd order approximation of
  Cn1)

18
4th Order ESDIRK4 Implicit RK Scheme

• If L is a linear operator then

19
Coefficients

• Gamma 1/4
• a21 1/4
• a31 8611/62500
• a32 -1743/31250
• a41 5012029/34652500
• a42 -654441/2922500
• a43 174375/388108
• a51 15267082809/155376265600
• a52 -71443401/120774400
• a53 730878875/902184768
• a54 2285395/8070912

b1 82889/524892 b2 0 b3
15625/83664 b4 69875/102672 b5 -2260/8211
bhat1 4586570599/29645900160 bhat2 0 bhat3
178811875/945068544 bhat4 814220225/1159782912
bhat5 -3700637/11593932 bhat6 61727/225920
20
Implementation Details

• As they are computed we are required to store
  Ctilde1, Ctilde2, , Ctilde6 until the end of the
  6th ESDIRK stage.
• Each stage 2,3,4,5,6 requires the solution of the
  same linear system, with different data on the
  right hand side.
• Chat6 may be computed without solving a linear
  system.
• The systems can be solved using an iterative (say
  conjugate gradient) method. i.e. the L matrix
  does not need to be computed and stored.
• We need to be able to compute (MdtgammaL)x
  for a given vector x.
• Ctilde6-Chat6 will give an estimate of the
  difference between a 3rd order and the 4th order
  approximation in time. This can be used to
  estimate the error made in this time step and can
  suggest a change in dt.

21
Solving the System

• We are going to use an iterative method to solve
  the following sequence of symmetric matrix
  problems

22
Summary of Temporal Implicit Schemes

• Backwards Euler is unconditionally stable for
  non-negative diffusion parameter D (i.e. any
  dtgt0) and first order in dt.
• Crank-Nicholson is unconditionally stable for
  non-negative diffusion parameter D (i.e. any
  dtgt0) and second order in dt.
• ESDIRK4 generalizes to fourth order in dt.

23
Homework 10

• Q1) Modify/use the umDIFFUSION.zip files to solve
  the diffusion equation with the following time
  integrators
• a) Backwards Euler
• b) Crank-Nicholsonc) ESDIRK4
• Q2) Create a domain and determine convergence
  rates in h for p1,3,5,7 for a small dt (i.e.
  make sure the convergence does not bottom out
  above 1e-10). Repeat for two of the above time
  integrators.
• Q3) Choose p7 and h small and verify rates of
  convergence in dt for integration to some fixed
  time T. Use the one time integrator you did not
  use in Q2.

24
Homework 10

• This homework can be completed in pairs or
  individually.
• This homework is due Monday 04/21/03
• Remember no more than 5 email questions will be
  replied to for this homework.

25
Driver for Backwards Euler Diffusion
Scheme(umDIFFUSIONdemo.m)
26
Time Loop (umDIFFUSIONrun.m)

• 1-12) Set parameters
• 14-16) Set initial conditions
• 19) Set solver parameters
• 22-40) Time stepping loop
• 28) Solves linear system using conjugate
  gradient
• 32-40) Plot solution and compute integral
  of C

27
umDIFFUSIONop.m

• Set up storage for computing (MC dtgammaD
  LC)

28
umDIFFUSIONop.m

• 15-25) compute qx,qy using Neumann bcs
• 27-40) compute div(qx,qy) using Dirichlet
  bcs
• 42-48) compute (MCgammadtDLC)

29
Demo Run
30
Discussion of Final Project

• Now is the time to ready your final projects.
• Submit a one page description of the PDE you
  intend to discretize with DG by 04/21/03
• This is a proposal and may be rejected
  requiring a resubmission or assignment of a set
  of PDEs.
• Include
• List of PDEs to be discretized
• List boundary conditions and initial conditions
  to be discretized
• Relevance to your own research (if any)
• List of group members (max 3, with one proposal
  per group)
• Interesting issues (application of PML, creation
  of specific PML, specific physical
  application/model)