I am working on a 1D drift-driffusion problem in a finite-difference (FD) approach.

I hade 3 equations per node ($3N$ in total): electron continuity $E_i$, Poisson $P_i$, hole continuity $H_i$. With corresponding variables: electron concentration $n_j$, potential $V_j$, hole concentration $p_j$.

By rearranging the problem has a banded matrix, the Jacobian $A$ takes the following form, with shape $(3N\times 3N)$

$$A = \begin{pmatrix} A_{1,1} & A_{1,2} & & & & & \\ A_{2,1} & A_{2,2} & A_{2,3} & & & & \\ & A_{3,2} & A_{3,3} & A_{3,4} & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & A_{N-2,N-3} & A_{N-2,N-2} & A_{N-2,N-1} & \\ & & & & A_{N-1,N-2} & A_{N-1,N-1} & A_{N-1,N} \\ & & & & & A_{N,N-1} & A_{N,N} \\ \end{pmatrix}$$


$$ A_{i, j} = \begin{pmatrix} {\partial E_i}/{\partial n_j} & {\partial E_i}/{\partial V_j} & {\partial E_i}/{\partial p_j} \\ {\partial P_i}/{\partial n_j} & {\partial P_i}/{\partial V_j} & {\partial P_i}/{\partial p_j} \\ {\partial H_i}/{\partial n_j} & {\partial H_i}/{\partial V_j} & {\partial H_i}/{\partial p_j} \end{pmatrix} $$

Owing to that data structure, I have been using scipy.linalg.solve_banded to solve for $X$ in

$$AX = b$$

where $b$ is a $(3N \times 1)$ shaped vector.

My problem comes from that I want to add an equation that controls the potential at one end (say on the left) based on the value of total fluxes of $n$ and $p$ at both ends.

As a result the additional equation depends on $(n_1, V_1, p_1, n_2, V_2, p_2)$ but also $(n_{N-1}, V_{N-1}, p_{N-1}, n_N, V_N, p_N)$.

Therefore, if I was to include this new equation in $A$ to form a $(3N+1 \times 3N+1)$ Jacobian $A'$, the bandwith would suddenly become much larger. Though I can optimize this by inserting the new equation in the middle row of $A$, the banded representation would remain very sparse.

My question is then: is there a way I can solve for $X'$ in $$A' X' = b'$$

while keeping the advantages of the banded structure of $A$?

EDIT addressing comments

Here is a visual example of the Jacobian before (left) and after reordering (right). The last column/line, highlighted in yellow, is the additional equation that fails to fit the banded structure. Elements in dark gray background are non-zero, elements in white are zeros, elements in light gray background are zero, but a future version of my physical could in principle make them non-zero so an ideal solution would not exclude that possibility.

visual sparsity representation (click to enlarge)

Disclaimer: I do not pretend to be the author of this ordering method. I have read it a few years ago in a software manual I believe ...

  • 1
    $\begingroup$ I think sparse direct solver should be able to handle efficiently. Try using scipy.sparse.linalg.spsolve? $\endgroup$ – piyush_sao Jan 23 '19 at 13:21
  • $\begingroup$ @zimbra314 Thank for your reply. In that case, does it matter that $A$ is reordered as shown or the natural (non banded) ordering will perform equally well? $\endgroup$ – Alexis Jan 23 '19 at 13:23
  • 5
    $\begingroup$ Sparse direct solver will reorder the matrix to reduce the amount of computation, before solving it. $\endgroup$ – piyush_sao Jan 23 '19 at 13:25
  • $\begingroup$ Can you show the sparsity structure of the matrix you have in the end? If you can decompose it into a 2x2 set of blocks where the large diagonal block is tridiagonal, then it may be worth computing the Schur complement. $\endgroup$ – Wolfgang Bangerth Jan 24 '19 at 16:55
  • 3
    $\begingroup$ Nice graphs of your matrix structure. But the bottom line is that the suggestions of @zimbra314 are the right approach. Sparse direct solvers are now mature and efficient to the point that band matrix solvers are largely relics of the 1960's and 1970s. I'm assuming your objective is to solve your equations rather that to experiment with matrix computations. $\endgroup$ – Bill Greene Jan 25 '19 at 0:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.