|
[Sponsors] |
June 7, 2010, 19:53 |
OpenFOAM Preconditioners
|
#1 |
Senior Member
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 621
Rep Power: 0 |
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 |
|
September 12, 2018, 12:52 |
|
#2 |
New Member
Thomas M
Join Date: Aug 2018
Posts: 20
Rep Power: 8 |
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 |
|
February 26, 2021, 09:59 |
|
#3 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22 |
Is there a paper or other documentation describing the "Simplified Diagonal-based Incomplete Cholesky preconditioner = DIC" implemented in OpenFOAM?
|
|
February 26, 2021, 14:00 |
|
#4 |
Super Moderator
Tobias Holzmann
Join Date: Oct 2010
Location: Bad Wörishofen
Posts: 2,711
Blog Entries: 6
Rep Power: 52 |
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 |
|
February 28, 2021, 14:40 |
|
#5 |
Senior Member
|
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} } |
|
March 2, 2021, 17:50 |
|
#6 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22 |
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 */ |
|
March 3, 2021, 06:26 |
|
#7 |
Senior Member
|
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 |
|
March 3, 2021, 07:34 |
|
#8 |
Super Moderator
Tobias Holzmann
Join Date: Oct 2010
Location: Bad Wörishofen
Posts: 2,711
Blog Entries: 6
Rep Power: 52 |
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 |
|
March 3, 2021, 17:00 |
|
#9 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22 |
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 } } // ************************************************************************* // |
|
December 3, 2023, 15:20 |
|
#10 | |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 281
Rep Power: 22 |
Quote:
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]]; } } |
||
Tags |
linear system, openfoam 1.5-dev, preconditioner |
|
|
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 |