UNsymmetric Banded Matrix Solver

[inkant]

Hello everyone,

I am working on a 2D CFD (Computational Fluid Dynamics) code. Let me give you a general view of the problem.

A general 2D discretization leads to an equation like this:

-a(i-1,j) -a(i+1,j) + a(i,j) -a(i,j-1) -a(i,j+1) = b (similar to a poisson eqn in 2D).

When solved just like this, this equation gives rise to a Penta-Diagonal matrix (a tridiagonal matrix with two more diagonals situated far away from the main diagonal), with large bandwidth. A way to solve this problem is to choose marching direction in either i or j.

-a(i-1,j) -a(i+1,j) + a(i,j) = b +a(i,j-1) +a(i,j+1) (both RHS a()'s taken from previous iteration, thus known).

The above equation is tridiagonal. which can be solved with the speed of LIGHT ! using TDMA.

Then same thing is done with the other direction:

a(i,j) -a(i,j-1) -a(i,j+1) = b + a(i-1,j) + a(i+1,j)

This process is iterative and takes many iterations. The order of complexity is is O(n). The process becomes non-converging when the equation set is stiff. (because of stiff chemistry).

My question to the community is if there are faster solvers on similar matrices?

1. LU factorization in my experience may not be helpful as there is only 1 RHS. (available in MKL)
2. I looked at xGByyy implementations in LAPACK, but I am not sure whether they would use this special nature of matrix (sparsity between the last diagonal)
3. Can Multigrid method be faster? The order of complexity for that is also O(n). (available in MKL ? no ?)
4. How will be Conjugate gradient methods on this type of system ? (also available in MKL)
5. I have also viewed CVODE and DASSL. But these are large solvers which may take me a long time to figure out how to implement in PDEs.

I would deeply appreciate any advice from the community members, since it is a research problem not really an MKL issue.

Best Regards,

Inkant.

[j_stults]

Quoting - inkant

The process becomes non-converging when the equation set is stiff. (because of stiff chemistry).
...
How will be Conjugate gradient methods on this type of system ? (also available in MKL)

I'd suggest conjugate gradient, it can be very fast, especially with an appropriate preconditioner (which it sounds like you need anyway because of the chemistry).

Here's the convergence of a couple iterative schemes on a 50x50x50 grid for Poisson's equation, the preconditioned conjugate gradient uses symmetric Gauss-Seidel as the preconditioner.

It isn't really an apples-to-apples comparison because each iteration of conjugate gradient requires more work than an SOR iteration, but for most problems it's worth the extra effort.

The cool thing about conjugate gradient is that you can use the line relaxation solver you already implemented as a preconditioner.

[inkant]

Quoting - j.stults

I'd suggest conjugate gradient, it can be very fast, especially with an appropriate preconditioner (which it sounds like you need anyway because of the chemistry).

Here's the convergence of a couple iterative schemes on a 50x50x50 grid for Poisson's equation, the preconditioned conjugate gradient uses symmetric Gauss-Seidel as the preconditioner.

It isn't really an apples-to-apples comparison because each iteration of conjugate gradient requires more work than an SOR iteration, but for most problems it's worth the extra effort.

The cool thing about conjugate gradient is that you can use the line relaxation solver you already implemented as a preconditioner.

Hi Stults,

Sorry for using the word symmetric(A terrible mistake). It does not happen to be symmetric.

THE MATRIX structure is like this tridiagonal near the main diagonal and two other diagonals (with lots of 0s in between) and it does not need to be symmetric.

I would suppose other family of Krylov methods like MINRES or GMRES or BiCG should work. Of course the only way to know is to try and give a shot but do these methods make use of the sparsity of the matrix (the 5 diagonal nature -The 0s inside the band, since by ADI- alternating direction Implicit scheme I can actually work with a tridiagonal matrix, by taking the diagonals far away to be known from prev. iteration).

Best Regards,

Inkant

[ArturGuzik]

Quoting - inkant

I would suppose other family of Krylov methods like MINRES or GMRES or BiCG should work. Of course the only way to know is to try and give a shot but do these methods make use of the sparsity of the matrix

Take a look at MKL (docs and examples) FGMRES solver (in Sparse Solvers Chapter). You can exploit the sparsity by using sparse matrix, matrix-vector procedures/tools you have available via, say, Sparse BLAS. The solver is using Reverse Communication (RCI) and that gives you a full control over solution process. Good preconditioner (ILU0, ILUT are already prepared for you in MKL) and you should be fine.

A.

[j_stults]

Yes, GMRES seems to be the most popular one for fluid problems if you can afford the memory hit. The Krylov methods just require the action of your operator on a vector, so there's nothing special about having a banded matrix as opposed to any other general sparse matrix (which you'd get with unstructured grids).

Even though plain conjugate gradient isn't guaranteed to converge for unsymmetric matrices, I've had good luck with it for toy problems like Burgers' equation that have unsymmetric Jacobians, YMMV.