[marluc]

Dear all,

I am implementing a new function object (for post-processing purposes) to write the cell center values and the selected field values of the internal mesh and the boundary patches to a file.
It works in serial (part of the code for single processor only) but I am facing problems for the parallel case.

Here attached you find the code of the templated function for writing the data to the file. It accepts as argument the name of the field in the function object dictionary and the name of the file to be written.

My main concern is when the "slave" processors pass the data to the master in order to write them on the file.
Indeed, while the cells (mesh.C().internalField()) and face center (mesh.boundaryMesh()[patchi].faceCentres()) data are correctly stored as "pointField", the same is not correct for the field data. "pointField" is a vector type and it is therefore not correct to store a generic field value that can be a scalar, a vector or a tensor.

My question is then which object type should I use to store the fields passed from the slave to the master instead of pointField? How should I modify that part of the code?

Thank you in advance for your help.
Luca

#include "writeFieldsToFile.H"
#include "volFields.H"

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

template<class Type>
void Foam::writeFieldsToFile::writeData
(
    const word& fieldName,
    const fileName& outputFile
)
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

    if (obr_.foundObject<fieldType>(fieldName))
    {

        const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
        const fvMesh& mesh = field.mesh();

        Info<< nl << "Writing components of " << field.name() 
            << " at cellCentres"
            << " to " << outputFile
            << endl;
        Info<< "    Number of cells written " << field.size() << endl;

        OFstream os(outputFile);

        label nCells = mesh.C().size();
        label nFields = field.size();

        reduce(nCells, sumOp<label>());
        Info << "Number of global cells: " << nCells << endl;
        reduce(nFields, sumOp<label>());
        Info << "Number of global fields: " << nFields << endl;

        // slaves send data to master, master writes data to file
        if (Pstream::parRun())
        {
            if (Pstream::myProcNo() != Pstream::masterNo())
            {
                OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                toMaster << mesh.C().internalField() << field.internalField();
            }

            else // master part
            {
                forAll(mesh.C(), cellI)
                {
                    os<< mesh.C()[cellI].component(0) << tab 
                      << mesh.C()[cellI].component(1) << tab 
                      << mesh.C()[cellI].component(2) << tab
                      << field.internalField()[cellI] << endl;
                }

                pointField bufferCellCenters;
                pointField bufferFieldValues;

                for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                {
                    IPstream fromSlave(Pstream::blocking, slave);
                    fromSlave >> bufferCellCenters >> bufferFieldValues;
                        
                    Info<< "BufferCellCenters size: " << bufferCellCenters.size() << endl;
                    
                    forAll(bufferCellCenters, cellI)
                    {
                        os<< bufferCellCenters[cellI].component(0) << tab 
                          << bufferCellCenters[cellI].component(1) << tab 
                          << bufferCellCenters[cellI].component(2) << tab
                          << bufferFieldValues[cellI] << endl;
                    }
                }
            }

            // write boundary data to file
            const fvPatchList& patches = mesh.boundary();

            forAll(patches, patchi)
            {
                const polyPatch& cPatch = mesh.boundaryMesh()[patchi];
                label nFaces = cPatch.size();
                reduce(nFaces, sumOp<label>());

                const vectorField& faceCentres = mesh.boundaryMesh()[patchi].faceCentres();
                label nFields = faceCentres.size();
                reduce(nFields, sumOp<label>());

                forAll(field.boundaryField()[patchi],facei)
                {
                    Info<< nl <<  field.boundaryField()[patchi][facei] << endl;
                }

                if (Pstream::myProcNo() != Pstream::masterNo())
                {
                    OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                    toMaster << faceCentres << field.boundaryField()[patchi];
                }

                else // master part
                {
                    forAll(faceCentres, faceI)
                    {
                        // get coordinate for face center
                        const vector& c = faceCentres[faceI];
                        os<< c[0] << tab << c[1] << tab << c[2] << tab 
                          << field.boundaryField()[patchi][faceI] << endl;
                    }

                    pointField bufferFaceCenters;
                    pointField bufferFieldValues;

                    for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                    {
                        IPstream fromSlave(Pstream::blocking, slave);
                        fromSlave >> bufferFaceCenters >> bufferFieldValues;
                            
                        Info<< "BufferFaceCenters size: " << bufferFaceCenters.size() << endl;
                        
                        forAll(bufferFaceCenters, faceI)
                        {
                            // get coordinate for face center
                            const vector& c = bufferFaceCenters[faceI];
                            os<< c[0] << tab << c[1] << tab << c[2] << tab 
                              << bufferFieldValues[faceI] << endl;
                        }
                    }
                }
            
            }

        }
        else // single processor only
        {
            forAll(mesh.C(), cellI)
            {
                os<< mesh.C()[cellI].component(0) << tab 
                  << mesh.C()[cellI].component(1) << tab 
                  << mesh.C()[cellI].component(2) << tab
                  << field.internalField()[cellI] << endl;
            }


            const fvPatchList& patches = mesh.boundary();

            forAll(patches, patchi)
            {
                const fvPatchVectorField& faceCentres = mesh.C().boundaryField()[patchi];

                forAll(faceCentres, faceI)
                {
                    // get coordinate for face center
                    const vector& c = faceCentres[faceI];
                    os<< c[0] << tab << c[1] << tab << c[2] << tab 
                      << field.boundaryField()[patchi][faceI] << endl;
                }
            }
        }
    }
}

