Trouble with boundary conditions in Newton-Raphson method for Coupled PDEs

[nabeem]

Hi all,

This pertains to a pretty common method of simulating semiconductor devices, but unfortunately I've looked through tons of sources that have been unable to answer my question:

I'm currently working on a 1D device simulator in MATLAB that uses a Newton-Raphson iteration to solve the coupled Poisson's equation, electron continuity equation, and hole continuity equation in discretized form. The device is divided into N nodes (N-2 interior nodes and 1 node each at the left and right boundaries). I'm trying to solve for the electrical potential, electron quasi-fermi level, and hole quasi-fermi level profiles along the device. So essentially I'm trying to solve for the variables:
$$\mathbf{V}(x) = ( V(1), V(2), ... , V(N) )$$
$$\mathbf{\phi_n}(x) = ( \phi_n(1), \phi_n(2), ..., \phi_n(N) )$$
$$\mathbf{\phi_p}(x) = ( \phi_p(1), \phi_p(2), ..., \phi_p(N) )$$

which are the electrical potential, the electron quasi-fermi level, and the hole quasi-fermi level, respectively. I've set up the poisson equation, the electron continuity equation, and the hole continuity equation in discretized, residual form, so I'm solving the system:

$$\mathbf{F} = ( \mathbf{F1(\mathbf{V},\mathbf{\phi_n},\mathbf{\phi_p}}), \mathbf{F2(\mathbf{V},\mathbf{\phi_n},\mathbf{\phi_p}}), \mathbf{F3(\mathbf{V},\mathbf{\phi_n},\mathbf{\phi_p}}) ) = 0$$.

where F1 is the Poisson equation, F2 is the electron continuity equation, and F3 is the hole continuity equation, at each point 1, ..., N. I've set up a Newton-Raphson iteration to solve these equations. The Jacobian is a 3N x 3N sparse, banded, block matrix that can be computed analytically. So the first row of the Jacobian consists of the blocks $$\frac{\partial{F1}}{\partial{V}}, \frac{\partial{F1}}{\partial{\phi_n}}, \frac{\partial{F1}}{\partial{\phi_p}}$$, the second row the same thing except with F2, you get the picture.

However, my problem is related to the implementation of the boundary conditions. I'm using simple Dirichlet boundary conditions with V(1), V(N), fn(1), fn(N), fp(1), and fp(N) set to fixed values. Naturally, I set up the first and last equations in each solution block to correspond to the fact that the correction values for V(1), V(N), etc are 0. So for example, in each block, I have the (1,1) and the (N,N) elements set to 1, with the rest of the elements in those two rows set to 0, and then set F1(1) and F1(N), etc. to 0. Unfortunately, since these rows get repeated in each block, I obviously end up having a singular Jacobian, and MATLAB spits out a correction vector of NaN's. When I try to use the method with just one variable, V(x), I don't have a singular matrix and the iteration converges and I get a solution and everything, but at the same time the solution still acts strangely at the boundaries (for example, staying flat through the device but then shooting up abruptly to the boundary values), so I'm thinking that even just implementing the boundary conditions in this manner doesn't/wouldn't give me an entirely accurate picture.

I've tried to simply eliminate the first and last rows and columns of each block altogether and use the Newton iteration to solve only for the INTERIOR points, but for whatever reason the iteration then takes a painfully slow time to converge (the norm of the correction vector decreases by something like 0.0001 per step, which leads me to think that it causes the iteration to just stop functioning properly altogether instead of actually having slow convergence rates- typically, the Newton iteration should take between 3 and 9 steps to converge in these situations). I've also tried to actually fill in meaningful values for the (1,1), (1,2) and (N,N-1) and (N,N) elements in each block, but this becomes a problem when trying to fill in corresponding values for F(1) and F(N) since they both rely on discretized second derivatives, and thus values of each variable that exist outside the solution domain.

Does anyone know what's going wrong/how I should go about fixing this? It's clear that people have been able to implement a Newton-Raphson method in this situation successfully before, but I've found nothing that has been able to point me to my error. I don't want to bother you guys with the code, since I don't think that's where the problem is, but if it's necessary to see it I can post it.

Thanks for the help.

 

[jvicens]

Hi there,

Thank you for sharing your question and background information with us. It seems like you have a solid understanding of the problem you are trying to solve and the methods you are using. From your description, it is difficult to pinpoint exactly where the issue may be occurring, but I have a few suggestions that may help you in your troubleshooting process.

First, have you tried using a different method for implementing the boundary conditions? Instead of setting the first and last rows and columns to 1 and the rest to 0, have you considered using a ghost node approach? This involves creating two additional nodes at the boundaries and using them to enforce the boundary conditions. This method may help to avoid the singularity issue you are encountering with the Jacobian matrix.

Second, have you checked your discretization scheme and boundary conditions to ensure they are consistent? In some cases, inconsistencies between the discretization and boundary conditions can lead to numerical instabilities and slow convergence rates. Double checking these aspects may help to identify any errors that may be causing issues in your simulation.

Finally, have you tried using a different solver or a different programming language? Sometimes, using a different solver or programming language can help to identify any issues in the implementation and lead to a successful solution.

I hope these suggestions are helpful to you. If you need further assistance, please don't hesitate to reach out. Keep up the good work and good luck with your simulation!