Parametrized maximum principle preserving flux limiters for high - - PowerPoint PPT Presentation

Parametrized maximum principle preserving flux limiters for high order schemes: Algorithm, Development, Application

Zhengfu Xu Michigan Technological University Collaborators: J.-M. Qiu, T. Xiong (U. Houston), Y. Jiang (ISU) , C. Liang (MTU), Y. Liu, Q. Tang, A. Christlieb (MSU) May 23, 2014

Outline

◮ Discrete maximum principle preserving (MPP). ◮ A sufficient condition for conservative schemes to be MPP. ◮ General MPP flux limiters ◮ Partial analysis results regarding CFL and accuracy. ◮ Applications and results. ◮ Conclusion with ongoing and future work.

Four major concerns for solving HCL

ut + f (u)x = 0 numerically:

• 1. Accuracy: High order reconstruction.
• 2. Stability: MPP provides a weak stability.
• 3. Conservation: This is the basic framework for discussion.
• 4. Efficiency: CFL number shall not be too small.

High order finite volume and finite difference method

High order finite volume formulation: integrate the PDE over Ij = [xj−1/2, xj+1/2] d¯ uj dt = −fj+1/2 − fj−1/2 ∆x (1) evolving ¯ uj =

1 ∆x

• Ij udx.
• 1. ¯

uj → u±

j+1/2 by polynomial reconstruction. For example, find

p(x) of order 2 such that ¯ uk = 1 ∆x

• Ik

p(x)dx, k = j − 1, j, j + 1. Then let uj+1/2 ≈ p(xj+1/2).

• 2. u±

j+1/2 → fj+1/2 by applying suitable approximate flux

function.

High order finite difference method

• 1. Finite difference formulation with a sliding function

1 ∆x x+ ∆x

2

x− ∆x

2

h(ξ)dξ = f (u(x, t)), f (u)x = 1 ∆x

• h(x + ∆x

2 ) − h(x − ∆x 2 )

• ,

duj dt = − 1 ∆x

• hj+/2 − hj−1/2
• .
• 2. Treat f (u) as the average volume of some function h, then

back to the finite volume reconstruction: f (uj) → h±

j+1/2.

• 3. Reconstruction: ENO (Essentially-Non-Oscillatory); WENO

(Weighted ENO); or others.

MPP high order scheme for hyperbolic equation

◮ One dimensional scalar hyperbolic problem

ut + f (u)x = 0, u(x, 0) = u0(x), (2) with boundary condition. MPP means: um ≤ u(x, t) ≤ uM if um = min

x u0(x) and uM = max x u0(x)

MPP high order scheme for hyperbolic equation

◮ One dimensional scalar hyperbolic problem

ut + f (u)x = 0, u(x, 0) = u0(x), (2) with boundary condition. MPP means: um ≤ u(x, t) ≤ uM if um = min

x u0(x) and uM = max x u0(x) ◮ A typical conservative scheme with Euler forward

un+1

j

= un

j − λ(ˆ

Hj+1/2 − ˆ Hj−1/2). (3)

MPP high order scheme for hyperbolic equation

◮ One dimensional scalar hyperbolic problem

ut + f (u)x = 0, u(x, 0) = u0(x), (2) with boundary condition. MPP means: um ≤ u(x, t) ≤ uM if um = min

x u0(x) and uM = max x u0(x) ◮ A typical conservative scheme with Euler forward

un+1

j

= un

j − λ(ˆ

Hj+1/2 − ˆ Hj−1/2). (3)

◮ The goal: To achieve MPP numerically in high order,

conservative manner: um ≤ un

j ≤ uM

MPP finite volume scheme by Zhang & Shu

◮ Numerical formulation:

¯ un+1

j

= ¯ un

j − λ{F(u− j+1/2, u+ j+1/2) − F(u− j−1/2, u+ j−1/2)}.

MPP finite volume scheme by Zhang & Shu

◮ Numerical formulation:

¯ un+1

j

= ¯ un

