Solving the Poisson equation in 3D using Python

Question

I'm implementing a Python code where I need to solve the following Poisson equation as one of the steps:

$$\nabla^2p=f(\mathbf r)$$

I am using a 3D rectangular grid with ~100 points on each direction. The finite difference equation I get is

$$p_{i, j, k}=\frac{k^2l^2(p_{i+1, j, k}+p_{i-1, j, k})+h^2l^2(p_{i, j+1, k}+p_{i, j-1, k})+h^2k^2(p_{i, j, k+1}+p_{i, j, k-1})+h^2k^2l^2 f_{i, j, k}}{2(h^2+k^2+l^2)}$$

where $h, k, l$ are the mesh spacings in the $x, y, z$ directions, respectively. This can obviously be expressed as a matrix equation and in principle I could just invert the matrix acting on the $p$s. However it would be a ~$10^6\times10^6$ matrix, which would require me to have ~$10^{13}$ bits (assuming I use 16 bit floating point numbers) of memory, which comes to about 1TB at the very least, which is ~80 times more than what I have installed in my computer. There must be a better way to solve this, but I'm not proficient in numerical methods for PDEs.

I saw this article by Profs. Pöplau and Portratz that deals with a generalized version of the code I'm trying to implement (they use a variable mesh spacing near the boundaries to simulate arbitrary shaped boundaries). However, the code that is presented in the article is incomplete, and I've been unable to reverse-engineer it to use it myself. I tried sending an email to Prof. Pöplau, who is the correspondence author, but the address no longer exists.

Could you point me towards literature about methods similar to Pöplau and Portratz's that I can use for my code? Or is there another already implemented library/package/code I could adapt or use to solve the Poisson equation numerically?

Note: The boundary conditions I'm using could be expressed as either Dirichlet or Von Neumann. However, periodic boundary conditions would not work, because the function I'm interested in is not necessarily periodic.

Answer 1

Too long for a comment.

1. It is a bad idea to denote one index and one mesh spacing, both, by $\color{red}{k}$. Better use $h_i,h_j,h_k$ for the mesh spacings.

2. From

\begin{align} \frac{\partial p}{\partial x}&\approx\frac{p_{i+1,j,k}-p_{i,j,k}}{h_i}\\[3mm] \frac{\partial^2 p}{\partial x^2}&\approx\frac{p_{i+1,j,k}-2p_{i,j,k}+p_{i-1,j,k}}{h_i^2}\\ \end{align} I get \begin{align} \nabla^2p&\approx\frac{p_{i+1,j,k}-2p_{i,j,k}+p_{i-1,j,k}}{h_i^2}\\ &+\frac{p_{i,j+1,k}-2p_{i,j,k}+p_{i,j-1,k}}{h_j^2}\\ &+\frac{p_{i,j,k+1}-2p_{i,j,k}+p_{i,j,k-1}}{h_k^2}\,. \end{align} What has this to do with your expression ?

3. When I imagine a cube with $100$ cells in each direction I see a total number of $10^6$ grid points, not a $10^6\times 10^6$ matrix.