AMS 691 Special Topics in Applied Mathematics Lecture 9

1. AMS 691Special Topics in Applied MathematicsLecture 9 James Glimm Department of Applied Mathematics and Statistics, Stony Brook University Brookhaven National Laboratory

2. Numerical Partial Differential Equations (PDEs) Elliptic equations: usually solved iteratively, unless the dimensions of the problem are small. Usually use software packages to solve these equations. Petsi is an example. Parabolic equations: d/dt U = elliptic U Time step is very small: delta t = O(delta x)^2 if solved explicitly. Usually solved implicitly: requires solving an elliptic sub problem. Use software packages.

3. Mixed equation types Often use operator splitting: To solve (1), alternately solve (2) and (3) each for a time step. Many problems are mixed, hyperbolic, parabolic, elliptic. Solution methods are specialized for each type, so use special method and Use operator splitting to separate types.

4. Hyperbolic Equations Linear problems: use high order methods and central differences Difficulty if coefficients are not smooth Nonlinear: often conservation laws. Use conservative methods Main ideas in 1D. Use operator splitting to extend to 3D

5. Reference for numericalconservation laws @Book{LeV92, author = "R. LeVeque", title = "Numerical Methods for Conservation Laws", publisher = "Birkh{\"a}user Verlag", address = "Basel--Boston--Berlin", year = "1992", }

6. Splitting to reduce the spatial dimension Operator splitting: (1) is the same as (2a-2c) each solved is succession. Difference is not zero but is second order in delta t. More complicated ordering of the substeps changes the error to (delta t)^2 or higher.

7. Numerical Conservation Laws in 1D Two problems central to any difference scheme: Stability issues related to solution discontinuities Convergence to a weak solution of equation #1: addressed by artificial viscosity or by limiters #2: addressed by conservative methods Illustrate difficulties associated with #2: Consider two equivalent forms of Burger’s equation. They have the same smooth solutions but different jump relations and different weak solutions. Difference schemes based on (b) will not converge to a solution of (a).

8. Conservative Difference Methods ujn = cell average ujn+1 = ujn +(delta t/delta x)[Fn+1/2j+1/2 –Fn+1/2j-1/2] F must be a function of cell averages. Many ways to do this. Many different conservative numerical methods Numerical flux Fn+1/2j+1/2 = F(unj-p,…,unj+q) for stencil size p+q+1. Consistent if F(u,…,u) = F(u).

9. Gudonov methodsbased on Riemann problems and characteristic variables First order Godunov (very simple, very diffusive): Fn+1/2j+1/2 = solution of Riemann problem between unj and unj+1, evaluated at x = xj+1/2 Higher order Godunov (MUSCL and PPM): Use neighbor u values (wider stencil) to fit polynomial approximation to discrete solution values (grid cell averages). Linear for MUSCL, piecewise parabolic for PPM. Use PDE to advance solution values to ½ time level n+1/2. Get left and right values for un+1/2j+1/2. Solve Riemann problem (or approximate Riemann problem), to get numerical flux Fn+1/2j+1/2. PDE solution uses characteristic coordinates, eigenvectors, wave curves

10. Primitive variable methodsLax-Friedrichs etc. Lax-Friedrichs is a first order method. It is very diffusive. But there are a number of very good higher order methods that use only primitive variables. A primitive variable is one of the (conserved) components of the solution U. Lax-Friedrichs is very simple:

11. Lax Wendroff Works well for smooth solutions. Second order accurate. Needs artificial viscosity to cure post shock oscillations if used for discontinuous flows. Not competitive with modern methods for such flows.

12. Nonlinear Stability Severe oscillations in the numerical solution following a shock wave. Similar to Gibbs phenomena, where by convergence of a Fourier series for discontinuous data is highly oscillatory near the discontinuity. Cure (a): add “artificial viscosity”, to dissipate the oscillations. (b) Slope limiters in the reconstruction step of higher order Godunov solvers. Consider MUSCL, which is a piecewise linear reconstruction, to define the solution U at every point, in terms of the cell averages over a stencil of grid cells adjacent to the cell in question. The linear reconstruction at the face between two adjacent cells will produce a jump discontinuity. We want this discontinuity to be interpolative not extrapolative, and we reduce the slope (slope limiter) in order to achieve this goal. In other words, the endpoints of the values for the solution within the cell must lie within the values given from the neighbor cells as extrapolated to their cell endpoints.

13. Refinements on slope limiters:ENO/WENO If the basic scheme uses 3 or 5 or j points in its stencil the basic idea of Eno/Weno is to slide the stencil to the left or right to choose an optimal stencil, one which has a smoother solution and a reduced requirement for limiting. In other words it is to limit the limiters! Eno considers a number of these translated stencils and chooses the best, ie the smoothest one with the least limiting. Weno considers all possible translated stencils and performs a weighted average of all, with the weights largest for the smoothest and smaller for the ones that need more limiting. Eno/Weno is a stencil method and can be used in combination with any solver, and either using primitive or characteristic methods

14. Comments on methods Many methods, others not yet mentioned (discontinuous Galerkin, compact schemes, hybrid combinations of above). Which to use? Problem dependent, some methods are faster, some more accurate, but accuracy is problem dependent. Usually any good method will do. Other issues: AMR = automatic mesh refinement Implicit methods: stable for large time steps. Useful if system has fast but not important waves. Many methods, others not yet mentioned (discontinuous Galerkin), compact schemes, hybrid combinations of

15. Other issues: entropy It is desirable to define a discrete entropy and to determine that the entropy inequality is maintained in time. This prevents rarefaction shocks from occuring.

16. Incompressible Navier-Stokes Equation (3D)

17. The Projection Method The role of the pressure is to ensure that the velocity is divergence free. The pressure is not a dynamic variable. It is derived from the divergence free condition on the velocity field. The essential steps in the solution of incompressible NS are elliptic steps: to solve the pressure equation and the diffusion equation.