High Order Asymptotic Preserving Schemes for Some Discrete-Velocity - - PowerPoint PPT Presentation

1

✬ ✫ ✩ ✪

High Order Asymptotic Preserving Schemes for Some Discrete-Velocity Kinetic Equations

Fengyan Li Department of Mathematical Sciences Rensselaer Polytechnic Institute Troy, NY, USA Joint work with: J. Jang (UC Riverside), J.-M. Qiu (Houston), T. Xiong (Houston)

lif@rpi.edu

2

✬ ✫ ✩ ✪ Outline

• Discrete-velocity kinetic equations
• High order asymptotic preserving (AP) methods
• Theoretical results
• Numerical examples

lif@rpi.edu

3

✬ ✫ ✩ ✪

Discrete-velocity Kinetic Equations

Consider the discrete-velocity model ∂tf + v · ∇xf = C(f) (1)

f = f(x, v, t): distribution function of particles dependent of time t > 0, position x, and velocity v ∈ {−1, 1}. C(f): collision operator (

• C(f)dv = 0).

lif@rpi.edu

4

✬ ✫ ✩ ✪ We are particularly interested in the discrete-velocity model in a diffusive scaling (under the scaling of t′ := ε2t and x′ := εx): ε∂tf + v · ∇xf = 1 εC(f) (2)

The parameter ε > 0 can be regarded as the mean free path of the particles, and it measures the distance of the system to the equilibrium

• state. The smaller ε is, the closer the system is to the equilibrium.

lif@rpi.edu

5

✬ ✫ ✩ ✪ Examples: Jin-Pareschi-Toscani (1998) dµ: discrete Lebesgue measure on {−1, 1} ·: f =

• fdµ = 1

2(f(x, v = −1, t) + f(x, v = 1, t))

C(f) when ε → 0, with ρ = f (a.1) f − f ∂tρ = ∂xxρ (a.2) f − f + Aεvf ∂tρ + A∂xρ = ∂xxρ (a.3) Cfm(f − f) ∂tρ =

C 1−m ∂xx(ρ1−m), m = 1

(a.4) f − f + Cε

• f2 − (f − f)2

v ∂tρ + C∂xρ2 = ∂xxρ Note:

• (a.1): telegraph equation
• ∂tρ =

1 1−m ∂xx(ρ1−m) with m < 0: porous medium equation

lif@rpi.edu

6

✬ ✫ ✩ ✪ In some applications, ε may differ in several orders of magnitude from the rarefied regime (ε = O(1)) to the hydrodynamic (diffusive) regime (ε << 1). It is desirable to design a class of numerical methods which work uniformly with respect to the parameter ε.

lif@rpi.edu

7

✬ ✫ ✩ ✪ Objective: to design high order asymptotic preserving methods for the discrete-velocity kinetic equation in the diffusive scaling Asymptotically preserving (AP) methods: Jin, Levermore, Naldi, Pareschi, Degond, Toscani, Klar, Filbet, Carrillo, Lemou, Mieussens, Hauck, Liu, · · ·

• uniformly stable with respect to ε ranging from O(1) to 0;
• When ε → 0, the methods are consistent for the limiting

equation on fixed mesh.

lif@rpi.edu

8

✬ ✫ ✩ ✪

Numerical challenges and considerations ε∂tf + v · ∇xf = 1 εC(f)

• Stiffness in both the convective and the collision terms
• The characteristic speed of the homogeneous hyperbolic part is 1

ε.

Standard schemes for hyperbolic problems with stiff relaxation (implicit treatment for the collision term + explicit treatment for the

convection term) may impose a restrictive stability condition: ∆t ≈ ε∆x (3)

• Not all stable schemes can capture the diffusive limit when ε → 0 on

under-resolved meshes (∆t, ∆x >> ε)

• Easy to solve

lif@rpi.edu

9

✬ ✫ ✩ ✪

High Order Asymptotic Preserving Methods

Three components

• Micro-macro decomposition of the equation
• Discontinuous Galerkin (DG) spatial discretization: numerical

