CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Understanding matrix_.Amul(..) in PCG.C (

michi May 26, 2009 11:37

Understanding matrix_.Amul(..) in PCG.C
Hallo everybody!

I am trying to understand the implementation of the linear solvers, for example the PCG solver. Most iterative methods have a matrix-vector multiplication which is performed in OpenFOAM by

matrix_.Amul(result, psi, interfaceBouCoeffs_, interfaces_, cmpt);

If we then have a further look in the implementation of Amul the following design can be found:
1. initMatrixInterfaces
2. Do the calculation: result = A*psi

3. updateMatrixInterfaces

I tried now to unterstand steps 1 and 3. As far as I see, steps 1 and 3 belong to some kind of correction of the final result. For example, if a case runs in parallel we have to transport some values of psi from one processor to the other. That are the values(cells) of psi which are close to the processor boundary. My question is now:
What exactly is the task of initMatrixInterfaces(....) and updateMatrixInterfaces(....)? Why is this correction of the result vector split into two steps and not just performed in one step?

It would be great if someone could help me!

Lukasz June 1, 2011 10:52


Have you received the answer for this question? I know OF uses zero halo layer approach but still have to figure out how does it refers to updatematrixinterfaces().

Let me know if you know the answer.


marupio June 1, 2011 14:26

I haven't delved too deeply into this, but I'm pretty sure these interfaces are used only for boundary conditions that contribute to off-diagonal matrix terms - cyclic, processor, possibly GGI (?). The regular boundary conditions are absorbed into the diagonal and source term of the matrix before calling an ldu solver... but the off-diagonal ones cannot be, as different ldu solvers handle them differently. Ultimately, that is why boundary coefficients and interfaces are given to the solver.

I'm guessing the Amul function is one of the ldu solver's methods for handling these off-diagonal contributions.

Amiin August 27, 2013 04:14

Hope to find my answer in this old thread!

My question is about the second part. (result = A*psi)
I know that OpenFOAM uses LDU method (lower, upper and diagonal) to store a matrix.
Is it a simple matrix multiplication as I show below?
[ApsiPtr](vector) = [lowerPtr + upperPtr + diagPtr](matrix) * [psiPtr](vector)

It seems a little bit different!
//Quote from lduMatrixATmul.C in OpenFOAM 2.2.0
for (register label cell=0; cell<nCells; cell++)
ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];

for (register label face=0; face<nFaces; face++)
ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];

Best regards,


luckycfd September 23, 2013 10:47

Dear Amiin;
Are you find your question?
Is the second part(result = A*psi) refers to the matrix multiplication?


All times are GMT -4. The time now is 21:54.