[Zeppo]

Originally Posted by marluc

Dear all,
My question is then which object type should I use to store the fields passed from the slave to the master instead of pointField? How should I modify that part of the code?

pointField

->

Field<T>

[marluc]

Dear Zeppo,
thanks a lot for your reply. I already tried that and actually it compiles. When you run the case you get the following error because the first keyword it finds is the type of the boundary field and not an integer number or a parentheses.

Do you have an idea how to fix it?
Thanks,
Luca

Writing components of U at cellCentres to "/home/marluc/OpenFOAM/marluc-2.4.0/run/tutorials/incompressible/pimpleFoam/channel395/processor0/../postProcessing/writeFieldsToFile/0.2/cellCenterData_U"
Number of cells written 60000
[0]
[0]
[0] --> FOAM FATAL IO ERROR:
[0] incorrect first token, expected <int> or '(', found on line 0 the word 'type'
[0]
[0] file: IOstream at line 0.
[0]
[0] From function operator>>(Istream&, List<T>&)
[0] in file /home/openfoam/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/lnInclude/ListIO.C at line 149.
[0]
FOAM parallel run exiting
[0]

EDIT:
Dear Sergei,

I think my answer to your post was not clear so I try to formulate it better here.
Here you find the new version of the code:

#include "writeFieldsToFile.H"
#include "volFields.H"

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

