Limiters and Riemann solvers – the basic ingredients of Finite Volume schemes

of Finite Volume schemes

Friedemann Kemm BTU Cottbus

kemm@math.tu-cottbus.de

-1 -0.5 0 0.5 1 1.5 2 2.5 3

-4 -2 0 2 4 6

TVD-region Sweby-region linear 3rd order Lax-Wendroff Beam-Warming Fromm

expansion

shock contact

We start with 1d linear advection, end with multidimensional nonlinear systems.

System

Multi Dimensional

Nonlinear

Wave equation

Full Gas Dynamics

Burgers Linear Advection

2d Advection

Shock tube

The simplest discretization of linear advection is unconditionally unstable.

q_iⁿ⁺¹ − q_iⁿ

∆t − aq_iⁿ₊₁ − aq_iⁿ₋₁

2∆x = 0

000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000

111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111

0000000000 0000000000 0000000000 0000000000 0000000000 0000000000

1111111111 1111111111 1111111111 1111111111 1111111111 1111111111

Eigenvalues Stability region

1

−i

−1

i

Stable schemes are obtained with

better time- or space discretizations.

Change time difference:

Leapfrog q_iⁿ⁺¹ − q_iⁿ⁻¹

2∆t

Lax-Friedrichs q_iⁿ⁺¹ − ¹₂(q_iⁿ₊₁ + q_iⁿ₋₁)

∆t

Rusanov/LLF q_iⁿ⁺¹ − ¹₂(νq_iⁿ₊₁ + 2(1 − ν)q_iⁿ + νq_iⁿ₋₁)

∆t Change space difference:

Upwind aq_iⁿ − aq_iⁿ₋₁

∆x

The different forms of the upwind scheme lead to different generalizations.

Difference form

q_iⁿ⁺¹ = q_iⁿ − ν∆q_iⁿ_−1/2

Fluctuation form

q_iⁿ⁺¹ = q_iⁿ + _∆x^∆t a (−∆w_iⁿ_−1/2) , w_iⁿ_−1/2 = ∆q_iⁿ_−1/2

Viscous form/Flux form

q_iⁿ⁺¹ = q_iⁿ + _2∆x^∆t [f (q_iⁿ₊₁) − f (q_iⁿ₋₁)] + _2∆x^∆t a [q_iⁿ₊₁ − 2q_iⁿ + q_iⁿ₋₁]

= q_iⁿ − _∆x^∆t [F_iⁿ_+1/2 − F_iⁿ_−1/2]

with F_i+1/2 = f (q_i) = ¹₂ [f (q_iⁿ₊₁) + f (q_iⁿ)] − ¹₂ |a|∆q_i+1/2

Lax-Wendroff can be written as a correction of first order upwind.

Lax Wendroff

advection

q_iⁿ⁺¹ = q_iⁿ −

ν + 1

2ν(1 − ν)

1

r_iⁿ − 1

∆q_iⁿ_−1/2

with

r = ∆q_up

∆q_{dow n} ν = distance covered in ∆t

∆x

Beam-Warming can be written as a correction of first order upwind.

Beam Warming

advection

q_iⁿ⁺¹ = q_iⁿ −

ν + 1

2ν(1 − ν)

r_iⁿ

r_iⁿ − r_iⁿ₋₁

∆q_iⁿ_−1/2

with

r = ∆q_up

∆q_{dow n} ν = distance covered in ∆t

∆x

Fromm scheme can be written as a correction of first order upwind.

Fromm

advection

q_iⁿ⁺¹ = q_iⁿ −

ν + 1

2ν(1 − ν)

(1 + r_iⁿ)/2

r_iⁿ − (1 + r_iⁿ₋₁)/2

∆q_iⁿ_−1/2

with

r = ∆q_up

∆q_{dow n} ν = distance covered in ∆t

∆x

TVD-schemes can be written as a correction of first order upwind.

General Case

advection

q_iⁿ⁺¹ = q_iⁿ −

ν + 1

2ν(1 − ν)

ϕ(r_iⁿ)

r_iⁿ − ϕ(r_iⁿ₋₁)

∆q_iⁿ_−1/2

with

r = ∆q_up

∆q_{dow n} ν = distance covered in ∆t

∆x

A mean function should be

at least consistent and homogeneous.

ϕ(r) = M(r,1)

Consistency M(a, a) = a ϕ(1) = 1

