[dohnie]

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:

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]

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.

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]

Thank you Niklas,
the reduce command should read

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.

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]

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)

    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]

Thank you, it seems to work fine!

[cfd@kgp]

this is a coded function added in controlDict,

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

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.

[Zeppo]

Originally Posted by cfd@kgp

reduce(volIntegral, maxOp<scalar>);

-->

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

[cfd@kgp]

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

parallel output at one time-step.

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

serial output at one time-step.

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

Thanks,

[Zeppo]

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)

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

// "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 }

// "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 }

// "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:

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

->

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

[cfd@kgp]

Thanks Zeppo for a informative reply!

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