CFD Online Discussion Forums

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)

Shibi October 7, 2021 07:31

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

Shibi October 7, 2021 13:44

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.


All times are GMT -4. The time now is 00:37.