A fast implicit solver for semiconductor models in one space dimension

TL;DR

This paper introduces a fast implicit solver for one-dimensional semiconductor models based on a simplified Boltzmann-Poisson system, utilizing a sweeping algorithm and acceleration techniques to improve computational efficiency.

Contribution

It develops and compares multiple iterative solvers with acceleration schemes for efficient implicit discretization of semiconductor models.

Findings

Accelerated iterative solvers significantly reduce convergence time.

The sweeping algorithm effectively solves fixed-point problems in this context.

Performance varies with electron mean-free-path, influencing solver choice.

Abstract

Several different approaches are proposed for solving fully implicit discretizations of a simplified Boltzmann-Poisson system with a linear relaxation-type collision kernel. This system models the evolution of free electrons in semiconductor devices under a low-density assumption. At each implicit time step, the discretized system is formulated as a fixed-point problem, which can then be solved with a variety of methods. A key algorithmic component in all the approaches considered here is a recently developed sweeping algorithm for Vlasov-Poisson systems. A synthetic acceleration scheme has been implemented to accelerate the convergence of iterative solvers by using the solution to a drift-diffusion equation as a preconditioner. The performance of four iterative solvers and their accelerated variants has been compared on problems modeling semiconductor devices with various electron mean-free-path.

Tables

Table 1. Table 1: Single-scale problem with ω ​ ( x ) = 1 𝜔 𝑥 1 \omega(x)=1 , δ = ϵ = 0.200 𝛿 italic-ϵ 0.200 \delta=\epsilon=0.200 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . Iteration counts and total computation time (sec) for the compared solvers and their accelerated variants.

Solver	SA	Iteration		Total	Solver	SA	Iteration			Total
		FP	SW	Time			FP	LS	SW	Time
Type-\@slowromancapi@	N	4.6	1.3	8.8	Type-\@slowromancapii@	N	3.2	3.0	1.3	21.4
PI	Y	4.5	1.3	9.7	AA/PI	Y	3.2	2.9	1.3	22.5
Type-\@slowromancapi@	N	3.0	1.4	7.9	Type-\@slowromancapii@	N	3.3	1.4	3.2	43.7
AA	Y	3.3	1.3	9.2	AA/GMRES	Y	7.3	2.1	3.3	97.8

Table 2. Table 2: Single-scale problem with ω ​ ( x ) = 1 𝜔 𝑥 1 \omega(x)=1 , δ = ϵ = 0.200 𝛿 italic-ϵ 0.200 \delta=\epsilon=0.200 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . The values of the efficiency indicator 𝒳 𝒳 \mathcal{X} ( 73 ) and the total number of sweeping iterations in an implicit step (Iter SW SW {}_{\text{SW}} ) are reported. This serves as a formal efficiency comparison between the proposed implicit scheme and the standard explicit and IMEX-AP schemes.

$1$ implicit step $\sim$ $10$ IMEX-AP steps $\sim$ $10$ explicit steps
Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$	Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$
Type-\@slowromancapi@	N	6.0	1.67	1.67	Type-\@slowromancapii@	N	12.5	0.80	0.80
PI	Y	5.9	1.69	1.69	AA/PI	Y	12.1	0.83	0.83
Type-\@slowromancapi@	N	4.2	2.38	2.38	Type-\@slowromancapii@	N	14.8	0.68	0.68
AA	Y	4.3	2.33	2.33	AA/GMRES	Y	50.6	0.20	0.20

Table 3. Table 3: Single-scale problem with ω ​ ( x ) = 1 𝜔 𝑥 1 \omega(x)=1 , δ = ϵ = 0.002 𝛿 italic-ϵ 0.002 \delta=\epsilon=0.002 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . Iteration counts and total computation time (sec) for the compared solvers and their accelerated variants.

Solver	SA	Iteration		Total	Solver	SA	Iteration			Total
		FP	SW	Time			FP	LS	SW	Time
Type-\@slowromancapi@	N	1399.3	1.0	2011.6	Type-\@slowromancapii@	N	2.5	440.1	1.0	2165.2
PI	Y	131.9	1.0	232.9	AA/PI	Y	6.1	84.4	1.0	918.0
Type-\@slowromancapi@	N	49.7	1.0	99.9	Type-\@slowromancapii@	N	6.1	22.7	2.1	392.2
AA	Y	29.0	1.1	67.8	AA/GMRES	Y	7.6	9.2	2.1	218.3

Table 4. Table 4: Single-scale problem with ω ​ ( x ) = 1 𝜔 𝑥 1 \omega(x)=1 , δ = ϵ = 0.002 𝛿 italic-ϵ 0.002 \delta=\epsilon=0.002 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . The values of the efficiency indicator 𝒳 𝒳 \mathcal{X} ( 73 ) and the total number of sweeping iterations in an implicit step (Iter SW SW {}_{\text{SW}} ) are reported. This serves as a formal efficiency comparison between the proposed implicit scheme and the standard explicit and IMEX-AP schemes.

$1$ implicit step $\sim$ $333$ IMEX-AP steps $\sim$ $1000$ explicit steps
Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$	Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$
Type-\@slowromancapi@	N	1399.3	0.24	0.71	Type-\@slowromancapii@	N	1100.3	0.30	0.91
PI	Y	131.9	2.53	7.58	AA/PI	Y	514.8	0.65	1.94
Type-\@slowromancapi@	N	49.7	6.71	20.12	Type-\@slowromancapii@	N	290.8	1.15	3.44
AA	Y	32.0	10.42	31.25	AA/GMRES	Y	146.8	2.27	6.81

Table 5. Table 5: Standard n + superscript 𝑛 n^{+} – n 𝑛 n – n + superscript 𝑛 n^{+} diode problem with ω ​ ( x ) 𝜔 𝑥 \omega(x) given in ( 76 ), δ = ϵ = 0.056 𝛿 italic-ϵ 0.056 \delta=\epsilon=0.056 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . Iteration counts and total computation time (sec) for the compared solvers and their accelerated variants.

Solver	SA	Iteration		Total	Solver	SA	Iteration			Total
		FP	SW	Time			FP	LS	SW	Time
Type-\@slowromancapi@	N	3.8	1.4	7.6	Type-\@slowromancapii@	N	2.6	4.8	1.4	28.6
PI	Y	4.6	1.6	10.1	AA/PI	Y	2.6	4.5	1.4	29.6
Type-\@slowromancapi@	N	2.6	1.4	7.3	Type-\@slowromancapii@	N	2.8	2.2	3.7	49.8
AA	Y	2.8	1.4	8.2	AA/GMRES	Y	6.0	3.5	3.6	110.7

Table 6. Table 6: Standard n + superscript 𝑛 n^{+} – n 𝑛 n – n + superscript 𝑛 n^{+} diode problem with ω ​ ( x ) 𝜔 𝑥 \omega(x) given in ( 76 ), δ = ϵ = 0.056 𝛿 italic-ϵ 0.056 \delta=\epsilon=0.056 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . The values of the efficiency indicator 𝒳 𝒳 \mathcal{X} ( 73 ) and the total number of sweeping iterations in an implicit step (Iter SW SW {}_{\text{SW}} ) are reported. This serves as a formal efficiency comparison between the proposed implicit scheme and the standard explicit and IMEX-AP schemes.

$1$ implicit step $\sim$ $35$ IMEX-AP steps $\sim$ $35$ explicit steps
Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$	Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$
Type-\@slowromancapi@	N	5.3	6.58	6.58	Type-\@slowromancapii@	N	17.5	2.00	2.00
PI	Y	7.4	4.76	4.76	AA/PI	Y	16.4	2.14	2.14
Type-\@slowromancapi@	N	3.6	9.62	9.62	Type-\@slowromancapii@	N	22.8	1.54	1.54
AA	Y	3.9	8.93	8.93	AA/GMRES	Y	75.6	0.46	0.46

Table 7. Table 7: Multiscale problem with ω ​ ( x ) 𝜔 𝑥 \omega(x) given in ( 77 ), δ = ϵ = 0.002 𝛿 italic-ϵ 0.002 \delta=\epsilon=0.002 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . Iteration counts and total computation time (sec) for the compared solvers and their accelerated variants. Here ∞ \infty denotes the case that diverging residual is observed.

Solver	SA	Iteration		Total	Solver	SA	Iteration			Total
		FP	SW	Time			FP	LS	SW	Time
Type-\@slowromancapi@	N	79.4	1.3	125.7	Type-\@slowromancapii@	N	1.5	107.2	1.3	395.0
PI	Y	$\infty$	$\infty$	$\infty$	AA/PI	Y	1.7	42.4	1.3	184.7
Type-\@slowromancapi@	N	14.2	1.3	26.0	Type-\@slowromancapii@	N	1.9	30.0	3.0	207.8
AA	Y	3.9	1.1	10.0	AA/GMRES	Y	2.2	9.9	2.8	92.0

Table 8. Table 8: Multiscale problem with ω ​ ( x ) 𝜔 𝑥 \omega(x) given in ( 77 ), δ = ϵ = 0.002 𝛿 italic-ϵ 0.002 \delta=\epsilon=0.002 , and Δ ​ t = Δ ​ x Δ 𝑡 Δ 𝑥 \Delta t=\Delta x . The values of the efficiency indicator 𝒳 𝒳 \mathcal{X} ( 73 ) and the total number of sweeping iterations in an implicit step (Iter SW SW {}_{\text{SW}} ) are reported. This serves as a formal efficiency comparison between the proposed implicit scheme and the standard explicit and IMEX-AP schemes.