template<class Type>
void Foam::writeFieldsToFile::writeData
(
    const word& fieldName,
    const fileName& outputFile
)
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

    if (obr_.foundObject<fieldType>(fieldName))
    {

        const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
        const fvMesh& mesh = field.mesh();
        const DimensionedField<Type,volMesh>& iField = field.dimensionedInternalField(); 

        OFstream os(outputFile);

        label nCells = mesh.C().size();
        label nFields = field.size();

        reduce(nCells, sumOp<label>());
        reduce(nFields, sumOp<label>());

        // slaves send data to master, master writes data to file
        if (Pstream::parRun())
        {
            if (Pstream::myProcNo() != Pstream::masterNo())
            {
                OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                toMaster << mesh.C().internalField() << field.internalField();
            }

            else // master part
            {
                Info<< nl << "Writing components of " << field.name() 
                << " at cellCentres"
                << " to " << outputFile
                << endl;
                Info<< "    Number of cells written " << nCells << endl;

                forAll(mesh.C(), cellI)
                {
                    os<< mesh.C()[cellI].component(0) << tab 
                      << mesh.C()[cellI].component(1) << tab 
                      << mesh.C()[cellI].component(2) << tab
                      << field.internalField()[cellI] << endl;
                }

                vectorField bufferCellCenters;
                Field<Type> bufferFieldValues;

                for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                {
                    IPstream fromSlave(Pstream::blocking, slave);
                    fromSlave >> bufferCellCenters >> bufferFieldValues;
                        
                    forAll(bufferCellCenters, cellI)
                    {
                        os<< bufferCellCenters[cellI].component(0) << tab 
                          << bufferCellCenters[cellI].component(1) << tab 
                          << bufferCellCenters[cellI].component(2) << tab
                          << bufferFieldValues[cellI] << endl;
                    }
                }
            }

            // write boundary data to file
            const fvPatchList& patches = mesh.boundary();

            forAll(patches, patchi)
            {
                const polyPatch& cPatch = mesh.boundaryMesh()[patchi];
                const vectorField& faceCentres = cPatch.faceCentres();

                if (Pstream::myProcNo() != Pstream::masterNo())
                {
                    OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                    toMaster << faceCentres << field.boundaryField()[patchi];
                }

                else // master part
                {
                    forAll(faceCentres, faceI)
                    {
                        // get coordinate for face center
                        const vector& c = faceCentres[faceI];
                        os<< c[0] << tab << c[1] << tab << c[2] << tab 
                          << field.boundaryField()[patchi][faceI] << endl;
                    }

                    pointField bufferFaceCenters;
                    Field<Type> bufferFieldValues;

                    for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                    {
                        IPstream fromSlave(Pstream::blocking, slave);
                        fromSlave >> bufferFaceCenters >> bufferFieldValues;
                            
                        forAll(bufferFaceCenters, faceI)
                        {
                            // get coordinate for face center
                            const vector& c = bufferFaceCenters[faceI];
                            os<< c[0] << tab << c[1] << tab << c[2] << tab
                              << bufferFieldValues[faceI] << endl;
                        }
                    }
                }
            
            }

        }
        else // single processor only
        {
            forAll(mesh.C(), cellI)
            {
                os<< mesh.C()[cellI].component(0) << tab 
                  << mesh.C()[cellI].component(1) << tab 
                  << mesh.C()[cellI].component(2) << tab
                  << field.internalField()[cellI] << endl;
            }


            const fvPatchList& patches = mesh.boundary();

            forAll(patches, patchi)
            {
                const fvPatchVectorField& faceCentres = mesh.C().boundaryField()[patchi];

                forAll(faceCentres, faceI)
                {
                    // get coordinate for face center
                    const vector& c = faceCentres[faceI];
                    os<< c[0] << tab << c[1] << tab << c[2] << tab 
                      << field.boundaryField()[patchi][faceI] << endl;
                }
            }
        }
    }
}

The code above compiles correctly but it returns the above error when running in parallel (in serial everything works fine). The reason is that "field.boundaryField()[patchi] does not fit with the type "Field<Type>" because it does not return a normal List but a sort of fvPatchField where the first token is the type of the boundary passed. The type Field<Type> expects as first token either an integer representing the length of the list or a scalar or a parentheses if a vector field of rank greater or equal 1 is passed. The same does not happen with the internal values that are passed with "field.internalField()" because these are in from of a list.

I hope the problem is now a little bit clearer. I would be greatful to anyone who can help me in fixing it.

Thank you in advance,
Luca

[Zeppo]

So the error comes from the second part of the code where you process boundary values? Slave processes send field.boundaryField()[patchi] which is of fvPatchField<T> type and a master process expects to receive Field<Type>. The data should be consistent, try

Field<Type> bufferFieldValues;

->

fvPatchField<T> bufferFieldValues;

[marluc]

Dear Sergei,
thanks a lot for your answer.
Actually I already tried this solution but without success. Here the first lines of the error message I get during compilation when using fvPatchField<Type> instead of Field<Type>.
If you have an idea how to fix it I would be very grateful.
Thanks,
Luca

