Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson - - PowerPoint PPT Presentation

Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson equation for semiconductors

Vittorio Rispoli

In collaboration with:

• G. Dimarco (Univ. of Toulouse, France)
• L. Pareschi (Univ. of Ferrara, Italy)

HYP2012

University of Padova ITALY

Padova, 28.06.2012

MOSFET modelling

Figure: Metal-Oxide-Semiconductor Field-Effect Transistor (MOSFET) device modelling from http://nanohub.org/topics/MOSFETLabPage

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 2 / 34

Outline

1 General framework 2 Application to kinetic equations for semiconductor

Diffusive limit Parity equations

3 Discretization

Time discretization: IMEX schemes Phase-space discretization

4 Results 5 Penalization technique for the collision term

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 3 / 34

Prototype example

A prototype hyperbolic system of conservation laws with diffusive relaxation that we will use to illustrate the subsequent theory is the following: ut + vx = 0, vt + 1 ε2 p(u)x = − 1 ε2 (v − q(u)), (1) where p(u)′ > 0. Characteristic speeds: ±

• p(u)′

ε

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 4 / 34

Prototype example

In the small relaxation limit ε → 0, the behavior of the solution is, at least formally, governed by the convection-diffusion equation ut + q(u)x = p(u)xx, v = q(u) − p(u)x. (2) Ensure that the schemes possess the correct zero-relaxation limit: the numerical scheme applied to system (1) should be a consistent and stable scheme for the limit system (2) as the parameter ε → 0. asymptotic preservation (AP)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 5 / 34

Problems

1 Characteristic speeds of the hyperbolic part is of order 1/ε, standard

IMEX R-K schemes developed for hyperbolic systems with stiff relaxation become useless in such parabolic scaling, because the CFL condition would require ∆t = O(ε∆x).

2 Most previous works on asymptotic preserving schemes, in the limit of

infinite stiffness become consistent explicit schemes for the diffusive limit equation: such explicit schemes clearly suffer from the usual stability restriction ∆t = O(∆x2).

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 6 / 34

IMEX Runge-Kutta schemes for diffusive relaxation

The starting point is to reformulate problem (1) as the equivalent system 1 ut = −(v + µp(u)x)x

• Explicit

+ µp(u)xx

Implicit

, ε2vt = −p(u)x − v + q(u))

• Implicit

. Here µ = µ(ε) ∈ [0, 1] is a free parameter such that µ(0) = 1.

1Boscarino, Pareschi, Russo (2011)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 7 / 34

IMEX Runge-Kutta schemes for diffusive relaxation

This system can be written in the form ut = f1(u, v)

Explicit

+ f2(u)

• Implicit

ε2vt = g(u, v)

Implicit

, to which we can apply and study the properties of an IMEX-RK scheme.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 8 / 34

IMEX Runge-Kutta schemes for diffusive relaxation

Applying an Implicit-Explicit (IMEX) Runge-Kutta scheme 2 we obtain: Uk = un + ∆t

k−1

• j=1
• akj f1
• Uj, V j

+ ∆t

k

• j=1

akj f2

• Uj

, ε2V k = ε2v n + ∆t

k

• j=1

akj g

• Uj, V j

for the internal stages for k = 1, . . . , ν and un+1 = un + ∆t

ν

• k=1
• wk f1
• Uk, Rk

+ ∆t

ν

• k=1

wkf2

• Uk

ε2v n+1 = ε2v n + ∆t

ν

• k=1

wk g

• Uk, V k

for the numerical solution .

2Pareschi, Russo (2010);

Asher, Ruuth, Spiteri (1997)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 9 / 34

IMEX Runge-Kutta schemes for diffusive relaxation

The ν × ν matrices A = ( aik) and A = (aik) and vectors b, b ∈ Rν characterize the scheme and can be represented by a double tableau in the usual Butcher notation:

• c
• A
• w t

c A w t

