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}$$