[tzaht]

I am trying to access the age of the particles in my simulation.

I have seen a part of code in the file: lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

template<class ParcelType>
inline Foam::scalar Foam::KinematicParcel<ParcelType>::age() const
{
    return age_;
}

In may time-marching loop, I tried to access this data with the following line of code:

#include "KinematicParcel.H"
Foam::KinematicParcel<Foam::particle>::age()

The error returned:

error: cannot call member function ‘Foam::scalar& Foam::KinematicParcel<ParcelType>::age() [with ParcelType = Foam:article; Foam::scalar = double]’ without object
Foam::KinematicParcel<Foam:article>::age()


This is somewhat unsurprising, however, I am not sure how to proceed. I feel like I need to look up a KinematicParcel object but I am unable to determine what it's name would be, or how to access it.

Any pointers would be appreciated.

Kind regards.

[raumpolizei]

Hey Tony,

you are absolutely right with your guess that you need a particle/parcel or something alike - object in order to call the function. The function is part of a class (- template) and since it is not a static function, you first need an existing object of that class in order to call it.
Generally it is not possible to track a specific parcel since they do not have a specific ID, like for instance the cells or faces of your mesh (practical reasons, the container storing the particle is a linked list). What you could do in your solver (assuming your solver is containing a sprayCloud object - or some cloud object which type is derived from lagrangian/basic/Cloud - named parcels, I am using of1612):

if(parcels.size())
{
auto particle = parcels.cbegin(); // const access to first parcel in cloud
Info << "position: " << particle().position() << nl;
}

Analogously, you could loop over the cloud and do something with all parcels (read specific properties, write the properties to a file etc.)

forAllConstIter(decltype(parcels),parcels,cParcelIter) // or forAllIter it you want to modify
{
    const auto pVelocity = cParcelIter().U(); // access data
    // do something with the data
}

Hope this helps,
Cheers
RP

[tzaht]

Hello raumpolizei,

Thank you for taking the time to reply to my thread.

In my createFields.H file, I have this part to create the lagrangian cloud.

basicKinematicCollidingCloud kinematicCloud
(
    kinematicCloudName,
    rhoInf,
    U,
    mu,
    g
);

The error I am receiving when trying to add on your suggestion:

 grabK.H:20:4: error: ‘parcels’ was not declared in this scope
 if(parcels.size())

I am using OF v1806.

I am nevertheless still a bit confused, I can see that there my code is creating an object called 'kinematicCloud'. I suggest this because in my main loop, we have the following line: kinematicCloud.evolve().

However, it is the case that the kinematicParcel contains the lagrangian particle data. And this itself cannot be access as an object. Is this correct?

[raumpolizei]

Hey there,
the reason why it is not working is hidden in the lines of code you posted. In OpenFOAM and C++, a lot is done via classes. Object may be variables with user defined types (which are defined by classes). kinematicCloud is an object of the specific type basicKinematicCollidingCloud. Now, if you want to use functions of that class (non-static, public functions), you need to put the name of the object + a dot + the name of the function you would like to call, as I have done it for the fictive object parcels in my previous post.

SomeCloud parcels; // create object parcels of type SomeCloud
parcels.doSomeRandomStuff(); // calls function doSomeRandomStuff on object parcels of type SomeCloud

In your code, there is no object called parcels so the code won't compile. You have to replace the name of the object with the name in your solver. I'll let you figure out which name that is .

Cheers

RP

[tzaht]

Hello raumpolizei,

Thanks for your help so far, it has been invaluable and got me going!

I have made some modifications so far and I am able to output the acceleration when I run OpenFOAM in serial.. but am having some issues in parallel.

I have been able to diagnose my issue.

First I will show you how I have modified the code:

I created a field 'Uold' in the file:
lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H and a few more files so that it is added to the KinematicParcel class:

template<class ParcelType>
inline Foam::KinematicParcel<ParcelType>::KinematicParcel
(
    const polyMesh& owner,
    const barycentric& coordinates,
    const label celli,
    const label tetFacei,
    const label tetPti
)
:
    ParcelType(owner, coordinates, celli, tetFacei, tetPti),
    active_(true),
    typeId_(-1),
    nParticle_(0),
    d_(0.0),
    dTarget_(0.0),
    U_(Zero),
    Uold_(Zero),
    acceleration_(Zero),
    rho_(0.0),
    age_(0.0),
    tTurb_(0.0),
    UTurb_(Zero)
{}

Next, in the following code before updating the particle velocity, I have stored the old velocity:

p.Uold() = p.U();

bool Foam::KinematicParcel<ParcelType>::move
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar trackTime
)
{
    typename TrackCloudType::parcelType& p =
        static_cast<typename TrackCloudType::parcelType&>(*this);
    typename TrackCloudType::parcelType::trackingData& ttd =
        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);

    ttd.switchProcessor = false;
    ttd.keepParticle = true;

    const cloudSolution& solution = cloud.solution();
    const scalarField& cellLengthScale = cloud.cellLengthScale();
    const scalar maxCo = solution.maxCo();

    p.Uold() = p.U();
    Info << "Old particle velocity" << p.Uold();
    while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)

This works fine when I am in serial but I am running into an annoying problem when I run in parallel.

For the case where particles enter from a different processor onto a given processor, for some reason the Uold and U field for that time-step are the same, and only the particles that came onto that processor store the old velocity in Uold and the updated velocity in U.

I can't work out why this is happening..

Any help would be appreciated, on even what I can try to debug this.

[raumpolizei]

Hey there,

