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 limiter:Monotonized central-difference limiter (MC limiter)

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

5. Flux Formulation with Piecewise Linear Reconstruction • Last time we showed how the ansatz of a piecewise linear reconstruction and Godunov’s method allowed us to compute the time averaged flux contribution at each end of the I’th cell Notice: we can obtain the i-1/2 flux by setting i->i-1 in the i+1/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 formulationwith piecewise reconstruction.

7. cont • So far we have assumed u>0 but we can generalize this for u<0 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 ( ): u>0 Limited downwindcell interface fluxcontribution Limited upwind cell interfaceflux contribution Upwind schemeflux contibution u<0

13. Harten’s 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/SINUM/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. Beam-Warming Fromm 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/0721062.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 (theta<0) to Fromm (at theta=1/3) then switches to a constant(at theta=3). 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 u>0 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.edu/~caam452/CodeSnippets/fluxlimiter.m With the initial condition supplied by: http://www.caam.rice.edu/~caam452/CodeSnippets/fluxlimiterexact.m

25. Homework 4 Q1) Using N=80,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(pi*x) ii. sin(pi*x) + (abs(x-.5)<.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.