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/)
-   -   Looping over a volScalarField in a parallel run (https://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!

cfd@kgp December 3, 2016 09:56

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.

Zeppo December 10, 2016 19:04

Quote:

Originally Posted by cfd@kgp (Post 628100)
Code:

reduce(volIntegral, maxOp<scalar>);

-->
Code:

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

cfd@kgp December 11, 2016 05:53

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,

Zeppo December 11, 2016 12:24

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>());


cfd@kgp December 13, 2016 06:32

Thanks Zeppo for a informative reply!

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


All times are GMT -4. The time now is 18:35.