Numerical Methods for Partial Differential Equations

1
Numerical Methods for Partial Differential
Equations

• CAAM 452
• Spring 2005
• Lecture 12
• Instructor Tim Warburton

2
Godunov Scheme Summary

• To complete this scheme we now specify how to
  compute the slopes.
• Standard formulae

3
With Limiting
Minmod slope limiterMonotonized
central-difference limiter (MC limiter)
4
Today

• More limiters
• Flux limiting function formulation.
• We will discuss Hartens sufficient conditions
  for a numerical method (including limiter) to be
  TVD
• Sweby TVD diagrams for flux limiting functions.
• Extension to systems of linear PDEs
• Extension to nonlinear PDEs

5
Flux Formulation with Piecewise Linear
Reconstruction

• Last time we showed how the ansatz of a piecewise
  linear reconstruction and Godunovs method
  allowed us to compute the time averaged flux
  contribution at each end of the Ith cell

Notice we can obtain the i-1/2 flux by setting
i-gti-1 in the i1/2 flux formula (i.e. the flux
formula is continuous at the cell boundary)
6
cont

• Using this notation the scheme becomes
• This is known as the flux formulation with
  piecewise reconstruction.

7
cont

• So far we have assumed ugt0 but we can generalize
  this for ult0 using the same approach as before
• To simplify this we write it as

8
cont

• By writing the time interval averaged flux
  function in this way
• We are philosophically moving away from a local
  cell reconstruction approach towards controlling
  the flux contribution from jumps in the averages
  between elements.

9
Flux Limiters

• The idea is limit the flux of q between cells
  and you will subsequently limit spurious growth
  in the cell averages near discontinuities
• A general approach is to multiply the jump in
  cell averages by a limiting function

10
cont

• The theta ratio can be thought of as a smoothness
  indicator near the cell interface at x_i-1/2.
• If the data is smooth we expect the ratio to be
  approximately 1 (except at extrema)
• Near a discontinuity we expect the ratio to be
  far away from 1.
• The flux limiting function, phi, will range
  between 0 and 2. The smaller it is, the more
  limiting is applied to a jump in cell averages.
  Above 1 it is being used to steepen the effective
  reconstruction.

11
cont

• Using this formulation we can recover the methods
  we have seen before and some new limiters

12
cont

• Using this notation we can write down the scheme
  in terms of the flux limiter function (
  )

ugt0
Limited downwindcell interface fluxcontribution
Limited upwind cell interfaceflux contribution
Upwind schemeflux contibution
ult0
13
Hartens Theorem

• Theorem Consider a general method of the form
• for one time step, where the coefficients C and
  D are arbitrary values (which in particular may
  depend on qbar in some way).
• Then provided
  that the following conditions are satisfied

14
Sweby Diagramshttp//locus.siam.org/fulltext/SINU
M/volume-21/0721062.pdf

• We need to express the flux limited scheme
• In the form
• And then satisfying the Harten conditions will
  guarantee the method is TVD.
• An appropriate choice (which we can work with) is

15
cont

• In this case since the D coefficients are zero
  and the Harten TVD conditions reduce to
• This will hold if
• We can summarize this in terms of the minmod
  function
• In addition we require
• See LeVeque p 116-118 for details

16
cont
i.e. any flux limiting function must
satisfy to be TVD. Graphically, the shaded
region is the TVD region Clearly non
of these linear limiters generate a TVD scheme.
Fromm
Beam-Warming
2
1
Lax-Wendroff
1
2
3
17
cont

• To guarantee second order accuracy and avoid
  excessive compression of solutions, Sweby
  suggested the following reduced portion of the
  TVD region as a suitable range for the flux
  limiting function

2
1
1
2
3
http//locus.siam.org/fulltext/SINUM/volume-21/072
1062.pdf
18
Minmod Flux Limiter on Sweby Diagram
2
1
1
2
3
It is apparent that the minmod flux limiter
applies the maximum possible limiting allowed
within the second order TVD region. (i.e. it will
be rather dissipative and smear out
discontinuities somewhat as seen on the right
hand side figure).
19
Superbee Flux Limiter on Sweby Diagram
2
1
1
2
3
The Superbee limiter applies the minimum limiting
and maximum steepeningpossible to remain TVD. It
is known to suffer from excessive sharpening
ofslopes as a result. On the right we show what
happens to a smooth sine wave after 20
periods. Notice the flattening of the peaks and
the steepening of the slopes.
20
MC Flux Limiter on Sweby Diagram
2
1
1
2
3
The MC limiter transitions from upwind (thetalt0)
to Fromm (at theta1/3) then switches to a
constant(at theta3). This is a compromise
between Superbee and minmod
21
van Leer Flux Limiter
The van Leer limiter charts a careful compromise
path throughthe Sweby TVD region.
22
Summary of Some Flux Limiting Functions
Linear non-TVD limiters
Nonlinear second orderTVD limiters
23
Implementation

• For ugt0 the scheme looks like
• We can easily achieve this in matlab

24
Matlab Version
This is a sample Matlab implementation of a
piecewise linear reconstructed Godunov approach
with a selection of flux limiters. Available
from the course home page http//www.caam.rice.e
du/caam452/CodeSnippets/fluxlimiter.m With the
initial condition supplied by http//www.caam.ri
ce.edu/caam452/CodeSnippets/fluxlimiterexact.m
25
Homework 4

• Q1) Using N80,160,320,640,1280 estimate the
  solution order of accuracy of the
• flux limited scheme
• with flux limiting functions
• i. Fromm
• ii. minmod
• iii. MC
• using initial conditions
• i. sin(pix)
• ii. sin(pix) (abs(x-.5)lt.25)
• on the periodic interval -1,1).
• Use the fluxlimiter.m Matlab code from the web
  page.
• You will also need to download fluxlimiterexact.m
  and minmod.m
• Measure error both using the maxmimum norm, l2
  norm and finally the
• maximum norm with data points near the
  discontinuity excluded.
• Use error plots and tables with discussion to
  describe your results.

26
Homework 4 cont

• Q2a) Invent your own 2nd order TVD flux limiter
  function (i.e. a function with range contained in
  the Sweby TVD region)
• Q2b) Modify sweby.m to plot your flux limiter
  function and compare with the limiter functions
  already used.
• Q2c) Estimate order of accuracy for a smooth
  initial condition to the advection equation
• Q2d) Estimate order of accuracy for a
  discontinuous initial condition to the advection
  equation
• Q2e) Compare results (with diagrams,results and
  comments) and discuss how your limiter differs
  from the other limiters we have seen.