[iordtila]

I am currently writing my thesis in accelerating the Gauss-Seidel smoother of the OpenFoam project. I have trouble understanding the OpenFoam implementation of the Gauss-Seidel method. The exact part I need to change is

double psii;
    int fStart;
    int fEnd = ownStartPtr[0];

    for (int celli=0; celli<nCells; celli++)
    {

        fStart = fEnd;
        fEnd = ownStartPtr[celli + 1];


        psii = bPrimePtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
        }


        psii /= diagPtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
        }

        psiPtr[celli] = psii;
    }
}

How can I change this part to the Jacobi method?

[tas38]

Originally Posted by iordtila

I am currently writing my thesis in accelerating the Gauss-Seidel smoother of the OpenFoam project. I have trouble understanding the OpenFoam implementation of the Gauss-Seidel method. The exact part I need to change is

double psii;
    int fStart;
    int fEnd = ownStartPtr[0];

    for (int celli=0; celli<nCells; celli++)
    {

        fStart = fEnd;
        fEnd = ownStartPtr[celli + 1];


        psii = bPrimePtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
        }


        psii /= diagPtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
        }

        psiPtr[celli] = psii;
    }
}

How can I change this part to the Jacobi method?

I think something like the following

double psii;
    int fStart;
    int fEnd = ownStartPtr[0];

    for (int celli=0; celli<nCells; celli++)
    {

        fStart = fEnd;
        fEnd = ownStartPtr[celli + 1];


        psii = bPrimePtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
            psii -= lowerPtr[facei]*psiPtr[lPtr[facei]];
        }

        psii /= diagPtr[celli];

        psiPtr[celli] = psii;
    }
}

where you have added something like...

const label* const __restrict__ lPtr =
         matrix_.lduAddr().lowerAddr().begin();

to the member function definition

[iordtila]

Originally Posted by tas38

I think something like the following

double psii;
    int fStart;
    int fEnd = ownStartPtr[0];

    for (int celli=0; celli<nCells; celli++)
    {

        fStart = fEnd;
        fEnd = ownStartPtr[celli + 1];


        psii = bPrimePtr[celli];


        for (int facei=fStart; facei<fEnd; facei++)
        {
            psii -= upperPtr[facei]*psiPtr[uPtr[facei]];
            psii -= lowerPtr[facei]*psiPtr[lPtr[facei]];
        }

        psii /= diagPtr[celli];

        psiPtr[celli] = psii;
    }
}

where you have added something like...

const label* const __restrict__ lPtr =
         matrix_.lduAddr().lowerAddr().begin();

to the member function definition

using potentialFoam it takes way too much time for an output (100 times more than gauss seidel method) and using simpleFoam I get a Floating point excepton... I tried updating the bPrimePtr[ ] as well the same way with:

for (int facei=fStart; facei<fEnd; facei++)
{
bPrimePtr[uPtr[facei]] -= psii;
}

but I keep getting Floating point exception in potentialFoam too now.

Thank you for helping if you have any more ideas I'd appreciate it!