[oswald]

Dear Foamers,

I have a problem with the lagrangian intermediate library. My goal is to have access to certain particle properties (e.g. d, U, ...) from within the calculation of a force I'm trying to implement. Unfortunately, the forces are called with code similar to this (from SphereDragForce.C):

template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const

Because p is defined as const I get an error like the following when I try to change for example p.d()

error: passing 'const parcelType {aka const Foam::KinematicParcel<Foam::particle>}' as 'this' argument of 'void Foam::KinematicParcel<ParcelType>::d() [with ParcelType = Foam::particle]' discards qualifiers [-fpermissive]
         p.d();

I played around with a definition as "mutable" - it said, I can't define the property as mutable. I tried to remove the const definition of p in the force, it compiled, but the force calculation then was not called when running a simulation.

So my question is: How can I access distinct parcel properties from within the force calculation?

Thanks in advance!

[tomislav_maric]

If the parcel is passed as the constant reference to a function, it's object should not be changed.

Why should a function that computes a force density from particles/parcels be able to change the particle/parcel properties?

Your question:

I have a problem with the lagrangian intermediate library. My goal is to have access to certain particle properties...

already implies that you have made the right design decision in having non-constant access in the function calcCoupled.

What is it you want to do?

Please do not use C++ terms and use Lagrangian terms in answering the question.

If you are trying to implement only force calculation, you should not modify parcel properties. If your function is computing the force and doing something else to the parcels (modifying them), you should consider refactoring the function into two of them: one for force calculatio, the other function for the property modification. Obviously, property modification will have to happen outside the force calculation function, that does not, and should not, change parcel properties.

[oswald]

Dear Tomislav,

thank you very much for your answer.

I'm trying to implement the Basset or History Force. For this purpose I need to store for example old values of this force and certain other values for every particle, which are not available outside the force calculation and which are also needed to be written out for a restart of calculations. The diameter in the original post was just an example, I'm not trying to change "usual" particle porperties as density, diameter, ...

Sometimes it's a bit strange... I tried on this for several days now, but within some minutes after my post I found out that adapting the force calls in forceSuSp.C/H and ParticleForce.C/H was the solution to my problem. But if you (Tomislav) or anyone else has an approach were I wouldn't have to change the const-definitions I would be glad to here about it!

[tomislav_maric]

Just a few ideas on how to attack this without breaking the function signatures/changing the intermediate library:

Originally Posted by oswald

I'm trying to implement the Basset or History Force. For this purpose I need to store for example old values of this force and certain other values for every particle, which are not available outside the force calculation and which are also needed to be written out for a restart of calculations.

I need to store for example old values of this force

The force is a geometric field, that allows you to store its old and old-old values, do you need more than the last two time steps? If so, you can store a list of fields in the force calculator, in that case, that is the point where you extend the library without modifications.

certain other values for every particle

To me this sounds like an extension of the particle/parcel. Each parcel has a nested class called TrackingData, you can extend that class with additional values.

This is just a rough sketch, but by checking the code, it is prepared for extension without tampering with the intermediate library...

[oswald]

Thank you again for your answer!

I wasn't really correct or specific enough in my last post. I don't only need the force values itself, but also some values that are calculated inside the force calculation. I want to calculate the Basset force according to
Hinsberg, Boonkamp, Clercx (2010): "An efficient, second order method for the approximation of the Basset history force".

In principal, this force needs an integral over the whole path of the particle, which is quite time and memory consumptive. The approach of Hinsberg et al is to split this in the calculation of a "window" term (e.g. the last five time steps), where the actual values are calculated and an approximated "tail". This tail is modeled with 10 exponential functions and each is split in a direct and a recursive part. The value of each of the 10 functions is then used in the next time step for the recursive part.

Do you see any chance to keep these values? It would also be okay if they are lost when the calculation stops.

[tomislav_maric]

You're welcome! I've just browsed through the paper you mentioned. As far as I understand it, it would be possible to compute something like this, with the code extensions I proposed in the previous post. You need a Basset parcel/particle that stores the data needed for the force calculation within a prescribed range of time steps, and a Basset force calculator that expects the same data to be provided by the parcel / particle.

Since the force calculator is templated:

template<class CloudType> 
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled (const typename CloudType:parcelType& p,

You can implement your own parcel type, instantiate the cloud with your parcel type on the solver level (I think it's in createFields.H usually) and then extend the force calculation.

In any case, the article defines your classes and extensions, exactly like your second post - you have a parcel with extended information and a new force calculator that expects that info.

[oswald]

Thank you again, Tomislav! At the moment it's working with the not-so-nice const-removal and I'm doing my simulations with it, because I'm in a hurry (deadline's coming near...). But I'll try to switch to your approach when I have the time for it.

[tomislav_maric]

Originally Posted by oswald

Thank you again, Tomislav! At the moment it's working with the not-so-nice const-removal and I'm doing my simulations with it, because I'm in a hurry (deadline's coming near...). But I'll try to switch to your approach when I have the time for it.

Your're welcome! Don't worry about breaking the interface a bit, you'll have time to refactor once the results are there. I usually open issues on bitbucket so that I don't forget what I planned for code cleanup. Deadlines are definitely more important than clean code, good luck!

[cojua8]

I'm working with the same algorithm, and so far I've managed to do this in a new BassetForce class (using the same structure as existing forces):

    // Private data

    // Constantes Basset
    const scalar N; //pasos basset
    // scalar C1_; //tiempo inicial (0)
    scalar C2_; //tiempo final (N)
    scalarList C3_; //tiempos intermedios (n)

and constructor:

template<class CloudType>
Foam::BassetForcePM<CloudType>::BassetForcePM
(
    CloudType& owner,
    const fvMesh& mesh,
    const dictionary& dict
)
:
    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
    N(readScalar(this->coeffs().lookup("N")))
{
    C2_ = (N-4/3)/((N-1)*sqrt(N-1)+(N-3/2)*sqrt(N));

    C3_.setSize(N-1);

    scalar k;

    forAll(C3_, i)
    {
        k=i+1.0;

        C3_[i] = (k+4.0/3.0)/((k+1.0)*sqrt(k+1.0)+(k+3.0/2.0)*sqrt(k))
               + (k-4.0/3.0)/((k-1.0)*sqrt(k-1.0)+(k-3.0/2.0)*sqrt(k));
    }
}

now I need to store all particle dt and U-Uc in a list of size N.
Should (and can I?) do this in the Basset code or I must modify the kinematicParcel code?