j − λ{F(u− j+1/2, u+ j+1/2) − F(u− j−1/2, u+ j−1/2)}. ◮ u− j+1/2, u+ j+1/2 are obtained by ENO/WENO reconstruction.

u±

j+1/2 = p± E/W (xj+1/2)

MPP finite volume scheme by Zhang & Shu

◮ Numerical formulation:

¯ un+1

j

= ¯ un

j − λ{F(u− j+1/2, u+ j+1/2) − F(u− j−1/2, u+ j−1/2)}. ◮ u− j+1/2, u+ j+1/2 are obtained by ENO/WENO reconstruction.

u±

j+1/2 = p± E/W (xj+1/2) ◮ A sufficient condition: rescaling p± E/W (x)

p±

E/W (x) = θ(p± E/W (x) − ¯

u) + ¯ u such that um ≤ p±

E/W (xα) ≤ uM

at quadratue points xα (include the end points).

Development

◮ The rescaling technique has been fitted into high order finite

volume and DG schemes on various mesh structures.

Development

◮ The rescaling technique has been fitted into high order finite

volume and DG schemes on various mesh structures.

◮ It has been applied to positivity preserving high order finite

volume and finite difference schemes for compressible Euler equations. Zhang & Shu 2010, 2011, 2012

Development

◮ The rescaling technique has been fitted into high order finite

volume and DG schemes on various mesh structures.

◮ It has been applied to positivity preserving high order finite

volume and finite difference schemes for compressible Euler equations. Zhang & Shu 2010, 2011, 2012

◮ The nonexistence of high order MPP finite difference scheme

is the motivation of this work.

The flow chart of high order method

◮ Reconstruction → Numerical flux → Evolution: Runge-Kutta

The flow chart of high order method

◮ Reconstruction → Numerical flux → Evolution: Runge-Kutta

↑ ↑ Zhang & Shu rescaling MPP

The flow chart of high order method

◮ Reconstruction → Numerical flux → Evolution: Runge-Kutta

↑ ↑ Zhang & Shu rescaling MPP

◮ Difficult to generalize to high order finite difference scheme:

no reference value of ˆ Hj+1/2 for rescaling p(x).

The flow chart of high order method

◮ Reconstruction → Numerical flux → Evolution: Runge-Kutta

↑ ↑ Zhang & Shu rescaling MPP

◮ Difficult to generalize to high order finite difference scheme:

no reference value of ˆ Hj+1/2 for rescaling p(x). An alternative:

◮ Reconstruction → Numerical flux → Evolution

↑ Parametrized flux limiters

Parametrized flux limiters

Looking for limiters of the type ˜ Hj+1/2 = θj+1/2(ˆ Hj+1/2 − ˆ hj+1/2) + ˆ hj+1/2 (4) such that um ≤ un

j − λ(˜

Hj+1/2 − ˜ Hj−1/2) ≤ uM. (5) ˆ hj+1/2 is the first order monotone flux.

Parametrized flux limiters

Looking for limiters of the type ˜ Hj+1/2 = θj+1/2(ˆ Hj+1/2 − ˆ hj+1/2) + ˆ hj+1/2 (4) such that um ≤ un

j − λ(˜

Hj+1/2 − ˜ Hj−1/2) ≤ uM. (5) ˆ hj+1/2 is the first order monotone flux.

◮ θj+1/2 = 0,

j = 1, 2, 3, ...: first order scheme with MPP property.

Parametrized flux limiters

Looking for limiters of the type ˜ Hj+1/2 = θj+1/2(ˆ Hj+1/2 − ˆ hj+1/2) + ˆ hj+1/2 (4) such that um ≤ un

j − λ(˜

Hj+1/2 − ˜ Hj−1/2) ≤ uM. (5) ˆ hj+1/2 is the first order monotone flux.

◮ θj+1/2 = 0,

j = 1, 2, 3, ...: first order scheme with MPP property.

◮ θj+1/2 = 1,

j = 1, 2, 3, ...: high order scheme most likely without MPP property.

Parametrized flux limiters

Looking for limiters of the type ˜ Hj+1/2 = θj+1/2(ˆ Hj+1/2 − ˆ hj+1/2) + ˆ hj+1/2 (4) such that um ≤ un