Matrix A is lower triangular, i.e. the implicit scheme is a Diagonally Implicit Runge-Kutta (DIRK ) scheme. This choice guarantees that implicit terms are, indeed, always explicitly evaluated.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 10 / 34

Boltzmann-Poisson equation for semiconductors

Let f (t, x, v) be the density distribution function for particles at time t ≥ 0, space point x ∈ Rd and that travel with velocity v ∈ Rd, where d = 1, 2, 3 is the dimension. Distribution function f solves a kinetic equation with diffusive scaling 3 ε ∂tf + v · ∇xf − q mE · ∇vf = 1 εQ(f ) . (3) In this formula: ε is the mean free path, q is the elementary charge and m the effective mass of the electron. E(t, x) = −∇xΦ(t, x) is the electric field. E is obtained from the electric potential Φ given by the solution of a Poisson equation.

3Poupaud (1991);

Klar (1998); Jin, Pareschi (2000)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 11 / 34

Collision term

The anisotropic collision term Q(f ) is defined by Q(f ) =

• σ(v, w)
• M(v)f (w) − M(w)f (v)
• dw,

where σ is the scattering kernel and M is the normalized Maxwellian at the temperature θ of the semiconductor M(v) = 1 (2πθ)d/2 exp

• −|v|2

2θ

• .

The collision frequency λ is defined by: λ(v) =

• σ(v, w)M(w) dw .
• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 12 / 34

Diffusive limit

Define the total mass ρ = ρ(t, x) as ρ =

• f (v) dv .

As ε → 0, one can show that f is approximated by f (t, x, v) ≈ ρ(t, x)M(v) (ε → 0) with ρ satisfying the drift-diffusion equation: ∂tρ = ∇x · (D∇xρ + ηρE) . (4) In this equation, D is the diffusion coefficient defined implicitly in terms of the cross section and constant η is the mobility given by the Einstein relation q D = η m θ.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 13 / 34

Even and Odd parities formalism

Define the even parity r and the odd parity j by: r(t, x, v) = f (t, x, v) + f (t, x, −v) 2 , j(t, x, v) = f (t, x, v) − f (t, x, −v) 2ε . Splitting eq. (3) into two equations, one for v and one for −v and adding and subtracting them we obtain: ∂tr + v · ∇xj − q mE · ∇vj = 1 ε2 Q(r), ∂tj + 1 ε2

• v · ∇xr − q

mE · ∇vr

• = − 1

ε2 λj , (5)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 14 / 34

Even and Odd parities formalism

In the fluid-dynamic limit, i.e. as ε → 0, we obtain: Q(r) = 0 , λj = −v · ∇xr + q mE · ∇vr , (ε = 0) from which r(t, x, v) = ρ(t, x)M(v) and j = 1 λ(v)

• −v · ∇xr + q

mE · ∇vr

• =

1 λ(v)

• −Mv · ∇xρ + q

mρE · ∇vM

• .

Substituting this relations into the first equation in (5) and integrating over v, we obtain the drift-diffusion equation (4).

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 15 / 34

Even and Odd parities formalism

Using parities formalism, our problem now reads: ∂tr + v · ∇xj − E · ∇vj = 1 ε2 Q(r), ∂tj + 1 ε2 (v · ∇xr − E · ∇vr) = − 1 ε2 λj , As was done before, we add and subtract in the first equation the term v · ∇x

• µ v

λ · ∇xr

• ,

where µ = µ(ε) is a real positive parameter, such that µ(0) = 1. The idea behind this choice is that we want to compute such term with the implicit solver when the equation is in the fluid-dynamic regime, in order to use the appropriate solver for the diffusion term.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 16 / 34

IMEX scheme for parity equations

The system we are going to solve now reads as follows: ∂tr + v · ∇x

• j + µ v

λ · ∇xr

• − E · ∇vj
• Explicit

= 1 ε2 Q(r) + v · ∇x

• µ v

λ · ∇xr

• Implicit

, ∂tj + 1 ε2 (v · ∇xr − E · ∇vr)

• Implicit

= − 1 ε2 λ j

Implicit

