Custom preconditioner implemantation approach?

[klausbu]

I'd like to implement an additional preconditioner for SPD/symmetric matrices i.e. to solve the pressure equation:

There are several candidates but they have many similarities from an computational perspective.

How can i implement linear algebra matrix operations involving the Identity Matrix (I) e.g. M = (I-LD^-1)(I-D^-1L^T) where D^-1 can be computed as in the diagonal preconditioner and L^T represents the transpose of L (should be U in OpenFOAM terms)?

[TonkomoLLC]

Hello,
I am not in a position to dig into this query in detail at the moment, but is here anything in the RapidCFD GAMG solver (RapidCFD-dev/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG) that may give you ideas of how to solve your problem? A quick grep showed matrix multiplication involving the identity matrix in the GAMG solver. If this is not helpful, my apologies: As mentioned, I cannot get too deep into this topic this week.

Good luck with your developments.

[klausbu]

No hurry.

I think, the current use of the term "identity" in RapidCFD sources refers in some cases to the identity of something like a processor number or to https://en.wikipedia.org/wiki/Identity_transform but not to operations involving the Identity Matrix.

[TonkomoLLC]

Thanks for the clarification. I am not sure I know how to help you with this matter offhand. Sorry about that.

[klausbu]

What's needed to compute "M = (I-LD^-1)(I-D^-1L^T)". Let's look at it step-by-step:

The basis could be diagonalPreconditioner.C / diagonalPreconditioner.H and as we look at symmetric matrices, only wA preconditioning will be needed.

1: The Identity Matrix of size "sol.matrix().diag().size()" needs to be created which is a diagonal matrix with ones only:

const scalargpuField& IdentityMatrix = ???

2: D^-1 is and should be stored temporarily as it's used twice:

const scalargpuField& Diag = solver_.matrix().diag();

thrust::transform
(
    Diag.begin(),
    Diag.end(),
    rD.begin(),
    divideOperatorSFFunctor<scalar,scalar,scalar>(1.0)
);

Unfortunately, I can't find operations supporting the following operations - maybe someone has an idea about how to approach these operations 3-6?:

3: Matrix * Matrix multiplication: Compute LD^-1

4: Transpose of a matrix (probably not needed if we can use U): Compute D^-1L^T where L^T is the transpose of L which should be "U" i.e. upper in OpenFOAM terms, I think

5: Subtraction of the Identity Matrix - Matrix: Compute the subtraction operations I-LD^-1 and I-D^-1L^T

6: Another matrix * matrix multiplication like in step 4 to complete the preconditioner

7: Apply the preconditioner i.e. compute wA, according to the diagonal preconditioner example

The AINV of RapidCFD is of no help as it uses Jacobi = diagonal preconditioning as the main building block according to the paper and from what I understand from the code.

[TonkomoLLC]

I probably can't get to thinking about your message in detail for a while because I am on holiday and I don't have the firsthand experience on this matter to reply quickly. Sorry about that.

[TonkomoLLC]

I haven't forgotten about your question, but I do not know the answer.

Perhaps you need to create functors to accomplish the required matrix operations? Do the matrix operations in lduMatrixSolverFunctors.H help?

Sorry I do not have a more direct answer to your question. I ported PBiCGStab to RapidCFD some years ago, but I am not sure this work will be helpful to you.

Perhaps there are others reading about this issue who can help.

[klausbu]

Yes, I think some functors need to be created or adapted. My problem is, that I need matrix x matrix operations instead of matrix x vector operations and those appear only in the AINV preconditioner implementation.

In the paper which is the basis for the AINV implementation in RapidCFD with the title "Algorithm for Sparse Approximate Inverse Preconditioners in the Conjugate Gradient Method" by Ilya B. Labutin and Irina V. Surodina, they describe and tested three implementations with different orders of approximations of the inverse to be used as preconditioner i.e. 1, 3, 7 orders using Jacobi as an initial approximation for the inverse to calculate D1 = J + J(I − AJ) where D1 is the first order approximation of the inverse of A which serves as preconditioner in one case.

There is a normalMult operation in the AINV implementation but I don't understand the code. I assume that's the matrix x matrix operation I am looking for.