j − λ(˜

Hj+1/2 − ˜ Hj−1/2) ≤ uM. (5) ˆ hj+1/2 is the first order monotone flux.

◮ θj+1/2 = 0,

j = 1, 2, 3, ...: first order scheme with MPP property.

◮ θj+1/2 = 1,

j = 1, 2, 3, ...: high order scheme most likely without MPP property.

◮ θj+1/2 exists, locally defined.

[Xu, Math Comp 2013], online. [Liang & Xu, JSC 2014].

Comments

◮ Early work of flux limiters (TVD stability): [P L Roe, 1982;

Van Leer, 1974; R. F. Warming AND R. M. Beam, 1976; P K Sweby, 1984]

◮ The parametrized limiters provide a sufficient condition for

high order conservative scheme to be maximum principle preserving.

◮ The accuracy of the high order scheme with MPP limiters will

be affected by CFL number.

◮ Question remains: Does it produce MPP solution with high

• rder accuracy when applied to FD ENO/WENO?

Numerical experiments

Third order finite volume ENO scheme with MPP parametrized limiters Consider f (u) = u, u(x, 0) = sin(2πx) on the interval [0, 1] with periodic BC. Let T be 1, one period of evolution. TVD third order Runge-Kutta time discretization is used in all the tests.

Table : Third order finite volume ENO scheme with parametrized limiters: advection equation

CFL=0.15 CFL=0.2 CFL=0.25 N L∞ error

• rder

L∞ error

• rder

L∞ error

• rder

40 2.18E-3 2.04E-3 2.13E-3 80 2.71E-4 3.00 2.56E-4 3.00 2.75E-4 2.95 160 3.37E-5 3.00 3.19E-5 3.00 4.15E-5 2.72 320 4.18E-6 3.01 3.99E-6 3.00 8.28E-6 2.32 640 5.18E-7 3.01 4.98E-7 3.00 2.07E-6 1.99

Numerical experiments

Third order finite volume scheme with MPP parametrized limiters for Burgers’ equation In the case of inviscid Burgers’ equation f (u) = u2

2 with initial

condition u(x, 0) = 0.5(1 + sin(2πx)) on the interval [0, 1], let the simulation run to T = 0.5

π .

Table : Third order finite volume scheme with parametrized limiters: Burgers’ equation

CFL=0.15 CFL=0.2 CFL=0.4 N L∞ error

• rder

L∞ error

• rder

L∞ error

• rder

40 1.90E-3 1.89E-3 1.87E-3 80 2.80E-4 2.75 2.80E-4 2.75 2.79E-4 2.75 160 3.71E-5 2.91 3.71E-5 2.91 3.69E-5 2.91 320 4.67E-6 2.98 4.67E-6 2.98 9.03E-6 2.03 640 5.84E-7 3.00 5.84E-7 3.00 2.30E-6 1.97

Numerical experiments

Third order finite difference ENO scheme with MPP parametrized limiters The linear advection problem.

Table : Third order finite difference ENO scheme with parametrized limiters: Advection equation

k h1.5 max|f ′(u)|=1 k hmax|f ′(u)|=0.1 k hmax|f ′(u)|=0.25

N L∞ error

• rder

L∞ error

• rder

L∞ error

• rder

40 2.16-3 2.15E-3 2.27-3 80 2.72E-4 2.98 2.94E-4 3.00 4.07E-4 2.48 160 3.39E-5 3.00 3.49E-5 2.94 9.80E-5 2.05 320 4.22E-6 3.00 5.15E-6 2.76 2.47E-5 1.98 640 5.17E-7 3.00 9.97E-7 2.36 6.17E-6 2.00

Numerical experiments

Third order finite difference scheme with MPP parametrized limiters The Burgers’ equation.

Table : Third order finite difference scheme with parametrized limiters: Burgers’ equation

k h1.5 max|(f ′(u))|=1 k hmax|f ′(u)|=0.2 k hmax|f ′(u)|=0.4

N L∞ error