flux

• Globally stiffly accurate implicit-explicit (IMEX) Runge-Kutta

temporal discretization: implicit-explicit strategy Major references:

• Jin-Pareschi-Toscani (1998)
• Lemou-Mieussens (2010), Liu-Mieussens (2010)
• Boscarino-Pareschi-Russo (2013)

lif@rpi.edu

10

✬ ✫ ✩ ✪ Main results: a family of high order methods are proposed for discrete-velocity kinetic equations in the diffusive scaling.

• For the telegraph equation, when a first order temporal

discretization is applied, uniform stability is established with respect to ε. Error estimates are also obtained for any ε.

• Formal asymptotic analysis shows that the proposed schemes

in the limit of ε → 0 provide explicit and consistent high order methods for the limiting equations.

lif@rpi.edu

11

✬ ✫ ✩ ✪

• 1. Micro-macro decomposition

Consider the Hilbert space L2(dµ) with the inner product ·, ·: f1, f2 =

• f1f2dµ,

and an orthogonal projection Π: Πf = f. Let f = Πf + (I − Π)f =: ρ + εg, then the micro-macro decomposition of (2) is ∂tρ + ∂xvg = 0 (4a) ∂tg + 1 ε2 v∂xρ + 1 ε{I − Π}(v∂xg) = 1 ε3 C(ρ + εg) (4b)

Note: Solving for f from (2) is equivalent to solving for ρ and g from (4).

lif@rpi.edu

12

✬ ✫ ✩ ✪

• 2. Discontinuous Galerkin (DG) spatial discretization

Mesh: Ii = [xi− 1

2 , xi+ 1 2 ],

hi = |Ii| Discrete space: Uh = U k

h = {u : u ∈ P k(Ii), ∀i}, with P k(Ii) consisting

• f polynomials of degree at most k

xi− 1

2

xi+ 1

2

u+

i− 1

2

u−

i+ 1

2

Ii Notations:

u± = u±(x) = lim

∆x→0± u(x + ∆x)

average : {u} = 1 2 (u− + u+) jump : [u] = u+ − u−

lif@rpi.edu

13

✬ ✫ ✩ ✪

Semi-discrete DG method: Look for ρh(·, t), gh(·, v, t) ∈ Uh, such that ∀φ, ψ ∈ Uh, and i,

• Ii

∂tρhφdx −

• Ii

vgh∂xφdx + vghi+ 1

2 φ−

i+ 1

2

− vghi− 1

2 φ+

i− 1

2

= 0, (5a)

• Ii

∂tghψdx + 1 ε

• Ii

(I − Π)Dh(gh; v) ψdx − 1 ε2

• Ii

vρh∂xψdx − v ρh,i+ 1

2 ψ−

i+ 1

2

+ v ρh,i− 1

2 ψ+

i− 1

2

• = 1

ε3

• Ii

C(ρh + εgh)ψdx. (5b)

Numerical flux:

Alternating: vg = vg−, ˆ ρ = ρ+, or vg = vg+, ˆ ρ = ρ− (6a) Central: vg = {vg}, ˆ ρ = {ρ} (6b)

lif@rpi.edu

14

✬ ✫ ✩ ✪ Dh(gh; v) ∈ Uh is an upwind approximation of v∂xg, defined by

• Dh(gh; v)ψdx =
• i
• −
• Ii

vgh∂xψdx + (vgh)i+ 1

2 ψ−

i+ 1

2

− (vgh)i− 1

2 ψ+

i− 1

2

• (7)

with

• vgh :=

   vg−

h ,

if v > 0 vg+

h ,

if v < 0 = v{gh} − |v| 2 [gh]. (8)

lif@rpi.edu

15

✬ ✫ ✩ ✪

Method in a compact form: Look for ρh(·, t), gh(·, v, t) ∈ Uh, such that ∀φ, ψ ∈ Uh, and i,

• ∂tρhφdx + ah(gh, φ) = 0
• ∂tghψdx + 1

