8/6/2019 000 Merita Citit Bine Scris
1/40
1Romain Teyssier 5th JETSET School
Numerical Solutions for Hyperbolic Systems of
Conservation Laws:
from Godunov Method to
Adaptive Mesh Refinement
Romain Teyssier
CEA Saclay
8/6/2019 000 Merita Citit Bine Scris
2/40
2Romain Teyssier 5th JETSET School
- Euler equations, MHD, waves, hyperbolic systems of conservationlaws, primitive form, conservative form, integral form
- Advection equation, exact solution, characteristic curve, Riemanninvariant, finite difference scheme, modified equation, Von Neumananalysis, upwind scheme, Courant condition, Second order scheme
- Finite volume scheme, Godunov method, Riemann problem,approximate Riemann solver, Second order scheme, Slope limiters,Characteristic tracing
- Multidimensional scheme, directional splitting, Godunov, Runge-Kutta, CTU, 3D slope limiting
- AMR, patch-based versus cell-based, octree structure, gradedoctree, flux correction, EMF correction, restriction and prolongation,
divB conserving interpolation- Parallel computing with the RAMSES code
Outline
8/6/2019 000 Merita Citit Bine Scris
3/40
3Romain Teyssier 5th JETSET School
Lagrangian form: fluid element of unit mass
- Mass conservation:
- Euler equation:
- First law of thermodynamics:
Chain rule:
from Lagrange derivative to Euler derivative
Volume expansion rate:
The Euler equations
8/6/2019 000 Merita Citit Bine Scris
4/40
4Romain Teyssier 5th JETSET School
Eulerian form: fixed volume element
- Mass conservation
- Momentum conservation
- Energy conservation
The Euler equations
Kinetic energy Internal energy Equation of state
8/6/2019 000 Merita Citit Bine Scris
5/40
8/6/2019 000 Merita Citit Bine Scris
6/40
6Romain Teyssier 5th JETSET School
The Euler equations
Primitive form of the Euler equations:
- Vector of primitive variables
- Quasi-linear form
8/6/2019 000 Merita Citit Bine Scris
7/40
7Romain Teyssier 5th JETSET School
The Euler equations
Hyperbolic system: J and A have real eigenvalues.
Eigenvectors are travelling waves.
Perturbation of an equilibrium state:
Wave equation:
Eigenvalues:A=
Sound speed: Eigenvectors:
8/6/2019 000 Merita Citit Bine Scris
8/40
8/6/2019 000 Merita Citit Bine Scris
9/40
9Romain Teyssier 5th JETSET School
Piecewise constant initial states:
Solution:
The Riemann problem
x
u
uL
uR
8/6/2019 000 Merita Citit Bine Scris
10/40
10Romain Teyssier 5th JETSET School
Finite difference scheme
Finite difference approximation of the advection equation
8/6/2019 000 Merita Citit Bine Scris
11/40
11Romain Teyssier 5th JETSET School
The Modified Equation
Taylor expansion in time up to second order
Taylor expansion in space up to second order
The advection equation becomes the advection-diffusion equation
Negative diffusion coefficient: the scheme is unconditionally unstable
8/6/2019 000 Merita Citit Bine Scris
12/40
12Romain Teyssier 5th JETSET School
The Upwind scheme
a>0: use only upwind values, discard downwind variables
Taylor expansion up to second order:
Upwind scheme is stable if C
8/6/2019 000 Merita Citit Bine Scris
13/40
13Romain Teyssier 5th JETSET School
Von Neumann analysis
Fourier transform the current solution:
Evaluate the amplification factor of the 2 schemes.
Fromm scheme:
Upwind scheme:
>1: the scheme is unconditionally unstable
8/6/2019 000 Merita Citit Bine Scris
14/40
14Romain Teyssier 5th JETSET School
The advection-diffusion equation
Finite difference approximation of the advection equation:
Central differencing unstable:
Upwind differencing is stable:
Smearing of initialdiscontinuity:
numerical diffusion
Thickness increases
as
8/6/2019 000 Merita Citit Bine Scris
15/40
15Romain Teyssier 5th JETSET School
Finite volume scheme
Finite volume approximation of the advection equation:
Use integral form of the conservation law:
Exact evolution of volume averaged quantities:
Time averaged flux function:
8/6/2019 000 Merita Citit Bine Scris
16/40
16Romain Teyssier 5th JETSET School
Sergei Konstantinovich Godunov
8/6/2019 000 Merita Citit Bine Scris
17/40
17Romain Teyssier 5th JETSET School
Godunov scheme for the advection equation
The time averaged flux function:
is computed using the solution of the Riemann problem defined
at cell interfaces with piecewise constant initial data.
x
u i
u i+1
For all t>0:
The Godunov scheme for the advection equation is identical tothe upwind finite difference scheme.
8/6/2019 000 Merita Citit Bine Scris
18/40
18Romain Teyssier 5th JETSET School
The time average flux function iscomputed using the self-similar solution of the inter-cell Riemannproblem:
Godunov scheme for hyperbolic systems
The system of conservation laws
is discretized using the followingintegral form:
This defines the Godunov flux:
Advection: 1 wave, Euler: 3 waves, MHD: 7 waves
Piecewise constantinitial data
8/6/2019 000 Merita Citit Bine Scris
19/40
19Romain Teyssier 5th JETSET School
Riemann solvers
Exact Riemann solution is costly: involves Raphson-Newtoniterations and complex non-linear functions.
Approximate Riemann solvers are more useful.
Two broad classes:
- Linear solvers
- HLL solvers
8/6/2019 000 Merita Citit Bine Scris
20/40
20Romain Teyssier 5th JETSET School
Linear Riemann solvers
Define a reference state as the arithmetic average or the Roe average
Evaluate the Jacobian matrix at this reference state.
Compute eignevalues and (left and right) eigenvectors
The flux function is given by the linear Riemann solution.
where
A simple example, the upwind Riemann solver:
8/6/2019 000 Merita Citit Bine Scris
21/40
21Romain Teyssier 5th JETSET School
Approximate the true Riemann fan by 2 waves and 1 intermediate state:
HLL Riemann solver
x
t
UL UR
U*
Compute U* using the integral form between S Lt and S Rt
Compute F* using the integral form between S Lt and 0.
8/6/2019 000 Merita Citit Bine Scris
22/40
22Romain Teyssier 5th JETSET School
Other HLL-type Riemann solvers
Lax-Friedrich Riemann solver:
HLLC Riemann solver: add a third wave for the contact (entropy) wave.
x
t
UL UR
U*L
S L S RS *
U*R
See Toro (1997) for details.
8/6/2019 000 Merita Citit Bine Scris
23/40
23Romain Teyssier 5th JETSET School
Higher Order Godunov schemes
Bram Van Leer
Godunov method is stable but very diffusive. It wasabandoned for two decades, until
8/6/2019 000 Merita Citit Bine Scris
24/40
24Romain Teyssier 5th JETSET School
Second Order Godunov scheme
x
u i
u i+1Piecewise linear approximation of the solution:
The linear profile introduces a length scale: theRiemann solution is not self-similar anymore:
The flux function is approximated using a predictor-corrector scheme:
The corrected Riemann solver has now predicted states as initial data:
8/6/2019 000 Merita Citit Bine Scris
25/40
8/6/2019 000 Merita Citit Bine Scris
26/40
26Romain Teyssier 5th JETSET School
Modified equation for the second order scheme
Taylor expansion in space and time up to third order:
We obtain a dispersive term and the scheme is stable for C
8/6/2019 000 Merita Citit Bine Scris
27/40
27Romain Teyssier 5th JETSET School
Summary: the MUSCL scheme for systems
Compute second order predicted states using a Taylor expansion:
Update conservative variables using corrected Godunov fluxes
8/6/2019 000 Merita Citit Bine Scris
28/40
28Romain Teyssier 5th JETSET School
Monotonicity preserving schemes
We use the central finite difference approximation for the slope:
In this case, the solution is oscillatory, and therefore non physical.
Second order linear scheme.
firstorder
secondorder
Oscillations are due to the non monotonicity of the numerical scheme.
A scheme is monotonicity preserving if:
- No new local extrema are created in the solution
- Local minimum (maximum) non decreasing (increasing) function of time.
Godunov theorem : only first order linear schemes are monotonicity preserving !
8/6/2019 000 Merita Citit Bine Scris
29/40
29Romain Teyssier 5th JETSET School
Slope limiters
Harten introduced the Total Variation of the numerical solution:
Hartens theorem : a Total Variation Diminishing (TVD) scheme ismonotonicity preserving.
Design non-linear TVD second order scheme using slope limiters:
where the slope limiter is a non-linear function satisfying:
8/6/2019 000 Merita Citit Bine Scris
30/40
30Romain Teyssier 5th JETSET School
No local extrema
We define 3 local slopes: left, right and central slopes
and
x
u i
u i+1
u i-1
New maximum !
For all slope limiters :
8/6/2019 000 Merita Citit Bine Scris
31/40
31Romain Teyssier 5th JETSET School
The minmod slope
x
u i
u i+1
u i-1
Linear reconstruction is monotone at time t n
Slope limiting is never truly second order !
8/6/2019 000 Merita Citit Bine Scris
32/40
8/6/2019 000 Merita Citit Bine Scris
33/40
33Romain Teyssier 5th JETSET School
The superbee slope
Predicted states dont overshoot initial average states.
TVD constraint is preserved by the Riemann solver.
The Courant factor now enters the slope definition.
8/6/2019 000 Merita Citit Bine Scris
34/40
34Romain Teyssier 5th JETSET School
The ultrabee slope
Use the final state to compute the slope limiter.
Upwind Total Variation constraint.
Strict Total Variation preserving limiter.
8/6/2019 000 Merita Citit Bine Scris
35/40
35Romain Teyssier 5th JETSET School
Summary: slope limiters
first order
minmod moncen
superbee ultrabee
8/6/2019 000 Merita Citit Bine Scris
36/40
36Romain Teyssier 5th JETSET School
Non linear systems: characteristics tracing.
x
u i
u i+1 Non-linear Riemann problems: waves speedsdepend on the input states.
TVD schemes are not necessary monotone.
Modify the predictor step according to the localRiemann solution: Piecewise Linear Method(PLM) and Piecewise Parabolic Method (PPM).
If else
If
else
8/6/2019 000 Merita Citit Bine Scris
37/40
37Romain Teyssier 5th JETSET School
Multidimensional Godunov schemes
2D Euler equations in integral (conservative) form
Flux functions are now time and space average.
2D Riemann problems interact along cell edges:
Even at first order, self-similarity does not apply to theflux functions anymore.
Predictor-corrector schemes ?
8/6/2019 000 Merita Citit Bine Scris
38/40
8/6/2019 000 Merita Citit Bine Scris
39/40
39Romain Teyssier 5th JETSET School
Unsplit schemes
Godunov schemeNo predictor step.
Flux functions computed using 1DRiemann problem at time t n in eachnormal direction.2 Riemann solves per step.Courant condition:
Runge-Kutta scheme
Predictor step using the Godunovscheme and t/2.Flux functions computed using 1DRiemann problem at time t n+1/2 ineach normal direction.4 Riemann solves per step.
Courant condition:
Corner Transport Upwind
Predictor step in transverse directiononly using the 1D Godunov scheme.Flux functions computed using 1DRiemann problem at time t n+1/2 ineach normal direction.4 Riemann solves per step.
Courant condition:
Second order schemes : multidimensional Taylor expansions and multidimensional slope limiters.
8/6/2019 000 Merita Citit Bine Scris
40/40
40Romain Teyssier 5th JETSET School
Beyond second order Godunov schemes ?
Smooth regions of the flowMore efficient to go to higher order.Spectral methods can show exponential convergence .More flexible approaches: use ultra-high-order shock-capturing schemes: WENO, discontinuous Galerkin anddiscontinuous element methods
Discontinuity in the flow
More efficient to refine the mesh, since higher order schemesdrop to first order.Adaptive Mesh Refinement is the most appealing approach.
What about the future ?Combine the 2 approaches.
Usually referred to as h-p adaptivity .
Top Related