I'm trying to simulate basic semiconductor models for pedagogical purposes--starting from the Drift-diffusion model. Although I don't want to use an off-the-shelf semiconductor simulator--I'll be learning other (common, recent or obscure) models, I do want to use an off-the-shelf PDE solver.

But even for the simple 1D case, the drift-diffusion model consists of a number of coupled nonlinear PDEs:

Current density equations $$J_n = q n(x) \mu_n E(x) + qD_n \nabla n$$ $$J_p = q p(x) \mu_p E(x) + qD_p \nabla p$$

Continuity equation $$\frac{\partial{n}}{\partial{t}} = \frac{1}{q} \nabla \cdot J_n + U_n $$ $$\frac{\partial{p}}{\partial{t}} = \frac{1}{q} \nabla \cdot J_p + U_p $$

Poisson equation $$\nabla \cdot (\epsilon \nabla V) = -(p - n + N_D^+ - N_A^-) $$

and a number of boundary conditions.

I have tried some python FEM solvers, FEniCS/Dolfin and SfePy, but with no luck, due to being unable to formulate them in the weak variational form with test functions.

There is of course the option of implementing the numerical solution from scratch but I haven't studied FEM/Numerical in depth yet, so I hope it's not my only option as I don't want to be overwhelmed with numerical issues.

So is there a package (pref. open source) that would take these equations, in that form, and solve them? Or perhaps is the variational form required by the tools is not as hard? In any case, what are my options?


Edit: Attempt of formulating the weak variational form for FEniCS/Dolfin or SfePy

Using three PDEs (Poisson + two continuity equations with J substituted), we are looking for V, n, and p. The Poisson equation (using a test function $u_V$) is straight forward. I'm having difficulty, however, with the continuity equations.

The second PDE (strong form) $$\frac{\partial{n}}{\partial{t}} = \nabla \cdot (C_1 n \nabla V +C_2\nabla n) + U $$ where $C_1, C_2$ are constants, $U, n, p, V$ are scalar functions

Let $f_n$ denote a test function for the second PDE. Then

$$\int_{\Omega}f_n\frac{n - n_1}{\Delta t}d\Omega - C_1 \int_{\Omega}f_n \nabla \cdot ( n \nabla V) d\Omega - C_2 \int_{\Omega}f_n \nabla^2 n d\Omega - \int_{\Omega}f_n U d\Omega$$

Especially worrying is the integral: $$C_1 \int_{\Omega} f_n \nabla \cdot ( n \nabla V) d\Omega$$

But $\mathbf{\nabla V}$ is a vector, and $V, u_n, n$ are scalars. Then using the identity $\nabla \cdot \phi \mathbf{A} = \mathbf{A} \cdot \mathbf{\nabla \phi} + \phi \nabla \cdot \mathbf{A}$

$$C_1 \int_{\Omega} f_n \nabla \cdot ( n \nabla V) d\Omega = C_1 \int_{\Omega} f_n (\nabla V \cdot \nabla n) + C_1 \int_{\Omega} f_n n \nabla \cdot \nabla V$$

Since V is solved by Poisson equation, can we use the recently computed value as allowed in software Dolfin/FEniCS and simplify how we treat V in this second coupled equation? These sort of techniques work in while discretizing (e.g. Gummel, ...), which I don't do in these ready solvers!

Also the boundary conditions are given in terms of $J_n$ not $n$, how do you implement this? Should I solved for the five variables $J_n, J_p, n, p, V$, even though $J_n$ is determined by V and n?

  • 1
    $\begingroup$ Why can't you write down their weak forms? $\endgroup$
    – Bill Barth
    Commented Dec 6, 2012 at 2:00
  • $\begingroup$ @BillBarth I edited my question, please have a look. Thanks. $\endgroup$
    – Weaam
    Commented Dec 7, 2012 at 7:18
  • 2
    $\begingroup$ Your integration by parts is wrong. Check the formula, there are signs missing, you have more derivatives on the right than on the left, and you forgot about the boundary integral. $\endgroup$ Commented Dec 7, 2012 at 10:01
  • $\begingroup$ Also, is there a reason you are using a dot product to represent multiplication by $u_n$? It's a scalar, right? $\endgroup$
    – Bill Barth
    Commented Dec 7, 2012 at 15:26
  • $\begingroup$ Yes, I should have been more careful. Please check my edit, especially my question regarding how we treat V since it should have been solved already by the previous PDE. Does this have any effect on the variational form? Thank you. $\endgroup$
    – Weaam
    Commented Dec 7, 2012 at 18:42

1 Answer 1


The Scharfetter-Gummel (SG) formulation is commonly used to solve the current density equations. This is a special formulation that overcomes difficulties in solving the nonlinear dependence between potential and current density.

A standard text discussing how these equations using box integration methods is in this book: Selberherr, S., Analysis and simulation of semiconductor devices. Springer-Verlag 1984

This type of simulation called Technology Computer Aided Design (TCAD). As opposed to the Finite Element Method (FEM), Finite Volume Method (FVM) is used to calculate the currents. This is since it fits into the SG formulation that has been shown (by practioners of this method) to work when solving the current density equations.

If you want to solve this using generalized PDE's, COMSOL has a Semiconductor Module which solves this problem using a hybrid FEM/FVM method.

In addition, commercial and open source TCAD simulators are listed here: http://www.tcadcentral.com

To my knowledge generalized-PDE TCAD solvers are DEVSIM, FLOOPS, PROPHET. The commercial tools tend to have most of the physical equations are hard coded in a compiled language such as C++.

  • $\begingroup$ I apologize for the extremely late reply. I realized such a direct application of DD (even with SG) was quite unstable (my implementation in Fenics at least), thus I abandoned it. In a later VLSI course, I indeed used Comsol and TCAD tools. Thank you for your comprehensive answer. $\endgroup$
    – Weaam
    Commented Sep 20, 2015 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.