0
$\begingroup$

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.

Improve this question
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$
  • $\begingroup$ I'm aware of why the DST works for the cuboid case. However, I don't think you can easily invert such a statement (inverse error fallacy). The DST would still be applicable if you found an efficient way to turn problem A (on a non-cuboid domain) into problem B (with cuboid domain). I actually already did something like this where this problem reduction step is very cheap. But my feeling is that for this current problem (with homogenious Dirichlet BC), that "problem reduction step" would be too costly and not worth the hassle. $\endgroup$ – sellibitze May 2 '19 at 8:45
  • $\begingroup$ Yes, I want to solve a linear equation system that you would get by using the 7-point stencil (common for finite difference discretization). But the interesting bit is that this is not an "approximation" to my original problem. It's exactly what I want to solve. :) $\endgroup$ – sellibitze May 2 '19 at 8:48
  • $\begingroup$ "tens of millions" should give you an idea of the size of the problem which is surely going to inform the choice of the solver, isn't it? $\endgroup$ – sellibitze May 2 '19 at 8:49
  • $\begingroup$ What would be your recommendation for an efficient solver? $\endgroup$ – sellibitze May 2 '19 at 8:49
  • $\begingroup$ My question about problem size was actually about WHY such a fine grid is required? You say this grid gives the "exact" solution to your problem. Why is an approximate solution from a coarser grid unacceptable? If you simply MUST solve this system of equations, since you have a regular grid, why don't you use standard multigrid? $\endgroup$ – Bill Greene May 2 '19 at 10:15

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.