εbh,v(gh, ψ) − v ε2 dh(ρh, ψ) = − 1 ε2 sh,v(ρh, gh, ψ) − s′

h,v(gh, ψ)

where ah(gh, φ) = −

• i
• Ii

vgh∂xφdx −

• i
• vghi− 1

2 [φ]i− 1 2 ,

bh,v(gh, ψ) =

• (I − Π)Dh(gh; v)ψdx =
• (Dh(gh; v) − Dh(gh; v))ψdx,

dh(ρh, ψ) =

• i
• Ii

ρh∂xψdx +

• i
• ρh,i− 1

2 [ψ]i− 1 2 ,

lif@rpi.edu

16

✬ ✫ ✩ ✪ and sh,v(ρh, gh, ψ) =              (gh, ψ) for (a.1) (gh − Avρh, ψ) for (a.2) (Cρm

h gh, ψ)

for (a.3) (gh − Cvρ2

h, , ψ)

for (a.4) s′

h,v(gh, ψ) =

             for (a.1) for (a.2) for (a.3) (Cvg2

h, ψ)

for (a.4)

lif@rpi.edu

17

✬ ✫ ✩ ✪

• 3. Globally stiffly accurate implicit-explicit (IMEX) temporal

discretization DG-IMEX1: ρn+1

h

− ρn

h

∆t φdx+ah(gn

h, φ) = 0,

(9a) gn+1

h

− gn

h

∆t ψdx+1 εbh,v(gn

h, ψ) − v

ε2 dh(ρn+1

h

, ψ) = − 1 ε2 sh,v(ρn+1

h

, gn+1

h

, ψ) − s′

h,v(gn h, ψ).

(9b)

Note: (1) The terms dh and sh,v are treated implicitly. (2) One first solves ρn+1

h

from (9a), and then gn+1

h

from (9b).

lif@rpi.edu

18

✬ ✫ ✩ ✪ Formal asymptotic analysis: As an example, apply the DG-IMEX1 to

the telegraph equation with C(f) = f − f = −εg:

• The equation in its micro-macro formulation:

∂tρ + ∂xvg = 0 ∂tg + 1 ε2 v∂xρ + 1 ε{I − Π}(v∂xg) = − 1 ε2 g

• DG-IMEX1 scheme: ∀φ, ψ ∈ Uh,

ρn+1

h

− ρn

h

∆t φdx+ah(gn

h, φ) = 0

gn+1

h

− gn

h

∆t ψdx+1 εbh,v(gn

h, ψ) − v

ε2 dh(ρn+1

h

, ψ) = − 1 ε2

• gn+1

h

ψdx lif@rpi.edu

19

✬ ✫ ✩ ✪ Formal asymptotic analysis: As an example, apply the DG-IMEX1 to

the telegraph equation with C(f) = f − f = −εg:

• The equation in its micro-macro formulation:

∂tρ + ∂xvg= 0 (12a) ∂tg + 1 ε2 v∂xρ + 1 ε{I − Π}(v∂xg) = − 1 ε2 g (12b)

• DG-IMEX1 scheme: ∀φ, ψ ∈ Uh,

ρn+1

h

− ρn

h

∆t φdx+ah(gn

h, φ) = 0

(13a) gn+1

h

− gn

h

∆t ψdx+1 εbh,v(gn

h, ψ) − v

ε2 dh(ρn+1

h

, ψ) = − 1 ε2

• gn+1

h

ψdx (13b) lif@rpi.edu

20

✬ ✫ ✩ ✪

When ε → 0,

• The limiting equation:

∂tρ + ∂xvg = 0 (14a) −v∂xρ = g (14b)

• DG-IMEX1 scheme in the limit: ∀φ, ψ ∈ Uh, ∀i
• Ii

ρn+1

h

− ρn

h

∆t φdx −

• Ii

vgn

h∂xφdx +

vgn

hi+ 1

2 φ−

i+ 1

2

− vgn

hi− 1

2 φ+

i− 1

2

= 0 v

• Ii

ρn

h∂xψdx −

ρn

hi+ 1

2 ψ−

i+ 1

