CFD Online Discussion Forums - Summation in Function Object for Multiple Processors

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

- OpenFOAM Post-Processing (https://www.cfd-online.com/Forums/openfoam-post-processing/)

- - Summation in Function Object for Multiple Processors (https://www.cfd-online.com/Forums/openfoam-post-processing/230887-summation-function-object-multiple-processors.html)

Summation in Function Object for Multiple Processors

This is my first time writing my own function object for OpenFOAM and am struggling to get it to work with mpirun. There are three quantities I want to calculate and write to a file in the postprocessing directory:

- Total volume of my mesh

sum(mesh_.V())

- Total volume of my liquid (this is multiphase)

sum(mesh_.V()*alpha.internalField())

- Sum of the position of the cell times the volume times the volume fraction (this should output a vector, not a scalar)

sum(mesh_.V()*mesh_.C().internalField()*alpha.inte rnalField())

My code seems to work perfectly when I just run interFoam, but gets hung on the summation part with I try to use it with mpirun. To test where it's getting hung, I separated out each of the calculation steps. It made it through until it hit the summation. Please see the code below which has both these test lines and, below that, the lines with correctly write the info to a file when run on a single processor.

I have tried both sum and gSum because it sounds like gSum is what should be used for multiple processors. Other forum posts have mentioned having good luck with gSum, but I can't figure out what I'm doing differently/incorrectly. Links to these posts are below.

I'm also concerned that even if I can get the summation to work, that it might only sum over the master processor, not all. Does anyone have any recommendations for how to get gSum to work in function objects with mpirun?

https://www.cfd-online.com/Forums/op...um-vs-sum.html
https://www.cfd-online.com/Forums/op...alarfield.html

Code:

/*---------------------------------------------------------------------------*\

  =========                 |

  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox

   \\    /   O peration     | Website:  https://openfoam.org

    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation

     \\/     M anipulation  |

-------------------------------------------------------------------------------

License

    This file is part of OpenFOAM.



    OpenFOAM is free software: you can redistribute it and/or modify it

    under the terms of the GNU General Public License as published by

    the Free Software Foundation, either version 3 of the License, or

    (at your option) any later version.



    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT

    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

    for more details.



    You should have received a copy of the GNU General Public License

    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.



\*---------------------------------------------------------------------------*/



#include "CM.H"

#include "Time.H"

#include "fvMesh.H"

#include "addToRunTimeSelectionTable.H"

#include "fvCFD.H"

#include "OSspecific.H"

#include <iostream>



// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //



namespace Foam

{

namespace functionObjects

{

    defineTypeNameAndDebug(CM, 0);

    addToRunTimeSelectionTable(functionObject, CM, dictionary);

}

}





// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //



Foam::functionObjects::CM::CM

(

    const word& name,

    const Time& runTime,

    const dictionary& dict

)

:

    fvMeshFunctionObject(name, runTime, dict),

    logFiles(obr_, name)//,

{

    read(dict);

    resetName(name);

}





// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //



Foam::functionObjects::CM::~CM()

{}





// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //



void Foam::functionObjects::CM::writeFileHeader(const label i) {



}



bool Foam::functionObjects::CM::read(const dictionary& dict)

{

    return true;

}





bool Foam::functionObjects::CM::execute()

{

    return true;

}





bool Foam::functionObjects::CM::end()

{

    return true;

}





bool Foam::functionObjects::CM::write()

{

  if (Pstream::master()) 

    {

     #include "volFields.H"

     #include "uniformDimensionedFields.H"

     const volScalarField& alpha = mesh_.lookupObject<volScalarField>("alpha.water");

     auto& runTime = fileObr_.time();





     logFiles::write();

     

     //This section just tests where the code is getting hung up. w, x, and y make it through fine. z is where it gets stuck when run in parallel

     auto w = alpha.internalField();

     auto x = alpha.internalField()*mesh_.V();

     auto y = mesh_.V()*mesh_.C().internalField()*alpha.internalField();

     auto z = gSum(mesh_.V()); //This is where it gets stuck



     //This prints out the results to a file- works on one processor, but not with mpirun

     file() << runTime.timeName() << " " << gSum(mesh_.V()) << " " << gSum(mesh_.V()*alpha.internalField()) << " " << gSum(mesh_.V()*mesh_.C().internalField()*alpha.internalField()) << endl;

     }



    return true;

}





// ************************************************************************* //

Just a guess: what happens if you use explicit OpenFOAM types, e.g. scalar for z?

Quote:

Originally Posted by jherb (Post 785097)

Just a guess: what happens if you use explicit OpenFOAM types, e.g. scalar for z?

It actually won't even get through wmake if I use either scalar or double, instead of auto. So maybe something is going wrong with dimensions when I use multiple processors? It just seems odd that it would be an mpi specific issue. This is the error message I get:

Code:

CM.C:121:30: error: cannot convert ‘Foam::dimensioned<double>’ to ‘Foam::scalar {aka double}’ in initialization

      scalar z = sum(mesh_.V());

I'm still reading more about dimensioned double vs. double, but it looks like one fix people found was that they need to put .value() at the end of their variable, so mesh_.V().value() for me. That doesn't work either and I get:

Code:

CM.C:121:26: error: ‘const class Foam::DimensionedField<double, Foam::volMesh>’ has no member named ‘value’

I will keep reading to see if I kind find other possible solutions, but any suggestions would be appreciated.

So I guess adding .value() after sum is what it needs to be a double. I no longer get the error during wmake when specifying that it should be a double/scalar. However, interFoam still gets hung on calculating z.

Alright, I have it figured out. I actually needed to do some calculations before moving to the master processor. Right now I have it set up to calculate each of my needed quantities on individual processors, gather lists, and then write the result on the master processor. I'm not actually convinced that it's even necessary to gather the list, so I'll post an update if I have a simplified code. But either way, this code works for what I wanted even if it isn't the prettiest.

Code:

/*---------------------------------------------------------------------------*\

  =========                 |

  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox

   \\    /   O peration     | Website:  https://openfoam.org

    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation

     \\/     M anipulation  |

-------------------------------------------------------------------------------

License

    This file is part of OpenFOAM.



    OpenFOAM is free software: you can redistribute it and/or modify it

    under the terms of the GNU General Public License as published by

    the Free Software Foundation, either version 3 of the License, or

    (at your option) any later version.



    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT

    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

    for more details.



    You should have received a copy of the GNU General Public License

    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.



\*---------------------------------------------------------------------------*/



#include "CM.H"

#include "Time.H"

#include "fvMesh.H"

#include "addToRunTimeSelectionTable.H"

#include "fvCFD.H"

#include "OSspecific.H"

#include <iostream>



// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //



namespace Foam

{

namespace functionObjects

{

    defineTypeNameAndDebug(CM, 0);

    addToRunTimeSelectionTable(functionObject, CM, dictionary);

}

}





// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //



Foam::functionObjects::CM::CM

(

    const word& name,

    const Time& runTime,

    const dictionary& dict

)

:

    fvMeshFunctionObject(name, runTime, dict),

    logFiles(obr_, name)



{

    read(dict);

    resetName(name);

}





// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //



Foam::functionObjects::CM::~CM()

{}





// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //



void Foam::functionObjects::CM::writeFileHeader(const label i) 

{

  //Print out a file header with what each number is

  file() << "Time(s)" << " " << "totalVol" << " " << "waterVol" << " " << "CMx" << " " << "CMy" << " " << "CMz" << endl;

}



bool Foam::functionObjects::CM::read(const dictionary& dict)

{

    return true;

}





bool Foam::functionObjects::CM::execute()

{

    return true;

}





bool Foam::functionObjects::CM::end()

{

    return true;

}





bool Foam::functionObjects::CM::write()

{

  #include "volFields.H"

  #include "uniformDimensionedFields.H"



  //Assign names to the quantities we want

  auto volume = mesh_.V(); //Volume of each cell in the mesh

  const volScalarField& alpha = mesh_.lookupObject<volScalarField>("alpha.water"); //Volume fraction of each cell in the mesh

  auto position = mesh_.C(); //Position of each cell in the mesh



  //Calculate the values we want on each processor

  double volsum = sum(volume).value(); //Volume of the mesh

  double watervol = sum(volume*alpha.internalField()).value(); //Volume of water

  auto   waterposit = sum(volume*alpha.internalField()*position).value(); //Sum of the volume*volume_fraction*cell_position

  double waterpositX = waterposit.component(0); //Extract x component of sum of the volume*volume_fraction*cell_position

  double waterpositY = waterposit.component(1); //Extract y component of sum of the volume*volume_fraction*cell_position

  double waterpositZ = waterposit.component(2); //Extract z component of sum of the volume*volume_fraction*cell_position



  //Create a variable proci which is the number of the current processor

  const label proci=Pstream::myProcNo();



  //Create lists which have the length equal to the number of processors

  List<double>volList(Pstream::nProcs());

  List<double>watervolList(Pstream::nProcs());

  List<double>waterXList(Pstream::nProcs());

  List<double>waterYList(Pstream::nProcs());

  List<double>waterZList(Pstream::nProcs());



  //Assign each variable a different position in the list based on their processor

  volList[proci]=volsum;

  watervolList[proci]=watervol;

  waterXList[proci]=waterpositX;

  waterYList[proci]=waterpositY;

  waterZList[proci]=waterpositZ;



  //Collapse the lists onto one processor

  Pstream::gatherList(volList);

  Pstream::gatherList(watervolList);

  Pstream::gatherList(waterXList);

  Pstream::gatherList(waterYList);

  Pstream::gatherList(waterZList);



  //On the master processor

  if (Pstream::master()) 

  {

  

     //Assign a name to the time step value

     auto& runTime = fileObr_.time();



     logFiles::write();

            

     //Calculate the sum of each value over all of the processors

     double totalVol = sum(volList)/(Pstream::nProcs());

     double waterVol = sum(watervolList)/(Pstream::nProcs());

     double volX = sum(waterXList)/(Pstream::nProcs());

     double volY = sum(waterYList)/(Pstream::nProcs());

     double volZ = sum(waterZList)/(Pstream::nProcs());



     //Center of mass is the sum(volume*position*volume_fraction)/sum(volume*volume_fraction)

     double CMx = volX/waterVol;

     double CMy = volY/waterVol;

     double CMz = volZ/waterVol;



     //Write the variables we want to the postprocessing log file

     file() << runTime.timeName() << " " << totalVol << " " << waterVol << " " << CMx << " " << CMy << " " << CMz << endl;

  }

  

  return true; 



}





// ************************************************************************* //

Center of Mass Function Object

In case anyone could benefit from it, I've attached the complete center of mass function object directory. You should be able to just run wmake in the attached directory, then include it in your controlDict using something like the following:

center
{
type CM;
libs ("libCMFunctionObject.so")
}

I saw this one a bit late. There are a number of things that are either more work that you need, or too inefficient. Here are some, with inline comments.

Quote:

Code:

//Assign names to the quantities we want

  auto volume = mesh_.V(); //Volume of each cell in the mesh

  auto position = mesh_.C(); //Position of each cell in the mesh

Those made full copies of the fields! Want this instead

Code:

const auto& volume = mesh_.V();

const auto& position = mesh_.C();

You need real care with the sum() function!

Quote:

Code:

double volsum = sum(volume).value();

double watervol = sum(volume*alpha.internalField()).value();

A sum() of volScalarField etc is actually a gSum() - ie a global sum, whereas a sum() of a primitive field is just processor-local sums.
You were right to doubt the need for all of those awkward processor lists. Instead of more comments, here are some simple tips

Code:

scalar someVal =...;   // processor local value

reduce(someVal, sumOp<scalar>());  // global sum



//  or 

scalar localVal =...; 

scalar globalVal = returnReduce(localVal, sumOp<scalar>());



// working with primitive fields

scalarField localField =... ;

scalar localSum = sum(localField);

scalar globSum1 = returnReduce(localSum, sumOp<scalar>());



// same thing

scalar globSum2 = gSum(localField);

average() and gAverage() etc.

If you start using these functions, you will be writing a lot less code, more quickly and easier to understand and maintain.

Quote:

Code:

     double CMx = volX/waterVol;

     double CMy = volY/waterVol;

     double CMz = volZ/waterVol;

IMO that is just plain ugly and complicated. Why three variables for x/y/z?? Use the vector container to organize things with x(), y(), z(). Same memory required, no overhead (all inlined), less typing.

Code:

// obvious filled directly elsewhere 

// with your volX, volY, volZ

vector volDirs(...) ; 

vector CMass = volDirs/waterVol;

If you make these changes, it would be interesting to know how many lines of code you can eliminate.

Quote:

Originally Posted by spaceplumber (Post 785170)

Code:

/*---------------------------------------------------------------------------*\

  =========                 |

  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox

   \\    /   O peration     | Website:  https://openfoam.org

    \\  /    A nd           | Copyright (C) 2020 OpenFOAM Foundation

     \\/     M anipulation  |

Is this a cut and paste error?
Or have you actually transferred your code copyrights to the OpenFOAM Foundation Ltd.?

Thank you so much for your feedback! It looks like you've got some excellent suggestions, so I'll go through them all as soon as I have the time!

Yes, the header is a copy and paste error. I used the OpenFOAM provided template, and was too worried about getting the code to work that I didn't check the header.

It may be a week or two before I'm able to get back to this, but I'll post a new code when I can.