[chegdan]

Hello All,

I am curious as to how the diagonal preconditioners work in openfoam. I know that a DIC preconditioner is a Simplified diagonal-based incomplete Cholesky preconditioner for symmetric matrices. I know that the diagonal is saved for later use, but is it just a diagonal matrix preconditioner? what is a good source covering these simplified diagonal preconditioners?

Dan

[tmik]

I think it would also be helpful if someone could shed some light on all the preconditioners in OF, how they work, their applications, etc.

OF only lists them, no description on any advantages or how they differ.
https://openfoam.com/documentation/u...fvSolution.php

[klausb]

Is there a paper or other documentation describing the "Simplified Diagonal-based Incomplete Cholesky preconditioner = DIC" implemented in OpenFOAM?

[Tobi]

You will always find the answer while checking the source code. Okay, c++ implementation is something different than the mathematics in general but if you go through the code, you will be able to figure out how it works.

[dlahaye]

Below is a list of references on the incomplete Cholesky factorization preconditioner.

In case that reverse engineering of the code results in good documentation, I will be happy to give that documentation a look.

Thanks, Domenico.

Original reference on IC preconditioner [1]

@article{meijerink1977iterative,
title={An iterative solution method for linear systems of which the coefficient matrix is a symmetric 𝑀-matrix},
author={Meijerink, J and Van Der Vorst, Henk A},
journal={Mathematics of computation},
volume={31},
number={137},
pages={148--162},
year={1977}
}

Modified IC preconditioner

@article{gustafsson1978class,
title={A class of first order factorization methods},
author={Gustafsson, Ivar},
journal={BIT Numerical Mathematics},
volume={18},
number={2},
pages={142--156},
year={1978},
publisher={Springer}
}

Relaxed IC preconditioner

@article{axelsson1986eigenvalue,
title={On the eigenvalue distribution of a class of preconditioning methods},
author={Axelsson, Owe and Lindskog, Gunhild},
journal={Numerische Mathematik},
volume={48},
number={5},
pages={479--498},
year={1986},
publisher={Springer}
}

Discussion on diagonal enties
@article{kershaw1978incomplete,
title={The incomplete Cholesky—conjugate gradient method for the iterative solution of systems of linear equations},
author={Kershaw, David S},
journal={Journal of computational physics},
volume={26},
number={1},
pages={43--65},
year={1978},
publisher={Elsevier}
}

[klausb]

Let's get started, for those who want to join in:

DIC is a "Simplified Diagonal-based Incomplete Cholesky preconditioner."

What are the simplifications compared to Incomplete Cholesky?

#include "DICPreconditioner.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(DICPreconditioner, 0);

    lduMatrix::preconditioner::
        addsymMatrixConstructorToTable<DICPreconditioner>
        addDICPreconditionerSymMatrixConstructorToTable_;
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::DICPreconditioner::DICPreconditioner
(
    const lduMatrix::solver& sol,
    const dictionary&
)
:
    lduMatrix::preconditioner(sol),
    rD_(sol.matrix().diag())
{
    calcReciprocalD(rD_, sol.matrix());
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::DICPreconditioner::calcReciprocalD
(
    scalarField& rD,
    const lduMatrix& matrix
)
{
    scalar* __restrict__ rDPtr = rD.begin();

    const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
    const scalar* const __restrict__ upperPtr = matrix.upper().begin();

    // Calculate the preconditioned diagonal
    const label nFaces = matrix.upper().size();
    for (label face=0; face<nFaces; face++) // iteration over all upper matrix coefficients 
    {
         rDPtr[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rDPtr[lPtr[face]]; 

        // What's the actual calculation? A(i,j) -= A(i,k)*A(j,k) limited to the diagonal?
     }


    // Calculate the reciprocal of the preconditioned diagonal
    const label nCells = rD.size();

    for (label cell=0; cell<nCells; cell++)
    {
        rDPtr[cell] = 1.0/rDPtr[cell];
    }
}


void Foam::DICPreconditioner::precondition
(
    scalarField& wA,
    const scalarField& rA,
    const direction
) const
{
    scalar* __restrict__ wAPtr = wA.begin();
    const scalar* __restrict__ rAPtr = rA.begin();
    const scalar* __restrict__ rDPtr = rD_.begin();

    const label* const __restrict__ uPtr =
        solver_.matrix().lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr =
        solver_.matrix().lduAddr().lowerAddr().begin();
    const scalar* const __restrict__ upperPtr =
        solver_.matrix().upper().begin();

    label nCells = wA.size();
    label nFaces = solver_.matrix().upper().size();
    label nFacesM1 = nFaces - 1;

    for (label cell=0; cell<nCells; cell++)
    {
    // here we multiply the new preconditioned diagonal rD "rDPtr[cell]" on residual rA "rAPtr[cell]" i.e. wA = preconditioned residual 
    // note: it is common practice to apply the preconditioner in each iteration to the residual rather than to the matrix
        wAPtr[cell] = rDPtr[cell]*rAPtr[cell]; // wA = rD * rA
    }

    for (label face=0; face<nFaces; face++)
     {
       
       // some sort of correction? What's the actual formula/calculation?

         wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]]; 
    }

    for (label face=nFacesM1; face>=0; face--)
     {
       
       // some sort of correction? What's the actual formula/calculation?

         wAPtr[lPtr[face]] -= rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]]; 

     }
}