2

+ ρn

hi− 1

2 ψ+

i− 1

2

• =
• Ii

gn

hψdx

Remark: In the limit of ε → 0, the proposed scheme becomes a consistent and

explicit discretization for the limiting equation on fixed mesh.

lif@rpi.edu

21

✬ ✫ ✩ ✪

For the two-velocity model, we can further write the limiting equation and scheme in terms of macro quantities ρ and j := vg

• The limiting equation:

∂tρ + ∂xj = 0 −∂xρ = j This is ∂tρ = ∂xxρ.

• DG-IMEX1 scheme in the limit: ∀φ, ψ ∈ Uh, ∀i
• Ii

ρn+1

h

− ρn

h

∆t φdx −

• Ii

jn

h ∂xφdx +

jn

h i+ 1

2 φ−

i+ 1

2

− jn

h i− 1

2 φ+

i− 1

2

= 0

• Ii

ρn

h∂xψdx −

ρn

hi+ 1

2 ψ−

i+ 1

2

+ ρn

hi− 1

2 ψ+

i− 1

2

=

• Ii

jn

h ψdx

This is a local DG method for ∂tρ = ∂xxρ. Cockburn-Shu (1998) lif@rpi.edu

22

✬ ✫ ✩ ✪

IMEX Runge-Kutta methods: in the double Butcher Tableau representation

below, the left is for the explicit term, and the right is for the implicit term

˜ c ˜ A ˜ bT c A bT ˜ A: lower triangular with zero diagonal entries A: lower triangular with nonzero diagonal entries (diagonally implicit RK

(DIRK))

1st order:

1 1 1 1 1 1

lif@rpi.edu

23

✬ ✫ ✩ ✪ Globally stiffly accurate: the last row of A is the same as b; the last row of ˜ A is the same as ˜ b; in addition, the last entries of c and ˜ c are 1. ˜ c ˜ A ˜ bT c A bT

Remark: The implicit-explicit strategy and the globally stiffly accurate property

• f the time discretizations ensure that numerical solutions from both inner

stages and discrete times will be relaxed to approximate the limiting problem when ε → 0.

lif@rpi.edu

24

✬ ✫ ✩ ✪

Boscarino-Pareschi-Russo (2013)

2nd order: ARS(2, 2, 2) (with γ = 1 −

1 √ 2 and δ = 1 − 1 2γ )

γ γ 1 δ 1 − δ δ 1 − δ γ γ 1 1 − γ γ 1 − γ γ

3rd order: ARS(4, 4, 3)

1/2 1/2 2/3 11/18 1/18 1/2 5/6 −5/6 1/2 1 1/4 7/4 3/4 −7/4 1/4 7/4 3/4 −7/4 1/2 1/2 2/3 1/6 1/2 1/2 −1/2 1/2 1/2 1 3/2 −3/2 1/2 1/2 3/2 −3/2 1/2 1/2

lif@rpi.edu

25

✬ ✫ ✩ ✪ Formal asymptotic analysis: the fully-discrete DG-IMEX scheme in the limit of ε → 0, becomes a consistent scheme for the limiting

• equation. This scheme involves high order DG discretization in

space where the numerical fluxes are naturally from the DG-IMEX schemes, and an explicit time discretization which is the explicit part of the IMEX scheme.

lif@rpi.edu

26

✬ ✫ ✩ ✪

Theoretical Results

Let ||φ|| = ||φ||L2(Ωx) and |||ψ||| = (||ψ||2)1/2. Without loss of generality, the mesh is assumed to be uniform with h = hi, ∀i. Inverse inequality: there exist constants Cinv = Cinv(k), ˆ Cinv = ˆ Cinv(k), s.t. ∀w ∈ P k([a, b]),

|w(y)|2(b − a) ≤ Cinv b

a

w(x)2dx, with y = a or b (b − a)2 b

a

|wx(x)|2dx ≤ ˆ Cinv b

a

w(x)2dx.

lif@rpi.edu

27

✬ ✫ ✩ ✪

