CFD Online Discussion Forums - IOStream Code

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

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - IOStream Code - dsmcParcelIO.C (https://www.cfd-online.com/Forums/openfoam-programming-development/126809-iostream-code-dsmcparcelio-c.html)

IOStream Code - dsmcParcelIO.C

Hi,

I have been developing the DSMC code to run with vibrational energy and this involves adding a scalarField of vibrational energies (one energy for each mode). My code runs well in serial, but I am getting segmentation faults when running in parallel:

[1] #3 cfree in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #3 cfree in "/lib/x86_64-linux-gnu/libc.so.6"
[1] #4 Foam::dsmcParcel::dsmcParcel(Foam::polyMesh const&, Foam::Istream&, bool) in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #4 Foam::dsmcParcel::dsmcParcel(Foam::polyMesh const&, Foam::Istream&, bool) in "/home/craig/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libdsmcPolyStrath.so"
[1] #5 in void Foam::Cloud<Foam::dsmcParcel>::move<Foam::dsmcParc el::trackingData>(Foam::dsmcParcel::trackingData&, double)"/home/craig/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libdsmcPolyStrath.so"
[0] #5 void Foam::Cloud<Foam::dsmcParcel>::move<Foam::dsmcParc el::trackingData>(Foam::dsmcParcel::trackingData&, double) in "/home/craig/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/li in "/home/craig/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lb/libdsmcPolyStrath.so"

I'm quite sure it's an issue with my code in dsmcParcelIO.C - I have created a IOField<scalarField> for the vibrational energy that successfully writes out the vibrational energies:

Code:

IOField<vector> U(c.fieldIOobject("U", IOobject::NO_READ), np);

    IOField<scalar> ERot(c.fieldIOobject("ERot", IOobject::NO_READ), np);

    IOField<scalarField> EVib(c.fieldIOobject("EVib", IOobject::NO_READ), np);

    IOField<label> typeId(c.fieldIOobject("typeId", IOobject::NO_READ), np);

    IOField<label> newParcel(c.fieldIOobject("newParcel", IOobject::NO_READ), np);

    IOField<label> classification(c.fieldIOobject("classification", IOobject::NO_READ), np);

    IOField<label> trackingNumber(c.fieldIOobject("trackingNumber", IOobject::NO_READ), np);



    label i = 0;

    forAllConstIter(dsmcCloud, c, iter)

    {

        const dsmcParcel& p = iter();



        U[i] = p.U();

        ERot[i] = p.ERot();

        EVib[i] = p.EVib();

        typeId[i] = p.typeId();

        newParcel[i] = p.newParcel();

        classification[i] = p.classification();

        trackingNumber[i] = p.trackingNumber();

        i++;

    }



    U.write();

    ERot.write();

    EVib.write();

    typeId.write();

    newParcel.write();

    classification.write();

    trackingNumber.write();

but I cannot get them to read back in and then run in parallel. I guess the question is does anyone know to successfully read back in a IOField<scalarField> and pass the values to the correct particles? I am currently trying this:

Code:

IOField<vector> U(c.fieldIOobject("U", IOobject::MUST_READ));

    c.checkFieldIOobject(c, U);



    IOField<scalar> ERot(c.fieldIOobject("ERot", IOobject::MUST_READ));

    c.checkFieldIOobject(c, ERot);

    

    IOField<scalarField> EVib(c.fieldIOobject("EVib", IOobject::MUST_READ));

    c.checkFieldIOobject(c, EVib);



    IOField<label> typeId(c.fieldIOobject("typeId", IOobject::MUST_READ));

    c.checkFieldIOobject(c, typeId);



    IOField<label> newParcel(c.fieldIOobject("newParcel", IOobject::MUST_READ));

    c.checkFieldIOobject(c, newParcel);

    

    IOField<label> classification(c.fieldIOobject("classification", IOobject::MUST_READ));

    c.checkFieldIOobject(c, classification);

    

    IOField<label> trackingNumber(c.fieldIOobject("trackingNumber", IOobject::MUST_READ));

    c.checkFieldIOobject(c, trackingNumber);



    label i = 0;

    forAllIter(dsmcCloud, c, iter)

    {

        dsmcParcel& p = iter();



        p.U_ = U[i];

        p.ERot_ = ERot[i];

        p.EVib_ = EVib[i];

        p.typeId_ = typeId[i];

        p.newParcel_= newParcel[i];

        p.classification_ = classification[i];

        p.trackingNumber_ = trackingNumber[i];

        i++;

    }

As I said, this works fine in serial, but crashes when running in parallel.

Any help would be appreciated here.

Thanks,

Craig

Could it be a problem with my Ostream operator function? The problem only seems to happen when a particle hits a boundary (even a processor patch).

Code:

Foam::Ostream& Foam::operator<<

(

    Ostream& os,

    const dsmcParcel& p

)

{    

    if (os.format() == IOstream::ASCII)

    {

        os  << static_cast<const particle&>(p) 

            << token::SPACE << p.U_

            << token::SPACE << p.ERot_

            << token::SPACE << p.EVib_

            << token::SPACE << p.typeId_

            << token::SPACE << p.newParcel_

            << token::SPACE << p.classification_

            << token::SPACE << p.trackingNumber_;

    }

    else

    {

        os  << static_cast<const particle&>(p);

        os.write

        (

            reinterpret_cast<const char*>(&p.U_),

            sizeof(p.U_)

            + sizeof(p.ERot_)

            + sizeof(p.typeId_)

            + sizeof(p.newParcel_)

            + sizeof(p.classification_)

            + sizeof(p.trackingNumber_)

        );

        os << p.EVib_;

    }



    // Check state of Ostream

    os.check

    (

        "Foam::Ostream& Foam::operator<<"

        "(Foam::Ostream&, const Foam::dsmcParcel&)"

    );



    return os;

}

p.EVib_ is a scalarField (all other variables here are vectors, scalars or labels) - I don't know how to deal with a scalarField in this context and I can't find any examples in OpenFOAM to copy, so this is merely my best guess at what to do...

I got this working, so I'll post my solution here in the hope it helps someone in the future.

I had to move the EVib_ entry to the end of the list of things to be initialised, which initially looked like:

Code:

inline Foam::dsmcParcel::dsmcParcel

(

    const polyMesh& mesh,

    const vector& position,

    const vector& U,

    const scalar ERot,

    const List<scalar>& EVib,

    const label cellI,

    const label tetFaceI,

    const label tetPtI,

    const label typeId,

    const label newParcel,

    const label classification,

    const label trackingNumber

)

:

    particle(mesh, position, cellI, tetFaceI, tetPtI),

    U_(U),

    ERot_(ERot),

    EVib_(EVib),

    typeId_(typeId),

    newParcel_(newParcel),

    classification_(classification),

    trackingNumber_(trackingNumber)

{}

and I changed it to (obviously also changing all the other entries throughout the code that needed changed due to this):

Code:

inline Foam::dsmcParcel::dsmcParcel

(

    const polyMesh& mesh,

    const vector& position,

    const vector& U,

    const scalar ERot,

    const List<scalar>& EVib,

    const label cellI,

    const label tetFaceI,

    const label tetPtI,

    const label typeId,

    const label newParcel,

    const label classification,

    const label trackingNumber

)

:

    particle(mesh, position, cellI, tetFaceI, tetPtI),

    U_(U),

    ERot_(ERot),

    typeId_(typeId),

    newParcel_(newParcel),

    classification_(classification),

    trackingNumber_(trackingNumber),

    EVib_(EVib)

{}

This allowed me to use the IO entries I posted above:

Code:

Foam::dsmcParcel::dsmcParcel

(

    const polyMesh& mesh,

    Istream& is,

    bool readFields

)

:

    particle(mesh, is, readFields),

    U_(vector::zero),

    ERot_(0.0),

    typeId_(-1),

    newParcel_(0),

    classification_(0),

    trackingNumber_(0),

    EVib_(0)

{    

    if (readFields)

    {           

        if (is.format() == IOstream::ASCII)

        {

            is >> U_;

            ERot_ = readScalar(is);

            typeId_ = readLabel(is);

            newParcel_ = readLabel(is);

            classification_ = readLabel(is);

            trackingNumber_ = readLabel(is);

            is >> EVib_;

        }

        else

        {            

            is.read

            (

                reinterpret_cast<char*>(&U_),

                sizeof(U_)

                + sizeof(ERot_)

                + sizeof(typeId_)

                + sizeof(newParcel_)

                + sizeof(classification_)

                + sizeof(trackingNumber_)

            );

            is >> EVib_;

        }

    }



    // Check state of Istream

    is.check

    (

        "Foam::dsmcParcel::dsmcParcel"

        "("

            "const Cloud<dsmcParcel>& cloud, "

            "Foam::Istream&, "

            "bool"

        ")"

    );

}

and

Code:

Foam::Ostream& Foam::operator<<

(

    Ostream& os,

    const dsmcParcel& p

)

{    

    if (os.format() == IOstream::ASCII)

    {

        os  << static_cast<const particle&>(p) 

            << token::SPACE << p.U()

            << token::SPACE << p.ERot()

            << token::SPACE << p.typeId()

            << token::SPACE << p.newParcel()

            << token::SPACE << p.classification()

            << token::SPACE << p.trackingNumber()

            << token::SPACE << p.EVib();

    }

    else

    {

        os  << static_cast<const particle&>(p);

        os.write

        (

            reinterpret_cast<const char*>(&p.U()),

            sizeof(p.U())

            + sizeof(p.ERot())

            + sizeof(p.typeId())

            + sizeof(p.newParcel())

            + sizeof(p.classification())

            + sizeof(p.trackingNumber())

        );

        os << p.EVib();

    }



    // Check state of Ostream

    os.check

    (

        "Foam::Ostream& Foam::operator<<"

        "(Foam::Ostream&, const Foam::dsmcParcel&)"

    );



    return os;

}

I do not entirely understand why I had to do this (and I only found it by trying anything that I could think of and lots of Pout and Info statements), but I guess it's something to do with the order that they are passed in the binary IO stream.