Grid adaptivity for systems of conservation laws M. Semplice 1 G. - - PowerPoint PPT Presentation

Grid adaptivity for systems of conservation laws

• M. Semplice1
• G. Puppo2

1Dipartimento di Matematica

Universit` a degli Studi di Torino

2Dipartimento di Scienze Matematiche

Politecnico di Torino

HYP2012

14th International Conference on Hyperbolic Problems: Theory, Numerics and Applications

Padova, 25-29 June 2012

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 1 / 34

A recipe for an adaptive scheme

Ingredients

(finite volume) grid numerical fluxes time marching method refine/coarsen strategy (and indicator)

Recipe

Mix all ingredients, stir until it compiles and runs smoothly, enjoy the movie show!

Tools

Object Oriented Programming and a library like DUNEa allows to code each part independently (e.g. the Runge-Kutta does not know the number

• f spatial dimensions)

aDune (Distributed and Unified Numerics Environment),

http://www.dune-project.org/index.html

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 2 / 34

1

Locally refined grids

2

Timestepping

3

Adaptivity

4

Local recomputation

5

Two space dimensions

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 3 / 34

Locally refined grids

2D grids and quad-trees

Data are attached only to coloured cells (leaves) The grid manager (dune-grid for us) allows to iterate easily

• n all the leaves
• n a given level

Note: similar structures holds cubes, triangles, tetrahedra, . . .

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 4 / 34

Timestepping

1

Locally refined grids

2

Timestepping

3

Adaptivity

4

Local recomputation

5

Two space dimensions

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 5 / 34

Timestepping

Semidiscrete scheme

∂tu + ∂xf (u) = 0 R-K: A b stages u(i)

j

= un

j − λ i−1

• k=1

aik

• F (k)

j+1/2 − F (k) j−1/2

• λ = ∆t

h fluxes F (i)

j+1/2 = F

• u(i),−

j+1/2, u(i),+ j+1/2

• step un+1

j

= un

j − λ ν

• i=0

bi

• F (i)

j+1/2 − F (i) j−1/2

• where u(i),±

j+1/2 are some (non-oscillatory) reconstructions.

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 6 / 34

Timestepping

Imposing the CFL globally

CFL = ⇒ ∀k : ∆t ≤ λhk = λh02−k = ⇒ ∆t = λh02−kmax t tn :

• i

biF (i) : cell averages tn+1 tn+2 tn+3 tn+4

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 7 / 34

Timestepping

Local timestepping

Choosing ∆t globally: easier to program wastes CPU resources for big cells, where the solution is well approximated (in theory) unneeded numerical diffusion if many cells have sub-optimal (local) mesh ratio

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 8 / 34

Timestepping

Local timestepping

Choosing ∆t globally: easier to program wastes CPU resources for big cells, where the solution is well approximated (in theory) unneeded numerical diffusion if many cells have sub-optimal (local) mesh ratio

Local timestepping

Each cell advances with different stepsizes ∆tk = λhk = λh02−k x t ∆t1 ∆t2 ∆t3 synchro

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 8 / 34

Timestepping

Review of local timestepping

Motivated by approximating (on uniform grids) conservation laws with strong spatial variations of the velocity: First order Osher-Sanders (1983) Second order, TVD property Dawson-Kirby (2001), Kirby (2003) multiresolution (wavelets) Muller-Stiriba(2007) More recent work has focused on AMR: Second order Runge-Kutta Constantinescu-Sandu (2007), Krivodonova (2010)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 9 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries everywhere

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries small

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries un

2 − ∆t h2 ( 1 4F3/2 − 1 4F1/2)

small

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries un

2 − ∆t h2

• 1

4F 1 3/2 + 1 4F 2 3/2

− 1

2F 3 1/2

• collect all

fluxes small, medium

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

First order scheme on nonuniform grid – local timestepping

∆tk = λhk ∆t = ∆tmacro x t tn tn+1 un

1

un

2

un

3 un 4

un

5

un

6

cell averages fluxes temporaries

un+1

1

un+1

2

un+1

3

un+1

4

un+1

5

un+1

6

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 10 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Second order scheme on nonuniform grid – global CFL

cell averages and stage values fluxes F (1)

j+1/2 and

F (2)

j+1/2

temporaries x t tn tn+1 automatically conservative! large reconstruction stencil ⇒ many temporaries

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 11 / 34

Timestepping

Let’s test the performance of nonuniform grids

1 build a nonuniform grid

ratio 1:2 ratio 1:4 ratio 1:4 ratio 1:4 ratio 1:4

2 leave it fixed in time 3 make different kinds of waves travel along them, crossing grid

discontinuities.

We expect:

poor behaviour on error (the grid is fixed!) check if grid discontinuities can deform waves

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 12 / 34

Timestepping

Smooth solution (linear transport)

u0(x) = [cos(2πx)]4

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 13 / 34

Timestepping

Shock (Burgers equation)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 14 / 34

Timestepping

Shock (Burgers equation)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 14 / 34

Timestepping

Errors on uniform/nonuniform grids

grid ratio 1:2

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 15 / 34

Timestepping

Smooth solution (Euler equations)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 16 / 34

Timestepping

Smooth solution (Euler equations)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 16 / 34

Adaptivity

1

Locally refined grids

2

Timestepping

3

Adaptivity

4

Local recomputation

5

Two space dimensions

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 17 / 34

Adaptivity

Semidiscrete scheme with numerical entropy production

∂tu + ∂xf (u) = 0 ∂tη(u) + ∂xψ(u) ≤ 0 R-K: A b stages u(i)

j

= un

j − λ i−1

• k=1

aik

• F (k)

j+1/2 − F (k) j−1/2

• λ = ∆t

h fluxes F (i)

j+1/2 = F

• u(i),−

j+1/2, u(i),+ j+1/2

• Ψ(i)

j+1/2 = Ψ

• u(i),−

j+1/2, u(i),+ j+1/2

• step un+1

j

= un

j − λ ν

• i=0

bi

• F (i)

j+1/2 − F (i) j−1/2

• entropy Sn

j = 1

∆t

• η
• un+1

j

• − η
• un

j

• + λ

ν

• i=1

bi

• Ψ(i)

j+1/2 − Ψ(i) j−1/2

• Note: most of the seminar is independent from this choice of error

indicator, although the numerical entropy production works well in general

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 18 / 34

Adaptivity

Numerical entropy production in the Sod test

With the second order scheme (Heun + KT flux + minmod), using periodic b.c. O(1) on contacts ⇒ threshold

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 19 / 34

Adaptivity

Choose your favourite indicator and loop until final time:

compute un+1

j

and indicator Sn

j

∃j :

• Sn

j > Sref

hj > hmin

?

coarsen to father if all sons have Sn

j < Scoa

locally refine (locally) recompute NO YES (next timestep)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 20 / 34

Adaptivity

How to set coarsen thresholds?

Coarsening will succede (in Rd) if 2d neighbouring cells have Sn

j < Scoa.

We expect to coarsen on smooth regions of the flow, so assume that at the next timestep the father cell might produce numerical entropy of size Sn+1

father ∼ 2dScoa

We want to avoid to immediately refine the cell again, so we need 2dScoa < ηSref for some η < 1 Practical guideline: η = 1

4

⇒ Scoa = Sref 2d+2

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 21 / 34

Adaptivity

Lax test (1)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 22 / 34

Adaptivity

Lax test (2)

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 23 / 34

Adaptivity

Shock-acoustic test problem (1)

Using minmod limiter

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 24 / 34

Adaptivity

Shock-acoustic test problem (2)

MM minmod limiter CT limiter from Cada, Torrilhon – J. Comput. Phys. (2009) ENO2 second order ENO AdLim minmod limiter applied only where error indicator is high

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 25 / 34

Local recomputation

1

Locally refined grids

2

Timestepping

3

Adaptivity

4

Local recomputation

5

Two space dimensions

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 26 / 34

Local recomputation

Need for local recomputation

The refine/recompute loop is highly inefficient if we recompute the solution over the whole grid! Using the numerical domain of dependence, it is possible to predict which cells need recomputing after a given cell is refined.

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 27 / 34

Local recomputation

Recomputation after refinement – uniform ∆t

• riginal

recomputed u(1) un+1 Example: 2-stage RK + linear reconstructions x t Note: the domain of dependence is enlarged by number of RK stages and reconstruction stencil

Local recomputation

Recomputation after refinement – uniform ∆t

• riginal

recomputed u(1) un+1 Example: 2-stage RK + linear reconstructions x t Note: the domain of dependence is enlarged by number of RK stages and reconstruction stencil

Local recomputation

Recomputation after refinement – uniform ∆t

• riginal

recomputed u(1) un+1 Example: 2-stage RK + linear reconstructions x t Note: the domain of dependence is enlarged by number of RK stages and reconstruction stencil

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 28 / 34

Local recomputation

Recomputation after refinement – local timestepping

Following the domain of dependence is hard, computationally expensive (and unneeded since we work with an adaptive framework) The error indicator was happy about the solution away from the first cell, so only recompute red quantitities in the diagram: x t tn tn+1 exact d. of dep. We can accumulate fluxes and don’t have to store them for each substep!

Local recomputation

Recomputation after refinement – local timestepping

Following the domain of dependence is hard, computationally expensive (and unneeded since we work with an adaptive framework) The error indicator was happy about the solution away from the first cell, so only recompute red quantitities in the diagram: x t tn tn+1 exact d. of dep. fixed buffer size conservative! We can accumulate fluxes and don’t have to store them for each substep!

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 29 / 34

Two space dimensions

1

Locally refined grids

2

Timestepping

3

Adaptivity

4

Local recomputation

5

Two space dimensions

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 30 / 34

Two space dimensions

Doswell “frontogenesis” problem

Scalar conservation law, mimicking the mixing of a cold and hot front: ut + v(x, y)∇ · u = 0

• n Ω = [−4, 4]2

with velocity v(x, y) = f (r)

.385

• − y

r , x r

• with r =
• x2 + y2 and f = tanh(r)/(cosh(r))2

and initial data u0(x, y) = − tanh(y/2).

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 31 / 34

Two space dimensions

Atmospheric instability: the solution

solution cell level

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 32 / 34

Two space dimensions

Euler equations: 2D Riemann problem

solution and cell size entropy indictor (log scale) initial

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 33 / 34

Two space dimensions

Euler equations: 2D Riemann problem

solution and cell size entropy indictor (log scale) middle

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 33 / 34

Two space dimensions

Euler equations: 2D Riemann problem

solution and cell size entropy indictor (log scale) final

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 33 / 34

Two space dimensions

Perspectives

third order CWENO reconstructions on nonuniform grids (with G. Russo, A. Coco) third order in time (done for global timestepping), coming for local timestepping . . . source terms and well balancing (with S. Noelle, . . . ) All tests were run using a finite volume library that I developed on top of the DUNE environment. The source code will be released under the GPL2 licence (by the end of 2012). Watch out on http://www.unito.it/dipmatematica

• M. Semplice et al. (Univ. Insubria)

Grid adaptivity HYP2012 34 / 34