Theorem (Stability) When the DG-IMEX1 method is applied to the telegraph equation with C(f) = f − f = −εg, given g0

h = 0, the

following stability result holds for the numerical solution,

||ρn+1

h

||2 + ε2|||gn

h|||2 ≤ ||ρn h||2 + ε2|||gn−1 h

|||2, ∀n (15)

under the CFL condition

∆t ≤ ∆tstab :=   

1 2 ˆ Cinv+8C2

inv

(h + 2Cinv min(ε, 2Cinv

ˆ Cinv h))h,

for k ≥ 1,

1 4C2

inv (h + 2Cinvε)h = 1

2εh + 1 4 h2,

for k = 0. (16) Note: (1) The stability result is uniform with respect to ε. (2) For k = 0, numerical stability requires ∆t = O(h2) in diffusive regime (ε << 1), and ∆t = O(εh) in rarefied regime (ε = O(1)). (3) For k ≥ 1, numerical stability requires ∆t = O(h2). It is conjectured that the restrictive timestep in the rarefied regime (ε = O(1)) can be improved and/or removed with higher order time discretizations. (supported by numerics)

lif@rpi.edu

28

✬ ✫ ✩ ✪

Theorem (Error estimates) When the DG-IMEX1 method is applied to the telegraph equation with C(f) = f − f = −εg, given g0

h = 0, the

following error estimates hold for the numerical solution,

||ρ(tn) − ρn

h||2 + ε2|||g(tn−1) − gn−1 h

|||2 ≤ C⋆   

• (1 + ε4)∆t2 + (1 + ε2)h2k+2 + εh2k+σ(k)

(alternating flux)

• (1 + ε4)∆t2 + ε2h2k+2 + h2k+σ(k)−1 + εh2k+σ(k)

(central flux)

for n : n∆t ≤ T under the CFL condition ∆t < ∆tstab. Here σ(k) =    1 for k ≥ 1, 2 for k = 0. (17) The positive constant C⋆ is independent of h, ∆t, and n. It depends on T, k, and some Sobolev norms of the exact solution (hence ε). lif@rpi.edu

29

✬ ✫ ✩ ✪ Spatial accuracy orders established by the error estimates

k ≥ 1, ε = 0 k ≥ 1, ε = 0 k = 0 alternative k + 1

2

k + 1 1 central k k

1 2

lif@rpi.edu

30

✬ ✫ ✩ ✪

Numerical Examples

DG(k+1)-IMEX(k+1): DG method with k-th order polynomials; (k + 1)-th order IMEX time discretization How to choose time step? ∆t = Cconvεh + Cdiffh2 (18) with

• DG1-IMEX1: Cconv = 0.5 and Cdiff = 0.25 (based on our analysis)
• DG2-IMEX2: Cconv = 0.5 and Cdiff = 0.01
• DG3-IMEX3: Cconv = 0.25 and Cdiff = 0.006

lif@rpi.edu

31

✬ ✫ ✩ ✪

Define j := vg

Example 1 (telegraph equation): C(f) = f − f. The exact solution is    ρ(x, t) = 1

r exp(rt) sin(x),

r =

−2 1+ √ 1−4ε2 ,

j(x, t) = exp(rt) cos(x) (19)

lif@rpi.edu

32

✬ ✫ ✩ ✪ Table 1: L1 error and order of ρ and j, T = 1.0, DG1-IMEX1 with the alternating flux.

N L1 error of ρ

• rder

L1 error of j

• rder

ε = 0.5 10 6.04E-02 – 7.46E-02 – 20 2.19E-02 1.46 3.38E-02 1.14 40 9.20E-03 1.25 1.60E-02 1.08 80 4.19E-03 1.14 7.81E-03 1.03 ε = 10−2 10 3.79E-02 – 8.05E-02 – 20 1.78E-02 1.09 3.77E-02 1.09 40 8.79E-03 1.02 1.85E-02 1.03 80 4.36E-03 1.01 9.22E-03 1.01 ε = 10−6 10 3.79E-02 – 8.03E-02 – 20 1.79E-02 1.08 3.76E-02 1.09 40 8.82E-03 1.02 1.85E-02 1.02 80 4.38E-03 1.01 9.21E-03 1.01

