CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

OpenFOAM Preconditioners

Register Blogs Community New Posts Updated Threads Search

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 7, 2010, 19:53
Default OpenFOAM Preconditioners
  #1
Senior Member
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 621
Rep Power: 0
chegdan will become famous soon enoughchegdan will become famous soon enough
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
chegdan is offline   Reply With Quote

Old   September 12, 2018, 12:52
Default
  #2
New Member
 
Thomas M
Join Date: Aug 2018
Posts: 20
Rep Power: 8
tmik is on a distinguished road
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
tmik is offline   Reply With Quote

Old   February 26, 2021, 09:59
Default
  #3
Senior Member
 
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22
klausb will become famous soon enough
Is there a paper or other documentation describing the "Simplified Diagonal-based Incomplete Cholesky preconditioner = DIC" implemented in OpenFOAM?
klausb is offline   Reply With Quote

Old   February 26, 2021, 14:00
Default
  #4
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Bad Wörishofen
Posts: 2,711
Blog Entries: 6
Rep Power: 52
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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.
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   February 28, 2021, 14:40
Default
  #5
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 783
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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}
}
dlahaye is offline   Reply With Quote

Old   March 2, 2021, 17:50
Default
  #6
Senior Member
 
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22
klausb will become famous soon enough
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?

Code:
#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

*/
klausb is offline   Reply With Quote

Old   March 3, 2021, 06:26
Default
  #7
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 783
Blog Entries: 1
Rep Power: 17
dlahaye is on a distinguished road
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
dlahaye is offline   Reply With Quote

Old   March 3, 2021, 07:34
Default
  #8
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Bad Wörishofen
Posts: 2,711
Blog Entries: 6
Rep Power: 52
Tobi has a spectacular aura aboutTobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to 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).
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   March 3, 2021, 17:00
Default
  #9
Senior Member
 
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22
klausb will become famous soon enough
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:


Code:
// * * * * * * * * * * * * * * * 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 is offline   Reply With Quote

Old   December 3, 2023, 15:20
Default
  #10
Senior Member
 
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22
klausb will become famous soon enough
Quote:
Originally Posted by dlahaye View Post
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?


Code:
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]];
    }
}
Code:
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]];
    }
}
klausb is offline   Reply With Quote

Reply

Tags
linear system, openfoam 1.5-dev, preconditioner


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36
Modified OpenFOAM Forum Structure and New Mailing-List pete Site News & Announcements 0 June 29, 2009 06:56
64bitrhel5 OF installation instructions mirko OpenFOAM Installation 2 August 12, 2008 19:07
Adventure of fisrst openfoam installation on Ubuntu 710 jussi OpenFOAM Installation 0 April 24, 2008 15:25
OpenFOAM Debian packaging current status problems and TODOs oseen OpenFOAM Installation 9 August 26, 2007 14:50


All times are GMT -4. The time now is 23:48.