A Brief Introduction to Godunov Methods

1. A Brief Introduction to Godunov Methods Tim Handy TexPoint fonts used in EMF. Read the TexPoint manual before you delete this box.: A

2. Outline • Basic Equations of Fluid Mechanics • Euler’s Equations • Divergence Refresher • Strong Forms • Weak Forms • Numerical Solution of Euler’s Equations • Terms • Discretization • Godunov’s Scheme

3. Euler’s Equations “Although to penetrate into the intimate mysteries of nature and thence to learn the true causes of phenomena is not allowed to us, nevertheless it can happen that a certain fictive hypothesis may suffice for explaining many phenomena.” – Leonhard Euler

4. Divergence Refresher • Divergence of a vector field at a point is the measure of how much the point acts like a source or sink • Defined as the limit of volume weighted flux as the volume goes to a point

5. Source: divF > 0

6. Sink: divF < 0

7. Incompresssible/Solenoidal: divF = 0

8. Euler’s Equations – Strong Form • General Equation • Statement of quantity conservation • At any point in space, the variation of the quantity in time is proportional to the amount of quantity flowing towards or away from the point (source/sink)

9. Euler’s Equations – Strong Form • Continuity Equation • Statement of mass conservation • At any point in space, the variation of the density in time is proportional to the amount of material flowing towards or away from the point (source/sink)

10. Euler’s Equations – Strong Form • Momentum Equation • Statement of momentum conservation • At any point in space, the variation of the momentum in time is proportional to the combination of advective momentum transfer, pressure gradient, and other body forces (gravity, electromagnetics [jxB], coordinate accelerations, etc.)

11. Euler’s Equations – Strong Form • Energy Equation • Statement of energy conservation For ideal gas situations, Not true for time dependent gravity situations

12. Euler Equations – Strong Form • In general, the Euler Equations are a system of nonlinear, hyperbolic PDE’s • General solutions are waves • Multi-valued solutions possible (Shocks; obey Rankine-Hugoniot; “weak” solutions) • Solutions travel along characteristics and come in pairs • If steady state • M<1: Elliptic • M=1: Parabolic

13. Euler Equations – Weak Form • General Formulation (Reynolds Transport Theorem) • (Rate of change of N) = (Accumulation in CV) + (Flux through CS due to CV moving) + (Flux through CS due to velocity field)

14. Euler Equations – Weak Form • Continuity (Stationary CV) • (Rate of change of mass) = (RoC of how much mass we have in CV) + (How much mass is crossing the boundary)

15. Euler Equations – Weak Form • Momentum (Stationary CV) • (Rate of change of momentum) = (RoC of how much momemtum we have in CV) + (How much momentum is crossing the boundary) = (Sum of forces on CV)

16. Euler Equations – Weak Form • Energy (Stationary CV) • (Rate of change of total energy) = (RoC of how much energy we have in CV) + (How much energy is crossing the boundary) = (Rate of heat transfer into CV) + (Rate of work done by system)

17. Equation Summary • Strong form: Finite Difference Methods, Spectral Methods • Weak form: Finite Volume Methods, Finite Element Methods, Boundary Element Methods

18. Numerical Solution of Euler’s Equations The Hard Part

19. Total Variation Diminishing • Total Variation: The total length in the y-direction that the ball moves traveling along the path • Total Variation Diminishing: At next timestep, the length the ball travels either stays the same or decreases • Wiggles do not grow, they are either conserved or decreased • Implies stability

20. Total Variation Diminishing • In 1D • A scheme is TVD if

21. Monotonicity • A scheme is monotonicity preserving if: • It does not create new local extrema within the spatial domain • The value of a local minima is non-decreasing and the value of a local maxima is non-increasing • Harten (1983) proved that: • A monotone scheme is TVD • A TVD scheme preserves monotonicity • Why do we want monotonicity? Fixes overshoots in approximation. For example, without it you may be constantly adding mass to your solution due to overshooting the density value in the continuity equation. • Monotonic schemes do not provide non-physical solutions

22. Godunov’s Theorem Monotonic, linear numerical schemes for solving partial differential equations can be at most first order accurate. • If we want a monotonic scheme, it can only be first-order accurate (piecewise constant) • Error very slowly decreases with mesh resolution

23. Consistency • Consistency: Truncation error goes to 0 as h→0 Truncation Error http://ct.gsfc.nasa.gov/amr.html

24. Convergent • A sequence {xn} is convergent to L if there exists an N such that |xn-L|<ε for all n>N • A simulation Sn is convergent to the true solution S if there exists a mesh spacing based on n such that ||Sn-S||<ε • Main point: Mesh spacing gets smaller, simulation becomes like true solution

25. Lax Equivalence Theorem Therefore, if we have a consistent method (truncation error goes to 0), and it is TVD (no extra wiggles), it should converge to the true solution as the mesh becomes increasingly more fine A consistent scheme is convergent if and only if it is stable.

26. Godunov’s Scheme • Reconstruct profiles for each variable • Constant, linear, parabolic, … • Force monotonicity by adjusting profile coefficients • Solve the local Riemann problem at each cell interface • Integrate in time and determine cell averages • Repeat

27. Riemann Problem • The Riemann problem describes the interaction of hyperbolic systems at an interface with a jump discontinuity PL uL ρL PR uR ρL What is the flux of our conserved quantities here?

28. Riemann Problem • Need to determine the average u/p state at the interface to determine fluxes • FLASH implements the two shock Riemann solver described by Colella and Glaz (1985). • The solver operates by decomposing A(U) into left, right, and center eigenvectors • These eigenvectors specify a set of characteristic equations that determine potential interface values • One is then able to integrate along these characteristics to determine the intersection point in the u-p plane

29. Discretization • Assume area is known, fluxes are constant w.r.t. CS (I2) using the average u/p at the interface, and we can determine the average amount of stuff in the CV (I1)

30. Determining New Cell Averages • Using the above information coupled with the reconstruction of density, we can determine the average fluxes • Now that we have the average fluxes, we can determine how much “stuff” should cross the boundary in ∆t time • Application of the discretization of a conservation law gives us the new cell average!

31. Shock Treatment • In order to accurately resolve shocks, the interpolated variable profiles need to be modified • In the vicinity of the shock front, all profiles are forced to be lower order fits (Remember Godunov’s Theorem!)

32. Courant-Friedrichs-Lewy (CFL) Condition • Information theory limits the admissible time step • If you move too far ahead in time, information (advective or acoustic) from one fluid packet may completely jump a cell • Therefore, must enforce that fluid packet signals may only move into the adjacent cell each time step