Jacobian Elements for Coupled Drift-Diffusion System using Vertex-Centered Finite Volume
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. 232
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? 65