On Composite Discontinuous Galerkin Method - arXiv

3 downloads 0 Views 297KB Size Report
Mar 4, 2016 - such a manner that E is a triangulation of Ω (figure 2). We will call ... i=1 Thi . For s > 0, we define the broken Sobolev spaces Hs(E) and ...... h,ϕh) + Jr(u∗ − u∗ ..... [1] Randolph E. Bank and Donald J. Rose. ... [19] Peter Wilkes.
arXiv:1603.01167v1 [math.NA] 3 Mar 2016

On Composite Discontinuous Galerkin Method for simulations of electric properties of semiconductor devices Konrad Sakowski∗ 1,3 , Leszek Marcinkowski3 , Pawel Strak1 , Pawel Kempisty1 , and Stanislaw Krukowski1,2 1

Institute of High Pressure Physics, Polish Academy of Sciences, ul. Sokolowska 29/37, 01-142 Warsaw, Poland 2 Interdisciplinary Centre for Materials Modeling, Warsaw University, ul. Pawińskiego 5a, 02-106 Warsaw, Poland 3 Department of Mathematics, Computer Science and Mechanics, Warsaw University, ul. Banacha 2, 02-097 Warsaw, Poland March 4, 2016

Abstract In this paper two variants of discretization of the van Roosbroeck’s equations in equilibrium state with the Composite Discontinuous Galerkin Method for rectangular domain are discussed. They base on Weakly Over-Penalized Symmertic Interior Penalty (WOPSIP) method and on Symmetric Interior Penalty Galerkin (SIPG) method. It is shown that the discrete problems are well-defined and that their solutions are unique. Error estimates are derived. Finally numerical simulations of gallium nitride semiconductor devices are presented.

1

Introduction

Numerical simulations are the important tool in development of semiconductor devices. Since our contemporary electronics relies on the semiconductors, there is strong demand on the progress in this domain. Examples of such devices are light emitting diodes, lasers, transistors, detectors, and many others. There are various approaches in simulations of such devices, depending on precision, efficiency and size of a simulated fragment. On the one hand there are so-called ab initio methods, which are used to investigate properties of elements composed of hundreds of thousands of atoms. These methods use fundamental laws of physics and they need days or weeks to to perform a single simulation on a computational cluster. On the other hand there is a drift-diffusion theory. ∗ [email protected]

1

In this case the model is much simpler and it allows to simulate whole semiconductor device on a standard desktop computer. This model describes two kinds of carriers, electrons and holes, which move in the electric field present in semiconductor devices. From the mathematical point of view, it consists of a system of three nonlinear elliptic differential equations, which are called van Roosbroeck equations [14]. Numerical modelling of semiconductor devices with the drift-diffusion model are performed since 1964, when Gummel [8] proposed a numerical algorithm based on the simple iteration method. Various methods were used for discretization of the van Roosbroeck equations, for example Finite Difference Method [17], Box method [1], Finite Element Method (FEM) [6]. Special variants of discretizations optimized for the so-called continuity equations were developed [11]. In our study, we are interested in laser and electroluminescent diodes based on gallium nitride. These devices have particular structure. They are composed of layers deposited one on another. These layers vary in material composition, doping level, dopant type or length, and they form a natural partition of a device. They are usually cuboids (see figure 1). In general, in these devices there are three important parts: two contacts and the active region. A potential can be applied to the contacts and then some current will flow through the device. Then, in the active region, the recombination occurs and fotons are generated. Due to such a structure, it is often sufficient to use one-dimensional model for simulations of the laser diodes and the electroluminescent diodes, because unless a device has some special design, the current flows more or less in the straight line between the contacts. However, to take into account more sophisticated designs, more spatial domensions have to be used. Motivated by this problem, in this paper we study the Composite Discontinuous Galerkin Method discretizations. In these methods, first we pick some partition {Ωi }N i=1 of a domain Ω. Then on every element Ωi we introduce some triangulation and a Finite Element space Xhi (Ωi ). These triangulations may be independent of each other, i.e. they do not need to agree on interfaces ∂Ωi ∩∂Ωj . Finally we search for a discrete solution, which is in the space Xhi (Ωi ) on every Ωi . On the interfaces we use a variant of the Discontinuous Galerkin Method (DGM) to “glue” functions from the different partition elements together. Our interest in application of these methods to modelling of the semiconductor luminescent devices is caused by the following reasons. First, the natural division of these devices described above is a good candidate for a partition of a computational domain Ω. It is so because the discontinuities of the coefficients are often on interfaces between layers due to variations of material parameters. Also it is often the case that one could easily predict on which layers the variations of the unknown functions are more frequent than on the others, due to physical properties of the device. Therefore it convenient to use separate meshes, as there is no problem in mesh refinement in arbitrary layers. We start with introduction of the differential problem in section 2. We propose two variants of Composite DGM discretization of this problem in section 3. Then we show existence and uniqueness of the introduced discrete problems in sections 4, 5. In section 6 we discuss interpolation properties of the discrete space. Then we pass to the error estimates in section 7. Finally we present results of numerical simulations in section 8 and we conclude in section 9.

2

p-Al0.3GaN (10 nm)

p-GaN

(50 nm)

p-GaN/p-Al0.16GaN (26 Å / 26 Å) * 80 SL

p-GaN

(70 nm)

In0.04GaN:Si (8 nm) n-GaN

In10%GaN /In4%GaN:Mg (45 Å / 80 Å) * 5 MQW

(140 nm)

n-GaN/Al0.16GaN (29 Å / 29Å )*110 SL

GaN:Si (530 nm)

bulk n-GaN

Figure 1: Example of a gallium nitride semiconductor laser structure.

2

Differential problem

The drift-diffusion model describes relationship between the electrostatic potential and the charge carrier concentrations: electrons and holes [19, 18]. Physical derivation of this model is beyond the scope of this work, therefore we will focus on the mathematical standpoint. We start with the domain Ω of our problem. Luminescent semiconductor devices are generally made of planar layers deposited one on another, which vary in composition of a semiconductor material or number of impurities (see figure 1). At opposite ends metal contacts are attached, where the current can be applied. If this is the case, it flows through the device perpendicular to the deposited layers. We assume that Ω is a rectangle. In this paper we deal with so-called equilibrium state. It is a physical state of a semiconductor device, where the device is disconnected from power sources, so the current does not flow through it. Then only the Poisson equation needs to be solved and the system reduce to the following problem: find u∗ ∈ H 1 (Ω), such that   ∗ ∗ ˆ = k1 , −∇ ε(x)∇u∗ + eu −ˆv − ew−u u∗ = uˆ on ΣD , ∇u · ν = 0 on ΣN .

(1)



where vˆ = w ˆ ≡ const. Since some results of this paper may be also applied to non-equilibrium case, we consider more general assumption that vˆ, w ˆ ∈ L∞ (Ω). Also we assume that ε, u ˆ ∈ H 1 (Ω) ∩ L∞ (Ω) and 0 < εm ≤ ε ≤ εM , εm , εM ∈ R. 3

7

8

9

10

11

12

1

2

3

4

5

6

Figure 2: An example of two-dimensional coarse grid of Ω. The following lemma is essential for the results presented in this paper. Its proof may be found in [9]. Lemma 2.1. Solution u∗ of problem (1) is bounded. A weak formulation of the differential problem (1) is as follows. Find u∗ ∈ u ˆ + H 1 (Ω), such that a(u∗ , ϕ) + b(u∗ , ϕ) = f (ϕ)

1 (Ω), ∀ϕ ∈ H0,Σ D

(2)

where a(u, ϕ) :=

Z

ε(x)∇u(x) · ∇ϕ(x) dx



b(u, ϕ) := f (ϕ) :=

Z 

ZΩ

 ˆ ϕ(x) dx. eu(x)−ˆv(x) − ew(x)−u(x)

(3)

k1 (x)ϕ(x) dx.



We use the following notation: ∞ (Ω) := {f ∈ C ∞ (Ω) : f |ΣD ≡ 0}, C0,Σ D   ∞ 1 (Ω) . (Ω) := clH 1 (Ω) C0,Σ H0,Σ D D

3 3.1

(4)

Discretization Discrete space

Let Ω ⊂ R2 be a rectangle, divided to disjoint subrectangles {Ωi }N i=1 =: E in such a manner that E is a triangulation of Ω (figure 2). We will call this division a coarse grid and we assume that each edge of any Ωi is contained either in ΣD or in ΣN . Let us define triangulations Thi := Ti,hi (Ωi ), where hi := max{diam(τ ) : τ ∈ Thi }. By Nhi we denote the nodes of the triangulation Thi . We assume that {Ti,hi (Ω)}hi is a regular uniform family of triangulations [5]. We will define 4

S s Th := N i=1 Thi . For s > 0, we define the broken Sobolev spaces H (E) and H s (Th ) as H s (E) := {v ∈ L2 (Ω) : ∀i ∈ {1, . . . , N } vi := v|Ωi ∈ H s (Ωi )} ⊂ L2 (Ω), H s (Th ) := {v ∈ L2 (Ω) : ∀τ ∈ Th v|τ ∈ H s (τ )} ⊂ L2 (Ω).

(5)

Then on every Ωi we define a discrete space Xhi (Ωi ) of piecewise linear functions on the triangulation Thi : n o Xhi := Xhi (Ωi ) := uh,i ∈ C(Ωi ) : ∀τ ∈ Thi uh,i τ ∈ P1 (τ ) (6) Finally we define Xh (Ω) as o n Xh (Ω) = (uh,1 , . . . , uh,N ) : uh,i ∈ Xhi (Ωi ), i ∈ {1, . . . , N } ⊂ L2 (Ω).

(7)

Note that Xh (Ω) 6⊂ H 1 (Ω) and Xh (Ω) 6⊂ H 2 (E), but Xh (Ω) ⊂ H 1 (E), H 1 (Ω) ⊂ H 1 (E) and Xh (Ω) ⊂ H 2 (Th ). By Γ we denote a set of all internal and boundary edges of E. Then Γ is a sum of disjoint sets ΓD , ΓN and ΓI , where ΓD := {e ∈ Γ : e ⊂ ΣD }, ΓN := {e ∈ Γ : e ⊂ ΣN }, ΓI := {e ∈ Γ : e ⊂ int(Ω)}.

(8)

Therefore ΓD (resp. ΓN ) contains edges lying on the boundary where Dirichlet (resp. Neumann) boundary conditions are imposed and in ΓI there are all internal edges, which we call interfaces, as they frequently correspond to the physical interfaces between different semiconductor materials. We also define ΓDI := ΓD ∪ ΓI ,

Γi := {e ∈ Γ : e ⊂ ∂Ωi }.

(9)

Let e ∈ Γ. Then two cases are possible. Either e ∈ ΓD ∪ ΓN , so there is an unique Ωi such that e is an edge of Ωi , or e ∈ ΓI and there are exactly two sets Ωi , Ωj ∈ E such that e is their common edge. We will often refer to these two cases. Also we define nb(Ωi ) := {Ωk ∈ E : Γi ∩ Γk 6= ∅}. For s > 1/2 we define operators [·] := [·]e : H s (E) → L2 (e), {·} := {·}e : s H (E) → L2 (e) as ( ui − uj if e ⊂ ΓI , e = ∂Ωi ∩ ∂Ωj , [u] := ui if e ⊂ ΓD ∪ ΓN , e = ∂Ωi ∩ ∂Ω, ( (10)  1 if e ⊂ ΓI , e = ∂Ωi ∩ ∂Ωj , ui + uj 2 {u} := ui if e ⊂ ΓD ∪ ΓN , e = ∂Ωi ∩ ∂Ω. For convenience, we will also use this notion for triangulation parameters, i.e.     1 1s + 1s n1o if e = ∂Ωi ∩ ∂Ωj , hj (11) := 21 hi {h−s } := s  s h if e = ∂Ωi ∩ ∂Ω. h i

5

For further analysis, we introduce so-called “broken” norm k · kh,σr in Xh (Ω) with Z N Z 2  X X ε ∇uh,i dx + ηr,e [uh ]2 ds. kuh k2h,σr := (12) i=1

Ωi

e∈ΓDI

e

Here ηr,e is a penalty coefficient for e. It depends on the triangulation parameters and penalty parameters σe > 0: ( −r 2σ e hi −r  e ∈ Γ D , e ⊂ Ωi , (13) ηr,e := 2σe {h } = −r −r e ∈ Γ I , e ⊂ Ωi ∩ Ωj . σe hi + hj

Depending on the method, we will use r = 1 or r = 2. To simplify the analysis, we assume that 0 < σ0 ≤ σe for all e ∈ ΓDI . Also we assume that 0 < hi < h0 ≤ 1 for all i ∈ {1, . . . , N }. The choice of σ0 and h0 will be discussed later in lemmata 3.2, 3.3 and 7.2. Lemma 3.1. For any uh ∈ Xh (Ω), Ωi ∈ E and e ∈ Γi , the following estimates hold −1/2

kuh kL2 (Ωi ) ,

(14)

−1/2

|uh |H 1 (Ωi ) .

(15)

kuh,i kL2 (e)



Chi

k∇uh,i · νkL2 (e)



Chi

Constant C does not depend on hi . Proof. These estimates are a consequence of the trace theorem applied to each edge of fine elements in Ωi coincident with e followed by a scalling argument.

3.2

Composite Discontinuous Galerkin variants

We propose two variants of the Composite Discontinuous Galerkin discretization. First one is based on Weakly Over-Penalized Symmetric Interior Penalty (WOPSIP) method (cf. [3] or [2]). Second formulation is derived from Symmetric Interior Penalty Galerkin (SIPG) method (cf. [13] or [12]). In each case we use the composite formulation (cf. [7]), i.e. inside every Ωi we use the Finite Element Method on the triangulation Thi , while on boundaries e ∈ ΓDI we use the respective variant of the Discontinuous Galerkin Method. Thus we call these methods Composite Weakly Over-Penalized Symmetric Interior Penalty (CWOPSIP) method and Composite Symmetric Interior Penalty Galerkin (CSIPG) method, respectively. 3.2.1

Composite Weakly Over-Penalized Symmetric Interior Penalty (CWOPSIP)

This discretization has a simpler formulation from the two methods we introduce. The discrete problem is defined as follows. Find u∗h ∈ Xh (Ω) such that a2,h (u∗h , ϕh ) + b(u∗h , ϕh ) = f2,h (ϕh ),

6

∀ϕh ∈ Xh (Ω),

(16)

where a2,h (uh , ϕh ) =

N Z X

Ωi

i=1

f2,h (ϕh ) =

Z

ε∇u∗h,i

· ∇ϕh,i dx +

X

η2,e

k1 ϕh dx +



η2,e

Z

[uh ] · [ϕh ] ds,

e

e∈ΓDI

X

Z

(17)

[ˆ u] · [ϕh ] ds,

e

e∈ΓD

and b is defined as in (3). 3.2.2

Composite Symmetric Interior Penalty Galerkin (CSIPG)

This problem is defined as follows. Find u∗h ∈ Xh (Ω) such that a1,h (u∗h , ϕh ) + b(u∗h , ϕh ) = f1,h (u∗h , ϕh ),

∀ϕh ∈ Xh (Ω),

(18)

where a1,h (uh , ϕh ) =

N Z X i=1

Ωi

ε∇u∗h,i

e∈ΓDI

X Z



e∈ΓDI

f1,h (ϕh ) =

Z

{ε∇ϕh · ν}[uh ] ds +

e



X Z

e∈ΓDI

X

η1,e

e∈ΓD

{ε∇uh · ν}[ϕh ] ds

e

X

η1,e

e∈ΓDI

k1 ϕh dx −

+

· ∇ϕh,i dx −

X Z

Z

Z

[uh ] · [ϕh ] ds,

e

(19)

{ε∇ϕh · ν}[ˆ u] ds

e

[ˆ u][ϕh ] ds,

e

and b is defined as in (3). Lemma 3.2. For any α ∈ (0, 1) there exist σ0 > 0, such that for every σe ≥ σ0 and uh ∈ Xh (Ω) X Z 2 {ε∇uh · ν}[uh ] ds ≤ αkuh k2h,σ1 , (20) e∈ΓDI

e

Proof. Let us take any e = Ωj ∩ Ωk . By the Schwarz inequality and a triangle inequality Z 

2 {ε∇uh · ν}[uh ] ds ≤ ε|Ωj ∇uh,j · ν L2 (e) e (21) 



[uh ] + ε|Ω ∇uh,k · ν k

L2 (e)

L2 (e)

Taking Ωi ∈ {Ωj , Ωk }, we use the Cauchy’s ǫ-inequality and lemma 3.1



[uh ]

ε|Ω ∇uh,i · ν i L2 (e) L2 (e)

2

2

1 1 ≤ ǫεM hi ∇uh,E L2 (e) + εM [uh ] L2 (e) 4ǫ hi

2

2 1 1 ≤ ǫεM ∇uh L2 (Ωi ) + εM [uh ] L2 (e) . 4ǫ hi 7

(22)

Therefore we get Z

2

2

1 2 {ε∇uh · ν}[uh ] ds ≤ 2εM ǫ ∇uh L2 (Ωj ) + εM {h−1 } [uh ] L2 (e) . 2ǫ e

(23)

On the other hand, if e ∈ ΓD then e ∈ ∂Ωi ∩ ∂Ω and by similar arguments we have Z

2

2 1 (24) 2 {ε∇uh · ν}[uh ] ds ≤ 2εM ǫ ∇uh L2 (Ωi ) + εM {h−1 } [uh ] L2 (e) ǫ e Summing these results up we get 2

N X Z X k∇uh k2L2 (Ωi ) {ε∇u · ν}[u ] ds ≤ ǫC h h e

e∈ΓDI

i=1

2 C X + {h−1 } [uh ] L2 (e) , ǫ

(25)

e∈ΓDI

where C depends on εM and the maximal number of edges of elements in coarse grid E. Finally for any α ∈ (0, 1) if we take ǫ := αεm /C and then taking σ0 (α) := C 2 /α2 εm we have 2

N X Z X εm k∇uh k2L2 (Ωi ) {ε∇uh · ν}[uh ] ds ≤ α

e∈ΓDI

e

i=1



X

e∈ΓDI

2 C2 {h−1 } [uh ] L2 (e) 2 εm α

(26)

≤ αkuh k2h,σ1 .

Lemma 3.3. There exist σ0 > 0 and c > 0, such that for every σe ≥ σ0 and uh ∈ Xh (Ω) ckuh k2h ≤ ah,1 (uh , uh ). (27) Proof. Due to definition of broken norm, we have X Z 2 ah,1 (uh , uh ) = kuh kh,σ1 − 2 {ε∇uh · ν}[uh ]ds. e∈ΓDI

Let us take α := 1/2 in lemma 3.2. Then we have X Z 2 ah,1 (uh , uh ) = kuh kh,σ1 − 2 {ε∇uh · ν}[uh ]ds e∈ΓDI



kuh k2h,σ1

e

1 1 − kuh k2h,σ1 ≥ kuh k2h,σ1 . 2 2

8

(28)

e

(29)

4

Existence

Let us denote the nodal base of Xh by {ϕ(j) }Ji=1 . Then we define P : Xh (Ω) → Xh∗ (Ω) as    (30) P (uh )ϕh := ar,h uh , ϕh + b uh , ϕh − fr,h ϕh .

We would like to use the following theorem [10]:

Theorem 4.1. (Brouwer) Let P : X → X ∗ be a continuous function on a finite-dimensional normed real vector space X, such that for suitable ρ > 0 we have P (x)x ≥ 0 ∀kxk ≥ ρ, (31) Then there exists x ∈ X, kxk ≤ ρ such that

4.1

P (x) = 0.

(32)

   P (uh )uh := a2,h uh , uh + b uh , uh − f2,h uh ,

(33)

 a2,h uh , uh = kuh k2h,σ2 .

(34)

Case of CWOPSIP

By definition (30)

Then we have that

On the other hand, by Schwarz inequality  u, k1 )kuh kL2 (Ω) ≥ −c(ˆ u, k1 , h)kuh kh,σ2 −f2,h uh ≥ −c(ˆ

(35)

Then let C := max{kˆ v kL∞ (Ω) , kwk ˆ L∞ (Ω) }. Then we may decompose Φ2 as Z    ˆ h uh dx euh −ˆv − ew−u b uh , uh = ZΩ   ˆ h uh χ{x∈Ω:|uh (x)|>C} dx = euh −ˆv − ew−u (36) Ω Z   ˆ h uh χ{x∈Ω:|uh (x)|≤C} dx + euh −ˆv − ew−u Ω

The first integral is non-negative, and the latter we can estimate from below Z   ˆ h (x) uh (x)χ{x∈Ω:|uh (x)|≤C} (x)dx ≥ −|Ω|2e2C C. (37) euh (x)−ˆv(x) − ew(x)−u Ω

In conclusion, we may use these estimations to obtain P (uh )uh ≥ kuh k2h,σ2 − c1 kuh kh,σ2 − c2

(38)

Note that constants ci in this inequality depend on h. It is therefore clear that for kuh kh,σ2 large enough, we have that P (uh )uh ≥ 0. Then by theorem 4.1 we have that there exists some u∗h , such that P (u∗h ) = 0.

9

4.2

Case of CSIPG

We proceed analogously to the CWOPSIP case. For b(uh , uh ) the argumentation is exactly the same. Then f1,h (uh ) has one additional term, which may be estimated using lemma 3.1 ant the trace inequality X Z X {ε∇uh · ν}[ˆ u] ds ≤ εM k∇uh kL2 (e) kˆ ukL2 (e) e∈ΓD

e

e∈ΓD

≤ cεM

N X

−1/2

hi

k∇uh kL2 (Ωi ) k[ˆ u]kL2 (Ωi )

(39)

i=1

ukL2 (Ω) , ≤ Ckuh kh,σ1 kˆ where C depends on εM and h. Then estimating a1,h (uh , uh ) by lemma (3.3), we have P (uh )uh ≥ c1 kuh k2h,σ1 − c2 kuh kh,σ1 − c3

(40)

The existence of u∗h is now proven.

5

Uniqueness

The uniqueness can be shown by contradiction for both cases. Here we present the CWOPSIP case, for CSIPG the proof is similar. Assume that there exist two solutions uh , vh ∈ Xh (Ω) of equation (18). Then by taking ϕh := uh − vh in (18) and subtracting these equations for uh and vh we obtain N Z    X e−ˆv evh − euh uh − vh dx a2,h (uh − vh , uh − vh ) = Ωi

i=1

+

N Z X i=1

Ωi

(41)

e

w ˆ



e

−uh

−e

−vh

  uh − vh dx.

By monotonicity of the exponential function, the right hand side is nonpositive. On the other hand 0 < kuh − vh k2h,σ2 = a2,h (uh − vh , uh − vh ).

(42)

Thus 0 < kuh − vh k2h,σ2 ≤ 0 since uh 6= vh , and we have a contradiction.

6

Interpolation operator

First let us take any Ωi ∈ E and let us define interpolation operator Πhi : H 2 (Ωi ) → Xhi ⊂ C 0 (Ωi ) as follows ∀x ∈ Nhi

Πhi u(x) = u(x).

(43)

Note that for s > 0 we have H 1+s (Ω) ⊂ C 0 (Ω) (see [13]), so this definition is well-posed. Then we define Πh : H 2 (E) → Xh by   ∀Ωi ∈ E Πh u := Πhi ui . (44) i

10

On any Ωi , we can use standard interpolation estimate for FEM [5] kui − Πhi ui kL2 (Ωi ) + hi kui − Πhi ui kH 1 (Ωi ) ≤ Ch2i |ui |H 2 (Ωi ) .

(45)

We remind that coarse grid E is independent of h. Thus ku − Πh uk2H 1 (Ω) =

N X

kui − Πhi ui k2H 1 (Ωi ) ≤ C

i=1

X

h2i |ui |2H 2 (Ωi ) .

(46)

Ωi ∈E

Let further uI := Πh u. Lemma 6.1. Let u ∈ H 2 (E), uI := Πh u. For any Ωi ∈ E and for any e ∈ Γi

ui − uI,i L2 (e) ui − uI,i 1 H (e)

3/2

(47)

1/2

(48)



Chi |ui |H 2 (Ωi ) ,



Chi |ui |H 2 (Ωi ) .

P Proof. For fixed e and Ωi , we have kui −uI,i k2L2 (e) = τ ∈Th ,e kui −uI,i k2L2 (e∩τ ) . i Note that on a single triangulation element τ we have that ui − uI,i ∈ H 2 (τ ), so using the trace inequality (see [13]) for H 2 (τ ) functions we have   1/2 −1/2 kui − uI,i kL2 (τ ) + hi ui − uI,i H1 (τ ) . (49) kui − uI,i kL2 (e∩τ ) = C hi

Then by (45) it follows that

−1/2

kui − uI,i kL2 (e) ≤ Chi

3/2

 

kui − uI,i kL2 (Ωi ) + hi ui − uI,i H 1 (Ωi )

(50)

≤ Chi |ui |H 2 (Ωi ) .

Proof of the latter estimate is analogous. Let us take any e ∈ ΓDI . For e ∈ ΓI we assume that e = Ωj ∩ Ωk for some Ωj , Ωk ∈ E and by a triangle inequality we have Z

e

[u − uI ]2 ds = k[u − uI ]k2L2 (e) ≤ 2

2 X

kui − uI,i k2L2 (e) ,

(51)

i∈{j,k}

while for e ∈ ΓD we have e ∈ Γi for some Ωi ∈ E and simply Z  Z 2 2 ui − uI,i ds = kui − uI,i k2L2 (e) [u − uI ] ds = e

(52)

e

Therefore it is sufficient to estimate kui − uI,i k2L2 (e) , for any e ∈ ΓDI , e ⊂ ∂Ωi , as in lemma 6.1 To proceed further, we need to take into account ηr,e . Therefore we will distinguish few cases.

11

In this case η1,e = σe {h−1 }. Let e ∈ ΓD . By (50) we have Z  2 2 η1,e ui − uI,i ds = σe h−1 i kui − uI,i kL2 (e)

CSIPG

e

(53)

3 2 2 2 ≤ Cσe h−1 i hi |ui |H 2 (Ωi ) = Cσe hi |ui |H 2 (Ωi ) .

On the other hand, if e ∈ ΓI then Z  2 −1 2 η1,e uj − uI,j ds = 0.5σe (h−1 j + hk )kuj − uI,j kL2 (e) e

Then if we sum up over e ∈ ΓDI X

e∈ΓDI

 h3j  |uj |2H 2 (Ωj ) . ≤ Cσe h2j + hk

Z  N 2 X C h2i + [u − uI ] ds ≤ η1,e e

i=1

X

Ωk ∈nb(Ωi )

! h3i |ui |2H 2 (Ωi ) . hk

(54)

(55)

Thus taking into account this estimate and (46) ku − uI k2h,σ1 ≤ C

N X

X

h2i +

i=1

Ωk ∈nb(Ωi )

! h3i |ui |2H 2 (Ωi ) . hk

(56)

If we increase destiny poportionally, i.e. hi := ci h, the result can be improved to ku − uI k2h,σ1 ≤ Ch2

N X

|ui |2H 2 (Ωi ) .

(57)

i=1

CWOPSIP - general case Here we have η2,e = σe {h−2 }. Proceeding as in previous case, we get for e ∈ ΓD Z  2 ui − uI,i ds ≤ Cσe hi |ui |2H 2 (Ωi ) , η2,e (58) e

and for e ∈ ΓI

η2,e

Z 

uj − uI,j

e

2

 h3j  ds ≤ Cσe hj + 2 |uj |2H 2 (Ωj ) . hk

(59)

So taking into account (46) ku −

uI k2h,σ2

≤C

N X

hi +

i=1

X

Ωk ∈nb(Ωi )

! h3i |ui |2H 2 (Ωi ) . h2k

(60)

Again assuming hi := ci h we obtain ku −

uI k2h,σ2

≤ Ch

N X

|ui |2H 2 (Ωi ) .

(61)

i=1

We must note that this estimate is generally worse than for CSIPG method. 12

Figure 3: An example of a mesh of Vh (Ω) space. Fine grid (mesh of Xh (Ω)) lines are in light gray, coarse grid (E) lines are in dark gray and the black color is used for mesh of Vh (Ω). Three levels of refinement are presented. CWOPSIP - special case Imposing some restrictions on the mesh and boundary conditions, we can improve the interpolation error estimate. Assume that for every space Xh (Ω) there is some continuous FEM space Vh (Ω) ⊂ C(Ω), such that Vh (Ω) ⊂ Xh (Ω). An example is presented on figure 3. Under such assumptions, we define the projection operator Πh : H 2 (E) → Vh (Ω) ⊂ Xh (Ω) to satisfy Πh u(x) = u(x) for all nodal points of discretization Vh . Then Πh u = uI is continuous and [uI ]|e ≡ 0 for every e ∈ ΓI . Also we have [u]|e ≡ 0 in this case as u ∈ H 1 (Ω) and Z  Z  2 2 X X η2,e [u − uI ] ds = η2,e (62) [ˆ u−u ˆI ] ds. e

e∈ΓDI

e

e∈ΓD

This integral do not vanish in general. However we additionally assume that u ˆ|e is constant on any e ∈ ΓD . In this paper, we deal with electrostatic potential, which do not vary within a given contact of a device. Thus by choosing the division in such a manner that boundary edges e ∈ ΓD do not cross the contact boundaries, we can satisfy this assumption. Then u ˆ|e = uI |e and simply by (46) ku − uI k2h,σ2 =

N Z X i=1

≤ Ch

2

Ωi

N 2  X h2i |ui |H 2 (Ωi ) εi ∇ u − uI dx ≤ C i=1

N X

(63)

|ui |H 2 (Ωi ) .

i=1

7

Error estimates

7.1

Definitions and additional assumptions

We start with auxiliary lemma. Lemma 7.1. Let u ∈ H s (E), s ≥ 1. Then kuk2L2 (Ω)

≤C

N Z hX i=1

Ωi

∇u

2

dx +

X

e∈ΓI

13

|e|

−1

Z

e

2

[u] ds +

X Z

e∈ΓD

e

i u2 ds .

(64)

Proof of lemma 7.1 may be found in [4]. Then we would like to have an analogue of a Poincare inequality for the H 1 (E) space. Lemma 7.2. Let u ∈ H s (E), s ≥ 1, r ∈ {1, 2}. Then there exists some h0 > 0, such that kukL2(Ω) ≤ ckukh,σr for 0 < h ≤ h0 , where c is independent of h. Proof. By definition of the broken norm (12), we have kuk2h,σr

:=

N Z X i=1



εi ∇u

Ωi

2

dx +

X

ηr,e

Z

[u]2 ds.

(65)

e

e∈ΓDI

Note that |e| do not depend on h and ηr,e → ∞ as h → 0. Thus we can find h0 > 0, such that ηr,e ≥ |e|−1 and ηr,e ≥ 1 for any 0 < h < h0 and then by lemma 7.1 N Z hX i X X 2 2 kukL2(Ω) ≤ C ∇u dx + |e|−1 k[u]k2L2 (e) + ku]k2L2(e) i=1

h

≤ C ε−1 m ≤

Ωi

e∈ΓI

N Z X i=1

εi ∇u

Ωi

C1 kuk2h,σr .

2

e∈ΓD

dx +

X

ηr,e

e∈ΓDI

Z

e

[u]2 ds

(66)

i

To prove error estimates of the proposed discretizaion, we would like to introduce the following assumptions:  u∗ ∈ H 1 (Ω) ∩ H 2 (E), ε ∈ v ∈ L∞ (Ω) : ∀i ∈ {1, . . . , n} v|Ωi ∈ C 1 (Ω) . (67)

7.2

Consistency

We start with an abstract result. Let u ∈ H 1 (Ω), f ∈ L2 (Ω). We pose two problems. The first is the following: find u ∈ H 1 (Ω) such that Z Z 1 (Ω), f ϕ dx ∀ϕ ∈ H0,Σ ε∇u · ∇ϕ dx = D (68) Ω Ω u = uˆ on ΣD . Second problem is posed in “broken” Sobolev space: find u ∈ H 1 (E), such that ∀ϕ ∈ H 1 (E) ∩ H 2 (Th ) N Z X i=1

Ωi

o X Z n ε∇u · ∇ϕ dx − ε∇u · ν [ϕ] ds e∈ΓDI

e

Z o X X Z n ε∇ϕ · ν [u] ds + +ξ ηr,e [u][ϕ] ds e∈ΓDI

=

N Z X i=1

Ωi

e

e∈ΓDI

e

Z o X Z n X f ϕ dx + ξ ε∇ϕ · ν [ˆ u] ds + ηr,e [ˆ u][ϕ] ds. e∈ΓD

e

e∈ΓD

e

(69)

14

We consider ξ ∈ {0, −1}, which will be used in CWOPSIP and CSIPG, respectively. We would like to prove the following result. Theorem 7.3. Assume that the solution u of problem (68) belongs to H 1 (Ω) ∩ H 2 (E) and ε∇u ∈ H 1 (E). Then u satisfies (69). Conversely, if u ∈ H 2 (E) ∩ H 1 (Ω) is a solution of (69) and ε∇u ∈ H 1 (E), then it is also a solution of (68). The proof is based on the standard approach in Discontinuous Galerkin Method [13]. 2 Lemma 7.4. Let u ∈ H 1 (Ω) ∩ H 2 (E), ε∇u ∈ H 1 (E) , 0 < εm ≤ ε ≤ εM and f ∈ L2 (Ω). The following statements are equivalent: • u satisfy:

Z

ε∇u · ∇ϕ =

Z

1 (Ω) ∀ϕ ∈ H0,Σ D

f ϕ,





(70)

• u satisfy: −

N Z X i=1

Ωi

Z   ∇ εi ∇ui ϕi = f ϕ,

∀ϕ ∈ L2 (Ω)



h

i ε∇u · ν = 0

∇u · ν = 0

on ΣI ,

(71)

on ΣN .

Proof. (71) ⇒ (70) follows simply from the Green formula. To prove (70) ⇒ 1 (71), take any ϕ ∈ C0∞ (Ω). Since C0∞ (Ω) ⊂ H0,Σ (Ω), then by (70) we have D Z Z f ϕ. (72) ε∇u · ∇ϕ = Ω



By the Green formula Z N Z X f ϕ dx = Ω

i=1

=−

ε∇u · ∇ϕ dx Ωi

N Z X i=1

Ωi

(73)   XZ ∇ ε∇u ϕ dx + [ε∇u · ν]ϕ ds e

e∈Γ

Since ϕ is zero on ∂Ω, we may rewrite last sum Z N Z   XZ X ∇ ε∇u ϕ dx + fϕ = − [ε∇u · ν]ϕ ds. Ω

i=1

Ωi

e∈ΓI

(74)

e

Then,h as a simple consequence of the Green formula, for every e ∈ ΓI ∪ ΓN we i have ε∇u · ν = 0, so −

N Z X i=1

Ωi

Z   f ϕ. ∇ ε∇u ϕ dx =

(75)



  This statement is true for ϕ ∈ C0∞ (Ω). Since ∇ ε∇u ∈ L2 (Ω) it is also true for any ϕ ∈ L2 (Ω) and first statement of (71) is shown. 15

Proof. (Theorem 7.3) First we prove (68) ⇒ (69). Assume that u is a solution of (68) and that it belongs to H 1 (Ω) ∩ H 2 (E). We have by definition Z Z 1 (Ω). ε∇u · ∇φ dx = f φ dx ∀φ ∈ H0,Σ (76) D Ω



We use lemma 7.4 and we obtain that for any φ ∈ L2 (Ω) Z Z   f φ dx. ∇ · ε∇u φ dx = −

(77)





Let us take any ϕ ∈ H 1 (E) ∩ H 2 (Th ) and substitute φ := ϕ. We may split integrals to N Z N Z   X X f ϕ dx. (78) ∇ · ε∇u ϕ dx = − i=1

Ωi

i=1

By the Green theorem, we have Z Z   ∇ · ε∇u ϕ dx =

ε∇u · ∇ϕ dx −

Z

ε∇u · νϕ dx.

(79)

∂Ωi

Ωi

Ωi

Ωi

Summing up these results in Ωi , we get N Z X

ε∇u · ∇ϕ dx −

Ωi

i=1

N Z X

ε∇u · νϕ dx =

∂Ωi

i=1

N Z X i=1

f ϕ dx.

(80)

Ωi

By lemma 7.4, we have that [ε∇u] = 0 on ΣI , thus {ε∇u · ν} = ε∇u · ν on any ∂Ωi and we have N Z X i=1

Ωi

N Z o XZ n X ε∇u · ∇ϕ dx − ε∇u · ν [ϕ] dx = e

e∈Γ

f ϕ dx.

(81)

Ωi

i=1

By homogeneous Neumann boundary condition (lemma 7.4) on e ∈ ΓN we have {ε∇u · ν} = 0 and N Z X i=1

ε∇u · ∇ϕ dx −

Ωi

N Z o X X Z n ε∇u · ν [ϕ] dx = e

e∈ΓDI

i=1

f ϕ dx.

(82)

Ωi

Since u ∈ H 1 (Ω), then [u] = 0 for any e ∈ ΓI and by assumption on e ∈ ΓD we have u = u ˆ so we have for any ϕ ∈ H 1 (E)(Ω) Z o X X Z n ηr,e [u][ϕ] ds + ξ ε∇ϕ · ν [u] ds e∈ΓDI

e

e∈ΓDI

=

X

e∈ΓD

e

ηr,e

Z

[ˆ u][ϕ] ds + ξ

e

X Z n

e∈ΓD

e

By adding this result side-by-side to (82) we obtain (69).

16

(83) o ε∇ϕ · ν [ˆ u] ds.

We proceed to (69) ⇒ (68). Assume (69) is true. First we recover the Dirichlet boundary conditions. Take any e ∈ ΓD , such that e ∈ ∂Ωi , and ϕ¯ ∈ C0∞ (e). Then let {ϕǫ }ǫ be a sequence of functions, such that ϕǫ ∈ C ∞ (Ω),

ϕǫ |e = ϕ, ¯ supp(ϕǫ ) ⊂ Ωi ∪ e, ∇ϕǫ · ν = 0, kϕǫ kL2 (Ω) −−−→ 0.

ϕǫ |∂Ωi \e ≡ 0,

ǫ→0

∂Ωi

Then ϕ ∈ H 1 (E) ∩ H 2 (Th ) and (69) becomes Z Z Z Z ε∇u·∇ϕǫ dx− ε∇u·ν ϕ ¯ ds+ηr,e uϕ¯ ds = e

Ωi

e

Ωi

By the Green theorem Z Z Z   ∇ ε∇u ϕǫ dx + ηr,e uϕ¯ ds = Ωi

f ϕǫ dx+ηr,e

e

f ϕǫ dx + ηr,e

Ωi

Z

uˆϕ¯ ds. (84)

e

Z

u ˆϕ¯ ds.

(85)

e

Passing to the limit ǫ → 0 ηr,e

Z

uϕ¯ ds = ηr,e

e

Z

u ˆϕ¯ ds.

(86)

e

Since ϕ¯ ∈ C0∞ (e) and e ∈ ΓD are arbitrary, we get ˆ|ΣD , u|ΣD = u

(87)

and the Dirichlet boundary conditions are satisfied. ∞ Then take any ϕ ∈ C0,Σ (Ω). Thus D X

ηr,e

e∈ΓDI

Z

[u][ϕ] ds =

e

X

ηr,e

Z

[ˆ u][ϕ] ds = 0.

(88)

e

e∈ΓD

∞ as [ϕ] = 0 for any e ∈ ΓI since ϕ ∈ C0,Σ (Ω) and on e ∈ ΓD we have [ϕ] = ϕ ≡ 0. D Analogously we see that o X Z n − ε∇u [ϕ] ds = 0. (89) e∈ΓDI

e

By assumptions of the theorem u ∈ H 1 (Ω), so [u] = 0 for any e ∈ ΓI while as we have already been shown u = u ˆ for e ∈ ΓD , so Z n o o X X Z n ξ ε∇ϕ · ν [u] ds = ξ ε∇ϕ · ν [ˆ u] ds. (90) e∈ΓDI

Thus we obtain

e

e

e∈ΓD

N Z X i=1

ε∇u · ∇ϕ dx =

Ωi

Z

f ϕ dx.

(91)



∞ Since this statement is true for any ϕ ∈ C0,Σ (Ω), then it is valid also for any D 1 ϕ ∈ H0,ΣD (Ω), so we regain the first statement of (68).

17

7.3

Auxiliary estimates

We define the following operators A(u, ϕ) :=

N Z X i=1

B(u, ϕ) :=

N Z X i=1

C(ϕ) :=

ε∇u · ∇ϕ dx,

Ωi

Ωi

N Z X i=1

 ˆ eu−ˆv − ew−u ϕ dx, k1 ϕ dx,

Ωi

X Z n ∂u o ε D(u, ϕ) := − [ϕ] ds, ∂ν e∈ΓDI e X Z n ∂ϕ o [u] ds, ε E(u, ϕ) := − ∂ν e∈ΓDI e X Z F (ϕ) := − {ε∇ϕ · n}[ˆ u] ds, e

e∈ΓD

Ir (ϕ) :=

X

ηr,e

e∈ΓD

Jr (u, ϕ) :=

X

(92)

Z

ηr,e

[ˆ u] · [ϕ] ds, e

Z

[u] · [ϕ] ds.

e

e∈ΓDI

Lemma 7.5. Let uh ∈ Xh (Ω) and ξ ∈ {0, −1}. Then A(uh , uh ) + J(uh , uh ) + ξD(uh , uh ) + ξE(uh , uh ) ≥ ckuh kh .

(93)

Proof. If ξ = 0, then simply A(uh , uh ) + J(uh , uh ) = kuh k2h .

(94)

Otherwise it is a simple consequence of lemma 3.3. Lemma 7.6. Let u, v ∈ L2 (Ω). Then B(u, u − v) − B(v, u − v) ≥ 0. Proof. Since the exponential function is monotone Z  e−ˆv eu − ev (u − v) dx B(u, u − v) − B(v, u − v) = ZΩ  + ewˆ e−v − e−u (u − v) dx ≥ 0.

(95)

Lemma 7.7. Let u, ϕ ∈ H 1 (E), r ∈ {1, 2}. Then A(u, ϕ) + Jr (u, ϕ) ≤ Ckukh,σr kϕkh,σr .

(96)



18

Proof. This is a simple consequence of the Schwarz inequality. Lemma 7.8. Let u, v, ϕ ∈ H 1 (E), r ∈ {1, 2} and α ≤ u, v ≤ β for some α, β ∈ R. Then (97) B(u, ϕ) − B(v, ϕ) ≤ Cku − vkh,σr kϕkh,σr , where C is a constant dependent on α, β, kˆ v kL∞ (Ω) and kwk ˆ L∞ (Ω) .

Proof. Note that the exponential function is locally Lipschitz-continuous, so since u, v are bounded keu − ev kL2 (Ω) ≤ Cku − vkL2 (Ω) .

(98)

The same is true for e−v − e−u . Thus using the Schwarz inequality and Poincare inequality for the “broken” norm (lemma 7.2) Z Z   v −ˆ v u ewˆ e−v − e−u ϕ dx e − e ϕ dx + e B(u, ϕ) − B(v, ϕ) = (99) Ω Ω ≤ Cku − vkL2 (Ω) kϕkL2 (Ω) ≤ Cku − vkh,σr kϕkh,σr . Lemma 7.9. Let u ∈ H 2 (E) and ϕh ∈ Xh (Ω). Then |D(u, ϕh )| ≤ Ch|u|H 1 (Ω) kϕh kh,σ2 .

(100)

Constant C depends on εM and σ0 . Proof. Using the Schwarz inequality X Z |{ε∇u · ν} [ϕh ]| D(u,ϕh ) ≤ e

e∈ΓDI





"

"

X Z

e∈ΓDI N X i=1

≤ Ch

"

e

−1 η2,e

h2i

{ε∇u · ν}

#1/2 "

X Z

e∈ΓDI

X

e∈ΓDI ∩Γi

X

2

Z

e

σe−1

{ε∇u · ν}

σe−1 {ε∇u · ν}

2

e∈ΓDI

#1/2

2

η2,e [ϕh ]

2

e

#1/2

#1/2

kϕh kh,σ2

(101)

kϕh kh,σ2

≤ Ch|u|H 1 (Ω) kϕh kh,σ2 , where the last inequality follows from the trace theorem. Constant C depends on εM , σ0 and the maximal number of edges of elements of the coarse triangulation E. Lemma 7.10. Let u ∈ H 2 (E), uI := Πh u (see section 6) and ϕh ∈ Xh (Ω). Then |D(u − uI , ϕh )| ≤ Ch

N X i=1

Constant C depends on εM and σ0 . 19

|ui |2H 2 (Ωi )

1/2

kϕh kh,σ1 .

(102)

Proof. We have o X Z n ε∇(u − uI ) · ν [ϕh ] ds

D(u − uI , ϕh ) = −

e∈ΓDI

(103)

e

Let us take any e ∈ ΓI , e ∈ ∂Ωj ∩ ∂Ωk . Then the Schwarz inequality yields that Z n o ε∇(u − uI ) · ν [ϕh ] ds ≤ εM k{∇(u − uI ) · ν}kL2 (e) k[ϕh ]kL2 (e) (104) e

Then by lemma 6.1 1/2

1/2

k{∇(u − uI ) · ν}kL2 (e) ≤ (hj |uj |H 2 (Ωj ) + hk |uk |H 2 (Ωk ) ) ≤ (hj + hk )1/2 (|uj |H 2 (Ωj ) + |uk |H 2 (Ωk ) ).

(105)

Therefore −1 η1,e k{∇(u − uI )}k2L2 (e) ≤ Cσe−1 hj hk (|uj |H 2 (Ωj ) + |uk |H 2 (Ωk ) )2

≤ Ch2 (|uj |H 2 (Ωj ) + |uk |H 2 (Ωk ) )2

(106)

If e ∈ ΓD , e ∈ ∂Ωi , then analogously −1 η1,e k{∇(u − uI ) · ν}k2L2 (e) ≤ Cσe−1 h2i |ui |2H 2 (Ωi ) ≤ Ch2 |ui |2H 2 (Ωi )

(107)

Therefore by Schwarz inequality and the inequalities derived above o X Z n ε∇(u − uI ) · ν [ϕh ] ds e∈ΓDI

≤C

e

 X

−1 η1,e k{∇(u − uI ) · ν}k2L2 (e)

e∈ΓDI

≤ Ch

N X i=1

1/2  X

η1,e k[ϕh ]k2L2 (e)

e∈ΓDI

|ui |2H 2 (Ωi )

1/2

1/2

(108)

kϕh kh,σ1 .

Constant C is independent of h. It depends on σ0 , εM and on the number of elements of ΓDI . Lemma 7.11. Let u ∈ H 2 (E), uI := Πh u (see section 6) and ϕh ∈ Xh (Ω). Then " N #1/2 X X h3i  2 2 |E(u − uI , ϕh )| ≤ Ckϕh kh,σ1 hi + . (109) |ui |H 2 (Ωi ) hk i=1 Ωk ∈nb(Ωi )

Proof. By the Schwarz inequality o X Z n |E(u − uI , ϕh )| ≤ ε∇ϕh · ν [u − uI ] ds, e∈ΓDI

≤ εM

e

X

k{∇ϕh · ν}kL2 (e) k[u − uI ]kL2 (e)

e∈ΓDI

20

(110)

Splitting this sum up



k{∇ϕh · ν}kL2 (e) ≤ ∇ϕh · ν

Ωj L2 (e)

Then using lemma 3.1

2



∇ϕh · ν

Ωi L2 (e)





+ ∇ϕh · ν



Ωk L2 (e)

.

(111)

−1 2 2 ≤ Ch−1 i k∇ϕh kL2 (Ωi ) ≤ Chi kϕh kh,σ1 .

(112)

On the other hand

k[u − uI ]kL2 (e) ≤ ku − uI Ωj kL2 (e) + ku − uI Ωk kL2 (e) .

By lemma 6.1 we have ku − uI Ωi k2L2 (e) ≤ Ch3i |u|2H 2 (Ωi ) , so |E(u − uI , ϕh )| ≤ Ckϕh kh,σ1

"

N  X

h2i

X

+

i=1

Ωk ∈nb(Ωi )

h3i  2 |ui |H 2 (Ωi ) hk

(113)

#1/2

(114)

The differential problem (69) satisfies: A(u∗ , ϕ) + B(u∗ , ϕ) + D(u∗ , ϕ) + ξE(u∗ , ϕ) + Jr (u∗ , ϕ) = C(ϕ) + ξF (ϕ) + Ir (ϕ)

∀ϕ ∈ H 1 (E) ∩ H 2 (Th ).

(115)

On the other hand, the family of discrete problems depending on parameter h is defined as A(u∗h , ϕh ) + B(u∗h , ϕh ) + ξD(u∗h , ϕh ) + ξE(u∗h , ϕh ) + Jr (u∗h , ϕh ) = C(ϕh ) + ξF (ϕh ) + Ir (ϕh )

∀ϕh ∈ Xh .

(116)

Here we use ξ = 0 for CWOPSIP (16) and ξ = 1 for CSIPG (18). We subtract these equations from each other with ϕ := ϕh and we obtain A(u∗ − u∗h , ϕh ) + B(u∗ , ϕh ) − B(u∗h , ϕh ) + D(u∗ − ξu∗h , ϕh ) + ξE(u∗ − u∗h , ϕh ) + Jr (u∗ − u∗h , ϕh ) = 0.

(117)

This is equivalent to LHS = RHS, where LHS := A(u∗I − u∗h , ϕh ) + B(u∗I , ϕh ) − B(u∗h , ϕh ) + ξD(u∗I − u∗h , ϕh ) + ξE(u∗I − u∗h , ϕh ) + Jr (u∗I − u∗h , ϕh ),

(118)

and RHS := A(u∗I − u∗ , ϕh ) + B(u∗I , ϕh ) − B(u∗ , ϕh ) + D(ξu∗I − u∗ , ϕh ) + ξE(u∗I − u∗ , ϕh ) + Jr (u∗I − u∗ , ϕh ).

21

(119)

7.4

Error estimate for CWOPSIP

In this case RHS := A(u∗I − u∗ , ϕh ) + B(u∗I , ϕh ) − B(u∗ , ϕh ) − D(u∗ , ϕh ) + J2 (u∗I − u∗ , ϕh ). (120) We take ϕh := u∗I − u∗h . By lemma 7.5 and lemma 7.6 we have that LHS ≥ cku∗I − u∗h kh,σ2 . Note that u∗ is bounded (lemma 2.1), thus u∗I is also bounded by the same constants, what allows us to use lemma 7.6. On the other hand, we may estimate RHS with lemmata 7.7, 7.8 and 7.9   RHS ≤Cku∗I − u∗h kh,σ2 ku∗ − u∗I kh,σ2 + h|u∗ |H 1 (Ω) . (121)

Therefore using LHS = RHS we obtain   ku∗I − u∗h k2h,σ2 ≤Cku∗I − u∗h kh,σ2 ku∗ − u∗I kh,σ2 + h|u∗ |H 1 (Ω) . Then diving both sides of this inequality by ku∗I − u∗h kh,σ2 > 0   ku∗I − u∗h kh,σ2 ≤C ku∗ − u∗I kh,σ2 + h|u∗ |H 1 (Ω) .

(122)

(123)

Thus by the triangle inequality and interpolation error estimate (60) we have ku∗ − u∗h kh,σ2 ≤ku∗ − u∗I kh,σ2 + ku∗I − u∗h kh,σ2 ≤C

N  X hi + i=1

X

Ωk ∈nb(Ωi )

h3i  ∗ 2 |ui |H 2 (Ωi ) h2k

!1/2

+ Ch|u∗ |H 1 (Ω) . (124)

With the additional assumption hi := ci h for every Ωi ∈ E we can simplify this expression to N  1/2 X ku∗ − u∗h kh,σ2 ≤ Ch1/2 |u∗ |2H 1 (Ω) + . |u∗i |2H 2 (Ωi )

(125)

i=1

If we use particular grids and Dirichlet boundary conditions are piecewise linear, as explained in section 6, we could get a better estimate ∗

ku −

u∗h kh,σ2

N 1/2  X ∗ 2 . |u∗i |2H 2 (Ωi ) ≤ Ch |u |H 1 (Ω) +

(126)

i=1

7.5

Error estimate for CSIPG

As for CWOPSIP, by lemma 7.5 and lemma 7.6 imply LHS ≥ cku∗I − u∗h kh,σ1 . Then we have RHS := A(u∗I − u∗ , ϕh ) + B(u∗I , ϕh ) − B(u∗ , ϕh ) + D(u∗I − u∗ , ϕh ) + E(u∗I − u∗ , ϕh ) + J1 (u∗I − u∗ , ϕh ).

22

(127)

p-GaN 100 nm

n-GaN 100 nm

Figure 4: Schema of the n-p diode used in numerical simulations. This time we may estimate RHS with lemmata 7.7, 7.8, 7.10 and 7.11 RHS ≤

Cku∗I



u∗h kh,σ1

ku∗I



− u kh,σ1

N  X h2i + +

"

i=1

X

Ωk ∈nb(Ωi )

h3i  ∗ 2 |ui |H 2 (Ωi ) hk

#1/2 !

(128)

Thus estimating LHS = RHS from below and above and dividing by ku∗I − u∗ kh,σ1 > 0 we obtain ku∗I



u∗h kh,σ1



ku∗I



− u kh,σ1 +

"

N  X

h2i

+

i=1

X

Ωk ∈nb(Ωi )

h3i  ∗ 2 |ui |H 2 (Ωi ) hk

#1/2 !

.

(129)

Thus by the triangle inequality and interpolation error estimate (56) we have ku∗ − u∗h kh,σ1 ≤ku∗ − u∗I kh,σ1 + ku∗I − u∗h kh,σ1 ≤C

N  X i=1

h2i

X

+

Ωk ∈nb(Ωi )

h3i  ∗ 2 |ui |H 2 (Ωi ) hk

!1/2

.

(130)

If hi := ci h for every Ωi ∈ E this estimate simplifies to ku∗ − u∗h kh,σ1 ≤Ch

N X i=1

8

|u∗i |2H 2 (Ωi )

1/2

.

(131)

Numerical experiments

We performed numerical simulations of a p-n junction (figure 4) to verify whether the theoretical estimates hold in practice. This device consists of two layers composed of the same material, but with different doping. We use homogeneous mesh inside layers of the device, which is nonmatching on the interface of the layers. It is formed by division of the layers to K parts in the longitude and then by dividing the first layer into 3K parts and dividing the second part into 2K parts in the perpendicular direction. We start with initial grid for K = 2, presented in figure 5. 23

.

Figure 5: Initial grid used in numerical simulations of pn junction.

p-GaN 110 nm

In0.1Ga0.9N 3 nm

n-GaN 120 nm

Figure 6: Schema of the quantum well structure used in numerical simulations.

Figure 7: Initial grid used in numerical simulations of the quantum well. To improve the convergence, p-GaN and n-GaN were splitted into two layers to thicken the grid near the quantum well.

24

Table 1: Relative errors of the potential u∗ . Simulation were performed for the pn junction in the equilibrium state. K 2 4 8 16 32 64 128 256

CSIPG L2 (Ω) H 1 (Ω) 1.8e-01 6.9e-01 4.5e-02 (4.0) 4.2e-01 (1.6) 2.3e-02 (1.9) 3.0e-01 (1.4) 9.7e-03 (2.4) 2.1e-01 (1.4) 2.8e-03 (3.4) 1.2e-01 (1.8) 6.1e-04 (4.6) 5.7e-02 (2.0) 1.5e-04 (4.2) 2.9e-02 (2.0) 3.5e-05 (4.2) 1.4e-02 (2.0)

CWOPSIP L2 (Ω) H 1 (Ω) 2.3e-01 9.8e-01 1.2e-01 (1.9) 5.9e-01 (1.7) 3.8e-02 (3.1) 3.2e-01 (1.9) 7.7e-03 (4.9) 2.1e-01 (1.5) 1.5e-03 (5.1) 1.2e-01 (1.8) 4.8e-04 (3.1) 5.7e-02 (2.0) 1.3e-04 (3.8) 2.9e-02 (2.0) 3.3e-05 (3.9) 1.4e-02 (2.0)

In numerical simulations, we use the Poisson equation scaled in SI units. Thus it is more elaborated than equation (1), as then physical constants and material parameters are present. Also it takes into account physical parameters and incomplete ionization of acceptors and donors, so k1 is actually a function of u. More details on the drift-diffusion equations, used by us in the simulations of realistic devices may be found in [16], while details of the linearization of the underlying discrete system may be found in [15]. Since the exact solution is not known, as a reference we use the result of one-dimensional simulation for K = 1024. For every function fK taken into account, we compute the relative errors defined as errorK,L2 (Ω) :=

kfK − fref kL2 (Ω) , kfref kL2 (Ω)

errorK,H 1 (Ω) :=

kfK − fref kH 1 (Ω) , (132) kfref kH 1 (Ω)

where fref is a numerical solution computed on a fine grid, as mentioned before. Errors are presented in function of K. Note that h = c/K for some constant c. Having in mind estimates of section 7, we expect errorK,H 1 (Ω) / error2K,H 1 (Ω) ≈ 2. Results of these simulations are presented in table 1. Both methods start slowly, for K ≤ 32 error ratio is < 2. From K = 64 we observe error ratio ≈ 2. The error norms for CSIPG and CWOPSIP are similar. Then we repeated our simulations for a single quantum well structure (figure 6). It is similar to the p-n junction, but it has a narrow layer between n region and p region, called the quantum well. In this case, we introduce five layers in our mesh, while two additional layers are used to improve the grid in n region and p region near the quantum well layer (figure 7). Results of this simulation is presented in table 2 and they generally agree with our previous simulation. Note that the slow start is absent in this case. This is due to additional layers introduced near the quantum well, where the function fluctuations are crucial.

9

Conclusions

We have presented two methods of composite Discontinuous Galerkin discretization of the drift-diffusion equations in the equilibrium state. These methods were derived from Symmetric Interior Penalty Galerkin method [13] and Weakly Over-Penalized Symmeric Interior Penalty method [3]. 25

Table 2: Relative errors of the potential u∗ . Simulation were performed for the single quantum well in the equilibrium state. K 2 4 8 16 32 64

CSIPG L2 (Ω) H 1 (Ω) 5.2e-02 3.8e-01 2.0e-02 (2.6) 2.3e-01 (1.7) 5.2e-03 (3.8) 1.1e-01 (2.0) 1.3e-03 (3.9) 5.8e-02 (1.9) 3.4e-04 (4.0) 2.9e-02 (2.0) 8.5e-05 (4.0) 1.5e-02 (2.0)

CWOPSIP L2 (Ω) H 1 (Ω) 7.5e-02 4.6e-01 2.4e-02 (3.1) 2.6e-01 (1.8) 7.3e-03 (3.3) 1.2e-01 (2.2) 1.8e-03 (4.0) 5.9e-02 (2.0) 4.5e-04 (4.0) 3.0e-02 (2.0) 1.1e-04 (4.0) 1.5e-02 (2.0)

Both discretizations are shown to be well-defined and the error is estimated. In case of uniform increase of grid density, the H 1 -norm of error of Composite Symmetric Interior Penalty Galerkin (CSIPG) method error is estimated by O(h), while for Composite Weakly Over-Penalized Symmeric Interior Penalty (WOPSIP) method we obtain only O(h1/2 ) estimate. In the latter case, the estimate may be improved to O(h) under additional assumptions on a mesh sequance and boundary conditions. In simulations of seminconductor luminescent devices, these assumptions are not excessively restrictive, as these structures’ design is often based on simple geometric shapes. Numerical simulations presented in this paper stay in a good agreement with these estimates.

References [1] Randolph E. Bank and Donald J. Rose. Some Error Estimates for the Box Method. SIAM Journal of Numerical Analysis, 24:777–787, 1987. [2] A. T. Barker, S. C. Brenner, E.-H. Park, and L.-Y. Sung. Two-Level Additive Schwarz Preconditioners for a Weakly Over-Penalized Symmetric Interior Penalty Method. Journal of Scientific Computing, 47(1):27, 2011. [3] S. C. Brenner, L. Owens, and L.-Y. Sung. A weakly over-penalized symmetric interior penalty method. Electronic Transactions on Numerical Analysis, 30:107–127, 2008. [4] Susanne C. Brenner. Poincare-Friedrichs Inequalities for Piecewise H1 Functions. SIAM Journal on Numerical Analysis, 41:306–324, 2004. [5] Susanne C. Brenner and L. Ridgway Scott. The Mathematical Theory of Finite Element Methods. Springer, New York, 2008. [6] R. K. Coomer and I. G. Graham. Massively Parallel Methods for Semiconductor Device Modelling. Computing, 56:1–27, 1996. [7] Maksymilian Dryja. On Discontinuous Galerkin Methods for Elliptic Problems with Discontinuous Coefficients. Computational Methods in Applied Mathematics, 3(1):76–85, 2003. [8] H. K. Gummel. A self-consistent iterative scheme for one-dimensional steady state transistor calculations. IEEE Trans. Elect. Dev., 11:455–464, 1964. 26

[9] Joseph W. Jerome. Consistency of Semiconductor Modeling: An Existence/Stability Analysis for the Stationary van Roosbroeck System. SIAM Journal of Applied Mathematics, 45(4):565–590, 1985. [10] J. L. Lions. Quelques m´ethodes de r´esolution des probl`emes aux limites non lin´eaires. Dunod/Gauthier-Villars, Paris, 1969. [11] J. J. H. Miller, W. H. A. Schilders, and S. Wang. Application of finite element methods to the simulation of semiconductor devices. Rep. Prog. Phys., 62:277–353, 1999. [12] Daniele A. Di Pietro and Alexandre Ern. Mathematical aspects of Discontinuous Galerkin Methods. Springer, Berlin, 2012. [13] Beatrice Riviere. Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation. Society for Industrial and Applied Mathematics, Philadelphia, 2008. [14] W. V. Van Roosbroeck. Theory of Flow of Electrons and Holes in Germanium and Other Semiconductors. The Bell System Technical Journal, 29:560–607, 1950. [15] Konrad Sakowski, Leszek Marcinkowski, and Stanislaw Krukowski. Modification of the Newton’s method for the simulations of gallium nitride semiconductor devices. Lecture Notes in Computer Science, 8385:551–560, 2014. [16] Konrad Sakowski, Leszek Marcinkowski, Stanislaw Krukowski, Szymon Grzanka, and Elzbieta Litwin-Staszewska. Simulation of trap-assisted tunneling effect on characteristics of gallium nitride diodes. Journal of Applied Physics, 111(12):123115, 2012. [17] Siegfried Selberherr. Analysis and Simulation of Semiconductor Devices. Springer-Verlag, Wien, 1984. [18] S.M. Sze and K.K. Ng. Interscience, Berlin, 2006.

Physics of Semiconductor Devices.

Wiley-

[19] Peter Wilkes. Solid State Theory in Metallurgy. Cambridge University Press, Cambridge, 1973.

27