I'm trying to solve the fully coupled drift-diffusion system using Newton's Method. Although I eventually plan to potentially use a Jacobian-Free Newton-Krylov approach, this is still something that I want to implement and look into.

A bit of notation for the equations presented below:

$d_{ij}$ and $l_{ij}$ are control area/volume constants, not physical entities. The $i$ and $j$ indices are for the central point being considered ($i$), and all neighbouring points that are connected ($j$'s). $A_{\Omega_i}$ is the control area/volume associated to each vertex/node, and is also constant.

All quantities that are explicitly required during calculation are at node points. In the equations $B(\Delta_{ij}) \equiv B(\phi_i - \phi_j)$, the Bernoulli function, used as a smoothing method for the differences. This difference should technically exist between node points, but that's for another time.

Here are the equations: $\phi$ is electrostatic potential, $n$ and $p$ are carrier concentrations. These are the independent variables to solve for. Everything else is either a constant, or a function of one of these variables.

$$\sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) = -\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i}$$

$$\frac{\partial{n_i}}{\partial t}A_{\Omega_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$

$$\frac{\partial{p_i}}{\partial t}A_{\Omega_i} = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i}$$

In residual form, and assuming equilibrium conditions so all time derivatives go to 0, these equations are:

$$\Phi = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}}(\phi_i - \phi_j) +\frac{q}{\epsilon}(p_i - n_i - N_{D_i}^+ - N_{A_i}^-)A_{\Omega_i} = 0$$

$$N = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{n_iB(\Delta_{ij}) - n_j B(-\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$

$$P = D_p \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \{p_iB(-\Delta_{ij}) - p_j B(\Delta_{ij})\}\right] + (G_i - R_i)A_{\Omega_i} = 0$$

Now, my question is whether I have these Jacobian elements correct. Here is the Jacobian element formula that I've tried to come up with. Only a few equations have been derived, since the remainining elements are just trivial variations of these 3 formulations. This post is mostly about getting more experienced eyes to look at these equations, and potentially see something that I haven't quite seen yet. Are these formulations for the Jacobian elements correct, or is there something wrong with them that I'm missing?

$$\frac{\partial\Phi}{\partial\phi_i} = \sum_{j\ne i} \frac{d_{ij}}{l_{ij}} +\frac{q}{\epsilon}(\frac{\partial p_i}{\partial \phi_i} - \frac{\partial n_i}{\partial \phi_i})A_{\Omega_i}$$

$$\frac{\partial{N}}{\partial \phi_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} \left\{\frac{\partial{n_i}}{\partial \phi_i}B(\Delta_{ij}) + n_i \frac{\partial B(\Delta_{ij})}{\partial \phi_i})\right\}\right] + \frac{\partial(G_i - R_i)}{\partial \phi_i}A_{\Omega_i}$$

$$\frac{\partial N}{\partial n_i} = D_n \sum_{j \ne i} \left[\frac{d_{ij}}{l_{ij}} B(\Delta_{ij})\right] + \frac{\partial (G_i - R_i)}{\partial n_i}A_{\Omega_i}$$


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.