In file included from writeFieldsToFile/writeFieldsToFile.H:240:0,
                 from writeFieldsToFile/writeFieldsToFile.C:26:
writeFieldsToFile/writeFieldsToFileTemplates.C: In instantiation of ‘void Foam::writeFieldsToFile::writeData(const Foam::word&, const Foam::fileName&) [with Type = double]’:
writeFieldsToFile/writeFieldsToFile.C:204:59:   required from here
writeFieldsToFile/writeFieldsToFileTemplates.C:220:40: error: no matching function for call to ‘Foam::fvPatchField<double>::fvPatchField()’
                     fvPatchField<Type> bufferFieldValues;
                                        ^
writeFieldsToFile/writeFieldsToFileTemplates.C:220:40: note: candidates are:
In file included from /opt/openfoam240/src/finiteVolume/lnInclude/fvPatchField.H:588:0,
                 from /opt/openfoam240/src/finiteVolume/lnInclude/volFields.H:40,
                 from writeFieldsToFile/writeFieldsToFileTemplates.C:27,
                 from writeFieldsToFile/writeFieldsToFile.H:240,
                 from writeFieldsToFile/writeFieldsToFile.C:26:
/opt/openfoam240/src/finiteVolume/lnInclude/fvPatchField.C:170:1: note: Foam::fvPatchField<Type>::fvPatchField(const Foam::fvPatchField<Type>&, const Foam::DimensionedField<Type, Foam::volMesh>&) [with Type = double]
 Foam::fvPatchField<Type>::fvPatchField
 ^

[Zeppo]

The error

writeFieldsToFile/writeFieldsToFileTemplates.C:220:40: error: no matching function for call to ‘Foam::fvPatchField<double>::fvPatchField()’

means there's no default constructor (a constructor which doesn't take any parameters) defined in type fvPatchField<T> and you have to pass some parameters to create an object of this type. The "minimal" constuctor is declared as fvPatchField (const fvPatch &, const DimensionedField< Type, volMesh > &). All this means that fvPatchField<T> bufferFieldValues won't work.

Also, fromSlave >> bufferFieldValues wouldn't work either even if we had an object of fvPatchField<T> because there is no Istream& operator>>(Istream&, fvPatchField<T>&) defined. When you write fromSlave >> bufferFieldValues the function template<class T> Istream& operator>>(Istream&, List<T>&) is invoked. This is possible since fvPatchField<T> is derived from List<T>. And it is reasonable that operator>>(Istream&, List<T>&) can only take a perfect list (check this: http://cpp.openfoam.org/v4/a07944_source.html)

At the same time it is perfectly ok to write toMaster << field.boundaryField()[patchi] as there's a specific operator to stream your fvPatchField<T> out:

template<class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const fvPatchField<Type>& ptf)
{
    ptf.write(os);
    os.check("Ostream& operator<<(Ostream&, const fvPatchField<Type>&"); 
    return os;
}

This operator writes to the stream not only (list)-values but also the type of a patch field.

Possible solutions:
1. If you are not interested in patch type information you can recast fvPatchField to List:

toMaster << field.boundaryField()[patchi];
->
toMaster << static_cast<const List<Type>&>( field.boundaryField()[patchi]);

and

List<Type> bufferFieldValues;
...
fromSlave >> bufferFieldValues;

2. Another option is more tedious. You parse fvPatchField into elemental parts (patch type, field values) and stream them in and out separately. I think you migth want to look at OpenFOAM-builtin reconstructPar utility's code. It does the same thing as you want your code to do.

[marluc]

Dear Sergei,
thanks a lot for your detailed explanation. Actually I tried to recast fvPatchField but not in the way you did. And your way is correct! Now it works.

I really learned something important and useful from you.

thanks again for the time spent on my question.
Best regards,
Luca

[Zeppo]

Luca, you are welcome.