. Remark: such scheme can be used for easy to invert collision terms Q; this is true, e.g., for the relaxed time approximation (RTA), which is the case when σ ≡ 1, which implies: Q(r) = Mρ − r.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 17 / 34

Hermite approximation of the velocity space

Because of the structure of the problem, it is convenient to write unknowns r and j as follows: r = φM and j = ψM, with φ = φ(t, x, v) and ψ = ψ(t, x, v). Set r = φM, j = ψM, with φ(v) =

N

• k=0

φk Hk(v) , ψ(v) =

N

• k=0

ψk Hk(v) , being the Hermite expansion.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 18 / 34

Hermite approximation of the velocity space

The system then becomes: ∂tφ + v∂x

• ψ + µv

λ∂xφ

• − E
• ∂vψ − 2θvψ
• =

= 1 ε2 Q(φ) + µv 2 λ ∂xxφ , for φ and, similarly for ψ: ∂tψ = − 1 ε2

• λψ + v∂xφ − E∂vφ + 2θvEφ
• .
• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 19 / 34

Hermite approximation of the velocity space

For the collision operator, we have Q(r)(v) = M(v)

N

• j=0

σ(v, vj)φ(vj)wj − λ(v)r(v) , and for the collision frequency λ(v) =

N

• j=0

σ(v, vj)wj , where (vj, wj) are the points and weights of the Gauss-Hermite quadrature rule.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 20 / 34

Hermite approximation of the velocity space

Derivatives with respect to v become easy to compute: ∂vr = M∂vφ − 2vMφ , ∂vj = M∂vψ − 2vMψ , which can be evaluated using ∂vφ =

N

• k=1

φk √ 2k Hk−1(v) =

N

• j=0

φ(vj)cj(v), ∂vψ =

N

• j=0

ψ(vj)cj(v) . for φ and for ψ. Coefficients cj(vi) = cij are given by cj(v) =

N

• k=1

√ 2k Hk(vj) Hk−1(v)wj . and are constant in time.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 21 / 34

Space discretization

General requirements for the space discretization:

1 correct diffusion limit: for an implicit approximation in the limit ε = 0,

require that µ(0) = 1;

2 compact stencil: it is important to obtain a scheme which uses a

compact stencil in the limit ε → 0: this is guaranteed by point (1) and by a suitable discretization for the diffusion term;

3 shock capturing: this is necessary for large values of ε and also for

convection-diffusion type limit equations with small diffusion;

4 avoid solving nonlinear algebraic equations: the implicit space-derivative

in the second equation is, indeed, explicitly evaluated. High-order finite volumes Weighted-Essentially Non Oscillatory (WENO) reconstructions are a good choice in order to achieve all these goals.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 22 / 34

Space discretization

Consider again the prototype system ut + vx = 0, vt + 1 ε2 ux = − 1 ε2 v. In a finite volume approximation, taking into account only the space-derivative discretization, this system reads: ¯ ut = − ˆ vi+1/2 − ˆ vi−1/2 ∆x , ¯ vt = −ˆ ui+1/2 − ˆ ui−1/2 ∆x − 1 ε2 ¯ v,

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 23 / 34

Numerical dissipation

Using, for example, Lax-Friedrichs type fluxes, we have 4: ˆ vi+1/2 = 1 2

• vi+1 + vi
• − 1

ε

• ui+1 − ui
• ,

ˆ ui+1/2 = 1 2

• ui+1 + ui
• − ε
• vi+1 − vi
• .

To be able to use our scheme also in the limit ε → 0, it is necessary to consider a modified flux. We consider modified artificial viscosity constants: ˆ vi+1/2 = 1 2

• vi+1 + vi
• − αv
• ui+1 − ui
• ,

ˆ ui+1/2 = 1 2

• ui+1 + ui
• − αu
• vi+1 − vi
• .

For example, we consider αu = αv = ε.

4Jin, Pareschi, Toscani (2000)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 24 / 34

Choice for µ