from my point of view, it looks like you are updating the velocity at the wrong place in the code. The piece of code updating Uold should be at the beginning of the function that is responsible for updating the velocity (I think the function name is something like calcVelocity). If you check the while statements,

while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)

you may see that if that switchProcessor flag is true, this very while loop is never entered, thus prohibiting a velocity update which is resulting in same values for Uold and U.
Cheers

RP

[tzaht]

Originally Posted by raumpolizei

Hey there,

from my point of view, it looks like you are updating the velocity at the wrong place in the code. The piece of code updating Uold should be at the beginning of the function that is responsible for updating the velocity (I think the function name is something like calcVelocity). If you check the while statements,

while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)

you may see that if that switchProcessor flag is true, this very while loop is never entered, thus prohibiting a velocity update which is resulting in same values for Uold and U.
Cheers

RP

Hey, I have taken your advice, now, I find that I am able to access the velocity in calcVelocity but when I assign it, I get the following error message:

lagrangian/intermediate/lnInclude/KinematicParcel.C:188:14: error: passing ‘const vector {aka const Foam::Vector<double>}’ as ‘this’ argument discards qualifiers [-fpermissive]
     p.Uold() = banana;

const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const
{
    const typename TrackCloudType::parcelType& p =
        static_cast<const typename TrackCloudType::parcelType&>(*this);
    typename TrackCloudType::parcelType::trackingData& ttd =
        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);

    const typename TrackCloudType::forceType& forces = cloud.forces();

    // Do a new old velocity
    Info << "Old particle velocity" << p.Uold() << endl;
    const vector banana  = p.U();
    p.Uold() = banana;

    Info << "banana" << banana;

Note, I have used the banana here just as a temporary variable to see if there was an issue in retrieving the data or in assigning it.

Any help would be kindly appreciated.

[einstein_zee]

Originally Posted by tzaht

Hey, I have taken your advice, now, I find that I am able to access the velocity in calcVelocity but when I assign it, I get the following error message:

lagrangian/intermediate/lnInclude/KinematicParcel.C:188:14: error: passing ‘const vector {aka const Foam::Vector<double>}’ as ‘this’ argument discards qualifiers [-fpermissive]
     p.Uold() = banana;

const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const
{
    const typename TrackCloudType::parcelType& p =
        static_cast<const typename TrackCloudType::parcelType&>(*this);
    typename TrackCloudType::parcelType::trackingData& ttd =
        static_cast<typename TrackCloudType::parcelType::trackingData&>(td);

    const typename TrackCloudType::forceType& forces = cloud.forces();

    // Do a new old velocity
    Info << "Old particle velocity" << p.Uold() << endl;
    const vector banana  = p.U();
    p.Uold() = banana;

    Info << "banana" << banana;

Note, I have used the banana here just as a temporary variable to see if there was an issue in retrieving the data or in assigning it.

Any help would be kindly appreciated.

Hey there,

I don't know what you are trying to do but I just saw the issue and wanted to reply if you need the answer in shorter time. The reason that the compiler is complaining is that you assign "banana" to a function call which will return a const access to the particle velocity (see p definition with "const" qualifier). And assignment for a const (if it is not at declaration) is not permitted. Take a look at how "const" works in C++ if you have doubts. Also, you may take a look at "https://github.com/OpenFOAM/OpenFOAM-dev/blob/aab3a22708d20fb3b31b5c897d151f831f9281d7/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H" to see there are two types of access functions "const" and "ref" ones.

[tzaht]

Originally Posted by einstein_zee

Hey there,

I don't know what you are trying to do but I just saw the issue and wanted to reply if you need the answer in shorter time. The reason that the compiler is complaining is that you assign "banana" to a function call which will return a const access to the particle velocity (see p definition with "const" qualifier). And assignment for a const (if it is not at declaration) is not permitted. Take a look at how "const" works in C++ if you have doubts. Also, you may take a look at "https://github.com/OpenFOAM/OpenFOAM-dev/blob/aab3a22708d20fb3b31b5c897d151f831f9281d7/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H" to see there are two types of access functions "const" and "ref" ones.

Thank you for your help, I worked in a different member function that wasn't 'const' and this fixed my problem.

However, one thing keeps bugging me:

const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const

I don't understand why this member function has const at the begging and the end? Why do we need to tell the compiler that we are dealing with a 'const' twice?

[olesen]

Originally Posted by raumpolizei

Hey Tony,

from lagrangian/basic/Cloud - named parcels, I am using of1612):

forAllConstIter(decltype(parcels),parcels,cParcelIter) // or forAllIter it you want to modify
{
    const auto pVelocity = cParcelIter().U(); // access data
    // do something with the data
}

In most recent OpenFOAM versions (perhaps 1712 and later), the following should work

for (const auto& p : parcels)
{
    const vector& vel = p.U(); // access data
    // do something with the data
}

IMO using 'auto' can be reasonable for the loop (since we probably don't really care what type of particle it is), but using 'vector' for the inner deferencing of the U() makes for somewhat better documentation of what is going on.

[einstein_zee]

Originally Posted by tzaht

Thank you for your help, I worked in a different member function that wasn't 'const' and this fixed my problem.

However, one thing keeps bugging me:

const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const

I don't understand why this member function has const at the begging and the end? Why do we need to tell the compiler that we are dealing with a 'const' twice?

Okay, may be you need to brush up your C++ stuff a bit . The first const means the output of the function will be a const. and the last const is only available while you are defining methods in your class. And literally it means you are promising the compiler you would not change any member variable in your class. But as in life promises can be broken .