How to solve a problem with structure similar to a finite difference discretization of the 2D Poisson equation, but with non-symetric coefficients?

Question

Recently, I've been asking about methods to solve a finite difference discretization of the 2D Poisson equation (see here and here) of the form:

$$U_{i-1,j} + U_{i+1,j} -4U_{i,j} + U_{i,j-1} + U_{i,j+1} = f_{i,j}$$

The left hand side of this system is often referred to as the 5-point discrete Laplacian. Suppose that left hand side of the linear system is not quite the traditional five-point discrete Laplacian, but instead has the form:

$$\sum_{k,l} a^{ij}_{kl} U_{kl} = a^{ij}_{i-1,j}U_{i-1,j} + a^{ij}_{i+1,j}U_{i+1,j} -a^{ij}_{i,j}U_{i,j} + a^{ij}_{i,j-1}U_{i,j-1} + a^{ij}_{i,j+1}U_{i,j+1} = f_{i,j}$$

with non-negative off-diagonal ($a^{ij}_{k,l} \ge 0$ for all off-diagonal entries $(i,j) \ne (k,l)$) and zero row sums, $$ \sum_{kl} a^{ij}_{kl} = 0 $$ for all $i,j$. From these conditions, the matrix $A = (a^{ij}_{kl})$ is the negative of an $M$-matrix.

The problem is sparse and diagonally dominant (I would characterize it as "weakly diagonally dominant" because the diagonal term is exactly equal to the sum of the off diagonal elements of the same row). Since it exhibits a similar structure to the five-point Laplacian (except perhaps the possibility of being symmetric?), I figured that there may be some similarities between methods to solve this system of equations and those of the five-point discrete Laplacian. Specifically, I was thinking that one of the following may be true:

1. We can compact the sparsity of the matrix into a thinly banded system (since there are only five variables per equation) and use a direct solver.
2. We can use some sort of Fast Fourier Transform to solve it. (I think this one is highly unlikely though.)
3. We can use some sort of multigrid approach to solve the system iteratively.

In the end, the systems will be very large and a parallel implementation will be necessary. Any feedback on fast methods to approach this problem in parallel would be greatly appreciated.

Answer 1

1. Please take the advice given by the answers to your several previous questions and stop suggesting that banded direct solvers would be appropriate or this problem. They are not, and sparse direct solvers do not attempt to find low-bandwidth orderings, they choose orderings that reduce fill (like nested dissection).

2. FFT is not appropriate for the case of general coefficients.

3. Geometric or (depending on problem structure) algebraic multigrid should work well for this problem, provided you handle a few possibilities correctly. Note that if the coefficients represent extreme anisotropy, then the continuum equations become similar to a bunch of decoupled 1D problems, so isotropic coarsening cannot be expected to work (at least without line smoothers). If the equations represent transport, then you may have to use upwinded interpolants and/or relaxation. Note that second order linear discretizations of transport produce oscillations, so you would likely be using the upwinded form already. If the coefficients are highly variable, you may need to construct low-energy interpolants.

In general, it is much better to say which equations you are solving in which regime with which discretization than it is to make superficial observations about the matrix structure.