lif@rpi.edu

33

✬ ✫ ✩ ✪ Table 2: L1 error and order of ρ and j, T = 1.0, DG2-IMEX2 with the alternating flux.

N L1 error of ρ

• rder

L1 error of j

• rder

ε = 0.5 10 1.35E-03 – 2.36E-03 – 20 3.00E-04 2.17 4.90E-04 2.27 40 7.23E-05 2.05 1.14E-04 2.10 80 1.79E-05 2.01 2.76E-05 2.04 ε = 10−2 10 4.83E-03 – 4.94E-03 – 20 1.19E-03 2.02 1.19E-03 2.06 40 2.96E-04 2.01 2.97E-04 2.00 80 7.40E-05 2.00 7.40E-05 2.00 ε = 10−6 10 4.82E-03 – 4.93E-03 – 20 1.19E-03 2.02 1.18E-03 2.06 40 2.96E-04 2.00 2.96E-04 2.00 80 7.40E-05 2.00 7.40E-05 2.00

lif@rpi.edu

34

✬ ✫ ✩ ✪ Table 3: L1 error and order of ρ and j, T = 1.0, DG3-IMEX3 with the alternating flux.

N L1 error of ρ

• rder

L1 error of j

• rder

ε = 0.5 10 6.33E-05 – 9.48E-05 – 20 7.54E-06 3.07 1.15E-05 3.04 40 9.31E-07 3.02 1.44E-06 3.00 80 1.16E-07 3.01 1.80E-07 3.00 ε = 10−2 10 2.53E-04 – 2.46E-04 – 20 3.11E-05 3.03 3.11E-05 2.98 40 3.89E-06 3.00 3.89E-06 3.00 80 4.87E-07 3.00 4.87E-07 3.00 ε = 10−6 10 2.53E-04 – 2.46E-04 – 20 3.11E-05 3.03 3.11E-05 2.98 40 3.89E-06 3.00 3.89E-06 3.00 80 4.87E-07 3.00 4.87E-07 3.00

lif@rpi.edu

35

✬ ✫ ✩ ✪ Let k be the polynomial degree of the discrete space. It is observed that

• with alternating flux: (k + 1)-th order
• with central flux: k-th order for odd k, and (k + 1)-th order for

even k

lif@rpi.edu

36

✬ ✫ ✩ ✪ Example 2 (telegraph equation): C(f) = f − f. Consider a Riemann problem with the initial condition    ρL = 2.0, jL = 0.0, −1 < x < 0; ρR = 1.0, jR = 0.0, 0 < x < 1. (20)

lif@rpi.edu

37

✬ ✫ ✩ ✪

x ρ

• 1
• 0.5

0.5 1 1 1.2 1.4 1.6 1.8 2

reference solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x ρ

• 1
• 0.5

0.5 1 1 1.2 1.4 1.6 1.8 2

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x j

• 1
• 0.5

0.5 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7

reference solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x j

• 1
• 0.5

0.5 1 0.5 1 1.5

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

Figure 1: Numerical solution of DG(k+1)-IMEX(k+1) (k = 0, 1, 2) with the alternating flux and h = 0.05. Left: ε = 0.7 at T = 0.25; Right: ε = 10−6 at T = 0.04. Top: ρ; Bottom: j. The reference solution is from DG3-IMEX3 with h = 0.004.

lif@rpi.edu

38

✬ ✫ ✩ ✪

x ρ

• 1
• 0.5

0.5 1 1 1.2 1.4 1.6 1.8 2

reference solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x ρ

• 1
• 0.5

0.5 1 1 1.2 1.4 1.6 1.8 2

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x j

• 1
• 0.5

0.5 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7

reference solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x j

• 1
• 0.5

0.5 1 0.5 1 1.5

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