Inclusion min{a, b} ≤ M(a, b) ≤ max{a, b} min{r,1} ≤ ϕ(r) ≤ max{r,1}

Homogeneity M(λa,λb) = λM(a, b)

Symmetry M(a, b) = M(b, a) ϕ(1/r) = ϕ(r)/r

Monotonicity M(·, b) , M(a,·) increasing ϕ(r) , r ϕ(1/r) increasing

The TVD-region is much larger than the Sweby region.

-1 -0.5 0 0.5 1 1.5 2 2.5 3

-6 -4 -2 0 2 4 6

CFL=0.1 TVD

Sweby

-1 0 1 2 3 4 5

-15 -10 -5 0 5 10 15

CFL=0.5 TVD

Sweby

0 5 10 15 20

-20 -10 0 10 20 30

CFL=0.9 TVD

Sweby

− 2

1 − |ν|≤ ϕ(r)

r − ϕ(R) ≤ 2

|ν|

−2≤ ϕ(r)

r − ϕ(R) ≤ 2

Third order schemes are upwind biased.

-1 -0.5 0 0.5 1 1.5 2 2.5 3

-6 -4 -2 0 2 4 6

CFL=0.1 3rd order

LW BW Fromm

-1 0 1 2 3 4 5

-15 -10 -5 0 5 10 15

CFL=0.5 3rd order

LW BW Fromm

0 5 10 15 20

-20 -10 0 10 20 30

CFL=0.9 3rd order

LW BW Fromm

ϕ₃(r) =

1 − 1 + |ν| 3

+ 1 + |ν|

3 r ϕ_LW(r) = 1 ϕ_BW(r) = r

Limiters might be constructed by sticking to third order as long as possible.

-1 -0.5 0 0.5 1 1.5 2 2.5 3

-6 -4 -2 0 2 4 6 CFL=0.1

θ=1.00 θ=0.75 θ=0.50 MC

-1 0 1 2 3 4 5

-15 -10 -5 0 5 10 15

CFL=0.5 θ=1.00

θ=0.75 θ=0.50 MC

0 5 10 15 20

-20 -10 0 10 20 30 CFL=0.9

θ=1.00 θ=0.75 θ=0.50 MC

ϕ_θ(r) = minn

maxn

−(1 − θ) 2

|ν|, ϕ₃(r)o ,

maxn

−(1 − θ) 2

1 − |ν|r, θ 2

|ν|ro

, θ 2 1 − |ν|

o

Superbee-type limiters provide a

good approximation of the amplitude.

0 0.2 0.4 0.6 0.8 1 1.2

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 CFL=0.1

Ultrabee β=2/3 Superbee

0 0.2 0.4 0.6 0.8 1 1.2

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 CFL=0.1

θ=0.50 θ=1.00 MC

Standard example with 200 cells after 10 full rounds (t = 20)

The squaring effect spoils the

convergence of Superbee type limiters.

0.01 0.1 1

100 1000 10000 100000

CFL=0.1

Superbee Ultrabee β=2/3 Superpower θ=1

0.01 0.1 1

100 1000 10000 100000

CFL=0.9

Superbee Ultrabee β=2/3 Superpower θ=1

L¹-norm of the error

There are more methods to obtain higher order.

PPM

Woodward Colella 1984

generalized to higher orders by Rider, Greenough and Kamm 2005

PHM, PRM, Double logarithmic limiter

Marquina, Xiao and Peng 2004, Artebrandt and Schroll, ˇCada and Torrilhon

ENO/WENO

Harten, Shu, Osher,. . .

measure oscillation of different stencils

2d transport can be realized by DCU or CTU.

v

CTU

DCU

One dimensional schemes can be used in multi-d by dimensional splitting.

Godunov

first order

Strang

preserves second order

Lin, Rood (Monthly Weather Review, 1996)

in addition resembles CTU

Geometric limiting can be done in all geometries.

Structured

direction wise multidimensional

Unstructured grids

Meshless grids (Sonar et al.)

CFL-independent limiters may be generalized for unstructured grids.

Minmod

Many variants

MC

Barth, Jespersen (1989)

Albada

Tu, Alliabadi (2005)

Multi-d limiters should be

at least consistent and homogeneous.

Swartz, 1999:

good stencil

bad stencil

Scalar methods can be applied to systems by characteristic decompositon.

