CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Looping over a volScalarField in a parallel run

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree5Likes
  • 2 Post By niklas
  • 2 Post By niklas
  • 1 Post By Zeppo

Reply
 
LinkBack Thread Tools Display Modes
Old   June 1, 2010, 08:46
Default Looping over a volScalarField in a parallel run
  #1
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 9
dohnie is on a distinguished road
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
dohnie is offline   Reply With Quote

Old   June 1, 2010, 09:55
Default
  #2
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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;
mm.abdollahzadeh and vitors like this.
niklas is offline   Reply With Quote

Old   June 2, 2010, 14:26
Default
  #3
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 9
dohnie is on a distinguished road
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
dohnie is offline   Reply With Quote

Old   June 4, 2010, 02:49
Default
  #4
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21
niklas will become famous soon enoughniklas will become famous soon enough
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;
    }
mm.abdollahzadeh and Traction like this.
niklas is offline   Reply With Quote

Old   June 16, 2010, 03:58
Default
  #5
Member
 
Florian Ettner
Join Date: Mar 2009
Location: Munich, Germany
Posts: 41
Rep Power: 9
dohnie is on a distinguished road
Thank you, it seems to work fine!
dohnie is offline   Reply With Quote

Old   December 3, 2016, 09:56
Default coded functions with forAll have similar problems with parallel execution
  #6
Member
 
a
Join Date: Oct 2014
Posts: 38
Rep Power: 4
cfd@kgp is on a distinguished road
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.... All suggestions are most welcome.

Thanks and regards.
cfd@kgp is offline   Reply With Quote

Old   December 10, 2016, 19:04
Default
  #7
Senior Member
 
Sergei
Join Date: Dec 2009
Posts: 142
Rep Power: 10
Zeppo is on a distinguished road
Quote:
Originally Posted by cfd@kgp View Post
Code:
reduce(volIntegral, maxOp<scalar>);
-->
Code:
reduce(volIntegral, maxOp<scalar>());
Zeppo is online now   Reply With Quote

Old   December 11, 2016, 05:53
Default
  #8
Member
 
a
Join Date: Oct 2014
Posts: 38
Rep Power: 4
cfd@kgp is on a distinguished road
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,
cfd@kgp is offline   Reply With Quote

Old   December 11, 2016, 12:24
Default
  #9
Senior Member
 
Sergei
Join Date: Dec 2009
Posts: 142
Rep Power: 10
Zeppo is on a distinguished road
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 likes this.
Zeppo is online now   Reply With Quote

Old   December 13, 2016, 06:32
Default
  #10
Member
 
a
Join Date: Oct 2014
Posts: 38
Rep Power: 4
cfd@kgp is on a distinguished road
Thanks Zeppo for a informative reply!

I will try this
Code:
reduce(volIntegral, sumOp<scalar>()); reduce(volIntegralV, sumOp<scalar>());
and get back to you!
cfd@kgp is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Unable to run OF in parallel on a multiple-node cluster quartzian OpenFOAM 3 November 24, 2009 14:37
Parallel run diverges, serial does not SammyB OpenFOAM Running, Solving & CFD 1 May 10, 2009 03:28
Time Parallel run Dr. FLow Squad FLUENT 5 October 11, 2007 03:18
Run in parallel a 2mesh case cosimobianchini OpenFOAM Running, Solving & CFD 2 January 11, 2007 07:33
How to run parallel in ICEM_CFD? Kiddo Main CFD Forum 2 January 24, 2005 09:53


All times are GMT -4. The time now is 13:59.