CFD Online Discussion Forums - Looping over a volScalarField in a parallel run

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

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

- - Looping over a volScalarField in a parallel run (https://www.cfd-online.com/Forums/openfoam-programming-development/76664-looping-over-volscalarfield-parallel-run.html)

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

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;

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

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;

    }

Thank you, it seems to work fine!

coded functions with forAll have similar problems with parallel execution

this is a coded function added in controlDict,

Quote:

functions
(

TIntegral
{
functionObjectLibs ("libutilityFunctionObjects.so");
type coded;
redirectType integral;
outputControl timeStep;
code
#{
const volScalarField& fl = mesh().lookupObject<volScalarField>("sMS1:alpha1") ;
// volScalarField magU(mag(U));

scalar volIntegral = 0;
scalar volIntegralV = 0;

forAll (fl, cellI)
{
volIntegral += fl[cellI]*mesh().V()[cellI];
volIntegralV += 1.0*mesh().V()[cellI];
}

//reduce(volIntegral, maxOp<scalar>);
//reduce(volIntegralV, maxOp<scalar>);
volIntegral=volIntegral/volIntegralV;
Info<<"GFL: " << volIntegral << endl;
//foo.write();

#};
}
);

In parallel this does not work, I tried to use the "educe(volIntegral, maxOp<scalar>);" command (which is commented above). But it was not compiled, it gives me a error

Quote:

Selecting radiationModel none
Courant Number mean: 0 max: 0

Starting time loop

Using dynamicCode for functionObject TIntegral at line 60 in ".TIntegral"
Creating new library in "dynamicCode/integral/platforms/linuxGccDPInt32Opt/lib/libintegral_b7227c78da17419d6b110e956c6a2cec82d452 d8.so"
Invoking "wmake -s libso /home/ojas/MY_works_VLP/OpenFoam_learning/practice/beta_-1/dynamicCode/integral"
wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file FilterFunctionObjectTemplate.C
Making dependency list for source file functionObjectTemplate.C
.TIntegral: In member function ‘virtual void Foam::integralFunctionObject::write()’:
.TIntegral:78:35: error: expected primary-expression before ‘)’ token
.TIntegral:79:35: error: expected primary-expression before ‘)’ token
make: *** [Make/linuxGccDPInt32Opt/functionObjectTemplate.o] Error 1

--> FOAM FATAL IO ERROR:
Failed wmake "dynamicCode/integral/platforms/linuxGccDPInt32Opt/lib/libintegral_b7227c78da17419d6b110e956c6a2cec82d452 d8.so"

file: .TIntegral from line 60 to line 65.

From function codedBase::createLibrary(..)
in file db/dynamicLibrary/codedBase/codedBase.C at line 213.

FOAM exiting

How to make use of "reduce" syntax here..

please help...:confused::confused::confused:. All suggestions are most welcome.

Thanks and regards.

Quote:

Originally Posted by cfd@kgp (Post 628100)

Code:

reduce(volIntegral, maxOp<scalar>);

-->

Code:

reduce(volIntegral, maxOp<scalar>());

Thanks Zeppo.
syntax wise (reduce(volIntegral, maxOp<scalar>());) works, but the answer of serial and parallel is different.

parallel output at one time-step.

Quote:

GFL: 0.012389639
Courant Number mean: 2.2797319e-07 max: 2.2945885e-05
Time = 0.36

serial output at one time-step.

Quote:

GFL: 0.0060198669
Courant Number mean: 2.2797193e-07 max: 2.2945764e-05
Time = 0.36

Thanks,

1. You might want to use a built-in feature instead of developing your own code. There's a function object volRegion (http://cpp.openfoam.org/v4/a06355_source.html)

Code:

pIntegral

{

    type            volRegion;

    libs              ("libfieldFunctionObjects.so");

    operation     volAverage;

    regionType   all;

    fields

    (

        p    

    );

}

This intergates p over all cells in a region.

2. Nevertheless, if you want to "reinvent the wheel" (or practice in OpenFOAM/C++ programming, if you like to think of it this way ;)) look at how this function object is implemented

Code:

// "src/functionObjects/field/fieldValues/volRegion/volRegion.H"

157 template<class Type>

158 bool Foam::functionObjects::fieldValues::volRegion::writeValues

159 (

160     const word& fieldName

161 )

162 {

163     const bool ok = validField<Type>(fieldName);

164 

165     if (ok)

166     {

167         Field<Type> values(setFieldValues<Type>(fieldName));

168         scalarField V(filterField(mesh().V()));

...

176         // Combine onto master

177         combineFields(values);

178         combineFields(V);

... 

181         if (Pstream::master())

182         {

183             Type result = processValues(values, V, weightField);

...

211         }

212     }

213 

214     return ok;

215 }

Code:

// "src/functionObjects/field/fieldValues/volRegion/volRegion.H"

 74 template<class Type>

 75 Type Foam::functionObjects::fieldValues::volRegion::processValues

 76 (

 77     const Field<Type>& values,

 78     const scalarField& V,

 79     const scalarField& weightField

 80 ) const

 81 {

 82     Type result = Zero;

 83     switch (operation_)

 84     {

...

105         case opVolAverage:

106         {

107             result = sum(V*values)/sum(V);

108             break;

109         }

...

149     }

150 

151     return result;

152 }

Code:

// "src/functionObjects/field/fieldValues/fieldValue/fieldValueTemplates.C"

32 template<class Type>

33 void Foam::functionObjects::fieldValue::combineFields(Field<Type>& field)

34 {

35     List<Field<Type>> allValues(Pstream::nProcs());

36 

37     allValues[Pstream::myProcNo()] = field;

38 

39     Pstream::gatherList(allValues);

40 

41     if (Pstream::master())

42     {

43         field =

44             ListListOps::combine<Field<Type>>

45             (

46                 allValues,

47                 accessOp<Field<Type>>()

48             );

49     }

50 }

3. Or just try this:

Code:

reduce(volIntegral, maxOp<scalar>());

reduce(volIntegralV, maxOp<scalar>());

Code:

reduce(volIntegral, sumOp<scalar>());

reduce(volIntegralV, sumOp<scalar>());

Thanks Zeppo for a informative reply!

I will try this
Code:
reduce(volIntegral, sumOp<scalar>()); reduce(volIntegralV, sumOp<scalar>());
and get back to you! :)