// ************************************************************************* //

/* For reference but not "simplified" as in OpenFOAM:

Description of IC by Huaguang Song following Gene H. Golub, Charles F. Van Loan [8]

Code: Algorithm 3: Incomplete Cholesky Factorization (IC)

for k = 1 : n
    A(k,k) = sqrt(A(k,k))
    for i = k + 1 : n
        if A(i,k) != 0
            A(i,k) = A(i,k)/A(k,k)
        end
    end
    for j = k + 1 : n
        for i = j : n
            if A(i,j) != 0
                A(i,j) = A(i,j) - A(i,k)A(j,k)
            end
        end
    end
end

Source: [8] Gene H. Golub, Charles F. Van Loan, “Matrix Computations”, JHU Press, Oct 15, 1996

or in the Octave scripting language (see: https://en.wikipedia.org/wiki/Incomp..._factorization)

function a = ichol(a)
    n = size(a,1);

    for k=1:n
        a(k,k) = sqrt(a(k,k));
        for i=(k+1):n
            if (a(i,k)!=0)
                a(i,k) = a(i,k)/a(k,k);            
            endif
        endfor
        for j=(k+1):n
            for i=j:n
                if (a(i,j)!=0)
                    a(i,j) = a(i,j)-a(i,k)*a(j,k);  
                endif
            endfor
        endfor
    endfor

    for i=1:n
        for j=i+1:n
            a(i,j) = 0;
        endfor
    endfor            
endfunction

*/

[dlahaye]

Wonderful to have lift-off on this discussion ;-)

Below is how I understand the code shared by Klaus.

In a next iteration we need to identify where exactly in src/OpenFoam the incomplete Cholesky factors are computed. This factorization in not include in the previous iteration by Klaus. It is unclear to me to what extend this code depends on the version of OpenFoam.

I Construction of the preconditioner

%..Stage (1/3): initialize reverse diagonal to diagonal A
%..this occurs in the constructor of the preconditioner
for i = 1:n
rD[i] = A[i,i]
end

%..Stage (2/3): modify the diagonal
%..this occurs in the loop over all faces
%..for each faces, owner (lPtr) and neighbor (uPtr) are identified
%..in a 2-by-2 grid, only entries 2, 3 and 4 are updated
%..unclear to Domenico in current iteration why only uPtr entries are modified
%..Observe that A is symmetric, thus A[i,j] = A[j,i]
for i = 1:n
for j = i+1:n %..upper triangular part only
rD[j] = rD[j] - A[i,j]*A[i,j]/rD[i];
end
end

%..Stage (3/3): reverse the diagonal
for i = 1:n
rD[i] = 1/rD[i]
end

II Application of the preconditioner: the preconditoner is applied to the residual in three stages
%..Stage (1/3): apply the diagonal part of the preconditioner to the residual vector
%..Stage (2/3): apply the lower triangular part of the preconditioner to the residual vector
%..Stage (3/3): apply the upper trt part of the preconditioner to the residual vector

[Tobi]

Just for reference. I do have the book "Numerical Precipices" that include code explanation actually and you can simply take it and use it. The algorithms mentioned here should be included in the book too (has around 1000 pages).

[klausb]

I had another look at it and compared DIC with FDIC in which factors are calculated differently. Differences, see comments in the DIC code section below:

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::DICPreconditioner::calcReciprocalD
(
    scalarField& rD,
    const lduMatrix& matrix
)
{
    scalar* __restrict__ rDPtr = rD.begin();

    const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
    const scalar* const __restrict__ upperPtr = matrix.upper().begin();

    const label nFaces = matrix.upper().size();
    
    // calculate the preconditioned diagonal DIC
    for (label face=0; face<nFaces; face++) 
    {
        rDPtr[uPtr[face]] -= upperPtr[face]*upperPtr[face]/rDPtr[lPtr[face]]; // DIC: What is actually calculated?
      //rDPtr[uPtr[face]] -= sqr(upperPtr[face])/rDPtr[lPtr[face]]; // FDIC is different
    }

    const label nCells = rD.size();

    // calculate the reciprocal of the preconditioned diagonal - this is what is cached/stored using DIC
    for (label cell=0; cell<nCells; cell++)
    {
        rDPtr[cell] = 1.0/rDPtr[cell];
    }
    
    // NOTE: In FDIC the upper coefficients divided by the diagonal are also calculated/cached/stored
}