• rder

L∞ error

• rder

L∞ error

• rder

40 2.43E-3 2.43E-3 2.43-3 80 3.21E-4 2.92 3.21E-4 2.92 3.16E-4 2.93 160 3.98E-5 3.01 3.98E-5 3.01 5.05E-5 2.64 320 4.84E-6 3.04 4.83E-6 3.04 1.33E-5 1.91 640 5.94E-7 3.02 9.22E-7 2.39 3.41E-6 1.96

Third order TVD Runge-Kutta

u(1) = un + ∆tL(un), u(2) = 3 4un + 1 4(u(1) + ∆tL(u(1))), (6) un+1 = 1 3un + 2 3(u(2) + ∆tL(u(2))). The previous approach requires un + ∆tL(un) ≤ uM, (7) u(1) + ∆tL(u(1)) ≤ uM, (8) u(2) + ∆tL(u(2)) ≤ uM. (9)

General MPP flux limiters:

◮ A rewriting of third order TVD RK FD WENO scheme:

un+1

j

= un

j − λ(ˆ

Hrk

j+ 1

2 − ˆ

Hrk

j− 1

2 ),

(10) where ˆ Hrk

j+ 1

2

. = 1 6 ˆ Hn

j+ 1

2 + 2

3 ˆ H(2)

j+ 1

2 + 1

6 ˆ H(1)

j+ 1

2 .

◮ Parametrized MPP flux limiters are applied to the final step

(the integral form along temporal direction) ˜ Hrk

j+ 1

2 = θj+ 1 2 (ˆ

Hrk

j+ 1

2 − ˆ

hj+ 1

2 ) + ˆ

hj+ 1

2

(11)

◮ Applied to incompressible flow problem. [Xiong, Qiu, Xu, JCP

2013]

Numerical tests: FD3RK3 for ut + ux = 0

Table : T = 0.5, CFL = 0.6, u0(x) = sin4(x), without limiters.

N L1 error

• rder

L∞ error

• rder

(uh)min 20 2.21e-02 – 4.43e-02 –

• 2.26E-02

40 3.49e-03 2.66 6.48e-03 2.77

• 3.69E-03

80 4.54e-04 2.94 8.77e-04 2.89

• 5.16E-04

160 5.76e-05 2.98 1.11e-04 2.98

• 6.68E-05

320 7.22e-06 3.00 1.40e-05 3.00

• 8.36E-06

Table : T = 0.5, CFL = 0.6, u0(x) = sin4(x), with limiters.

N L1 error

• rder

L∞ error

• rder

(uh)min 20 1.83e-02 – 4.43e-02 – 3.55E-14 40 3.24e-03 2.50 6.48e-03 2.77 1.23E-14 80 4.57e-04 2.82 8.77e-04 2.89 6.38E-23 160 5.75e-05 2.99 1.23e-04 2.83 1.72E-16 320 7.22e-06 2.99 1.71e-05 2.85 9.61E-22

Where is the proof?

◮ FD: For 1-D general nonlinear problem, the general MPP flux

limiters does not affect the third order accuracy when ˆ h is local Lax-Friedrich flux.

◮ FD: When ˆ

h is global Lax-Friedrich flux, MPP and third order accuracy are obtained when CFL is less than 0.886. [Xiong, Qiu & Xu, JCP, 2013]. O FV: When applied to FV WENO with arbitrary order solving problems ut + f (u)x = 0, f ′(u) > 0, there is no extra CFL requirement when ˆ h is a first order up-winding flux.

◮ FV: When applied to FV WENO solving general

ut + f (u)x = 0 with LF fluxes (Local or Global), MPP and third roder accuracy is obtained without extra CFL

• requirement. [Xiong, Qiu & Xu, Submitted to SIAM NUMER

ANALYSIS].

Key element used to prove case [O]

The entropy solution satisfies the maximum principle in the form of um ≤ ¯ uj(tn) − λ(ˇ fj+ 1

2 − ˇ

fj− 1

2 ) ≤ uM.

(12) where λ = ∆t/∆x and ¯ uj(t) = 1 ∆x xj+1/2

