CFD Online Discussion Forums - OpenFOAM sum row with LDU addressing

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - OpenFOAM sum row with LDU addressing (https://www.cfd-online.com/Forums/openfoam-programming-development/238860-openfoam-sum-row-ldu-addressing.html)

OpenFOAM sum row with LDU addressing

Hello to all,

I would like to know how to sum the value of a row in a OpenFOAM matrix through the LDU addressing.

Is there any already built in function to do this?

I have found for the upperTriangular part, the following:

Code:

const scalarField& diag = matrix_.diag();

const scalarField& lower = matrix_.lower();

const scalarField& upper = matrix_.upper();

 

const labelList& lowerAddr = matrix_.lduAddr().lowerAddr();

const labelList& upperAddr = matrix_.lduAddr().upperAddr();

const labelList& ownStartAddr = matrix_.lduAddr().ownerStartAddr();

 

label fStart, fEnd;

 

sumRowElements = 0;

fStart = ownStartAddr[rowI];

fEnd = ownStartAddr[rowI + 1];

 

for (label j= fStart; j < fEnd; j++)

{

    sumRowElements  += upper[j]*x[upperAddr[j]];

}

Know I need to add the terms of the lower triangular part.

https://i.ibb.co/DQwdMDP/Untitled.png

Can anyone give me a hand with this?

Best Regards

So, I ended up checking Gauss Seidel, and it appears to send the lower adressing contributions to the right hand side vector.

So I did:

Code:

const scalarField& diag = matrix_.diag();

 const scalarField& lower = matrix_.lower();

const scalarField& upper = matrix_.upper();

 

const labelList& lowerAddr = matrix_.lduAddr().lowerAddr();

const labelList& upperAddr = matrix_.lduAddr().upperAddr();

const labelList& ownStartAddr = matrix_.lduAddr().ownerStartAddr();

 

const label nRows = x.size();

 

scalarField store(x.size(), 0);

 

label fStart, fEnd;

 

for (label row = 0; row < nRows; row++)

{

    scalar sumRowElements = 0;

 

    // Start and end of this row

    fStart = ownStartAddr[row];

    fEnd = ownStartAddr[row + 1];

 

    // Sum terms in upper addressing



    for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)

    {

        sumRowElements  += upper[curCoeff]*x[upperAddr[curCoeff]];

    }

 

// Get row summation without diagonal element

    sumRowElements -= store[row];

 

    // Add lower addressing contributions to store vector for the next matrix row 

    for (label curCoeff = fStart; curCoeff < fEnd; curCoeff++)

    {

        store[upperAddr[curCoeff]] -= lower[curCoeff]*x[upperAddr[curCoeff]];

    }

}

Hopefully, this is correct.