[DSMC123]

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:olyMesh const&, Foam::Istream&, bool) in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #4 Foam::dsmcParcel::dsmcParcel(Foam:olyMesh 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:

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:

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

[DSMC123]

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

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

[DSMC123]

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:

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):

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:

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

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.