Numerical Methods for Partial Differential Equations

Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 11 Instructor: Tim Warburton

Today • Flux functions for finite-volume methods • Higher order reconstruction • Limiters • Some nice notes at: • http://www.astro.uiuc.edu/classes/archive/astr496/s03_cac/graphics_0214.pdf

The Basic FV Equation • Suppose the amount of q passing through a given x position per unit time, in the positive x-direction, is given by the function f which depends on the quantity. • Consequently, the conservation law can be written as: • In the case where f or q is discontinuous (i.e. multivalued at some set of x values) we consider the evaluations on the right hand side as the limits from inside the interval [a,b] x=b x=a

cont • We define the i’th cell average of q as: • i.e. the average value of q in the interval • Assume the line is divided into uniform cells of width dx then on each cell:

cont • The x axis is subdivided into cells (also known as elements) • The cells are centered at the x with integer sub indices • The cell boundaries are at the midpoints between the cell centers.

Recall: FV Approximation Issues • The evolution equation:only gives information how to evolve the cell mean of q at each cell. • To close this equation we need to be able to compute approximate pointwise values of the flux function at the ends of each cell. • Time stepping…

Exact Time-Integration A traditional first step is to integrate the evolution equation over one time step:

Godunov’s Method The mean evolution equation without approximation: Godunov’s approach: 0) At the start of each time step we only have cell average values for q • Use the cell averages to reconstruct a polynomial approximation to q defined on each cell • Form the exact right hand side time integrals, but for the new propagation problem with the reconstructed solution as initial condition • Compute the cell averages for the new time level.

1) Reconstruct • Given the cell averages of q we will compute a local representation (qtilde) of q for each cell. • We could take the easy way out: • i.e. assume that the reconstructed polynomial is a 0th order polynomial agreeing set to the mean value at each cell. • We could instead assume a linear (in x) behavior for the reconstructed qtilde : • For some approximation to the slope of q given by:

1) Local Reconstruction cont • At each cell we will establish an approximation to the slope. • The reconstructed solution will pass through the cell average at the cell center and have the computed slope (shown as blue lines):

2) Advect Reconstructed Solution Exactly • For the simple case of the advection equation it is obvious how to compute the flux integral. • We know how to analytically evolve the reconstructed solution in time, starting at t^n • The solution shifts to the positive x direction with speed u • Thus we can use a method of characteristics argument to evaluate the flux integral (for the reconstructed solution as it propagates through the end points of each cell). • In the more general case, a fairly standard approach is to linearize the flux function about some appropriate mean state, decompose the solution into characteristic variables of the linearized system and perform the same kind of linear wave propogation analysis as we are about to describe.

2) Advect Exactly -- outflow • In the space-time plane we can plot the characteristics for the advection equation. • For the outflow – we backtrack into the i’th cell at the n’th time level

2) Time Step Restriction • If the stencil for the slopes only includes the i’th and i+1’th cells then the time step must satisfy u*dt <=dx • Otherwise we would have to use the reconstruction from the i-1’th cell. dx dt u*dt

2) Advect Exactly -- inflow • In the space-time plane we can plot the characteristics for the advection equation. • For the inflow – we backtrack into the i-1’th cell

cont • In summary – we will use values of the polynomial reconstructed q at the n’th time level to evaluate the flux integral at the inflow and outflow boundaries.

2) Advect Exactly cont • Given a piecewise linear reconstructed solution we know how to exactly evaluate each of the integrals: • We back tracked along the characteristics to find the value of the qtilde field at

3) Update Averages (p113 LeVeque) • Repeating for the inflow integral (using the neighbor cell reconstruction), we can express approximations for the right hand integrals in terms of cell averages (and the as yet unspecified slopes):

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

Accuracy • We can establish method order accuracy by Taylor series expansions. • The latter 3 methods (Fromm, Beam-Warning, Lax-Wendroff) are 2nd order accurate for sufficiently smooth solutions. • This is one order better that the straight upwind version (assuming zero slope). • We will consider the slope computed in the presence of a discontinuity in the solution.

Lax-Wendroff Reconstruction(down wind slope) • We consider a piecewise constant solution • The blue line is the reconstructed version of the approx. solution • The red dots are the original cell averages,

Exact Advection dt=dx/(2u) • We perform the required exact propagation of the reconstructed solution with dt=dx/2 • The blue line is the reconstructed version of the approx. solution advected u*dt to the right. • The red dots are the original cell averages,

Updated Averages • The new averages immediately to the left and right of the discontinuity is larger than the correct average. • These are known as overshoots. • If left unchecked these can grow at each time step.

Limiting & Total Variation • Clearly, close to the discontinuity the computed slope was too large. • To avoid these overshoots we “limit” the slope. • We measure the amount of oscillation in a solution by defining the total variation of a function: • The more the cell means vary between cells, the larger the total variation  the more oscillatory the function.

Total Variation Diminishing Method • A one step method is said to be total variation diminishing (TVD) if for any set of cell averages: • Physically, the amount of oscillation of a solution for a TVD scheme does not increase over one time-step for any starting data.

Monotone Initial Data • If the data at the n’th time step is monotone, i.e. • Then the solution remains monotone for all future time steps if the scheme is TVD. Definition: A method is called monotonicity-preserving if:Note: any TVD scheme is monotonicity-preserving.

Interpretation • If initial data is monotone then a TVD scheme cannot possibly make the solution oscillatory. • Discontinuities can however be smeared out (smoothed).

Min-mod Slope • Our goal is to choose a reduced slope near a discontinuity (here we will choose between downwind, & upwind slopes). Blue – downwind slope Green – upwind slope It is apparent that the upwind slope is a good choice in this case

Min-mod Slope cont • We introduce the min-mod slope: • Where the minmod function satisfies:

cont • Thus: • Chooses the slope with the smallest magnitude if both are the same sign. • If the slopes are of different size it sets the slope to zero (indicates local turning point for data). • This choice of slopes gives a TVD scheme.

Monotonized Central-Difference (MC) Limiter • We can also compare the upwind, downwind & central slopes. • In this case we define the slope by: • This compares the central difference slope of Fromm’s method with twice the one-sided slope to either side. • In smooth regions this reduces to the central difference formula. • The slope doubling reduces means that the one-sided slopes will only be chosen if they are rather

MC Slope Limiter • We choose between central, 2*downwind, & 2*upwind slopes. Blue – downwind slope Green – upwind slope Red – central slope

Comparisons • First – no limiter

After One Advective Period Beam-Warming (upwind slope)Big oscillations – in front of discontinuity Upwind: Very dissipative Lax-Wendroff (downwind slope): Big oscillations – post discontinuity Fromm (central slope):Oscillations either side of discontinuity

With Limiting Minmod slope limiter:Monotonized central-difference limiter (MC limiter)

With Limiting After One Period • The limiter has functioned as advertised – there are now no overshoots or undershoots either side of the discontinuity. • The minmod limiter has clearly been activated at the turning point on the smooth part of the solution – which is not desirable. • The MC limiter functions very well.

Implementation • Here I wrote one piece of code which performs any of the given methods – dictated by an input variable. • See: • fvsolver.m • fvexact.m • minmod.m

Next Time • We will discuss Harten’s sufficient conditions for a numerical method (including limiter) to be TVD • More limiters • Discussion on stability. • Extension to systems of linear PDE’s • Extension to nonlinear PDE’s

Numerical Methods for Partial Differential Equations