xj−1/2

u(x, t)dx, ˇ fj−1/2 = 1 ∆t tn+1

tn

f (u(xj−1/2, t))dt. (13)

Theorem

Assuming f ′(u) > 0 and λ maxu |f ′(u)| ≤ 1, we have ¯ uj(tn) − λ(ˇ fj+ 1

2 − f (¯

uj−1(tn))) ≤ uM (14) if u(x, t) is the entropy solution to (2) subject to initial data u0(x).

Proof of the theorem

Consider the problem (2) with an initial condition at time level tn, ˜ u(x, tn) =

• u(x, tn)

x ≥ xj− 1

2 ,

¯ uj−1(tn) x < xj− 1

2 ,

(15) here u(x, tn) is the exact solution of (2) at time level tn. Let ˜ u(x, t) be the entropy solution with initial data ˜ u(x, tn), we have ¯ ˜ uj(tn) = ¯ uj(tn). (16) Since f ′(u) > 0, we have f (˜ u(xj− 1

2 , t)) = f (¯

uj−1(tn)), (17) for t ∈ [tn, tn+1]. Since λ maxu |f ′(u)| ≤ 1, the characteristic starting from xj− 1

2 would not hit the side xj+ 1 2 , therefore

˜ u(xj+ 1

2 , t) = u(xj+ 1 2 , t)

(18) for t ∈ [tn, tn+1].

End of the proof

Since ˜ u satisfies the maximum principle um ≤ ˜ u ≤ uM, we have um ≤ ¯ ˜ un+1

j

= ¯ ˜ un

j − λ(ˇ

˜ fj+ 1

2 − ˇ

˜ fj− 1

2 ) ≤ uM,

where ˇ ˜ fj−1/2 = 1 ∆t tn+1

tn

f (˜ u(xj−1/2, t))dt. (19) Substituting (16), (17) and (18) into the above inequality, it follows that um ≤ ¯ uj(tn) − λ(ˇ fj+ 1

2 − f (¯

uj−1(tn)) ≤ uM.

Results and development

MPP flux limiters are generalized to convection-dominated diffusion equation [Jiang & Xu, SIAM JSC 2013]. ∂u ∂t + ∂f (u) ∂x = ∂2A(u) ∂x2 , (20) where A′(u) ≥ 0. The porous medium equation ut = (um)xx, x ∈ R, t > 0, (21) Barenblatt solution Bm(x, t) = t−k[(1 − k(m − 1) 2m |x|2 t2k )+]1/(m−1), (22)

6 4 2 2 4 6 0. 0.2 0.4 0.6 0.8 1. X U 6 4 2 2 4 6 0. 0.2 0.4 0.6 0.8 1. X U

Figure : m=8

Results and development

Generalized to Euler system for positive density, pressure and internal energy [Xiong, Qiu & Xu, Revision submitted to JCP]. High Mach number astrophysical jets: Two high Mach number astrophysical jets without the radiative cooling

Figure : Top: density of Mach 80 at T = 0.07; Bottom: density of Mach 2000 at T = 0.001.

Results and development

Vlasov equation simulation: Ion-acoustic turbulence [Xiong, Qiu & Xu, accepted by JCP]. ∂tfe + v∂xfe − E(t, x)∂vfe = 0, (23) ∂tfi + v∂xfi + E(t, x) Mr ∂vfi = 0, (24) E(t, x) = −∇φ(t, x), −∆φ(t, x) = ρ(t, x). (25)

Conclusion and ongoing projects

• 1. A sufficient condition for conservative high order schemes to

be MPP is provided. Partial results are provided as for the choice of CFL and its effect on the designed order of accuracy.

• 2. Application to MHD problem. With H. Yang, F.-Y. Li.
• 3. Application to high order conservative high order FD WENO

method solving ideal MHD equation (CT framework). With

• Y. Liu, Q. Tang, A. Christlieb.
• 4. Develop MPP flux limiters for high order methods with IMEX

RK time integration solving general convection diffusion

• equation. With Y. Jiang.