void Foam::DICPreconditioner::precondition
(
    scalarField& wA,
    const scalarField& rA,
    const direction
) const
{
    scalar* __restrict__ wAPtr = wA.begin();
    const scalar* __restrict__ rAPtr = rA.begin();
    const scalar* __restrict__ rDPtr = rD_.begin();

    const label* const __restrict__ uPtr =
        solver_.matrix().lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr =
        solver_.matrix().lduAddr().lowerAddr().begin();
    const scalar* const __restrict__ upperPtr =
        solver_.matrix().upper().begin();

    label nCells = wA.size();
    label nFaces = solver_.matrix().upper().size();
    label nFacesM1 = nFaces - 1;

    for (label cell=0; cell<nCells; cell++)
    {
    // here we multiply the new preconditioned diagonal rD "rDPtr[cell]" on residual rA "rAPtr[cell]" i.e. wA = preconditioned residual 
    // note: it is common practice to apply the preconditioner in each iteration to the residual rather than to the matrix
        wAPtr[cell] = rDPtr[cell]*rAPtr[cell]; // wA = rD * rA
    }

    for (label face=0; face<nFaces; face++)
    {
        wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]]; // here we calculate and apply some of the upper coefficients
    }

    for (label face=nFacesM1; face>=0; face--)
    {
        wAPtr[lPtr[face]] -= rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]]; // here we calculate and apply the rest of the upper coefficients
    }
}


// ************************************************************************* //

[klausb]

Originally Posted by dlahaye

Wonderful to have lift-off on this discussion ;-)
...

II Application of the preconditioner: the preconditoner is applied to the residual in three stages
%..Stage (1/3): apply the diagonal part of the preconditioner to the residual vector
%..Stage (2/3): apply the lower triangular part of the preconditioner to the residual vector
%..Stage (3/3): apply the upper trt part of the preconditioner to the residual vector

I am wondering why the implementation of the application of the preconditioner differs e.g. between FDIC and DIC?

Isn't the application of the preconditioner simply a matrix (the preconditioner) x vector (residual = wA) multiplication where the matrix is split into a lower, diagonal and upper part due to the LDU storage format used within OpenFOAM hence we have to apply those three parts in three steps but why differ the algorithms?

FDIC:

void Foam::FDICPreconditioner::precondition
(
    solveScalarField& wA,
    const solveScalarField& rA,
    const direction
) const
{
    solveScalar* __restrict__ wAPtr = wA.begin();
    const solveScalar* __restrict__ rAPtr = rA.begin();
    const solveScalar* __restrict__ rDPtr = rD_.begin();
 
    const label* const __restrict__ uPtr =
        solver_.matrix().lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr =
        solver_.matrix().lduAddr().lowerAddr().begin();
 
    const solveScalar* const __restrict__ rDuUpperPtr = rDuUpper_.begin();
    const solveScalar* const __restrict__ rDlUpperPtr = rDlUpper_.begin();
 
    const label nCells = wA.size();
    const label nFaces = solver_.matrix().upper().size();
    const label nFacesM1 = nFaces - 1;
 
    for (label cell=0; cell<nCells; cell++)
    {
        wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
    }
 
    for (label face=0; face<nFaces; face++)
    {
        wAPtr[uPtr[face]] -= rDuUpperPtr[face]*wAPtr[lPtr[face]];
    }
 
    for (label face=nFacesM1; face>=0; face--)
    {
        wAPtr[lPtr[face]] -= rDlUpperPtr[face]*wAPtr[uPtr[face]];
    }
}

DIC:

void Foam::DICPreconditioner::precondition
(
    solveScalarField& wA,
    const solveScalarField& rA,
    const direction
) const
{
    solveScalar* __restrict__ wAPtr = wA.begin();
    const solveScalar* __restrict__ rAPtr = rA.begin();
    const solveScalar* __restrict__ rDPtr = rD_.begin();
 
    const label* const __restrict__ uPtr =
        solver_.matrix().lduAddr().upperAddr().begin();
    const label* const __restrict__ lPtr =
        solver_.matrix().lduAddr().lowerAddr().begin();
    const scalar* const __restrict__ upperPtr =
        solver_.matrix().upper().begin();
 
    const label nCells = wA.size();
    const label nFaces = solver_.matrix().upper().size();
    const label nFacesM1 = nFaces - 1;
 
    for (label cell=0; cell<nCells; cell++)
    {
        wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
    }
 
    for (label face=0; face<nFaces; face++)
    {
        wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
    }
 
    for (label face=nFacesM1; face>=0; face--)
    {
        wAPtr[lPtr[face]] -= rDPtr[lPtr[face]]*upperPtr[face]*wAPtr[uPtr[face]];
    }
}