Monotone iterative methods for the adaptive $nite element
solution of semiconductor equations
R.-C. Chen, Jinn-Liang Liu
∗Department of Applied Mathematics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan Received 15 March 2002; received in revised form 19 April 2003
Abstract
Picard, Gauss–Seidel, and Jacobi monotone iterative methods are presented and analyzed for the adap-tive $nite element solution of semiconductor equations in terms of the Slotboom variables. The adapadap-tive meshes are generated by the 1-irregular mesh re$nement scheme. Based on these unstructured meshes and a corresponding modi$cation of the Scharfetter–Gummel discretization scheme, it is shown that the resulting $nite element sti4ness matrix is an M-matrix which together with the Shockley–Read–Hall model for the generation–recombination rate leads to an existence–uniqueness–comparison theorem with simple upper and lower solutions as initial iterates. Numerical results of simulations on a MOSFET device model are given to illustrate the accuracy and e<ciency of the adaptive and monotone properties of the present methods.
c
2003 Elsevier B.V. All rights reserved.
Keywords: Monotone iteration; Drift-di4usion model; Adaptive $nite element
1. Introduction
Computer-aided simulation is one of the important processes in developing semiconductor devices. Numerical methods for the fundamental semiconductor equations play a signi$cant role in this de-velopment. For most practical device structures, the electrostatic potential and carrier concentrations exhibit extreme layers [26], particularly in the neighborhood of p–n junctions and the oxide [9]. The presence of layers implies that adaptive mesh generation of unstructured grids is inevitable if an accurate and e<cient device simulation platform is required [9].
To obtain numerical solutions of semiconductor equations, one must solve a system of nonlinear algebraic equations resulting from a discretization by, for example, the $nite element method. The
This work was supported by NSC under Grant 90-2115-M-009-028, Taiwan.
∗Corresponding author.
E-mail address:jinnliu@math.nctu.edu.tw(J.-L. Liu).
0377-0427/03/$ - see front matter c 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0377-0427(03)00538-7
standard method for the solution of the system is Newton’s method or its variation. Newton’s method is a local method that, in general, is very sensitive to the initial guess and converges quadratically in a su<ciently small neighborhood of the exact solution [28].
The method of monotone iterations is a classical tool for the study of the existence of the solution of semilinear PDEs of certain types [1,14,29,33]. It is also useful for the numerical solution of these types of problems approximated, for instance, by the $nite di4erence [12,19,30], $nite element [16], or boundary element [6,11,32] method. It is a constructive method that depends essentially on only one parameter, called the monotone parameter herein, which determines the convergence behavior of the iterative process. Based on adaptive 1-irregular $nite element meshes, we extend this classical method to device simulation by exploiting a very special nonlinear property of the drift-di4usion model that the carriers satisfy Maxwell–Boltzmann statistical laws by which the model can be expressed in terms of the electrostatic potential and the Slotboom variables [5,17,35,39].
Embedded in the widely used Gummel’s decoupling algorithm [13,34], three monotone iterative methods, namely, Picard, Gauss–Seidel, and Jacobi methods, are presented and analyzed in this pa-per. By extending the Scharfetter–Gummel discretization scheme proposed in [41] to the adaptive $nite element approximation, it is shown that the resulting sti4ness matrix is an M-matrix which together with the Shockley–Read–Hall model for the generation– recombination rate leads to an existence–uniqueness–comparison theorem with simple upper or lower solution as an initial guess. These methods also yield positivity of carrier concentrations under conditions of strong recombi-nation similar to that of the method presented in [35] for the drift-di4usion model in one-space dimension.
In the next section, we state the model from the Van Roosbroeck system to the corresponding Slotboom-variable formulation and Gummel’s decoupling algorithm. The model is subject to Dirich-let and Neumann types of conditions on various parts of the boundary of a real-life device domain. In Section 3, we $rst analyze the matrix properties of the resulting adaptive $nite element sys-tems for the Poisson equation, which then lead to the M-matrix properties for the semiconductor equations. Starting with the upper and lower solutions as initial guesses, it is shown in Section 4
that maximal and minimal sequences generated by Picard, Gauss–Seidel, and Jacobi iterations all converge monotonically from above and belowto the unique solution of the resulting nonlinear sys-tem. We then summarize in Section 5 our implementation procedures into two algorithms, namely, monotone-Gummel and adaptive algorithms which combine Gummel’s decoupling, monotone itera-tive, and adaptive methods. Section 6 represents a part of our extensive numerical experiments on various n-MOSFET device models to demonstrate the accuracy and e<ciency of adaptive and mono-tone properties of the proposed methods. Moreover, numerical results of the Jacobi and Gauss–Seidel monotone iterations are also given to verify the theoretical results.
2. The drift-diusion model
The steady-state Van Roosbroeck system [43], usually referred to as the drift-di4usion model of semiconductor devices, is
J =q
s(n − p + N
−
1
q∇ · Jn= R(; n; p); (2.2)
1
q∇ · Jp= −R(; n; p); (2.3)
where is the electrostatic potential, n and p are the electron and hole concentrations, q is the elementary charge, s is the permittivity constant of the semiconductor, NA− and ND+ are densities of ionized impurities, Jn and Jp are the electron and hole current densities, and R(; n; p) is a function describing the balance of generation and recombination of electrons and holes. The current densities Jn and Jp are de$ned as follows:
Jn= −q nn∇ + qDn∇n; (2.4)
Jp= −q pp∇ − qDp∇p; (2.5)
where n and p are the $eld-dependent electron and hole mobilities and the di4usion coe<cients of electrons and holes are expressed in the Einstein relations
Dn= VT n; Dp= VT p (2.6)
with VT = kT=q being the thermal voltage, k Boltzmann’s constant, and T a constant temperature. Based on Maxwell–Boltzmann statistical laws [17,39], we may write system (2.1)–(2.3) as
J = F(; u; v); (2.7) ∇ · Dnniexp VT ∇u = R(; u; v); (2.8) ∇ · Dpniexp − VT ∇v = R(; u; v); (2.9) where u = exp −’n VT ; (2.10) v = exp ’p VT (2.11) are the Slotboom variables in which the quasi-Fermi potentials ’n and ’p are expressed as
n = niexp − ’n VT ; p = niexp ’p− VT ;
and F(; u; v) =qn i s u exp VT − v exp − VT +q(NA−− ND+) s ; (2.12) R(; u; v) = n2i(uv − 1) 0
n(niv exp(−=VT) + pT) + 0p(niu exp(=VT) + nT): (2.13) Here we consider particularly the Shockley–Read–Hall generation–recombination model with ni being the intrinsic carrier concentration, 0
n and 0p the electron and hole lifetimes, and pT and nT the electron and hole densities associated with energy levels of the traps.
In device simulation, most numerical methods have been developed for the approximation of system (2.1)–(2.3) with the primal state variables (; n; p) [3,36]. Nevertheless, there are also some methods developed for the Slotboom-variable formulation (2.7)–(2.9) with the state variables (; u; v) [35,39]. It has been observed in [3] that the choice of the variables u and v de$ned in (2.10) and (2.11) is preferable since they lead to self-adjoint and positive de$nite matrices and that the resulting matrices are better scaled than those of ’n and ’p. The theory and solution methods for systems of
self-adjoint partial di4erential equations have reached a very high standard such that a solution of the static semiconductor equations in (; u; v) can be carried out very e<ciently. However, the major drawback of the variables u and v lies in the enormous dynamic range required for real number representation in actual computations. By recalling de$nitions (2.10) and (2.11), for example, we $nd that the exponential terms vary more than 32 orders of magnitude for ∈ [ − 1; 1]V . It is, therefore, obvious that computations are limited to low-voltage applications. Although Newton’s method can be successfully applied to (2.7)–(2.9) [3,25], it is very sensitive to the initial guess of those variables due to its local convergence property. In practical simulation, the device terminal characteristics of I–V curves (i.e., I–V points) is usually of interest. The conventional approach to obtain these curves is a homotopy process from lower biases to higher biases by Newton’s method, which can be very costly in terms of computing time and human work load associated with the convergence problems of the method. On the other hand, with formulation (2.7)–(2.9), the monotone iterative method presented in this paper is a global method for which the initial guess can be chosen in a more arbitrary way, see Theorem 4.1 and Section 6 below. This then allows us to have a simultaneous (parallel) computing of multiple I–V points with various biasing conditions and with independent constant initial guesses for each I–V point calculation. The computational e4ort can thus be dramatically reduced [22].
System (2.7)–(2.9) is subject to some appropriate conditions on the boundary of a rectangular region denoted by ⊂ R2 and shown in Fig. 1. The domain is bounded by the six segments
AB= AB, BC= BC, CD= CD, DE= DE, EF= EF, and FA= FA, i.e., its boundary
9 = AB∪ BC∪ CD∪ DE∪ EF∪ FA:
By assuming the charge neutrality condition and the mass-action law[40], the boundary conditions of the system in terms of the variables , u, and v are described as follows: The Dirichlet part of boundary conditions = VO+ Vb on 9D; (2.14) u = exp −VO VT on 9D; (2.15)
Fig. 1. Geometry of an n-MOSFET device. v = exp VO VT on 9D; (2.16)
where 9D= BC∪ DE∪ FA; VO= VS; VD, or VB, are the source, drain, or substrate voltage applied on the ohmic contact parts BC; DE, or FA, respectively. Vb represents the built-in potential [40]. On the boundary AB and EF, we assume that the normal components of the electric $eld E =−∇
and current densities are zero,
· E = 0; · Jn= 0 and · Jp= 0:
These conditions lead to the Neumann boundary conditions 9 9 = 0 on 9N= AB∪ EF; (2.17) 9u 9 = 0 on 9N; (2.18) 9v 9 = 0 on 9N; (2.19)
where is the outward normal on 9, and to the mixed and Neumann boundary conditions
−= +; s9y−= d9y+ on 9R= CD; (2.20)
9u
9 = 0 on 9R; (2.21)
9v
where VG is the voltage applied on the gate, tox is the gate oxide thickness, d is the permittivity constant of the gate oxide, and the + and − signs refer to as the limits from the oxide and the semiconductor regions, respectively, to the interface. The oxide region O (the rectangular region CIJD) is assumed to be free of charges, i.e., the Laplace equation J=0 holds there and IJ=Vb+VG. We nowdescribe Gummel’s decoupling algorithm for the DD model. With a given initial guess ((0); u(0); v(0)) and for each Gummel’s iteration index g ; g = 0; 1; : : : ; (g+1) is computed by solving the nonlinear Poisson equation in and the Laplace equation in O
J(g+1)= F((g+1); u(g); v(g)) in ; (g+1)= V O+ Vb on 9D; 9(g+1) 9 = 0 on 9N; J(g+1)= 0 in O; (g+1)= Vb+ VG on IJ; 9(g+1) 9 = 0 on CI∪ DJ; (g+1)− = (g+1)+ ; s9y(g+1)− = d9y(g+1)+ on 9R (2.23) and then u(g+1) is computed by solving the electron current continuity equation, with now the known functions (g+1) and v(g), ∇ · D(g+1) n niexp (g+1) VT ∇u(g+1)= R((g+1); u(g+1); v(g)) in ; u(g+1)= exp −VO VT on 9D; 9u(g+1) 9 = 0 on 9N∪ 9R; (2.24)
and $nally v(g+1) is computed by solving the hole current continuity equation, with both (g+1) and u(g+1) known, ∇ · D(g+1) p niexp −(g+1) VT ∇v(g+1) = R((g+1); u(g+1); v(g+1)) in ; v(g+1)= expVO VT on 9D; 9v(g+1) 9 = 0 on 9N∪ 9R (2.25)
This decoupling algorithm is widely used in practical simulations of semiconductor devices. Im-portant analyses of the algorithm pertaining to the existence, stability, convergence, e4ectiveness, etc. have been thoroughly studied by Jerome and Kerkhoven [17,18,20]. The algorithm de$nes an outer iteration in a simulation procedure. The monotone method proposed here is applied individually to each one of the nonlinear algebraic systems resulting from the discretization of di4erential equations (2.23)–(2.25).
Before going into the discrete systems, we describe the essence of the monotone iterative method by using a nonlinear Poisson model problem [16,29,30]. Our description is primarily based on [30]. Consider the semilinear elliptic PDEs
J = f(x; y; ) in
a99 + b = g(x; y) on 9; (2.26)
where a ≡ a(x; y) and b ≡ b(x; y) are nonnegative functions on 9 with a + b ¿ 0, and f and g are prescribed nonlinear and linear functions in their respective domains.
Applying the $nite element method to (2.26) on a certain partition of the domain, we obtain a system of nonlinear algebraic equations in a compact form
AZ = −F(Z); (2.27)
where A is an N×N matrix, Z ≡ (z1; : : : ; zN)Tis an unknown vector, and F(Z) ≡ (F1(Z); : : : ; FN(Z))T is a vector associated with both functions f and g. We denote by Nh the set of grid points associated with the partition of N = ∪ 9, i.e.,
Nh= {(xi; yi) ∈ ∪ 9 : i = 1; 2; : : : ; N};
where N is the total number of regular nodes (see Section 3 for the de$nition) of the $nite el-ement partition. The sets of grid nodes in ; 9; 9D; 9N, and 9R are similarly denoted by
h; 9h; 9h
D; 9hN, and 9hR, respectively.
Starting with a given initial vector Z(0) for problem (2.27), the monotone iterative method generates
a sequence of iterates {Z(m)}; m = 0; 1; : : : ; by solving the equation
AZ(m+1)+ *Z(m+1)= *Z(m)− F(Z(m)); (2.28)
where * is a nonnegative diagonal matrix in which its entries +k; k = 1; : : : ; N, are parameters that are determined by the nonlinear function f. Under various conditions on the matrices A and *; or equivalently on the discretization and the function f, it has been shown in [30] that, for the $nite di4erence approximation, the sequence {Z(m)} generated by Eq. (2.28) converges monotonically to a
solution of (2.27). Obviously, the convergence behavior of the monotone process (2.28) is essentially dedicated by *.
There are some variants of the Picard method (2.28), such as Jacobi, Gauss–Seidel, and block iterative monotone methods [30]. We shall discuss the Gauss–Seidel method and the Jacobi method below.
3. Matrix properties of the drift-diusion model
For each Gummel’s iteration and after the discretization, each one of the nonlinear problems (2.23)–(2.25) will result in a system of nonlinear algebraic equations similar to (2.27) with which our main concern nowis the property of the resulting sti4ness matrix A. The discretization considered here is particularly based on the adaptive $nite element method using the 1-irregular mesh re$nement scheme [10,24].
Let T be a $nite element partition of the domain such that T = {j; j = 1; : : : ; M; N =Mj=1 Nj},
and Sh(T) denote a $nite element space on T. Let Na be a set of N indices that are assigned to
active degrees of freedom (i.e., regular nodes) and Nc assigned to constrained degrees of freedom
(irregular nodes). By an active degree of freedom, we mean one that de$nes a parameter associated with the global sti4ness matrix whereas a constrained degree of freedom is a linear combination of active degrees of freedom that are associated with the constrained node by element connectivity. For each i ∈ Nc, there exists a set A(i) ⊂ Na of corresponding active degrees of freedom such that the
resulting $nite element space Sh(T) consists of continuous functions. Let h be an arbitrary function
in Sh(T). If the element is a square, then h is of the following form:
h= i∈Na iˆbi+ j∈Nc jˆbj = i∈Na iˆbi+ j∈Nc k∈A(j) 1 2kˆbj;
where i are scalars and ˆbi are unconstrained bilinear bases which can be constructed via the
following four shape functions: s1= (1 − 0)(1 − 1)=4;
s2= (1 + 0)(1 − 1)=4;
s3= (1 + 0)(1 + 1)=4;
s4= (1 − 0)(1 + 1)=4
de$ned on the reference element ˆ = {(0; 1): |0| 6 1; |1| 6 1}. For every i ∈ Na, let
C(i) = {j ∈ Nc| i ∈ A(j)}:
We rewrite h in the form
h= i∈Na iˆbi+ k∈Na j∈C(k) 1 2kˆbj = i∈Na i ˆbi+ j∈C(i) 1 2ˆbj :
Thus, the functions bi= ˆbi+
j∈C(i)
1
2ˆbj ∀i ∈ Na (3.1)
form constrained bilinear bases.
Let (xi; yi) ∈ N be a mesh point in T. For each i ∈ Na and using the standard notation i ≈ (xi; yi), there exists a set V (i) ⊂ Na; i ∈ V (i), of active degrees of freedom such that the $nite element approximation of problem (2.23) results in a system
0ii− k∈V (i) 0kk = −Fi(i) + Fi∗; (3.2) where 0k= −Bh(bi; bk) ≡ − ∈T ∇bi· ∇bkdx dy; 0i= Bh(bi; bi); Fi(i) ≡ ∈T
F(i; ui; vi)bidx dy and F∗
i is associated with the boundary conditions in (2.23) if (xi; yi) ∈ 9 and Fi∗=0 if (xi; yi) ∈ .
Since the partition consists of rectangular elements, the following result can be easily proved (see, e.g., [2]) with each type of the 1-irregular elements as given in [24].
Theorem 3.1. The matrix induced by (3.2) is diagonally dominant, i.e., 0i¿
k∈V (i)
0k; (3.3)
0k¿ 0 ∀k ∈ V (i): (3.4)
Furthermore, the strict inequality in (3.3) holds for at least one i ∈ Na.
However, for continuity equations (2.24) and (2.25), it is well known that the Scharfetter–Gummel discretization induces nonphysical di4usion in the direction normal to drift velocity for multidimen-sional problems, which has led to various modi$cations of the method [3,8,21,27,34,37,38,41,42]. In order to obtain the same matrix property as that of Theorem 3.1, we extend in particular the method proposed in [41] to the 1-irregular mesh re$nement scheme. By analogy, it su<ces to consider only the electron continuity equation.
For each j ∈ Na with u
j ≈ u(xj; yj), (2.24) is approximated by
L[uj] ≡ 1juj−
k∈V (j)
1kuk= −Rj(uj) + R∗j (3.5)
or in the more compact matrix form
where 1j= k∈V (j) 1k; (3.7) 1k= 0kdk; (3.8) dk= Dn|(k;j)niB (g+1)j − (g+1)k VT exp (g+1)j VT ; (3.9) 0k= −Bh(bj; bk); (3.10) Rj(uj) = R((g+1)j ; u(g+1)j ; vj(g)); Dn|(k;j)= (Dn(xk; yk) + Dn(xj; yj))=2; U = (u1; : : : ; uN); R(U) = (R1(u1); : : : ; RN(uN))
and B(t) = t=(et− 1) is the Bernoulli function for any real number t and R∗
j is associated with the boundary conditions in (2.24) if (xj; yj) ∈ 9 and R∗j = 0 if (xj; yj) ∈ . Note that, by the de$nition
of the Bernoulli function and of the di4usion coe<cient, the factors dk in (3.9) are positive. We
thus conclude the following result.
Theorem 3.2. The matrix A in (3.6) is diagonally dominant, i.e., 1j¿
k∈V (j)
1k;
1k¿ 0 ∀k ∈ V (j): (3.11)
Furthermore, the strict inequality in (3.11) holds for at least one j ∈ Na.
4. Monotone convergence results
The diagonal dominance of the resulting matrices (i.e., M-matrices) of the model problems (2.23)– (2.25) provides not only stability of numerical solutions (i.e., no nonphysical oscillations) but also convergence of iterative procedures when the special properties of the nonlinearity in these prob-lems are taken into account. Moreover, the existence and uniqueness of the solutions can also be guaranteed by means of the construction of lower and upper solutions which are de$ned as follows: De!nition 4.1. A vector ˜U ≡ ( ˜u1; : : : ; ˜uN) ∈ RN is called an upper solution of (3.6) if it satis$es
the following inequality: 1i˜ui−
k∈V (i)
and ˆU ≡ ( ˆu1; : : : ; ˆuN) ∈ RN is called a lower solution if 1iˆui−
k∈V (i)
1kˆuk6 − Ri( ˆui) + R∗i (4.2)
for i ∈ Na.
As in the previous section, we only consider the monotone convergence for the electron continuity equation (2.24) since the nonlinear functionals in (2.23)–(2.25) are at right-hand sides and these equations are all associated with self-adjoint operators. It is obvious that every solution of (3.6) is an upper solution as well as a lower solution. We say that ˆU and ˜U are ordered if ˆU 6 ˜U. Given any ordered lower and upper solutions ˆU and ˜U, we de$ne
ˆU; ˜U ≡ {U ∈ RN; ˆU 6 U 6 ˜U}; (4.3)
ˆui; ˜ui ≡ {wi∈ R; ˆui6 wi6 ˜ui}: (4.4)
We only consider positive solutions, i.e., 0 ¡ ˆU 6 U 6 ˜U, due to physical relation (2.10). Choose the nonnegative scalars +i such that
+i≡ max 9Ri(wi) 9ui ; wi∈ ˆui; ˜ui or in matrix form * ≡ diag(+i)
for i ∈ Nd. Then by adding the term +iui on both sides of (3.5) we obtain the equivalent system
L[ui] + +iui= +iui− Ri(ui) + R∗i: (4.5)
It is easily seen from the de$nition of +i that +i¿ 0 and
+iui− Ri(ui) + R∗i ¿ +ivi− Ri(vi) + R∗i when ˜ui¿ ui¿ vi¿ ˆui: (4.6) Let NU(0)= ˜U be an initial iterate. We construct a sequence { NU(m+1)} by solving the linear system
1iNu(m+1)i −
k∈V (i)
1kNu(m)i + +iNu(m+1)i = +iNu(m)i − Ri( Nu(m)i ) + R∗
i (4.7)
for m = 0; 1; 2; : : : and i ∈ Na. Similarly, by using U(0) = ˆU as another initial iterate, we obtain a sequence {U(m+1)} from the linear system
1iu(m+1)i −
k∈V (i)
1ku(m)i + +iu(m+1)i = +iu(m)i − Ri(u(m)i ) + R∗i (4.8) for m=0; 1; 2; : : : and i ∈ Na. We refer to { NU(m)} and {U(m)} as the maximal and minimal sequences. We nowshowthat these two sequences are monotone and converge to a solution of (3.6).
Lemma 4.1. The maximal and minimal sequences { NU(m)} and {U(m)} given by (4.7) and (4.8)
with NU(0)= ˜U and U(0)= ˆU possess the monotone property
ˆ
U 6 U(m)6 U(m+1)6 NU(m+1)6 NU(m)6 ˜U; m = 0; 1; 2; : : : : (4.9) Moreover for each m, NU(m) and U(m) are ordered upper and lower solutions.
Proof. Let w(0)i = Nu(0)i − Nu(1)i = ˜ui− Nu(1)i . By (4.7) (1i+ +i)w(0)i = (1i+ +i) ˜ui− k∈V (i) 1kNu(0)k + +iNu(0)i − Ri( Nu(0)i ) + R∗i = 1i˜ui− k∈V (i) 1k˜uk− [ − Ri( ˜ui) + R∗i] ¿ 0:
In viewof Theorem 3.2, w(0)i ¿ 0 for all i ∈ Na. This leads to Nu(0)
i ¿ Nu(1)i . A similar argument using
relations (4.8) and (4.2) gives u(0)i 6 u(1)i . Let wi(1)= Nu(1)i − u(1)i . By (4.7) and (4.8), we have (1i+ +i)w(1)i =
k∈V (i)
1k( Nu(0)k − u(0)k ) + +i( Nu(0)i − u(0)i ) − [Ri( Nu(0)i ) − Ri(u(0)i )]:
It then follows from the relation Nu(0)i ¿ u(0)i , the nonnegativity of 1k and (4.6) that
(1i+ +i)w(1)i ¿ 0;
which again leads to w(1)i ¿ 0 and hence u(0)i 6 u(1)i 6 Nu(1)i 6 Nu(0)i for all i ∈ Na. Assume, by induction,
that u(mi −1)6 u(m)i 6 Nu(m)i 6 Nu(mi −1) for some m ¿ 1. By (4.7), w(m)i = Nu(m)i − Nu(m+1)i satis$es (1i+ +i)w(m)i =
k∈V (i)
1k( Nu(mk−1)− Nu(m)k ) + +i( Nu(mi −1)− Nu(m)k ) − [Ri( Nu(mi −1)) − Ri( Nu(m)k )]:
It follows again from (4.6) and Theorem 3.2 that (1i+ +i)w(m)i ¿ 0:
This yields wi(m)¿ 0 which shows that Nu(m)i ¿ Nu(m+1)i . A similar argument gives u(m)i 6 u(m+1)i and Nu(m+1)i ¿ u(m+1)i . The monotone property (4.9) thus follows by induction.
To showthat NU(m) and U(m) are upper and lower solutions for each m, we observe from (4.7)
that
1iNu(m)i =
k∈V (i)
1kNu(mk −1)+ +i( Nu(mi −1)− Nu(m)i ) − Ri( Nu(mi −1)) + R∗i:
By (4.6), (4.9) and Theorem 3.2, we have 1iNu(m)i ¿
k∈V (i)
1kNu(m)k − Ri( Nu(m)i ) + R∗i:
This shows that NU(m) is an upper solution. The proof for the lower solution U(m) is similar.
Based on the monotone property of Lemma 4.1, we have the following monotone convergence result.
Theorem 4.1. Let ˜U; ˆU be a pair of ordered upper and lower solutions of (3.6). Then the sequences { NU(m)} and {U(m)} generated by solving (4.7) and (4.8) with NU(0) = ˜U and U(0)= ˆU converge
monotonically to the solutions NU and U of (3.6), respectively. Moreover ˆ
U 6 U(m)6 U(m+1)6 U 6 NU 6 NU(m+1)6 NU(m)6 ˜U; m = 1; 2; : : : (4.10) and if U∗ is any solution of (3.6) in ˆU; ˜U then U 6 U∗6 NU.
Proof. By Lemma 4.1, the limits
lim NU(m)= NU and lim U(m)= U as m → ∞
exist and satisfy relation (4.10). Letting m → ∞ in (4.7) and (4.8) shows that NU and U are solutions of (4.5). The equivalence between (3.6) and (4.5) ensures that NU and U are solutions of (3.6). Now if U∗∈ ˆU; ˜U is a solution of (3.6) then ˜U and U∗ are ordered upper and lower solutions. Using
N
U(0)= ˜U and U(0)= U∗, Lemma 4.1 implies that NU(m)¿ U∗ for every m. Letting m → ∞ gives
N
U ¿ U∗. A similar argument using U∗ and ˆU as ordered upper and lower solutions yields U∗¿ U.
This proves the theorem.
In viewof the relation NU ¿ U∗¿ U for any solution U∗ in ˆU; ˜U, NU and U are often called
maximal and minimal solutions in ˆU; ˜U, respectively. In general, these two solutions are not necessarily the same. Nevertheless, they are the same for model problems (2.23)–(2.25).
Theorem 4.2. If the conditions in Theorem 4.1 hold, then NU = U is the unique solution of (3.6). Proof. Let W = NU − U. Then W ¿ 0 and by (3.6)
AW = −R( NU) + R(U) 6 9( NU − U) = 9W; where 9 ≡ max −9Ri(wi) 9ui ; wi∈ ˆui; ˜ui; i ∈ N a : (4.11)
Equivalently, we have (A−9I)W 6 0 where I is the identity matrix. Hence, the inverse (A−9I)−1
exists and is nonnegative since 9 6 0 and A is symmetric due to (3.9) and (3.10). This implies that W 6 0 which leads to W = 0, i.e., in ˆU; ˜U, NU = U = U∗ is the unique solution of (3.6).
Finally, we give some comparison results for the monotone sequences obtained by the three basic iterative methods of Picard, Jacobi, and Gauss–Seidel. The iterative methods are based on system (3.6) with A written in the split form A=D−L−U, where D; L and U are the diagonal, lower-o4 diagonal and upper-o4 diagonal matrices of A, respectively. By Theorem 3.2, the elements of D are positive and those of L and U are nonnegative.
Using the split form of A the three iterative schemes are given as follows: (a) Picard method:
(A + *)U(m+1)= *U(m)− R(U(m)) + R∗: (4.12)
(b) Gauss–Seidel method:
(D − L + *)U(m+1)= *U(m)+ UU(m)− R(U(m)) + R∗: (4.13)
(c) Jacobi method:
(D + *)U(m+1)= *U(m)+ (L + U)U(m)− R(U(m)) + R∗: (4.14)
It is clear that, following Theorem 3.2, the inverses of the matrices (A + *), (D − L + *) and (D + *) all exist and are nonnegative.
Using either ˜U or ˆU as initial iterates, we can construct a sequence from each one of the iterative schemes (4.12)–(4.14). Denote the respective sequences by { NU(m)P }, { NUG(m)} and { NU(m)J } with NU(0)= ˜U
and by {U(m)P }, {U(m)G } and {U(m)J } with U(0)= ˆU, and refer to them again as maximal and minimal
sequences, respectively. The convergence of these sequences and the uniqueness of the limiting solutions can be shown in a similar way as in Theorems 4.1 and 4.2. Furthermore, these three minimal or maximal sequences exhibit an ordered behavior as given in the following theorem for which a proof can be found in [30].
Theorem 4.3. Let ˜U; ˆU be a pair of ordered upper and lower solutions of (3.6). Then the three pairs of sequences ({ NU(m)P }; {U(m)P }), ({ NU(m)G }; {U(m)G }) and ({ NU(m)J }; {U(m)J }) generated by solving (4.12), (4.13), and (4.14), respectively, with NU(0)= ˜U and U(0)= ˆU converge monotonically to the
solution of (3.6). Moreover N
U(m)P 6 NU(m)G 6 NU(m)J and U(m)P ¿ U(m)G ¿ U(m)J for m = 1; 2; : : : : (4.15) 5. Implementation algorithms
We brieRy summarize our implementation procedures for the proposed methods into two algo-rithms. The $rst algorithm is an adaptive process based on the general framework of the weak residual error estimation developed in [23] and on the object-oriented programming (OOP) pro-totype proposed in [24]. The data structure of the prototype is designed to combine 1-irregular mesh re$nement scheme, various types of PDEs, and various numerical methods in an integrated computational platform.
An adaptive algorithm.
Step 1: Initial mesh: Generate a coarse and structured mesh for which the number of nodes can be chosen as small as possible.
Step 2: Junction re8nement: Since the initial mesh is usually very coarse, a preprocessing re-$nement to capture irregularities caused by the junction layers of the doping pro$le and by the inversion layer due to the applied voltages proves to be an essential step for more e4ective re-$nement and faster convergence in the subsequent computations. For the junction layers, Poisson’s problem in (2.23) with low biasing conditions is pre-solved few times with the same procedures of error estimation and re$nement as that of Steps 5 and 6 below.
Step 3: Interface re8nement: For the inversion layer, several re$nement iterations are performed speci$cally along the interface boundary. The re$nement procedure is the same as that of Step 6. We nowhave a nonuniform mesh with a better resolution in the vicinity of the interface and the junction.
Step 4: Solution: On the current mesh, each one of nonlinear problems (2.23)–(2.25) is then approximated by FEM. The resulting nonlinear algebraic equations are then solved by some mono-tone iterative method. Note that this solution procedure consists of an outer loop associated with Gummel’s iteration solving (2.23)–(2.25) consecutively and an inner loop associated with the mono-tone iteration. Moreover, the assembly of global sti4ness and mass matrices of the resulting ap-proximation is not required, that is, the solution of the discretized nonlinear systems is performed
on a node-by-node (regular node) basis. This step is described in more details in the following monotone-Gummel algorithm.
Step 5: Error estimation: All adaptive methods require more or less a posteriori information about the computed solution for optimizing overall computational e4ort in the sense that the methods deliver a given level of accuracy with a minimum of degrees of freedom. In essence, the a posteriori error estimation can be regarded as the driving force of adaptive mechanism. Several methods for error estimation in semiconductor simulations have been proposed and incorporated with adaptive grid re$nement. PISCES [31], for example, limits the variation of the quasi-Fermi potential over one element; the element is re$ned if this limit is exceeded. Bank and Weiser [4] proposed a method to estimate the error based on the computation of the norm of the local residual of the elliptic equation and the jump in the normal derivative of the computed solution at inter-element boundaries. An estimation based on computation of the error in the current continuity equation is developed in [9]. We found that the variation of the approximate potential or electron concentration is much easier and inexpensive to obtain and yet su<ciently e4ective (since the variations di4er drastically from lowto high concentration) for adaptive mesh generation. Variations are calculated with respect to every two nodes in each element, from which the largest one is chosen as an error indicator for the element. Error indicators are obtained in an element-by-element manner according to the hierarchical tree structure of the elements in the OOP database (see [24] for more details). A set of criteria on such as global error (maximum) norms of approximated solutions in inner iteration and outer iteration will also be veri$ed (see the next section). If none of the stopping criteria is satis$ed, the adaptive process will continue to Step 6, otherwise it will go to Step 7 for postprocessing the computed solutions.
Step 6: Re8nement: For each element, if its error indicator is larger than a preset error tolerance, it is divided into four subelements according to the rules of the 1-irregular mesh re$nement scheme as given in [24]. We then move on to Step 4.
Step 7: Postprocessing: All computed solutions are then postprocessed for further analysis of the physical phenomena.
The second algorithm illustrates howthe monotone iterative methods are embedded into the Gum-mel decoupling algorithm. Here we use the notation g as GumGum-mel’s (outer) iteration index and m as the monotone (inner) iteration index.
A monotone-Gummel algorithm Step 0: Set g := 0.
Step 1: Solve the Poisson and Laplace equations in (2.23). Step 1.1: Set m := 0 and set the initial guess
(m)j = ˜j or ˆj if g = 0; (g)j otherwise
for all (xj; yj) ∈ Nh, where ˜j and ˆj are constant values that can be easily veri$ed to be an upper and lower solution of , respectively.
Step 1.2: If g = 0, set u(g) and v(g) by the charge neutrality condition.
Step 1.3: Compute (m+1)j (an unknown real value) by solving the discrete potential system 0j(m+1)j + +j()(m+1)j = k∈V (j) 0k(m)k − Fj((m)j ; u(g)j ; vj(g)) + +j()(m)j ∀(xj; yj) ∈ h; (m+1)j = VO+ Vb ∀(xj; yj) ∈ 9h D; 9(m+1)j 9 = 0 ∀(xj; yj) ∈ 9hN; 0j(m+1)j = k∈V (j) 0k(m)k ∀(xj; yj) ∈ hO; (m+1)j = VG+ Vb ∀(xj; yj) ∈ IJ; 9(m+1)j 9 = 0 ∀(xj; yj) ∈ CI∪ DJ; s9y(m+1)− = d9y(m+1)+ ∀(xj; yj) ∈ 9hR (5.1) with +j() = max qni sVT u(g)j exp j VT + vj(g)exp − j VT ; ˆj6 j6 ˜j :
Step 1.4: Set (m)j := (m+1)j ∀j and m := m + 1. Go to Step 1.3 until the stopping criteria
of the inner iteration are satis$ed. Step 1.5: Set (g+1)j := (m+1)j ∀j.
Step 2: Solve the electron continuity equation in (2.24). Step 2.1: Set m := 0 and set the initial guess
u(m)j =
˜uj or ˆuj if g = 0;
u(g)j otherwise
for all (xj; yj) ∈ Nh, where ˜uj and ˆuj are constant values for all (xj; yj) ∈ Nh that can be easily veri$ed to be an upper and lower solution of u, respectively.
Step 2.2: Compute u(m+1)j (an unknown real value) by solving the discrete electron system 1ju(m+1)j + +j(u)u(m+1)j = k∈V (j) 1ku(m)k − Rj((g+1)j ; u(m)j ; vj(g)) + +j(u)u(m)j ∀(xj; yj) ∈ h; u(m+1)j = exp −VO VT ∀(xj; yj) ∈ 9h D; 9u(m+1)j 9 = 0 ∀(xj; yj) ∈ 9hN∪ 9hR (5.2)
with +j(u) = max 0 n(nivj(g)exp(−(g+1)j =VT) + pT)n2ivj(g)+ p0nTn2ivj(g)+ n3i0pexp((g+1)j =VT) (0n(nivj(g)exp(−(g+1)j =VT) + pT) + 0p(niujexp((g+1)j =VT) + nT))2 : Step 2.3: Set u(m)j := u(m+1)j ∀j and m := m + 1. Go to Step 2.2 until the stopping criteria of
the inner iteration are satis$ed.
Step 2.4: Set u(g+1)j := u(m+1)j ∀j.
Step 3: The hole continuity equation in (2.25) is solved analogously as done in Step 2.
Step 4: Set g := g + 1 and go to Step 1 until the stopping criteria of the outer iteration are satis$ed.
Note that the solution procedure in the above algorithm is completely similar to that of the standard Gummel algorithm except that the initial iterates for the Gummel iterations in Steps 1.1 and 2.1 are chosen by the lower and upper solutions of the corresponding semilinear PDEs. Moreover, it can be seen from (5.1) and (5.2) that the unknown scalars are calculated in a node-by-node manner without explicitly computing the Jacobian matrix as required by Newton’s method. As shown in Theorem 4.1, we have rather large freedom to choose the values of the initial iterates. For elliptic PDEs, the upper or lower solution is readily determined by the boundary conditions. This in turn suggests that the initial guesses in the above algorithm can be deduced from applied voltages and built-in potentials. For example, we can choose ˜j= VD+ Vn and ˆj= VB+ Vp, ˜uj= exp(−VB=VT),
ˆuj=exp(−VD=VT), ˜vj=exp(VD=VT) and ˆvj=exp(VB=VT), where Vnis the built-in potential in n-region
and Vp is the built-in potential in p-region.
6. Numerical example
We consider an n-MOSFET device with a channel length of 0:34 m and with the gate oxide thickness of 7 nm. The doping pro$le for the test device is shown in Fig. 2, which is an elliptical Gaussian distribution with the concentration 1020cm−3 in the source and drain regions and 1016 cm−3
in the p-substrate region. The shallowimplantation is needed to obtain a ‘normal-o4’ device with positive threshold voltage and the deep implantation is necessary to avoid punchthrough. The junction depth is 0:2 m and the lateral di4usion under the gate is 0:08 m.
With the drain voltage VD = 1:0 V and the gate voltage VG = 1:0 V, Fig. 3 illustrates a typical
adaptive $nal mesh which was generated from an initial mesh of 16 elements. The inversion and junction layers are e4ectively captured by the adaptive process. Corresponding to the $nal mesh, these computed solutions were obtained by the Gauss–Seidel monotone iteration with, for exam-ple, an initial guess of the upper solution shown in Fig. 4 for the electron continuity equation. The $nal results of the potential and electron concentration are shown in Figs. 5 and 6, respec-tively.
In order to observe more clearly the monotone convergence behavior, we present some numerical results of Jacobi and Gauss–Seidel monotone iterations for the electron continuity equation. Since exp(−VD=VT) 1:7e − 17, we have chosen ˜u = 1 and ˆu = 1:0e − 18 as an ordered pair of upper
and lower solutions for the iterations. The approximate maximal and minimal solutions Nui and ui
Fig. 2. Doping concentration.
Fig. 3. The $nal adaptive mesh.
respectively. The last rowin these tables gives the number of iterations for the maximal and minimal sequences { Nu(m)}, {u(m)}. It is clearly seen from these tables that Nui is always larger than ui at every mesh point i. Moreover, the number of iterations for the Gauss–Seidel iteration is smaller than that of the Jacobi iteration. The stopping criterion for these iterations is u(m)− u(m−1) 6 0:001VT
where · is the maximum error norm. Note that, for the hole current continuity equation, the initial values are determined by exp(VD=VT). Numerical results for this case are similar to those of Tables 1 and 2.
Fig. 4. The upper solution for the Gauss–Seidel iteration.
Fig. 5. Electrostatic potential.
We make a concluding remark on some numerical aspects in connection with the model formu-lation. Obviously, the dynamic range of the Slotboom values of u and v in (2.10) and (2.11) is enormously large in computations. For this, various scaling strategies have been proposed in the literature to avoid catastrophic roundo4 e4ect [3]. The worst case of the numerics for the variables u and v that we have experienced during the course of the development of our code is about of order 10100 on our computing systems (Unix on DEC workstations and Linux on Pentium III) with
the machine number of order 10300. The ranges of applied voltages that have been tested with our
code are −10 V (the reverse bias) to 10 V (the forward bias) for a p–n diode, 0–5 V (the drain bias) and 0–4 V (the gate bias) for a MOSFET.
Fig. 6. Electron concentration.
Table 1
Upper and lower solutions (in log scale) at nodes (xi; yi) by the Jacobi monotone iteration (No.: Number of iterations)
yi\ xi 0.05 0.5 0.7 0.75 0.8 0.95 Nui 0.95 −4:78E − 10 −0:2423 −16:790 −16:794 −16:797 −16:799 ui −2:07E − 9 −0:4330 −16:799 −16:799 −16:799 −16:799 Nui 0.85 −1:33E − 9 −0:281 −16:253 −16:795 −16:796 −16:798 ui −5:80E − 9 −0:310 −16:278 −16:799 −16:799 −16:799 Nui 0.8 −2:17E − 6 −0:247 −6:168 −11:796 −13:217 −13:746 ui −2:18E − 6 −0:259 −6:169 −11:797 −13:218 −13:746 Nui 0.75 −1:38E − 2 −0:189 −1:040 −1:594 −1:9250 −2:2226 ui −1:39E − 2 −0:193 −1:041 −1:595 −1:9254 −2:2229 Nui 0.5 −6:69E − 2 −0:1133 −0:1782 −0:1945 −0:1945 −0:1945 ui −6:72E − 2 −0:1136 −0:1785 −0:1947 −0:1947 −0:1947 Nui 0.05 −8:16E − 3 −0:0111 −0:0139 −0:0146 −0:0146 −0:0146 ui −8:19E − 3 −0:0112 −0:0139 −0:0146 −0:0146 −0:0146 No. 1619 ( Nui) 1750 (ui)
To quantitatively discuss the issue of matrix conditioning [15] associated with the implementation of these Slotboom variables, we present the conditioning numbers of the sti4ness matrices of (5.1) and (5.2) in Table 3 for several bias conditions.
The second rowfrom the bottom of the table clearly shows that matrix conditioning deteriorates dramatically for large biases due to the term exp(=VT) in (3.9). To improve the conditioning, we
Table 2
Upper and lower solutions by the Gauss–Seidel iteration
yi\ xi 0.05 0.5 0.7 0.75 0.8 0.95 Nui 0.95 −6:46E − 11 −0:1736 −16:794 −16:796 −16:798 −16:799 ui −2:00E − 9 −0:4274 −16:799 −16:799 −16:799 −16:799 Nui 0.85 −1:80E − 10 −0:2681 −16:245 −16:797 −16:797 −16:799 ui −5:61E − 9 −0:3099 −16:277 −16:799 −16:799 −16:799 Nui 0.8 −2:17E − 6 −0:241 −6:167 −11:796 −13:217 −13:746 ui −2:18E − 6 −0:259 −6:169 −11:797 −13:218 −13:746 Nui 0.75 −1:38E − 2 −0:188 −1:040 −1:594 −1:9249 −2:2224 ui −1:39E − 2 −0:193 −1:041 −1:595 −1:9254 −2:2229 Nui 0.5 −6:68E − 2 −0:1131 −0:1781 −0:1943 −0:1943 −0:1943 ui −6:71E − 2 −0:1136 −0:1785 −0:1947 −0:1947 −0:1947 Nui 0.05 −8:15E − 3 −0:0111 −0:0139 −0:0146 −0:0146 −0:0146 ui −8:19E − 3 −0:0112 −0:0139 −0:0146 −0:0146 −0:0146 No. 826 ( Nui) 1627 (ui) Table 3
The matrix conditioning associated with the Slotboom variables (VS= 0 and VB= 0)
Bias conditions VD= 0 VD= 1 VD= 3 VD= 3 VD= 5
VG= 0 VG= 1 VG= 0 VG= 2 VG= 3
Degrees of freedom 749 749 3452 3452 4393
Matrix in (5.1) 2.346E+2 2.346E+2 2.443E+2 2.443E+2 8.646E+2 u Matrix in (5.2) 3.487E+17 1.587E+39 4.419E+187 7.274E+121 1.396E+188 u Matrix in (6.1) 1.515E+4 5.062E+4 1.133E+5 1.106E+5 1.706E+5
can divide (5.2) by this term and ni to obtain
ˆ1ju(m+1)j + ˆ+j(u)u(m+1)j = k∈V (j) ˆ1ku(m)k − R((g+1)j ; u(m)j ; vj(g)) niexp (g+1)j VT + ˆ+j(u)u(m)j ; (6.1) where ˆ1j= exp(−(g+1)j =VT)1j=ni, ˆ1k= exp(−(g+1)j =VT)1k=ni and ˆ+j= exp(−(g+1)j =VT)+j=ni. Before solving (6.1), the division is performed node-by-node for all the known terms in the equation. The conditioning is indeed improved as shown in the last row of the table. Note that the sti4ness matrix is an M-matrix and Theorem 4.1 still holds with this scaling technique. Another way to improve the conditioning is to perform (at the discrete level) the change of the variables n=niue=VT
with ˆ1k= 0kDn|(k;j)B − j+ k VT ; (6.2) ˆ1j= (m) k∈V (j) 0kDn|(k;j)B j− k VT (6.3) is then no longer an M-matrix. For example, for some node i, if B((j− k)=VT) = 1 for k = i and j¿ i then B((j − i)=VT) ¡ B((−j + i)=VT) and ˆ1j¡(m)k∈V (j) ˆ1k. Nevertheless, the mixed or hybrid methods proposed in [7] can be used to recover the M-matrix property and furthermore to have the current conservation property. However, the implementation of these methods is more complicated than that of (6.1) since the discrete system is enlarged by these methods and the matrix reduction by means of static condensation requires an element-wise inversion of the block-diagonal matrix associated with the auxiliary variable. Moreover, a suitable numerical integration formula for the local and global matrices and for the right-hand side vector is required (see [7] for more details). The monotone parameters +j(u) in (5.2) will also be more involved with these methods.
Acknowledgements
The authors would like to thank the referees for many valuable and helpful comments and sug-gestions on the manuscript.
References
[1] H. Amann, Supersolution, monotone iteration and stability, J. Di4erential Equations 21 (1976) 367–377.
[2] O. Axelsson, V.A. Barker, Finite Element Solution of Boundary Value Problems, Academic Press, NewYork, 1984. [3] R.E. Bank, D.J. Rose, W. Fichtner, Numerical methods for semiconductor device simulation, IEEE Trans. Electron
Dev. ED-30 (1983) 1031–1041.
[4] R.E. Bank, A. Weiser, Some a posteriori estimators for elliptic partial di4erential equations, Math. Comp. 44 (1985) 283–301.
[5] M. Berger, Nonlinear mathematical phenomena associated with semiconductor devices, AMS Lectures Appl. Math. 25 (1990) 99–106.
[6] C.A. Brebbia, S. Walker, Boundary Element Techniques in Engineering, Newnes-Butterworths, London, 1980. [7] F. Brezzi, L.D. Marini, P. Pietra, Two-dimensional exponential $tting and applications to drift-di4usion models,
SIAM J. Numer. Anal. 26 (1989) 1342–1355.
[8] J.F. BUurgler, R.E. Bank, W. Fichtner, R.K. Smith, A newdiscretization scheme for the semiconductor current continuity equation, IEEE Trans. Comput. Aided Des. 85 (1989) 479–489.
[9] J.F. BUurgler, W.M. Coughran Jr., W. Fichtner, An adaptive grid re$nement strategy for the drift-di4usion equations, IEEE Trans. Comput.-Aided Des. 10 (1991) 1251–1258.
[10] L. Demkowicz, J.T. Oden, W. Rachowicz, O. Hardy, Toward a universal h–p adaptive $nite element strategy. Part 1. Constrained approximation and data structure, Comput. Methods Appl. Mech. Eng. 77 (1989) 79–112.
[11] Y. Deng, G. Chen, W.M. Ni, J. Zhou, Boundary element monotone iteration scheme for semilinear elliptic partial di4erential equations, Math. Comp. 65 (1996) 943–982.
[12] D. Greenspan, S.V. Parter, Mildly nonlinear elliptic partial di4erential equations and their numerical solution II, Numer. Math. 7 (1965) 129–146.
[13] H.K. Gummel, A self-consistent iterative scheme for one-dimensional steady-state transistor calculations, IEEE Trans. Electron Dev. ED-11 (1969) 455–465.
[14] S. Heikkila, V. Lakshmikantham, Monotone Iterative Techniques for Discontinuous Nonlinear Di4erential Equations, Marcel Dekker, NewYork, 1994.
[15] N.J. Higham, Fortran codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation, ACM Trans. Math. Software 14 (1988) 381–396.
[16] K. Ishihara, Monotone explicit iterations of the $nite element approximations for the nonlinear boundary value problem, Numer. Math. 43 (1984) 419–437.
[17] J.W. Jerome, Consistency of semiconductor modeling: an existence/stability analysis for the stationary van Roosbroeck system, SIAM J. Appl. Math. 45 (1985) 565–590.
[18] J.W. Jerome, Analysis of Charge Transport, Springer, Berlin, Heidelberg, 1996.
[19] R. Kannan, M.B. Ray, Monotone iterative methods for nonlinear equations involving noninvertible linear part, Numer. Math. 45 (1984) 219–225.
[20] T. Kerkhoven, On the e4ectiveness of Gummel’s method, SIAM J. Sci. Statist. Comput. 9 (1988) 48–60.
[21] J.P. Kreskovsky, A hybrid central di4erence scheme for solid state device simulation, IEEE Trans. Electron Dev. ED-345 (1987) 1128–1133.
[22] Y. Li, J.-L. Liu, S.M. Sze, T.-S. Chao, A newparallel adaptive $nite volume method for the numerical simulation of semiconductor devices, Comput. Phys. Comm. 142 (2001) 285–298.
[23] J.-L. Liu, On weak residual error estimation, SIAM J. Sci. Comput. (1996) 1249–1268.
[24] J.-L. Liu, I.-J. Lin, M.-Z. Shih, R.-C. Chen, M.-C. Hsieh, Object oriented programming of adaptive $nite element and $nite volume methods, Appl. Numer. Math. 21 (1996) 439–467.
[25] P.A. Markowich, The Stationary Semiconductor Device Equations, Springer, Vienna, 1986.
[26] P.A. Markowich, C.A. Ringhofer, S. Selberherr, M. Lentini, A singular perturbation approach for the analysis of the fundamental semiconductor equations, IEEE Trans. Electron Dev. ED-30 (1983) 1165–1180.
[27] C.C. McAndrew, K. Singhal, E.L. Heasell, A consistent nonisothermal extension of the Scharfetter–Gummel stable di4erence approximation, IEEE Trans. Electron Dev. Lett. EDL-6 (1985) 446–447.
[28] J. Ortega, W.C. Rheinboldt, Iterative Solutions of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.
[29] C.V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, NewYork, 1992.
[30] C.V. Pao, Block monotone iterative methods for numerical solutions of nonlinear elliptic equations, Numer. Math. 72 (1995) 239–262.
[31] M.R. Pinto, C.S. Ra4erty, R.W. Dutton, PISCES II: Poisson and Continuity Equation Solver, Stanford University, Stanford, CA, 1984.
[32] M. Sakahihara, An iterative boundary integral equation method for mildly nonlinear elliptic partial di4erential equations, in: C.A. Brebbia, G. Maier (Eds.), Boundary Elements IX, Vol. II, Springer, Berlin, Heidelberg, 1985, pp. 13.49–13.58.
[33] D. Sattinger, Topics in Stability and Bifurcation Theory, in: Lecture Notes in Mathematics, Vol. 309, Springer, Berlin, Heidelberg, NewYork, 1973.
[34] D.L. Scharfetter, H.K. Gummel, Large-signal analysis of a silicon Read diode oscillator, IEEE Trans. Electron Dev. ED-16 (1969) 64–77.
[35] T.I. Seidman, S.C. Choo, Iterative scheme for computer simulation of semiconductor devices, Solid-State Electron. 15 (1972) 1229–1235.
[36] S. Selberherr, Analysis and Simulation of Semiconductor Devices, Springer, NewYork, 1984.
[37] S. Selberherr, A. Schutz, H.W. Potzl, MINIMOS—a two-dimensional MOS transistor analyzer, IEEE Trans. Electron Dev. ED-278 (1980) 64–77.
[38] M. Sharma, G.F. Carey, Semiconductor device simulation using adaptive re$nement and Rux upwinding, IEEE Trans. Comput. Aided Des. 8 (1989) 590–598.
[39] J.W. Slotboom, Computer-aided two-dimensional analysis of bipolar transistor, IEEE Trans. Electron Dev. ED-20 (1973) 669–679.
[40] S.M. Sze, Physics of Semiconductor Devices, 2nd Edition, Wiley Interscience, NewYork, 1981.
[41] G.-L. Tan, X.-L. Yuan, Q.-M. Zhang, W.H. Ku, A.-J. Shey, Two-dimensional semiconductor device analysis based on new$nite-element discretization employing the S-G scheme, IEEE Trans. Comput. Aided Des. 8 (1989) 468– 478.
[42] T.W. Tang, Extension of the Scharfetter–Gummel algorithm to the energy balance equation, IEEE Trans. Electron Dev. ED-31 (1984) 64–77.
[43] W.V. Van Roosbroeck, Theory of Rowof electrons and holes in germanium and other semiconductors, Bell Systems Tech. J. 29 (1950) 560–607.
