[bagbagwan]

Hi all

I am trying to make continuous random walk model for Lagrangian particle tracking.

The problem I faced with is below.
Whenever I try, number generated from random generation dose not change. Always same number came out.

It seems that seed for random function is the problem but I am not familiar with C++ and all method I found turned out to be fail.

Could you check my code and give me advice to fix it.

I am using openFoam v7.

xi is scalar variable I made for saving random number for each time step.

PHP Code:

on kinematicParcel.C
.....

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class ParcelType>
template<class TrackCloudType>
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 scalarField& cellLengthScale = cloud.cellLengthScale();
    const scalar maxCo = cloud.solution().maxCo();

    while (ttd.keepParticle && !ttd.switchProcessor && p.stepFraction() < 1)
    {
        // Cache the current position, cell and step-fraction
        const point start = p.position();
        const scalar sfrac = p.stepFraction();

        // Total displacement over the time-step
        const vector s = trackTime*U_;

        // Cell length scale
        const scalar l = cellLengthScale[p.cell()];

        // Deviation from the mesh centre for reduced-D cases
        const vector d = p.deviationFromMeshCentre();

        // Fraction of the displacement to track in this loop. This is limited
        // to ensure that the both the time and distance tracked is less than
        // maxCo times the total value.
        scalar f = 1 - p.stepFraction();
        f = min(f, maxCo);
        f = min(f, maxCo*l/max(small*l, mag(s)));
        if (p.active())
        {
            // Track to the next face
            p.trackToFace(f*s - d, f);
        }
        else
        {
            // At present the only thing that sets active_ to false is a stick
            // wall interaction. We want the position of the particle to remain
            // the same relative to the face that it is on. The local
            // coordinates therefore do not change. We still advance in time and
            // perform the relevant interactions with the fixed particle.
            p.stepFraction() += f;
        }

        const scalar dt = (p.stepFraction() - sfrac)*trackTime;

        // Avoid problems with extremely small timesteps
        if (dt > rootVSmall)
       {
            // Update cell based properties
            p.setCellValues(cloud, ttd);

            // Generate new random function

//          Random& randObj = cloud.rndGen(); 
            Random randObj(p.origId()+p.age()); 
            scalar oldxi = p.xi();
            scalar dxi = randObj.scalarNormal() - oldxi;
            p.xi() = randObj.scalarNormal() + dxi*0.0;

            // update Dispersion modol
//
            p.calcDispersion(cloud, ttd, dt);
//          p.calcDispersion(cloud, ttd, dt, dxi);

            if (cloud.solution().cellValueSourceCorrection())
            {
                p.cellValueSourceCorrection(cloud, ttd, dt);
            }

            p.calc(cloud, ttd, dt);
        }

        p.age() += dt;

        if (p.active() && p.onFace())
        {
            cloud.functions().postFace(p, ttd.keepParticle);
        }

        cloud.functions().postMove(p, dt, start, ttd.keepParticle);

        if (p.active() && p.onFace() && ttd.keepParticle)
        {
            p.hitFace(f*s - d, f, cloud, ttd);
        }
    }

    return ttd.keepParticle;
}
..... 


https://drive.google.com/file/d/13Fv...ew?usp=sharing

you can install my code by

PHP Code:

cd ./myKinematicParcelFoam /lagrangian/intermediate/
wclean
wmake

cd ../../../myKinematicParcelFoam
wclean
wmake

cd ../particeTest
myKinematicParcelFoam 


[crubio.abujas]

I think it may be related with the seed provided. You always feed it with the same numbers so it returns the same numbers. I checked any OF code using Random and found that in brownianMotionForce a random source is included. In this case it uses the Random defined on the inherited class DSCMCloud as:

Random& rnd = this->owner().rndGen();

I see that you have a slightly different version of this code commented upwards, so maybe tried it and didn't work for you. Maybe you can try that form, and see if that fix your problem.

Anyway, if you keep pulling the thread you find that owner().rndGen() is calling a private method called rndGen_, defined inside the DSMCCloud. The code looks like this:

// DSMCCloud.H
...
Random rndGen_;
...


// DSMCCloud.C
...
rndGen_(label(971501) + 1526*Pstream::myProcNo()),
...

I am not familiar with the nature of the pseudo-random algorithm itself, but it seem that it is taking some big number as an offset and afterwards adding some scaled variable external to the simulation (Process ID, I think) and using it to create an unique seed. That seed will change on each simulation, as the PID will change as well, so you can try to mimic this type of initialization.

[bagbagwan]

Yes, it worked. Thanks Crubio.

For other's information. Main reason to my problem was that after Random generated, I loaded random number with specified seed from could object.

[hbulus]

You have already solved your problem, but i share a perfect way of generating random numbers in C++ that i learned recently:
The solution is nearly perfect, the only issue left is running in parallel. This might seem a non-issue when we just want to implement random numbers for an application we only will use in serial. However, the trick is
rather easy.
We use the current time as seed value and add the PID. This will ensure, that when multiple processes are spawned at the same time, when starting a parallel run, each process has its unique seed value thanks to the
contribution of the PID.

1 // random stuff
2 # include " Random .H"
3 # include " clock .H"
4 Random ranGen ( clock :: getTime ()+pid ());
5
6 for (int j = 0; j < 20; j ++)
7 {
8 Info << ranGen . integer (1, 100) << endl ;
9 }

[crubio.abujas]

Hi hbulus,

I haven't check it, but I don't know if the parallel might be an issue with the approach I've suggested. I've been checking the Pstream::myProcNo() and it seems to return the process number (0, 1, 2, ...), so the seed shall be different for each of the processors.

What I do found problematic is that the process number is not randomly distributed, so the seed are going to be the same for each processor (although different among them). The usage of clock and pid numbers you proposed seems to make a much more randomly seed. I find your approach more convenient so I'm wondering why the implementation of the cloud uses this fix seed approach.