CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Looping over a volScalarField in a parallel run (http://www.cfd-online.com/Forums/openfoam-programming-development/76664-looping-over-volscalarfield-parallel-run.html)

dohnie June 1, 2010 08:46

Looping over a volScalarField in a parallel run
 
Dear all,
I'd like to do some on-the-run monitoring of my computation where I am interested in a temperature front. Actually, of all the cells in my domain where the temperature exceeds 1000 K, I'd like to determine the one which has the highest x-coordinate. My code looks like this:

Code:

label maxindex=0;

forAll(mesh.C(),i)
{
    if(T[i]>1000.0)
    {
    if(mesh.C()[i].component(0) >  mesh.C()[maxindex].component(0))
          maxindex=i;
    }
}

Info<< "front      @ " << mesh.C()[maxindex] << endl;

Not very elegant, but it works well - as long as I run on one CPU. However, in a parallel run, the forAll loop is done over one CPU only. Can anybody suggest a solution how to loop over the whole domain, please?

Regards,
Florian

niklas June 1, 2010 09:55

ooo, a well formulated problem with a simple answer, my favourite :)

you cant use the index cause it only makes sense on the own processor.
I would do it like this.

Code:

scalar  xMax = -1.0e+10;

forAll(mesh.C(),i)
{
    if(T[i]>1000.0)
    {
        if(mesh.C()[i].component(0) >  xMax)
        {
            xMax = mesh.C()[i].component(0);
        }
    }
}

reduce(xMax, maxOp<scalar>);
Info<< "front      @ " << xMax << endl;


dohnie June 2, 2010 14:26

Thank you Niklas,
the reduce command should read
Code:

reduce(xMax, maxOp<scalar>());
then it works fine.

A more tricky problem: I also want to know the y and z coordinate of the location, where the highest x coordinate is (not the maximum y and z coordinate). So a separate reduce on the scalars won't work.

Code:

scalar  xMax = -1.0e+10;
vector maxVector(-1e10,0.0,0.0);

forAll(mesh.C(),i)
{
    if(T[i]>1000.0)
    {
        if(mesh.C()[i].component(0) >  xMax)
        {
            xMax = mesh.C()[i].component(0);
            yMax = mesh.C()[i].component(1);
            zMax = mesh.C()[i].component(2);

            maxVector = mesh.C()[i];
        }
    }
}

reduce(xMax, maxOp<scalar>);
Info<< "front      @ " << xMax << endl;

Can I apply the reduce/maxOp operators to a vector, with the first component as argument??


Thanks,
Florian

niklas June 4, 2010 02:49

That was a bit trickier...
I dont know if there's such an operation already implemented, but I would transfer the location vector for each coordinate to the main processor and do the comparison there.

Here's a piece of code that generates a random vector on each processor and sends it to the 0-processor, it should be straightforward to apply it to your case.
You might wanna put it inside an if(Pstream:: parRun()) if you want it to work for both serial and parallell runs. (edit:: actually its not necessary)

Code:

    Random rnd(0);

    // get the number of processors
    label n=Pstream::nProcs();

    // generate a random vector, since the seed is the same for all procs they will be the same
    // if we only do it once
    vector localV = rnd.vector01();
    for(label i=0; i<Pstream::myProcNo(); i++)
    {
        localV = rnd.vector01();
    }

    // print the vector on the processor (so we can compare it later when we print it on the main proc)
    Sout << "[" << Pstream::myProcNo() << "] V = " << localV << endl;

    if (Pstream::myProcNo() == 0)
    {
        List<vector> allV(n, vector::zero);
        allV[0] = localV;
        for(label i=1; i<n; i++)
        {
            // create the input stream from processor i
            IPstream vStream(Pstream::blocking, i);
            vStream >> allV[i];
        }
        // print the list of all vectors on the main proc
        Info << allV << endl;
    }
    else
    {
        // create the stream to send to the main proc
        OPstream vectorStream
        (
            Pstream::blocking, 0
        );
        vectorStream << localV;
    }


dohnie June 16, 2010 03:58

Thank you, it seems to work fine!


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