Fast Poisson solver (with Dirichlet BC zero) on a *truncated* Cartesian 3D grid

Question

I find myself in the position of having to solve

$-\Delta u = f$ on a subset of Cartesian grid points that don't necessarily form a cuboid domain subject to a homogenious Dirichlet boundary condition ($u = 0$). The Laplacian filter I require is the standard 7-point stencil for 3D.

Now, for a cuboid domain, I know this can be efficiently solved using a type-I discrete sine transform (DST) since the DST diagonalizes the equation system.

I wonder, whether it's possible to leverage the fast DST for domains that are not cuboid. It's still a Cartesian grid, but only a subset which is not necessarily convex and might include holes.

Are there fast algorithms that can exploit the Cartesian grid structure? Or do I have to use a more general solver (like an algebraic multigrid preconditioned conjugate gradient method or something like this)?

Maybe there is an efficient method to reduce my problem to another with an emcompassing cuboid domain so I still get to use the DST? As far as I can tell, I could introduce source terms in $f$ on my actual boundary so that the DST-based approach will reconstruct a boundary of zero. But the problem with this is that determining those source terms seems to require solving a big and "ugly" linear equation system (dense matrix which is not trivially diagonalizable).

What would be your recommended solver for this problem? The number of unknowns are in the tens of millions.

Answer 1

The short answer is, no; you can't use the DST approach for a case with general geometry or boundary conditions.

The best way to understand this is to consider WHY the DST approach works for the "rectangular" case. For this case, we happen to know that the complete set of eigenvectors of the differential equation are products of sin functions. We can then conveniently express the solution as a sum of the eigenvectors. The DST is a fast way to compute the eigenvectors in this case.

For a general geometry, the eigenvectors are not simple sin functions and we have no cheap way to compute them; that makes an eigenvector basis unattractive.

Since you say your grid is Cartesian, it seems like algebraic multigrid would be an attractive option and is theoretically faster than using FFT. As shown on page 2 of these notes by Demmel, FFT has performance proportional to $N log N$ while multigrid is order $N$, where $N$ is the number of equations.