For large value of ε, we want to avoid adding and subtracting terms which may cause loss of stability. A first choice is based on the simple remark that when ε ≥ ∆x, we want to keep the explicit scheme; this leads to: µ = µ(ε, ∆x) = 1 if ε < ∆x, if ε ≥ ∆x,

• r some smoothed version of it (e.g. µ = exp(−ε/∆x)).

Another possibility is to relate µ to the time step, with the goal to avoid adding backward diffusion up to O(ε2). For the IMEX Euler scheme, this leads to µ = µ(ε, ∆t) = ∆t ε2 + ∆t . More accurate choices can be done by taking µ = µ(ε, ∆x, ∆t).

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 25 / 34

Boundary conditions

For x ∈ (xL, xR) and v ∈ R, problem (1) is complemented with boundary conditions: for v > 0 f (t, xL, v) = FL(v), f (t, xR, −v) = FR(v), where FL and FR are assigned non-negative functions. To get boundary conditions for r and j, it is possible to use relations 5: for v > 0, r + εj

• x=xL = FL ,

r − εj

• x=xR = FR .

5Jin, Pareschi (2000)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 26 / 34

Boundary conditions

When ε ≪ 1, a reasonable approximation for j is λj = −v∂xr + E∂vr. Applying this to the boundaries: r − ε λ(v∂xr − E∂vr)

• x=xL = FL ,

r + ε λ(v∂xr − E∂vr)

• x=xR = FR .

Remark: from r ± εj = FL,R we have ∂vr = ∂vFL,R + O(ε). Then, up to O(ε2): r − ε λ(v∂xr − E∂vFL)

• x=xL = FL ,

r + ε λ(v∂xr − E∂vFR)

• x=xR = FR .
• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 27 / 34

Properties of the IMEX scheme

General requirements for the time discretization:

1 AP; 2 Globally stiffly accurate.

Examples of schemes which uniformly work for our problem are: second order: SSP2-(3,3,2) scheme and ARS-(2,2,2) scheme; third order: ARS-(4,4,3) scheme and BPR-(3,5,3) scheme 6.

6Boscarino, Pareschi, Russo 2011

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 28 / 34

Numerical results

Figure: Test1 at Tf = 0.03, ε = 0.002

∆tJP = 5 ∗ 10−5, ∆tDPR = 5 ∗ 10−4.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 29 / 34

Numerical results

Figure: Electric field given by a potential well

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 30 / 34

Numerical results

Figure: Test 2 at time Tf = 0.4 ε = 0.001

∆tJP = 10−4, ∆tDPR = 2 ∗ 10−3 ⇒ ∆tJP/∆tDPR = ∆x.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 31 / 34

Penalization technique for the collision term

How do we treat more realistic collision terms? We substitute the full collision term Q adding and subtracting its “linearized” version 7: Q(r)

• Implicit

→

• Q(r) − L(r)
• Explicit

+ L(r)

• Implicit

. For example, if Q(r) = 0 ⇒ r = ρM, then such strategy can be effectively applied using the RTA approximation: L(r) := ρM − r. Remark: the linear part is computed implicitly to stabilize the non-linear collision operator, without changing the asymptotic behavior of the solution.

7Jin, Filbet (2010);

Dimarco, Pareschi (2011)

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 32 / 34

Penalization technique for the collision term

Using such technique, the IMEX scheme for diffusive relaxation can be used also for very general collision operators: ∂tr + v · ∇x

• j + µ v

λ · ∇xr

• − E · ∇vj − 1

ε2

• Q(r) − L(r)
• Explicit

= = 1 ε2 L(r) + v · ∇x

• µ v

λ · ∇xr

• Implicit

, ∂tj + 1 ε2 (v · ∇xr − E · ∇vr)

• Implicit

= − 1 ε2 λ j

Implicit

. Remark: the asymptotic behavior is the same as in the “linear” case.

• V. Rispoli (University of Ferrara)

IMEX-RK schemes for semiconductors Padova, 28.06.2012 33 / 34

THANK YOU!