$1$ implicit step $\sim$ $1000$ IMEX-AP steps $\sim$ $1000$ explicit steps
Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$	Solver	SA	Iter${}_{\text{SW}}$	${𝒳}_{\text{IMEX}}$	${𝒳}_{\text{EX}}$
Type-\@slowromancapi@	N	103.2	9.69	9.69	Type-\@slowromancapii@	N	209.0	4.78	4.78
PI	Y	$\infty$	—	—	AA/PI	Y	93.7	10.67	10.67
Type-\@slowromancapi@	N	18.5	54.05	54.05	Type-\@slowromancapii@	N	171.0	5.85	5.85
AA	Y	4.3	232.56	232.56	AA/GMRES	Y	61.0	16.39	16.39

Equations

$${\partial }_{t}f+v{\partial }_{x}f+\frac{q_e}{m}E{\partial }_{v}f=\omega (\rho M_{\Theta }-f),$$

$$E={\partial }_x\Phi ,{\partial }_x^2\Phi =\frac{q_e}{{\epsilon }_p}(\rho -D).$$

$$M_{\Theta }(v):=(2\pi \Theta {)}^{-\frac{1}{2}}{e}^{-v^2/2\Theta },$$

$$\hat{f}(\hat{t},\hat{x},\hat{v})=\frac{f(t,x,v)}{f_0},\hat{\omega }(\hat{x})=\frac{\omega (x)}{{\omega }_0},\hat{D}(\hat{x})=\frac{D(x)}{D_0},\hat{\rho }(\hat{t},\hat{x})=\frac{\rho (t,x)}{{\rho }_0}.$$

$$\delta {\partial }_{t}f+$$

$$E$$

$${\Theta }^{\text{nicefrac}12}={\left(\frac{k_{\text{B}}}{m}T\right)}^{\text{nicefrac}12},B_0:={\left(\frac{q_e[\Phi ]}{m}\right)}^{\text{nicefrac}12},\text{mbox}and{C}_0:=x_0{\omega }_{\text{pe}}:=x_0{\left(\frac{q_e^2{D}_0}{{\epsilon }_{p}m}\right)}^{\text{nicefrac}12}.$$

$$f(t,x,v)=\varrho (t,x)M_{{\alpha }^2}(v)+\mathcal{O}(ϵ).$$

$$\xi {\partial }_t\varrho -$$

$$E={\partial }_x\Phi ,{\partial }_x^2\Phi =\frac{{\gamma }^2}{{\beta }^2}(\varrho -D),$$

$$v{\partial }_x{f}^{n+1}+{\beta }^2{E}^{n+1}{\partial }_v{f}^{n+1}+\left(\frac{\delta }{\Delta t}+\frac{\omega }{ϵ}\right)f^{n+1}=\frac{\omega }{ϵ}{\rho }^{n+1}{M}_{{\alpha }^2}+\frac{\delta }{\Delta t}{f}^n,$$

$$v{\partial }_{x}f+\eta E{\partial }_{v}f+\sigma f=q,$$

$$\partial {\mathcal{D}}^{\mathrm{in}}:=\{(x,v):(v,E)\cdot n(x,v)\le 0\}\text{mbox}and\partial {\mathcal{D}}^{\mathrm{out}}:=\{(x,v):(v,E)\cdot n(x,v)>0\}.$$

