[spaceplumber]

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?

gSum( ) vs sum()
How to sum up a volScalarField

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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;
}


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

[jherb]

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

[spaceplumber]

Originally Posted by jherb

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:

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:

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.

[spaceplumber]

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.

[spaceplumber]

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.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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; 

}


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

[spaceplumber]

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

[olesen]

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.

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

const auto& volume = mesh_.V();
const auto& position = mesh_.C();

You need real care with the sum() function!

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

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.

     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.

// 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.

[olesen]

Originally Posted by spaceplumber

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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.?

[spaceplumber]

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.