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?
Thanks
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?