$$\{\begin{array}{ll}v{\partial }_{x}f+\eta E{\partial }_{v}f+\sigma f=q& \\ f(x,v)=(\Lambda f)(x,v),& (x,v)\in \partial {\mathcal{D}}^{\Lambda }\\ f(x,v)=w(x,v),& (x,v)\in \partial {\mathcal{D}}^{\mathrm{Data}}\end{array}$$

$$\mathcal{F}={\mathcal{F}}^{\Lambda }\cup {\mathcal{F}}^{\mathrm{Data}}\cup {\mathcal{F}}_x\cup {\mathcal{F}}_v,$$

$$g^{h,\pm }(x,v)=\underset{ϵ\to 0^{+}}{\lim }{g}^h(x\pm ϵ,v\pm ϵ),\text{for all }(x,v)\in e.$$

$${\hat{g}}^h(x,v)=\{\begin{array}{ll}{g}^{h,-}(x,v),& {\bar{v}}_{e_x}>0,\\ g^{h,+}(x,v),& {\bar{v}}_{e_x}<0,\end{array}{\hat{g}}^h(x,v)=\{\begin{array}{ll}{g}^{h,-}(x,v),& {\bar{E}}_{e_v}>0,\\ g^{h,+}(x,v),& {\bar{E}}_{e_v}<0.\end{array}$$

$$g^{h,\partial }(x,v)=\underset{ϵ\to 0^{+}}{\lim }{g}^h(x-ϵn_x,v-ϵn_v),\text{for all }(x,v)\in e,$$

$$\mathcal{A}(f^h,g^h)=\mathcal{Q}(g^h),\forall g^h\in {\mathcal{Z}}^h,$$

$$\mathcal{A}(f^h,g^h):=$$

$$-\sum _{e\in {\mathcal{F}}_x}{\bar{v}}_e{\int }_e{\hat{f}}^h(g^{h,+}-g^{h,-})\text{d}v-\sum _{e\in {\mathcal{F}}_v}\eta {\bar{E}}_e{\int }_e{\hat{f}}^h(g^{h,+}-g^{h,-})\text{d}x$$

$$-\sum _{e\in {\mathcal{F}}_x^{\Lambda }}|{\bar{v}}_e|{\int }_e(\Lambda f^h)g^{h,\partial }\text{d}v-\sum _{e\in {\mathcal{F}}_v^{\Lambda }}\eta |{\bar{E}}_e|{\int }_e(\Lambda f^h)g^{h,\partial }\text{d}x$$

$$\mathcal{Q}(g^h):=\sum _{C\in \mathcal{T}}{\iint }_{C}q{g}^h\text{d}x\text{d}v+\sum _{e\in {\mathcal{F}}_x^{\mathrm{Data}}}|{\bar{v}}_e|{\int }_{e}w{g}^{h,\partial }\text{d}v+\sum _{e\in {\mathcal{F}}_v^{\mathrm{Data}}}\eta |{\bar{E}}_e|{\int }_{e}w{g}^{h,\partial }\text{d}x,$$

$$g_1^h(x,v)=\{\begin{array}{ll}{g}^h(x,v),& (x,v)\in {\mathcal{D}}_1,\\ 0,& \text{otherwise},\end{array}{g}_2^h(x,v)=\{\begin{array}{ll}{g}^h(x,v),& (x,v)\in {\mathcal{D}}_2,\\ 0,& \text{otherwise}.\end{array}$$

$$\mathcal{A}(f^h,g^h)=\mathcal{A}(f_1^h,g_1^h)+\mathcal{A}(f_2^h,g_2^h)+\mathcal{A}(f_1^h,g_2^h)+\mathcal{A}(f_2^h,g_1^h).$$

$$\mathcal{A}(f_1^h,g_2^h)$$

$$\mathcal{A}(f_2^h,g_1^h)$$

$${\hat{f}}_1^{h,\ast }={\mathcal{P}}_1(f_1^h):=\{\begin{array}{ll}{\hat{f}}_1^h,& \text{on }{\mathcal{F}}_0^{-},\\ 0,& \text{on }{\mathcal{F}}_0^{+},\end{array}\text{mbox}and{\hat{f}}_2^{h,\ast }={\mathcal{P}}_2(f_2^h):=\{\begin{array}{ll}{\hat{f}}_2^h,& \text{on }{\mathcal{F}}_0^{+},\\ 0,& \text{on }{\mathcal{F}}_0^{-}.\end{array}$$

$$\mathcal{A}(f_1^h,g_1^h)$$

$$\mathcal{A}(f_2^h,g_2^h)$$

$${\hat{f}}_1^{h,\ast }$$

$${\hat{f}}_2^{h,\ast }$$

$$A\mathbf{f}$$

$${\hat{\mathbf{f}}}^{\ast }$$

$$(I+P{A}^{-1}B){\hat{\mathbf{f}}}^{\ast }=P{A}^{-1}\mathbf{q},$$

Full text

1 Introduction

The Boltzmann-Poisson system is considered an accurate kinetic model of electron transport in semiconductor devices [43]. This system describes the evolution of an electron distribution function using a semi-classical Boltzmann kinetic equation and generates a self-consistent electric field by coupling the Boltzmann equation to a Poisson equation that is driven by the electron density. Numerical simulation of the Boltzmann-Poisson system is known to be difficult for several reasons, including the nonlinear coupling between equations, the nonlinear collision operator that describes electron-electron and electron-background interactions, and the dimension of the computational domain. Indeed, simulating a three-dimensional device requires the solution to a six-dimensional Boltzmann equation.

Under a low-density assumption, electron-electron interactions become negligible and electrons can be treated as classical particles interacting with a material background. In such cases, the nonlinear collision operator can be replaced by a linear, relaxation time approximation [43, 23, 13] when the steady-state solution to the Boltzmann-Poisson system is of interest. In the case that the electron transport is through channels that are parallel to the electric field, the semiconductor device is effectively one-dimensional [13, 43] and can therefore be simulated using the approximated Boltzmann-Poisson system in one space dimension.

In addition to traditional Direct Simulation Monte Carlo methods [32], many deterministic numerical schemes have been developed for solving the Boltzmann-Poisson system and its simplified variants. The deterministic schemes considered in previous works discretize the position-velocity phase space of the Boltzmann equation and its simplified variants using the weighted essentially non-oscillatory (WENO) finite difference method [11, 10, 12], the discontinuous Galerkin (DG) method [17, 16, 18], and the spectral-difference method [48]. These schemes either consider a steady-state Boltzmann-Poisson system or use explicit time-stepping schemes to capture transient behavior. To guarantee stability, explicit time-stepping schemes usually require the size of the time steps to be proportional to the mean-free-path of particles in the system. Such a restriction can be computational prohibitive for highly collisional problems, where the mean-free-path is small. An explicit asymptotic preserving scheme was introduced to address this issue in [34], where the stability is guaranteed under a parabolic time-step restriction (independent of the mean-free-path) for highly collisional problems. Later, an implicit-explicit (IMEX) asymptotic preserving scheme was developed in [21] to relax the parabolic time-step restriction. However, these and similar approaches do not allow for large-scale variations in the mean-free-path that are common in multiscale problems.

In this paper, we consider a fully implicit numerical scheme for solving the simplified Boltzmann-Poisson system under the low density assumption in one space dimension. The fully implicit time-stepping method allows for larger time steps that are independent of the mean-free-path, regardless of the collisionality of the problem. Such stability comes at the cost of solving large, possibly ill-conditioned, linear and nonlinear algebraic equations. Hence efficient numerical solvers are needed to update the numerical solution at each time step. In [24], a fast, fully implicit solver was proposed for the nonlinear Vlasov-Poisson system, which is the collisionless variant of the simplified Boltzmann-Poisson system considered in this paper. At each time step, this solver applies a special decomposition of the phase space to allow for the use of the sweeping technique that are commonly used to accelerate the solution of radiation transport problems [2, 37, 39]. To utilize the fast solver [24] in the collisional case, we consider the scattering term as a source and formulate the simplified Boltzmann-Poisson system as a nonlinear fixed-point problem at each implicit time step. At each fixed-point iteration, a collisionless problem with source is solved and the electron distribution that solves the collisionless problem is used to update the collision term, which becomes the source at the next fixed-point iteration. The fixed-point problem reaches a solution when the collisionless problem gives an electron distribution that is consistent with the collision term.

Two types of fixed-point formulations for the simplified Boltzmann-Poisson systems are considered and compared in this paper. The main difference between the two formulations is in the treatment of the scattering source in the relaxation operator: in the first case, both the electric field and the scattering source are lagged; in the second, only the electric field is lagged. As a result, problems in the first formulation are solved in a single iteration loop, while the ones in the second formulation require two nested loops. We apply various iterative solvers on these problems and compare their performance. To solve problems in the first formulation, we consider Picard iteration (see, e.g., [27, Section I.8]) and also Anderson acceleration [5, 55]. For problems in the second formulation, we solve the nonlinear outer loop via Anderson acceleration and solve the linear inner loop using either Picard iteration or the generalized minimal residual method (GMRES) [49]. We do not apply Picard iteration on the outer loop since preliminary numerical results suggest that this approach is not competitive.

We also consider accelerated/preconditioned variants of the solvers described above, based on the idea of synthetic acceleration (SA) [38, 35, 2, 39], an approach that accelerates convergence of iterative solvers by applying a correction term in between each iteration. This correction term is obtained by solving coarse, cheap, or low-order approximate equations to the error equation of the base iterative solver. Thus, many of the SA schemes can be viewed as two-level multigrid algorithms [36] or preconditioned iterative solvers [2, 20, 9]. For neutron transport problems, the correction terms can be computed by solving the transport equation on a coarse mesh [42, 9] or solving a diffusion equation [3, 4, 2] that is a low-order approximation to the transport equation near the collision limit where the mean-free-path is small. In this paper, we compute the correction term by solving a drift-diffusion equation that approximates the simplified Boltzmann-Poisson system in the highly collisional, low-field regime.

The various strategies are tested on a one-dimensional silicon $n^{+}$–$n$–$n^{+}$ diode problem [11, 14, 17, 13, 31]. The algebraic equations to be solved are derived via a backward Euler discretization in time and a discontinuous Galerkin (DG) discretization in the position-velocity phase space. The low-order time discretization is chosen primary to simplify the presentation; however the DG discretization is important for capturing the drift-diffusion limit. While other discretizations are possible, our focus here is on the efficiency of the solver strategy. Thus for various levels of collisionality, the computation time and iteration count of each fixed-point formulation and each iterative solver are compared. These results provide a guideline on the selection of iterative solvers for problems with different material profiles.

The remainder of the paper is organized as follows. In Section 2, the simplified Boltzmann-Poisson system, the drift-diffusion equation, and the time, space, and velocity discretization for solving them are described. Section 3 provides details of fixed-point formulations for the discretized equations as well as the iterative solvers for these fixed-point problems. In Section 4, the implementation details and numerical results for the various iterative solvers are reported for the $n^{+}$–$n$–$n^{+}$ diode problem. Conclusions and discussion are given in Section 5.

2 Preliminaries

2.1 Semiconductor models

We consider the kinetic model

[TABLE]

Here (1a) describes the evolution of the electron distribution $f=f(t,x,v)$, which is a function of position $x\in[0,L]$, velocity $v\in\mathbb{R}$, and time $t\geq 0$; the electric field $E=E(t,x)$ in (1b) is the spatial gradient111There is a sign difference between the electric field defined in (1) and the usual physics definition convention. We make this choice to match the definition used in the fast sweeping algorithm in [24] for solving Vlasov-Poisson systems. In particular, the sign choice in (1) implies that sign of $E$ determines the direction of flow in the velocity variable. The sweeping algorithm is used extensively in this paper. See Section 2.2.3 for details. of a potential $\Phi=\Phi(t,x)$ that satisfies a Poisson equation with a source due to the balance between a given doping profile $D=D(x)$ and the particle (electron) concentration $\rho=\rho(t,x)=\int_{\mathbb{R}}f(t,x,v)\text{d}v$. The constants $q_{e}$, $m$, and $\varepsilon_{p}$ denote, respectively, the magnitude of the electron charge, the effective electron mass, and the electric permittivity of the material. The collision frequency $\omega=\omega(x)$ takes the form $\omega=\frac{q_{e}}{m\mu}$ with electron mobility $\mu=\mu(x)$, and the absolute Maxwellian

[TABLE]

where the background temperature $\Theta:=\frac{k_{\textup{B}}}{m}T$ with $k_{\textup{B}}$ the Boltzmann constant and $T$ the lattice temperature. For the detail derivations of this model, we refer to the reader to [43, 51].

2.1.1 Scaled semiconductor models

Since the qualitative behavior of solutions to (1) largely depends on the scales of the system, we introduce non-dimensional variables $\hat{t}=\frac{t}{t_{0}}$, $\hat{x}=\frac{x}{x_{0}}$, $\hat{v}=\frac{v}{v_{0}}$ and express (1) in terms of the scaled variables

[TABLE]

Let $\Phi_{0}$ be a nominal value of the potential $\Phi$ and let $\Delta\Phi(t,x):=\Phi(t,x)-\Phi_{0}$. Because the solution of (1) is independent of $\Phi_{0}$, we consider the scaled potential $\hat{\Phi}=\Phi_{0}+\frac{\Delta\Phi}{[\Phi]}$. By assuming $\hat{\rho}=\int_{\mathbb{R}}\hat{f}\text{d}\hat{v}$, the scaled system takes the form

[TABLE]

where the hats on the variables dropped for simplicity. In (4), the kinetic Strouhal and Knudsen numbers [26, 50] are given by $\delta=\frac{x_{0}}{v_{0}t_{0}}$ and $\epsilon=\frac{v_{0}}{x_{0}\omega_{0}}$, respectively, and the ratios are defined as $\zeta:=\frac{\rho_{0}}{D_{0}}$, $\alpha:=\frac{\Theta^{\nicefrac{{1}}{{2}}}}{v_{0}}$, $\beta:=\frac{B_{0}}{v_{0}}$, and $\gamma:=\frac{C_{0}}{v_{0}}$ with

[TABLE]

Here $\Theta^{\nicefrac{{1}}{{2}}}$ is the thermal velocity, $B_{0}$ is the ballistic velocity, and in plasma physics, $\omega_{\textup{pe}}$ is known as the plasma frequency [7, 30].

2.1.2 The drift-diffusion limit

To describe semiconductors with different characteristics, there exist various scaling of the semiconductor model (4), such as the low-field scaling [43, 46], the high-field scaling [52, 45, 13], and the ballistic scaling [11]. In this paper, we use the solution of a drift-diffusion equation as a preconditioner to accelerate the solution procedure for (4) in the low-field, highly collisional regime. While this preconditioner is expected to work well only in this regime, the discretizations of (4) and the formulation of solvers for the resulting algebraic equations do not rely on any particular scaling.

The low-field scaling of (4) assumes that the ratio $\beta$ is an $\mathcal{O}(1)$ quantity and that the ratio $\zeta=1$, i.e., the scaling of the particle concentration $\rho$ and the doping profile $D$ is identical. Under these assumptions, when $\epsilon$ is small, i.e., when the electron mean-free-path $\lambda:=\frac{v_{0}}{\omega_{0}}$ is much smaller than the spatial scale $x_{0}$, the collision term on the right-hand side of (4a) is much larger than the drift term $\beta^{2}E\partial_{v}f$; it thus becomes the dominant term. In this situation, it is necessary to choose $\delta\approx\epsilon$ in order to observe the nontrivial dynamics in the long time scale, in which case (4) can be approximated by a drift-diffusion-Poisson model.

It is shown in [46] that when the potential $\Phi$ is known and sufficiently smooth, a standard drift-diffusion model can be derived from (4a) in the low-field, collision limit $(\delta\approx\epsilon\to 0)$ by expanding the distribution function $f$ via the Hilbert expansion as

[TABLE]

This result has been extended in [1] and [44] respectively to the one-dimensional and multi-dimensional Boltzmann-Poisson systems with a self-consistent potential via Poisson coupling as in (4b). The resulting drift-diffusion-Poisson model takes the form

[TABLE]

where (7a) is a drift-diffusion equation coupled with the Poisson system (7b), and $\xi>0$ is the ratio between $\delta$ and $\epsilon$, i.e., $\delta=\xi\epsilon$. We refer to this low-field, collision limit as the “drift-diffusion limit”.

In this paper, we use the solution to the drift-diffusion equation as a preconditioner for solving the scaled semiconductor kinetic equation (4a). The numerical method we used to solve the drift-diffusion equation is discussed in Section 2.4, and the drift-diffusion preconditioning approach is introduced in Section 3.3.

2.2 Solving the kinetic equation

In this paper, implicit time discretization of the kinetic equation (4a) is considered. In the simplified case when the electric field $E$ and the particle concentration $\rho$ are known a priori, implicit time discretization of (4a) leads to linear systems that can be solved efficiently via the fast sweeping algorithm proposed in [24]. The implicit time discretization, the position-velocity phase space discretization, and the fast sweeping approach for solving (4a) in the simplified case are discussed in Sections 2.2.1, 2.2.2, and 2.2.3, respectively. Based on the method presented in this section, we propose several iterative solvers in Section 3 for solving (4a) in the self-consistent case that $E$ and $\rho$ are coupled via the Poisson equation (4b).

2.2.1 Time discretization

In the temporal domain $[0,\,t_{\textup{final}}]$, we apply a uniform discretization with time step size $\Delta t$ and denote $f^{n}\approx f(t^{n},\cdot,\cdot)$, where $t^{n}=n\Delta t$. To simplify the presentation, we consider the backward Euler scheme in this paper. Although this scheme is only first-order accurate, it can be used as a building block for higher-order implicit schemes, such as the singly diagonally implicit Runge-Kutta method (SDIRK) [27, 28]. Applying the backward Euler scheme to (4a) leads to

[TABLE]

where $\rho^{n+1}=\int_{\mathbb{R}}f^{n+1}\text{d}v$ and $E^{n+1}$ is coupled via the Poisson equation (4b) with $\rho=\rho^{n+1}$. For the remainder of Section 2.2, we assume that $\rho^{n+1}$ and $E^{n+1}$ are known a priori at time $t^{n}$. In this simplified case, (8) becomes a linear steady-state Vlasov problem with a source:

[TABLE]

where, in an abuse of notation, $f=f^{n+1}(x,v)$ denotes the steady-state unknown, $E=E^{n+1}(x)$ denotes a given electric field, $\eta$ and $\sigma$ are positive constants, and $q=q(x,v)$ denotes a general source. This simplified steady-state problem (9) can be solved efficiently using the fast sweeping algorithm proposed in [24]. We give the details of this algorithm in Section 2.2.3.

2.2.2 Phase space discretization

For the position-velocity ($x$-$v$) phase space, we truncate the velocity domain from $\mathbb{R}$, the entire real line, to a finite interval $[a_{v},\,b_{v}]$. The position-velocity computational domain is then $\mathcal{D}:=[0,\,L]\times[a_{v},\,b_{v}]$. Let $\partial\mathcal{D}$ be the boundary of $\mathcal{D}$ and $n(x,v)\in\mathbb{R}^{2}$ be the outward normal at $(x,v)\in\partial\mathcal{D}$. As usual, we decompose the boundary into two disjoint pieces: $\partial\mathcal{D}=\partial\mathcal{D}^{\rm{in}}\cup\partial\mathcal{D}^{\rm{out}}$ where

[TABLE]

We further decompose $\partial\mathcal{D}^{\rm{in}}$ into pieces: $\partial\mathcal{D}^{\rm{in}}=\partial\mathcal{D}^{\Lambda}\cup\partial\mathcal{D}^{\rm{Data}}$, where $\partial\mathcal{D}^{\Lambda}$ is the portion of the inflow boundary that depends on the interior solution, e.g., periodic or reflecting boundary, and $\partial\mathcal{D}^{\rm{Data}}$ is the portion upon which data is given. With these notations, the steady state equation (9) takes the form

[TABLE]

where $w$ is known and the abstract linear operator $\Lambda$, which maps functions on $\partial\mathcal{D}^{\text{out}}$ to functions on $\partial\mathcal{D}^{\Lambda}$, can be used to describe periodic or reflecting boundary conditions.

The computational domain $\mathcal{D}$ is discretized into $N_{x}\times N_{v}$ rectangular cells of uniform size $\Delta x\times\Delta v$. For $i=1,\dots,N_{x}$ and $j=1,\dots,N_{v}$, the cell $C_{i,j}$ is centered at $(x_{i},v_{j}):=((i-\frac{1}{2})\Delta x,a_{v}+(j-\frac{1}{2})\Delta v)$. We denote the set of all cells by $\mathcal{T}$ and the set of all edges by $\mathcal{F}$. The set $\mathcal{F}$ is then decomposed into disjoint sets

[TABLE]

where $\mathcal{F}^{\Lambda}$ contains cell edges on the boundary component $\partial\mathcal{D}^{\Lambda}$, $\mathcal{F}^{\rm{Data}}$ contains cell edges on the boundary component $\partial\mathcal{D}^{\rm{Data}}$, and $\mathcal{F}_{x}$ and $\mathcal{F}_{v}$ contains the remaining cell edges that are perpendicular to the $x$ and $v$ axes, respectively. We further decompose $\mathcal{F}^{\Lambda}=\mathcal{F}^{\Lambda}_{x}\cup\mathcal{F}^{\Lambda}_{v}$ and $\mathcal{F}^{\rm{Data}}=\mathcal{F}^{\rm{Data}}_{x}\cup\mathcal{F}^{\rm{Data}}_{v}$, where the subscripts denote the axis to which the edges are perpendicular.

Let $\mathcal{Z}=\{g\in L^{2}(\mathbb{R}\times\mathbb{R})\colon g|_{\mathcal{D}^{c}}=0\}$ and $\mathcal{Z}^{h}:=\{g^{h}\in\mathcal{Z}\colon g^{h}|_{C}\in\mathbb{P}^{1}(C)\>,\,\forall C\in\mathcal{T}\}$, where $\mathbb{P}^{1}(C)$ denotes the space of polynomials up to degree one on the cell $C$. For $g^{h}\in\mathcal{Z}^{h}$, the traces on the two sides of an edge $e\in\mathcal{F}_{x}\cup\mathcal{F}_{v}$ are defined as

[TABLE]

For these edges, the numerical trace of $g^{h}$ is defined via upwinding. Specifically, let $\bar{v}_{e_{x}}$ denote the value of $v$ at the center of edge $e_{x}\in\mathcal{F}_{x}$, and let $\bar{E}_{e_{v}}$ denote the value of $E$ at the center of edge $e_{v}\in\mathcal{F}_{v}$, the numerical traces on $e_{x}$ and $e_{v}$ are respectively defined as

[TABLE]

This definition guarantees a constant upwind direction on each edge. For a test function $g^{h}\in\mathcal{Z}^{h}$, the traces on a boundary edge $e\in\mathcal{F}^{\Lambda}\cup\mathcal{F}^{\rm{Data}}$ are defined as

[TABLE]

where $n_{x}$ and $n_{v}$ are the first and second components of the outward normal $n$, respectively.

With these definitions, the discontinuous Galerkin method solves for $f^{h}\in\mathcal{Z}^{h}$ that satisfies

[TABLE]

with the bilinear operator

[TABLE]

and the source

[TABLE]

where $\bar{v}_{C}$ and $\bar{E}_{C}$ are the values of $v$ and $E$ at the center of cell $C$, respectively.

In the case of neutral particles ($E=0$), the upwind definition of the numerical traces in (14) allows (16) to be solved with an explicit sweeping procedure that moves through the computational domain in a direction determined by the sign of $v$. Such sweeping procedures are commonly used for solving radiation transfer problems [2, 39, 37]. However, in the case of charged particles, the procedure no longer applies when both $v$ and $E$ are allowed to change sign over the phase space domain. This is because, in such cases, changes in the upwind direction may create cyclic dependencies in the elements, see, e.g., [24, Figure 3.1]. To address this challenge, a domain decomposition approach was introduced in [24] to break these dependencies. We briefly discuss this approach and the associated sweeping method in the next subsection.

2.2.3 Domain decomposition and fast sweeping

The domain decomposition method separates the phase space into subdomains along the line $\{v=0\}$. Let $\mathcal{D}_{1}:=[0,\,L]\times(0,\,b_{v}]$ and $\mathcal{D}_{2}:=[0,\,L]\times[a_{v},\,0)$ be the subdomains, and let $\Gamma:=[0,\,L]\times\{0\}$. We assume that there exists some index $j_{0}$ such that $C_{i,j_{0}}\subseteq\mathcal{D}_{2}$ and $C_{i,j_{0}+1}\subseteq\mathcal{D}_{1}$, and denote the set of cell edges in $\Gamma$ as $\mathcal{F}_{0}$. We further decompose $\mathcal{F}_{0}$ into disjoint sets $\mathcal{F}_{0}^{+}$ and $\mathcal{F}_{0}^{-}$ based on the sign of the electric field at the edge center. For $g^{h}\in\mathcal{Z}^{h}$, we define

[TABLE]

The bilinear form in (16) then can be expanded as

[TABLE]

By the definition of $\mathcal{A}$, it is straightforward to verify that

[TABLE]

where the edges in $\mathcal{F}_{0}^{+}$ do not appear in (21a) since $f_{1}^{h}$ does not contribute to the numerical traces $\hat{f}^{h}$ on these edges due to upwinding. Similarly, the edges in $\mathcal{F}_{0}^{-}$ are not included in (21b). We then define for all $\mathcal{F}_{0}^{\pm}\subset\Gamma$, the edge values

[TABLE]

The system (16) is then equivalent to the coupled system

[TABLE]

where $\mathcal{B}(\hat{f}^{h,*}_{1},g^{h}_{2}):=\mathcal{A}(f^{h}_{1},g^{h}_{2})$, $\mathcal{B}(\hat{f}^{h,*}_{2},g^{h}_{1}):=\mathcal{A}(f^{h}_{2},g^{h}_{1})$, and the equations in (23a) and (23b) are coupled only through the projections in (23c) and (23d). Suppose that $\hat{f}^{h,*}_{1}$ and $\hat{f}^{h,*}_{2}$ are known; then (23a) and (23b) are fully decoupled. In each subdomain ($\mathcal{D}_{1}$ or $\mathcal{D}_{2}$) of the phase space, the sign of $v$ is fixed and only $E$ is allowed to change sign. Thus, there is no cyclic dependency in these subdomains, and the decoupled systems (23a) and (23b) can be solved independently via explicit sweeping approach in $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, respectively.

Let $\mathbf{f}$ be the expansion coefficients of $f^{h}:=f^{h}_{1}+f^{h}_{2}$ on an orthogonal basis of $\mathbb{P}^{1}(C)$ for all $C\in\mathcal{T}$, and let $\hat{\mathbf{f}}^{*}$ be the expansion coefficients of $\hat{f}^{h,*}:=\hat{f}^{h,*}_{1}+\hat{f}^{h,*}_{2}$ on an orthogonal basis of $\mathbb{P}^{1}(e)$ for all $e\in\mathcal{F}_{0}$. Then (23) can be expressed a linear system in $\mathbf{f}$ and $\hat{\mathbf{f}}^{*}$ as

[TABLE]

where $\mathbf{q}$ is the vector of expansion coefficients of the source term $\mathcal{Q}(g^{h})$, and the matrices $A$, $B$, and $P$ are defined based on the operators $\mathcal{A}$, $\mathcal{B}$, $\mathcal{P}_{1}$, and $\mathcal{P}_{2}$. See [24] for the detailed definitions. Applying $A^{-1}$ from the left on both sides of (24a) and plugging the resulting equation into (24b) leads to a much smaller linear system

[TABLE]

where the operation $A^{-1}$ can be performed efficiently via the sweeping approach. Then $\hat{\mathbf{f}}^{*}$ can be computed by solving (25) with a Krylov solver, such as the generalized minimal residual method (GMRES) [49]. After obtaining $\hat{\mathbf{f}}^{*}$, the full expansion coefficient $\mathbf{f}$ is then computed by a final sweeping procedure

[TABLE]

To summarize, (25)–(26) defines a mapping from the discretized source $\mathbf{q}$ to the vector $\mathbf{f}$ which solves the discretized form of the steady-state equation (11).

2.3 Solving the Poisson equation

We solve the Poisson equation (4b) with a continuous Galerkin method with $\mathbb{Q}^{1}$ elements on the same spatial mesh as given in Section 2.2.2. Because the method is standard, we omit the details and refer the reader to, for example, [8, 22] for complete presentation. For a given particle concentration $\rho$, this method maps the Galerkin discretization of $\rho$ to a discretized potential $\Phi$, and, since $\mathbb{Q}^{1}$ elements are used, the discretized electric field $E$ can be directly calculated from $\Phi$.

2.4 Solving the drift-diffusion equation

In the drift-diffusion limit, numerically solving the scaled semiconductor kinetic equation (4a) is difficult since the system is stiff. As discussed in Section 2.1.2, the drift-diffusion equation (7a) serves as a good approximation to (4a) near the drift-diffusion limit. To accelerate the solution procedure of (4a) near this limit, we apply a synthetic acceleration [2, 37] and use the solution to (7a) as a preconditioner when solving (4a). The detailed discussion of this acceleration technique is given in Section 3.3. Here we focus on the discretization of (7a). Discretizing in time with backward Euler gives

[TABLE]

If $E^{n+1}$ is known a priori at time $t^{n}$, then (27) takes the steady-state form

[TABLE]

where $E$ is a given electric field and $Q$ is a general source.

We solve (28) with the direct Discontinuous Galerkin method with interface correction (DDG-IC). The original DDG scheme [40] is derived based on the weak formulation of (28) with numerical fluxes approximating derivatives of the solution at element boundaries. The interface correction for DDG was introduced later in [41] to obtain the optimal $(k+1)$-th order of accuracy for polynomial approximations of degree $k$. It is shown in [15] that, with proper choices of numerical flux and limiters, the DDG-IC method satisfies the maximum principle with accuracy up to third order.

Let the spatial domain $[0,L]$ be divided into $N_{x}$ cells $\{I_{j}\}_{j=1}^{N_{x}}$, where $I_{j}=[x_{j-1},x_{j}]$ with $x_{j}=j\Delta x$, and let $\mathcal{V}^{h}=\{\varphi^{h}\in L^{2}(\mathbb{R})\colon$ $\varphi^{h}|_{[0,L]^{c}}=0\>,\,\varphi^{h}|_{I_{j}}\in\mathbb{P}^{1}(I_{j})\>,\,\forall j\ =1,\dots,N_{x}\}$ denote the numerical solution space. For $\varphi^{h}\in\mathcal{V}^{h}$, we define the numerical trace of $\varphi^{h}$ at cell interface $x_{j}$ as $\varphi^{h,\pm}_{j}:=\lim_{\epsilon\to 0^{+}}\varphi^{h}(x_{j}\pm\epsilon)$ for $j=0,\dots,N_{x}$. The jump and average of $\varphi^{h}$ at $x_{j}$ are defined respectively as

[TABLE]

The DDG-IC scheme then solves (28) by finding the solution $\varrho^{h}\in\mathcal{V}^{h}$ such that, for any test function $\varphi^{h}\in\mathcal{V}^{h}$ and on any cell $I_{j}$,

[TABLE]

where $\tau^{h}$ is the $L^{2}$ orthogonal projection of $\omega^{-1}$ onto $\mathcal{V}^{h}$, and the fourth and fifth terms in (30) are the interface correction terms. Here the numerical flux $\widehat{\tau^{h}\partial_{x}{\varrho^{h}}}$ at $x_{j}$ is defined as

[TABLE]

and the Lax-Friedrich flux is used for $\widehat{\tau^{h}E\varrho^{h}}$, i.e.,

[TABLE]

where $E_{j}$ and $E_{j+1}$ are the values of $E$ at the cell centers of $I_{j}$ and $I_{j+1}$, respectively. Here the term $\omega^{-1}$ is not involved in the numerical fluxes since the collision frequency $\omega$ is assumed to be known on the entire spatial domain, including the cell boundaries.

3 Nonlinear solution strategies

In this section, we propose several strategies for solving (4). To simplify the discussion, we first introduce a concise operator notation for the fast sweeping method discussed in Section 2.2.3. Specifically, we write (8) as

[TABLE]

where

[TABLE]

We refer to $\mathcal{S}\rho$ and $s$ as the scattering source and the general source, respectively. If $E^{n+1}=\widetilde{E}$ and $\rho^{n+1}=\widetilde{\rho}$, where $\widetilde{E}$ and $\widetilde{\rho}$ are known, then $f^{n+1}$ satisfies a steady-state problem of the form (9):

[TABLE]

We let

[TABLE]

denote the discretization of (35) as described in Section 2.2.2. Here the operators $\mathsf{L}$ and $\mathsf{S}$ are discretized versions of $\mathcal{L}$ and $\mathcal{S}$ while $\mathsf{f}$, $\widetilde{\mathsf{E}}$, $\widetilde{\uprho}$, and $\mathsf{s}$ denote the discretizations of $f$, $\widetilde{E}$, $\widetilde{\rho}$, and $s$, respectively. The solver discussed in Section 2.2.3 then computes $\mathsf{f}^{n+1}$ by solving the linear problem

[TABLE]

where the operation $\mathsf{L}_{\widetilde{\mathsf{E}}}^{-1}$ is performed using the sweeping algorithm from [24] that is summarized in Section 2.2.3.

In the self-consistent setting, $\rho^{n+1}:=\int_{\mathbb{R}}f^{n+1}\textrm{d}v$, and $E^{n+1}$ is coupled to $\rho$ via the Poisson equation (4b). Thus, instead of the linear problem (37), we solve

[TABLE]

where $\mathsf{P}$ denotes integration over the computational velocity domain $[a_{v},\,b_{v}]$ and $\mathsf{F}$, which maps a given particle concentration to an electric field, denotes the solution procedure of the Poisson equation (4b) using the continuous Galerkin method in Section 2.3. This problem is nonlinear since $\mathsf{E}$ depends on $\mathsf{f}$.

The problem (38) can be solved via nonlinear fixed-point iterative solvers. However, the cost of solving (38) could be prohibitive in the multi-dimensional setting due to the high dimensionality of $\mathsf{f}$. To reduce the problem dimension, a common trick (see, e.g., [39, 2]) is to integrate the first equation in (38) with respect to $v$ and solve the resulting fixed-point problem for $\uprho$ and $\mathsf{E}$:

[TABLE]

The solution to (39) gives $\uprho^{n+1}$ and $\mathsf{E}^{n+1}$. Thus $\mathsf{f}^{n+1}$ can be computed by a final sweeping procedure by setting $\widetilde{\mathsf{E}}=\mathsf{E}^{n+1}$ and $\widetilde{\uprho}=\uprho^{n+1}$ in (37). In the remainder of this section, we consider two different formulations of (39).

3.1 Type-@slowromancapi@ formulation

The type-@slowromancapi@ approach formulates (39) as a nonlinear fixed-point problem on $\uprho$, i.e.,

[TABLE]

where $\mathsf{G}_{1}$ is nonlinear due to the coupling between $\mathsf{E}(=\mathsf{F}(\uprho))$ and $\uprho$. To solve (40), two iterative solvers are considered: standard Picard iteration (see, e.g., [27, Section I.8]) and Anderson acceleration [5, 55].

3.1.1 Type-@slowromancapi@ – Picard iteration

With an initial guess $\uprho^{(0)}$, Picard iteration lags both the scattering source and the electric field terms in $\mathsf{G}_{1}$ and updates the electron concentration by evaluating $\mathsf{G}_{1}$. Specifically, the Picard iteration update at iteration $k+1$ is given by

[TABLE]

It is well-known that Picard iteration converges when $\mathsf{G}_{1}$ is a contraction mapping, and the rate of convergence depends on the spectral radius of the Jacobian of $\mathsf{G}_{1}$.

3.1.2 Type-@slowromancapi@ – Anderson acceleration

Anderson acceleration was first proposed in [5] as an acceleration method based on nonlinear Krylov solvers for fixed-point problems. Here we adopt the variant given in [55] for solving (40). At iteration $k+1$, Anderson acceleration first computes the residual

[TABLE]

then solves the least-squares problem

[TABLE]

with $m_{k}:=\min\{m,k\}$, and finally updates

[TABLE]

Here the truncation parameter $m$ is a nonnegative integer that indicates the maximum number of residuals maintained in memory. When $m=0$, Anderson acceleration reduces to standard Picard iteration. For $m>0$, Anderson acceleration updates $\uprho$ with a convex combination of the previous $m_{k}$ iterates that leads to the minimum residual. It is proved in [53] that Anderson acceleration converges if Picard iteration converges. Since Anderson acceleration utilizes information from the previous iterations, it is expected to converge faster than Picard iteration in practice, but at a cost of additional memory usage. When the problem is linear, Anderson acceleration has been shown to be equivalent to GMRES under some mild assumptions [55].

3.2 Type-@slowromancapii@ formulation

The type-@slowromancapii@ fixed-point formulation of (39) aims to reduce the nonlinearity of the type-@slowromancapi@ formulation (40) by relaxing the coupling between the electric field $\mathsf{E}$ and electron concentration $\uprho$ in the iterative procedure. In other words, the type-@slowromancapii@ formulation gives fixed-point problems that can be solved by a nested iterative procedure that consists of a nonlinear outer loop on $\mathsf{E}$ and a linear inner loop on $\uprho$. The intent is that the inner loop will provide a fast, accurate update of $\uprho$ to feed into the outer loop, thereby improving the overall efficiency. To derive the type-@slowromancapii@ formulation, we write (39) as

[TABLE]

where $\mathsf{I}$ is the identity operator. Thus the nonlinear fixed-point problem on $\mathsf{E}$ takes the form

[TABLE]

where the right-hand side involves solving a linear system. To formulate a fixed-point on $\rho$ as the type-@slowromancapi@ problem (40), the type-@slowromancapii@ formulation takes the following equivalent form of (46):

[TABLE]

Here $\mathsf{G}_{2}(\uprho)$ depends on $\uprho$ only through the electric field $\mathsf{F}(\uprho)$, and evaluating $\mathsf{G}_{2}(\uprho)$ also requires solving the linear system. Specifically, to compute $\mathsf{G}_{2}(\widetilde{\uprho})$ for a given $\widetilde{\uprho}$, we solve

[TABLE]

Nonlinear iterative solvers are used in the nested procedure to solve (47) in the outer loop, while linear iterative solvers are considered for solving (48) in the inner loop. Both Picard iteration and Anderson acceleration can serve as the nonlinear solver in the outer loop. However, we only consider Anderson acceleration in this paper, since preliminary numerical results indicate that using Picard iteration here is not competitive in terms of computation time. For clarity, Anderson acceleration for solving (47) is stated here. At iteration $k+1$ in the outer loop, Anderson acceleration first computes the residual

[TABLE]

then solves the least-squares problem

[TABLE]

with $m_{k}:=\min\{m,k\}$, and finally updates

[TABLE]

To evaluate $\mathsf{G}_{2}(\uprho^{(k)})$ in (49), we replace $\widetilde{\uprho}$ with $\uprho^{(k)}$ in (48) and solve the resulting linear system using either Picard iteration or GMRES [49]. These iterative solvers form the inner loop of the nested procedure. In the following subsections, we discuss the details of the application of Picard iteration and GMRES in the inner loop.

3.2.1 Type-@slowromancapii@ – Anderson acceleration with Picard iteration

Let $\uprho^{(k)}$ denote the $k$-th iterate in the outer loop. To evaluate $\mathsf{G}_{2}(\uprho^{(k)})$, we apply Picard iteration on an equivalent fixed-point formulation of (48). Specifically, at iteration $k+1$ in the outer loop, Picard iteration updates

[TABLE]

where $\uprho^{(k,0)}=\uprho^{(k)}$. Let $\uprho^{(k,*)}:=\lim_{\ell\to\infty}\uprho^{(k,\ell)}$ denote the limit point of iterates generated by (52), then it follows that $\mathsf{G}_{2}(\uprho^{(k)})=\uprho^{(k,*)}$.

3.2.2 Type-@slowromancapii@ – Anderson acceleration with GMRES

Let $\uprho^{(k)}$ still denote the $k$-th iterate in the outer loop. $\mathsf{G}_{2}(\uprho^{(k)})$ can also be evaluated by solving the linear system

[TABLE]

using GMRES. In general, the GMRES solver is expected to converge in fewer iterations than the Picard iteration in (52). We will verify this in the numerical results reported in Section 4.

3.3 Synthetic acceleration for semiconductor equations

Synthetic acceleration (SA) schemes were first developed in [38, 35] to improve efficiency of iterative solvers for transport equations. The basic idea of these schemes is to compute a coarse, cheap, or low-order correction term from residuals of the base iterative solver, and apply this correction to the current iterate to accelerate convergence of the base solver. As noted in [36], many of the synthetic acceleration schemes can be formulated as two-level multigrid algorithms.

In Section 3.3.1, we derive an SA scheme for the semiconductor equation (4a), and then apply this scheme to both type-@slowromancapi@ and type-@slowromancapii@ Picard iteration solvers considered in Sections 3.1.1 and 3.2.1. In Section 3.3.2, we follow the same approach as in [2, 20, 9] to formulate the synthetic acceleration as a preconditioner. We then accelerate the type-@slowromancapi@ and type-@slowromancapii@ Krylov solvers in Sections 3.1.2 and 3.2.2 by applying these solvers on the preconditioned problems.

3.3.1 SA scheme on Picard iteration

The derivation of SA schemes for steady-state linear transport equations with neutral particles can be found in, for example, [2, Section II.B.] and [39, Section 2-3]. These equations are well approximated in collisional by diffusion equations, which are used to compute cheap corrections to iterates in a solver, resulting in the diffusion synthetic acceleration (DSA) scheme [3, 4, 2, 39].

For the semiconductor equations, it is known [43, 46, 44] that drift-diffusion equations serve as proper low-order approximations to the semiconductor equations in the drift-diffusion limit, as discussed in Section 2.1.2. In this section, we derive an SA scheme for semiconductor equations with the correction term computed by solving a drift-diffusion equation. The derivation is mostly a straightforward extension of the derivation for the DSA scheme in [2], and we include it here for completeness.

To derive the SA scheme, we first rewrite the semi-discrete semiconductor equation (33) as a steady-state equation

[TABLE]

Here $\mathcal{P}$ denotes integration over the velocity domain and $\mathcal{F}$ denotes the mapping from the particle concentration $\rho=\mathcal{P}f$ to the electric field $E=\mathcal{F}(\mathcal{P}f)$ via the Poisson equation (4b). Applying $\mathcal{L}_{\mathcal{F}(\mathcal{P}f)}^{-1}$ on both sides of (54) and solving the resulting equation with Picard iteration leads to

[TABLE]

where index of the update is now $k+\frac{1}{2}$ instead of $k+1$. Integrating (55) with respect to $v$ gives

[TABLE]

which is a continuous version of the type-@slowromancapi@ Picard iteration (41). To derive a correction for (56), we write (55) as

[TABLE]

and subtract (57) from (54). By adding and subtracting terms in the resulting equation, we have

[TABLE]

where $\psi:=f-f^{(k+\nicefrac{{1}}{{2}})}$ denotes the error. If $\psi$ can be computed by solving (58), then it can be used as a correction to Picard iterates by taking $f^{(k+1)}=f^{(k+\nicefrac{{1}}{{2}})}+\psi$. However, solving (58) is equivalent to finding $f$.

The SA scheme considered here computes corrections to Picard iterates by solving a reduced-order equation that approximates (58). To obtain the reduced-order equation, we apply $\mathcal{P}$ to both sides of (58), which leads to

[TABLE]

where $\phi:=\mathcal{P}\psi$ is the integral of $\psi$ over the velocity domain. Here the first term in (58) vanishes since, for any electron distribution $g$ and any electric fields $E_{1}$ and $E_{2}$, it follows from (34) that $\mathcal{P}(\mathcal{L}_{E_{1}}-\mathcal{L}_{E_{2}})g=\beta^{2}(E_{1}-E_{2})\int_{\mathbb{R}}\partial_{v}g\textrm{d}v=0$, provided that $g$ goes to zero as $v\to\pm\infty$. Motivated by the drift-diffusion limit, we approximate the operators on the left-hand side of (59) by a drift-diffusion operator $\mathcal{D}_{E}$ defined as

[TABLE]

Thus, (59) is approximated by

[TABLE]

We correct the Picard iterate in (56) with $\widetilde{\phi}$, a solution to (61), by taking $\rho^{(k+1)}=\rho^{(k+\nicefrac{{1}}{{2}})}+\widetilde{\phi}$, which results in the SA scheme

[TABLE]

where (62b) follows from the definition of $\widetilde{\phi}$ in (61). The SA scheme based on the type-@slowromancapi@ Picard iteration (41) is then given by a discretized version of (62)

[TABLE]

where $\mathsf{D}_{\mathsf{F}(\uprho^{(k)})}^{-1}$ denotes the solution procedure of the drift-diffusion equation (61) using the DDG-IC solver presented in Section 2.4.

An analogous SA scheme based on the type-@slowromancapii@ Picard iteration (52) can be derived by repeating the analysis above with minor modification in the indices in (55). The resulting SA takes the form

[TABLE]

where $\mathsf{D}_{\mathsf{F}(\uprho^{(k)})}^{-1}$ is still performed using the DDG-IC solver in Section 2.4.

3.3.2 Preconditioner form of the SA scheme

The SA scheme in Section 3.3.1 is derived specifically for Picard iteration. It is well-known [9, 20, 2] that many SA schemes can be formulated as preconditioners. Following this approach, we derive the SA scheme in the preconditioner forms for the type-@slowromancapi@ Anderson acceleration and the type-@slowromancapii@ GMRES solver considered in Sections 3.1.2 and 3.2.2, respectively.

We first consider the type-@slowromancapi@ Anderson acceleration in Sections 3.1.2. In this case, the SA scheme applies Anderson acceleration (42)–(44) to a preconditioned version of fixed-point problem (40). This preconditioned problem is derived from the SA scheme based on type-@slowromancapi@ Picard iteration (63) by first rewriting the correction process (63b) in the residual form as

[TABLE]

Plugging (63a) into (65) then gives

[TABLE]

which is equivalent to a standard Picard iteration update on the (preconditioned) fixed-point problem

[TABLE]

Here (67) is a preconditioned version of (40) with preconditioner $(\mathsf{I}+\mathsf{D}_{\mathsf{F}(\uprho)}^{-1}\mathsf{P}\mathsf{S})$. We obtain the SA scheme based on type-@slowromancapi@ Anderson acceleration by replacing each $\mathsf{G}_{1}$ in (42)–(44) with $\widetilde{\mathsf{G}}_{1}$.

For the type-@slowromancapii@ GMRES solver in Section 3.2.2, we follow a similar approach and use (64) to derive a preconditioned version of the linear system (53):

[TABLE]

Therefore, the SA scheme based on type-@slowromancapii@ GMRES solver computes $\mathsf{G}_{2}(\uprho^{(k)})$ by solving the preconditioned system (68).

4 Numerical results

The iterative solvers from the previous section are tested and compared on the one-dimensional silicon $n^{+}$–$n$–$n^{+}$ diode problem [11, 14, 17, 13, 31] with different collision frequencies. In Section 4.1, we describe the silicon diode problem and state the implementation details. In Section 4.2, we consider “single-scale” problems where the collision frequency is assumed to be constant throughout the spatial domain. Results from these single-scale tests illustrate characteristics of the different solvers. In Section 4.3, we consider more realistic “multiscale” problems with collision frequencies varying in the spatial domain. We first model the diode using the collision frequency specified in [14] which depends on the doping profile. We then consider a more challenging problem with collision frequency that changes more drastically over the spatial domain.

4.1 Silicon diode problem setup and implementation details

In the silicon (Si) $n^{+}$–$n$–$n^{+}$ diode problem, we simulate the electron movement in an one-dimension Si diode of length $L=0.6$\mathrm{\SIUnitSymbolMicro}\mathrm{m} with a bias voltage $V_{\textup{bias}}=1$\mathrm{V}. Here the electron charge, the effective electron mass, and the electric permittivity of the Si material are respectively given by $q_{e}=1.602\times 10^{-19}$\mathrm{C}, $m=2.368\times 10^{-31}$\mathrm{k}\mathrm{g}, $\varepsilon_{p}=1.034\times 10^{-10}$\nicefrac{{F}}{{m}}. The Boltzmann constant $k_{\textup{B}}=1.38\times 10^{-23}$\nicefrac{{J}}{{K}}, and the lattice temperature $T=300$\mathrm{K}$$. As in [14], for $x\in[0,0.6]($\mathrm{\SIUnitSymbolMicro}\mathrm{m}$)$, the doping profile is

[TABLE]

With these data, we apply the following scaling: $x_{0}=10^{-6}$\mathrm{m}, $[\Phi]=1$\mathrm{V}, and $D_{0}=10^{21}$\mathrm{m}^{-3}$$. We set $\zeta=1$ in (4b) (low-field scaling). The velocities in (5) take values

[TABLE]

We set the reference velocity to be $v_{0}=\max\{\Theta^{\nicefrac{{1}}{{2}}},B_{0},C_{0}\}=C_{0}$. Thus, the ratios in (4) are $\alpha=0.129$, $\beta=0.803$, and $\gamma=1$. The nondimensional doping profile on $x\in[0,0.6]$ is

[TABLE]

In the numerical simulations, we impose smooth transitions into the doping profile as in [31]. These smooth transitions are constructed by cubic splines, and the transition regions are of width 0.04, centered at 0.1 and 0.5.

The initial condition and the incoming boundary data for the kinetic equation (4a) are respectively given by $f(t_{0},x,v)=D(x)M_{\alpha^{2}}(v)$ and

[TABLE]

Without loss of generality, we let $\Phi_{0}=\Phi(t,0)=0$\mathrm{V} in the Poisson equation ([4b](#S2.E4.2)). The boundary data for ([4b](#S2.E4.2)) then become $\Phi(t,0)=0$\mathrm{V} and $\Phi(t,L)=V_{\textup{bias}}$. Since $[\Phi]=1$\mathrm{V}$$, the scaled boundary data are $\Phi(t,0)=0$ and $\Phi(t,L)=1$.

The computation is performed on a truncated domain $\mathcal{D}:=[0,0.6]\times[-2,2]$. The velocity space is truncated from $\mathbb{R}$ to $[-2,2]$, i.e., $v_{\max}=2$, since the value of Maxwellian $M_{\alpha^{2}}(v)$ is smaller than the machine precision when $|v|>2$. We discretize $\mathcal{D}$ into $200\times 50$ uniform rectangular elements of size $\Delta x\times\Delta v$, and solve the semiconductor model (4) from initial time $t_{0}=0$ to final time $t_{\textup{f}}=0.5$, at which point the system is essentially in steady state. The kinetic equation (4a) is solved using the fast sweeping algorithm detailed in Section 2.2, the Poisson equation (4b) is solved via the continuous Galerkin method discussed in Section 2.3, and the drift-diffusion equation (61), which is involved in the synthetic acceleration procedure, is solved by the direct discontinuous Galerkin method in Section 2.4. The Poisson equation and the drift-diffusion equation are both solved on a uniform mesh with $200$ elements on $[0,0.6]$. We choose the parameter $m=3$ in Anderson acceleration (42)–(44). The relative tolerances for the type-@slowromancapi@ and type-@slowromancapii@ iterative solvers are set to $10^{-8}$, while the relative tolerance in the GMRES solver for solving (25) in the fast sweeping algorithm is set to $10^{-10}$. The lower tolerance on the fast sweeping GMRES solver is due to the fact that it is used to evaluate $\mathsf{L}^{-1}_{\mathsf{E}}$, which is a fundamental building block in both type-@slowromancapi@ and type-@slowromancapii@ solvers. The iterative solvers are terminated once the relative residual is below the set tolerance, or when the solvers reach the maximum allowed number of iterations, which is set to $10,000$ for all solvers. This number is set to be large solely for studying the behavior of different iterative solvers. For practical applications, the maximum allowed number of iterations should be set much lower.

In the following sections, we consider the diode problem described in this subsection with various collision frequencies and compare the performance of the eight iterative solvers introduced in Section 3. We also make a formal efficiency comparison of the proposed implicit scheme to standard explicit schemes and standard implicit-explicit (IMEX) asymptotic preserving (AP) schemes. In all tests, the implicit time step is chosen to be $\Delta t=\Delta x$. We note that the CFL condition for standard explicit schemes takes the form $\Delta t\leq\min\{C_{1}\delta{\Delta x}/{{v}_{\max}},C_{2}{\delta\epsilon}/{{\omega}_{\max}}\}$ and the CFL condition for standard IMEX-AP schemes takes the form $\Delta t\leq\max\{C_{3}\delta{\Delta x}/{{v}_{\max}},C_{4}\delta{{\omega}_{\min}}{\Delta x^{2}}/\epsilon\}$ [33], where ${v}_{\max}$, ${v}_{\min}$, ${\omega}_{\max}$ and ${\omega}_{\min}$ are the maximum and minimum values of the velocity and collision frequency over the spatial domain, respectively. For simplicity, we assume that the $O(1)$ constants $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ are all equal to one. Here we do not consider the time step restrictions associated to $\Delta v$, since $\Delta x$ is much smaller than $\Delta v$. Also, the condition $\Delta t\leq{\delta\epsilon}/{{\omega}_{\max}}$ for standard explicit schemes is never active in problems tested in this section due to the small value of $\Delta x$. In the formal comparison, we assume that each explicit or IMEX step roughly requires the same computation time as each sweeping iteration in an implicit time steps. Thus for each iterative solver, we compute the ratio

[TABLE]

as an efficiency indicator of the proposed scheme. We also note that parallelization of the explicit or IMEX updates is often possible, while the sweeping procedure in the implicit scheme requires serial implementation. We do not take this fact into account in the comparison.

4.2 Single-scale test

In this section, we test the iterative solvers on problems with constant collision frequency on the spatial domain. We first consider the low collision case, where the electron mobility is approximated (see [14]) by

[TABLE]

with $D(x)=2\times 10^{21}$\mathrm{m}^{-3}, which determines the collision frequency $\omega(x)=\frac{q_{e}}{m\mu_{\textup{Si}}(x)}$. We scale the collision frequency by $\omega_{0}=5.114\times 10^{12}$\mathrm{s}^{-1}, and the nondimensional collision frequency is $\omega(x)=1$. The Knudsen number is $\epsilon=\frac{v_{0}}{x_{0}\omega_{0}}=0.200$ and we choose $\delta=\epsilon=0.200$.

Table 1 reports the iteration counts and computation time for each iterative solver. Here the column “FP” gives the iteration counts for solving the fixed-point problems (40) and (47), respectively. For type-@slowromancapii@ methods, the column “LS” reports the iteration counts for solving the linear system (48) when evaluating $\mathsf{G}_{2}$ in (47). The column “SW” gives the iteration counts for the GMRES solver for solving (25) in the fast sweeping algorithm when computing $\mathsf{L}^{-1}_{\mathsf{E}}$. As for the solvers, “PI” and “AA” stand for Picard iteration and Anderson acceleration, respectively.

From the results reported in Table 1, we first observe that the accelerated solvers are slower than the unaccelerated ones. This result is to be expected, since the system is away from the drift-diffusion limit due to the relatively small collision frequency. We also observe that the type-@slowromancapi@ solvers are faster than the type-@slowromancapii@ solvers. We conclude that for the type-@slowromancapii@ solvers, the additional computation cost of solving the linear system (48) when evaluating $\mathsf{G}_{2}$ outweighs the benefit of the more accurate updates for the fixed-point problem As we expected, the Krylov-type solvers (Anderson acceleration and GMRES) converge in fewer iterations than Picard iteration does. For type-@slowromancapi@ solvers, this results in less computation time, while for type-@slowromancapii@ solvers, the higher computation time per iteration of GMRES makes it slower than Picard iteration.

For this problem, the explicit and IMEX-AP CFL conditions both take the form $\Delta t\leq\delta\Delta x/v_{\max}$, which results in time steps that are 10x smaller than the implicit time step $\Delta t=\Delta x$ used in the test. For each tested iterative solver, the value of the efficient indicator $\mathcal{X}$ defined in (73) is reported in Table 2. Here Iter${}_{\text{SW}}$ denotes the total number of sweeping iterations required in one implicit step. From Table 2, the proposed scheme with type-@slowromancapi@ solvers results in $\mathcal{X}>1$ when comparing to both the explicit and IMEX-AP schemes. When using type-@slowromancapii@ solvers, $\mathcal{X}<1$ as the higher iteration counts outweigh the benefit of larger implicit time steps.

We next test the iterative solvers on problems with large collision frequency that is 100x of the one in the previous problem, and we scale the collision frequency by $\omega_{0}=5.114\times 10^{14}$\mathrm{s}^{-1}$$ so that the nondimensional collision frequency is still $\omega(x)=1$. The Knudsen number is then $\epsilon=\frac{v_{0}}{x_{0}\omega_{0}}=0.002$ and again we choose $\delta=\epsilon=0.002$. We still choose the time step to be $\Delta t=\Delta x$. For this problem, the system is close to the drift-diffusion limit due to the large collision frequency. Thus we expect that the drift-diffusion based synthetic acceleration would provide sufficient accurate corrections and result in faster convergence for the iterative solvers. The iteration counts and computation time for this problem are reported in Table 3. From these results, we observe that SA schemes indeed require fewer number of iteration and speed up the base iterative solvers from 1.5x to 8.7x. We also note that the computation time per iteration remains roughly the same, which implies that the time spent on computing the drift-diffusion correction term is essentially negligible.

For this problem, the standard explicit time step satisfies $\Delta t\leq\delta{\Delta x}/{v_{\max}}$, which is 1000x smaller than the implicit time step $\Delta t=\Delta x$. The standard IMEX-AP time step satisfies $\Delta t\leq\delta{\omega_{\min}}{\Delta x^{2}}/\epsilon$, which is 333x smaller than the implicit time step. From Table 4, in most cases, the proposed scheme results in $\mathcal{X}>1$ when comparing to both the explicit and IMEX-AP schemes. The exceptions are the unaccelerated type-@slowromancapi@ PI solver and the type-@slowromancapii@ AA/PI solvers.

4.3 Multiscale test

In this section, we test the solvers on multiscale problems with collision frequencies varying in the spatial domain. The first multiscale problem is the “standard” silicon diode problem from [14, 31]. Here the collision frequency is determined by the approximated electron mobility based on the doping profile. Specifically, with the approximate formula $\mu_{\textup{Si}}(x)=0.0088+\frac{0.1793\times 10^{23}}{1.4320\times 10^{23}+D(x)}$ at $T=300$\mathrm{K}$$ and the doping profile $D(x)$ in (69), the electron mobility is given by

[TABLE]

The nondimensional collision frequency on $x\in[0,0.6]$ is then

[TABLE]

with the scaling $\omega_{0}=1.843\times 10^{13}$\mathrm{s}^{-1}$$. Here the Knudsen number $\epsilon=\frac{v_{0}}{x_{0}\omega_{0}}=0.056$, and we choose $\delta=\epsilon=0.056$. As mentioned in Section 4.1, the doping profile used in the numerical tests includes the smooth transitions as in [31]. The collision frequency at these transition regions are computed from the smooth doping profile using the approximate formula given above.

Figure 1 shows the scaled doping profile $D(x)$ with smooth transitions, Maxwellian $M_{\alpha^{2}}(v)$, and collision frequency $\omega(x)$ in this silicon diode problem. Figure 2 illustrate the electron concentration $\rho$, the electric field $E$, and the full electron distribution $f$ at the final, steady-state time $t_{\textup{f}}=0.5$. As expected, all tested solvers give identical results up to the set numerical tolerance, and these results agree with the those reported in other places in the literatures, such as [11, 14].

Table 5 reports the iteration counts and computation time for each of the iterative solvers. The system is away from the drift-diffusion limit with the given collision frequency. Hence the results in Table 5 are similar to the ones in Table 1 and the drift-diffusion synthetic acceleration does not accelerate the convergence of the iterative solvers. The type-@slowromancapi@ solvers still outperform the type-@slowromancapii@ solvers for this test.

Here the implicit time step $\Delta t=\Delta x$ is roughly 35x larger than the standard explicit and IMEX-AP time steps $\Delta t\leq\delta{\Delta x}/{v_{\max}}$. We observe in Table 6 that the implicit scheme generally gives $\mathcal{X}>1$, except when the accelerated type-@slowromancapii@ AA/GMRES solver is used.

We next consider another multiscale problem with stronger variation in the collision frequency. Specifically, the scaled collision frequency considered in this problem is

[TABLE]

with $\omega_{0}=5.114\times 10^{14}$\mathrm{s}^{-1}$$ and $\delta=\epsilon=\frac{v_{0}}{x_{0}\omega_{0}}=0.002$.

The iteration counts and computation time for this problem are reported in Table 7. These results show that for this multiscale problem, synthetic acceleration with the drift-diffusion model speeds up the convergence of most iterative solvers by roughly a factor of two. However, we observe diverging residuals (denoted as $\infty$) when using the accelerated type-@slowromancapi@ PI scheme, which indicates that applying the drift-diffusion based synthetic acceleration on multiscale problems may lead to unstable schemes that give divergent results. This observation is related to a known deficiency of diffusion synthetic acceleration (DSA). Specifically, it has been reported in [25] and analyzed in [47] that DSA becomes unstable and gives divergent results when applied to highly collisional problems. To guarantee convergence, the spatial discretization of the diffusion equation has to be consistent to the discretization of the original transport equation. We refer the reader to [2] and [37] for a complete discussion.

Finally, we note that the standard explicit and IMEX-AP time steps both satisfy $\Delta t\leq\delta\Delta x/v_{\max}$, which are 1000x smaller than the implicit time step $\Delta t=\Delta x$. As shown in Table 8, the efficient indicator $\mathcal{X}>1$ for all iterative solvers that lead ti a convergent implicit scheme. Among these iterative solvers, the most efficient one is the accelerated type-@slowromancapi@ AA solver.

5 Conclusions and discussion

We have proposed a fully implicit numerical scheme for solving Boltzmann-Poisson systems with an approximate collision operator that describes linear relaxation. At each implicit time step, the updated solution comes from a nonlinear fixed-point problem. We have formulated Boltzmann-Poisson systems as two types of fixed-point problems: the type-@slowromancapi@ problems that are solved in a single strongly coupled iterative loop and the type-@slowromancapii@ problems that require two nested iterative loops. We have applied Picard iteration and Anderson acceleration to solve the type-@slowromancapi@ problems. For the type-@slowromancapii@ problems, the outer loop of the nested iterative procedure is performed by Anderson acceleration while the inner loop uses either Picard iteration or GMRES. The performance of these iterative solvers and their synthetically accelerated (SA) variants is compared on several scaled versions of a standard silicon diode problem. Numerical results show that (i) solving the type-@slowromancapi@ fixed-point problems requires fewer iterations than solving the type-@slowromancapii@ problems and (ii) Anderson acceleration is more efficient than Picard iteration on type-@slowromancapi@ problems, in terms of both iteration counts and computation time. These results also confirm that SA schemes with a drift-diffusion model converge faster than the standard iterative solvers when the system is near the drift-diffusion limit (when collision frequency is large). For systems away from this limit, there is no observed advantage in using these accelerated schemes. To address this issue, some potential approaches to be investigated in the future include (i) modifying the drift-diffusion SA schemes to incorporate the boundary conditions as proposed in [54] for a neutron transport problem and (ii) deriving SA schemes based on low-order/low-cost approximations other than the drift-diffusion equation. For (ii), possible candidates of such approximation include the $S_{2}$-SA scheme considered in [42, 9] and the more general, two-level multigrid algorithms as the one found in [6, Section 6.2].

Another potential future work is to apply the hybrid schemes [19, 29] proposed for linear transport equations. These schemes decompose the transport solution into collisional and non-collisional components. At each time step, the hybrid schemes use cheap, low-resolution approximations for the collisional component and reserve expensive, high resolution approximations for the the non-collisional component. Since the primary difficulty on solving the Boltzmann-Poisson system is the nonlinear coupling between the electric field and collision term, we expect that the hybrid approach would significantly lower the computation cost while producing quality solutions that are comparable to the ones given by an uniformly high-resolution solver.