Figure 2: Numerical solution of DG(k+1)-IMEX(k+1) (k = 0, 1, 2) with the central flux and h = 0.05. Left: ε = 0.7 at T = 0.25; Right: ε = 10−6 at T = 0.04. Top: ρ; Bottom: j. The reference solution is from DG3-IMEX3 with h = 0.004.

lif@rpi.edu

39

✬ ✫ ✩ ✪ Example 3 (viscous Burgers’ equation): The initial condition is given as two local Maxwellian    ρL = 2.0, −10 < x < 0 ρR = 1.0, 0 < x < 10 (21) with j = ρ2/(1 +

• 1 + ρ2ε2).

lif@rpi.edu

40

✬ ✫ ✩ ✪

x ρ

• 10
• 5

5 10 1 1.2 1.4 1.6 1.8 2

x ρ

• 10
• 5

5 10 1 1.2 1.4 1.6 1.8 2

x j

• 10
• 5

5 10 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

x j

• 10
• 5

5 10 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Figure 3: Numerical solution of DG3-IMEX3 with the alternating flux with h = 0.5 at T = 2.0. Left: ε = 0.4; Right: ε = 10−6. Top: ρ; Bottom: j. Symbol: numerical solution; Solid line: reference solution.

lif@rpi.edu

41

✬ ✫ ✩ ✪ Example 4 (porous medium equation): C(f) = fm(f − f) = −ερmg with m = −1. The initial condition is chosen to be the same as the Barenblatt solution for the limiting equation      ρ(x, t) =

1 R(t)

• 1 −
• x

R(t)

2 , j(x, t) = ρ

4x R(t)3 ,

|x| < R(t), ρ(x, t) = 0, j(x, t) = 0, |x| > R(t), where R(t) = [12(t + 1)]1/3, t ≥ 0.

Note:

• Nodal DG method is used.
• Alternating flux needs to be suitably chosen around the interface of ρ = 0.

lif@rpi.edu

42

✬ ✫ ✩ ✪

x ρ

• 10
• 5

5 10 0.05 0.1 0.15 0.2 0.25 0.3

exact solution DG2-IMEX2 DG3-IMEX3

x j

• 10
• 5

5 10

• 0.04
• 0.02

0.02 0.04

exact solution DG2-IMEX2 DG3-IMEX3

Figure 4: Numerical solution of DG(k+1)-IMEX(k+1) (k = 1, 2) with the alternating flux (Right half domain: vg = vg−, ˆ ρ = ρ+, left half domain vg = vg+, ˆ ρ = ρ−); ε = 10−6 with h = 0.5 at T = 3.0 compared with the exact solution. Left: ρ; Right: j.

Note: with k = 0, the scheme in the center element is inconsistent.

lif@rpi.edu

43

✬ ✫ ✩ ✪

x ρ

• 10
• 5

5 10 0.05 0.1 0.15 0.2 0.25 0.3

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

x j

• 10
• 5

5 10

• 0.04
• 0.02

0.02 0.04

exact solution DG1-IMEX1 DG2-IMEX2 DG3-IMEX3

Figure 5: Numerical solution of DG(k+1)-IMEX(k+1) (k = 0, 1, 2) with the central flux. ε = 10−6 with h = 0.5 at T = 3.0 compared with the exact solution. Left: ρ; Right: j.

lif@rpi.edu

44

✬ ✫ ✩ ✪

Summary: a family of high order methods are proposed for discrete-velocity kinetic equations in the diffusive scaling.

• For the telegraph equation, stability is established for DG-IMEX1

uniformly with respect to ε. Error estimates are also obtained for any ε.

• Formal asymptotic analysis shows that the proposed schemes in the

limit of ε → 0 provide explicit and consistent high order methods for the limiting equations.

• Numerical experiments demonstrate the performance of the methods

with higher order temporal accuracy and other collision operators.

• The proposed methods and most analysis can be naturally extended

to kinetic models with v ∈ {v1, · · · , vN}, or v ∈ [−v∗, v∗]. Current and future effort: kinetic equations in other scalings and (or) with general collision operators; boundary conditions; analysis lif@rpi.edu