q_t + Aq_x = 0

A = RΛL w = Lq

w_t + Λw_x = 0

The numerical flux for linear advection has many nonlinear counterparts.

Godunov: Numerical flux from exact solution of RP Flux Splitting:

F_i+1/2 = f ⁺(q_i) + f ⁻(q_i+1) Enquist Osher:

f ⁺(q) =

Z q

0

max{f ⁰(z),0} d z

f ⁻(q) =

Z q

0

min{f ⁰(z),0} d z

F_i+1/2 = 1

2(f (q_i+1) + f (q_i)) − 1 2

Z q_i+1

q_i

|f ⁰(z)| d z Linearization (Roe)

For nonlinear equations, numerical fluxes can be computed by Roe-linearization.

f (q_i+1) − f (q_i) =

Z 1 0

f ⁰(θq_i+1 + (1 − θ)q_i) d θ

| {z }

=:˜a(q_i+1,q_i)

(q_i+1 − q_i)

More general:

q = q(w) , b = d q

d w f (q) = F(w) , c = d F d w

∆q_i+1/2 = ˜b(q_i+1, q_i)∆w_i+1/2 ∆F_i+1/2 = c˜(q_i+1, q_i)∆w_i+1/2 leads to

˜

a(q_i+1, q_i) = c˜(q_i+1, q_i)

˜b(q_i+1, q_i)

The scalar concepts can be directly generalized to nonlinear systems.

Godunov Numerical flux from exact solution of RP

Osher-Solomon

F_i+1/2 = 1

2(f(q_i+1)+f(q_i))−1 2

Z

Γ

|f⁰(z)|dz

other FVS Steger Warming, Vijayasundaram, van Leer, AUSM, Eberle etc.

Roe Linearization parameter vectors

HLL-type solvers are Roe-type solvers and vice versa.

V_Roe = 1

2 |A(q˜ _r,q_l)|

V_HLL = 1 2

S_R + S_L S_R − S_L

A(q˜ _r, q_l) − S_RS_L S_R − S_LI

Rusanov/LLF

S_R = −S_L = |λ_max| fastest wave speed in RP

Lax-Friedrichs

S_R = −S_L = ∆x

∆t maximal allowed wave speed on grid

Smooth limiters are a

good choice for nonlinear waves.

3 3.5 4 4.5

0 0.5 1 1.5 2 2.5 3

reference solution mixed Superpower and Ultrabee

3 3.5 4 4.5

0 0.5 1 1.5 2 2.5 3

reference solution mixed β=2/3 and Ultrabee

Detail of Shu-Osher problem, t = 1.8, 400 cells

Characteristic CFL-dependent limiting enhances the quality of TVD-schemes.

0 1 2 3 4 5 6 7

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35

Superbee, primitive Minmod, characteristic Superbee, characteristic Ultrabee, characteristic

Detail of Toro’s test case 3

Godunov scheme is

not always the best choice.

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60

5 0 10 15 20 30 25

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60

5 0 10 15 20 30 25

0 10 20 30 40 50 60 70

Godunov HLLE

Standard carbuncle fix is based on nonlocal data.

flux to compute

strong shock?

strong shock?

The carbuncle can be addressed from within the Riemann solver.

By entropy consistency Roe 2008

Bouchut 2003/2004

By adjusting viscosity on shear waves Kemm 2008

Nishikawa and Kitamura 2008

0 10 20 30 40 50 60 70

0 10 20 30 40 50 60

0 5 10 15 20 25 30

0 10 20 30 40 50 60 70

HLLEMCC

JST-schemes choose viscosity terms by indicator functions.

differentiable

mainly for steady states and implicit time discretizations

stabilized by implicit time schemes

naive use yields bad results

viscosity via second and fourth derivatives

New schemes provide

multidimensional flux calculations.

FVEG

Morton, Warnecke, Lukaˇcova,. . .

2d HLL/HLLC/LFC

Wendroff 1999

Nonstaggered central schemes

Tadmor et al.

A careful choice of limiters and

Riemann solvers yields fine results.

3 3.5 4 4.5

0 0.5 1 1.5 2 2.5 3

CFL, characteristic, Roe

3 3.5 4 4.5

0 0.5 1 1.5 2 2.5 3

Minmod, primitive, HLLE

Detail of Shu-Osher problem, t = 1.8, 400 cells