Inverting a pressure matrix for fluid simulation

Question

I am implementing a fluid simulator as my numerical methods course project and I have to compute pressure at each simulation step. Basically, that means solving an equation $Ap = b$, where $A$ is a sparse matrix: all its entries are zero except:

• Its main diagonal
• Diagonals before and after main diagonal
• For some $k$, the $k$th diagonal before and the $k$th diagonal after main

Somewhat like this:

$$\begin{bmatrix} a_{11} & a_{12} & 0 & \dots & a_{1k} & 0 & 0 & \dots & 0 \\ a_{21} & a_{22} & a_{23} & \dots & 0 & a_{2(k+1)} & 0 & \dots & 0 \\ 0 & a_{32} & a_{33} & \dots & 0 & 0 & a_{3(k+2)} & \dots & 0 \\ 0 & 0 & a_{43} & \dots & 0 & 0 & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ a_{k1} & 0 & 0 & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & a_{(k + 1)2} & 0 & \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & a_{(k + 2)3} & \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots\\ 0 & 0 & 0 & \dots & 0 & 0 & a_{(n-1)(n-2)} & a_{(n-1)(n-1)} & a_{(n - 1) n} \\ 0 & 0 & 0 & \dots & 0 & 0 & 0 & a_{n(n-1)} & a_{nn} \end{bmatrix}$$

Since I sumulate the fluid without free surface, the matrix is not changing during the simulation, so I want to invert it before the simulation and then just compute $A^{-1} b$ at each step. However, I believe there are faster algorithms for doing it than Gaussian elimination which could exploit sparsity and the structure of such a matrix.

So, my questions are:

1. Is there a commonly used name for such a matrix?
2. Are there any effective algorithms for inverting it? (something similar to tridiagonal matrix algorithm, maybe?)
3. If 2 is positive, are there any Java libraries doing that?

Answer 1

1. I suppose you are solving the so called Pressure Poisson Equation and that your $A$ is a discrete representative of a Laplacian.

2. Typically, you don't want to have the inverse computed and stored, cf. the comments to this question. The main reason is that the inverse will be a dense matrix. I rather suggest you compute a sparse LU factorization of $A$, so that $Ax=LUx=b$ can be solved for $x$ very fast. For the 2D Laplacian the factors $L$ and $U$ are known to be computable with little fill in.

3. Routines for the sparse LU decomposition are available here or as a part of high performant libraries like umfpack. However, I wasn't able to find a java implementation within 3 minutes of googling.

Finally, as your $A$ is typically positive definite, you can go with the CG algorithm which takes only a few lines of code to implement.

Answer 2

Your matrix is not "banded" in the sense of "banded" direct solvers, so don't bother with those. Multigrid is absolutely the best way to solve these "pressure Poisson" problems. There are lots of libraries and it's not difficult to implement a simple multigrid algorithm for structured grids. CG and direct solvers are fundamentally non-scalable, so the cost will increase superlinearly as you increase resolution.