CFD Online Discussion Forums - Particle tracking error

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

- OpenFOAM Bugs (https://www.cfd-online.com/Forums/openfoam-bugs/)

- - Particle tracking error (https://www.cfd-online.com/Forums/openfoam-bugs/174933-particle-tracking-error.html)

Particle tracking error

This bug is possibly related to this http://www.cfd-online.com/Forums/ope...an-fields.html

OpenFOAM sometimes forgot to change particle cell when particle move from one mesh cell to another.

I was carrying out some particle simulations on 1 core machine and when I tried to decompose it so I could continue calculation on 8 core machine I received following error:

Code:

Time = 0.000655

Identified lagrangian data set: "kinematicCloud"



--> FOAM FATAL ERROR: 

position (-0.03446333 -0.03248561 0.06133807)

    for requested cell 4132

    If this is a restart or reconstruction/decomposition etc. it is likely that the write precision is not sufficient.

    Either increase 'writePrecision' or set 'writeFormat' to 'binary'



    From function void Foam::particle::initCellFacePt()

    in file /home/alchem/OpenFOAMRep/OpenFOAM-dev/src/lagrangian/basic/lnInclude/particleI.H at line 718.

That was simple icoUncoupledKinematicParcelFoam simulation, just particles moving in a constant velocity field. Write precision was more than enough for this case and I tried to look at this particle through paraview. At 0.000655 time paraview throws similar error but on previous time step 0.00065 I managed to see particle positions. Actually there was no particle at (-0.03446333 -0.03248561 0.06133807) coordinates and it could not move too far away within time step. When I looked at its coordinates in 0.00065/lagrangian/kinematicCloud/positions file I found that there is a particle really close to the coordinates which is not shown by paraview.

After some time of incomprehension and study of this impossible problem I saw that this particle moved normally for some time, then slowed and its coordinates on paraview screen and in positions file started to diverge. Here is the particle coordinates from paraview on 0.00065 time (-0.0346717 -0.0330058 0.0613266) and here is the "positions" file coordinates (-0.0344722 -0.03249445 0.06133771). Then I looked at cell id of particle and everything fall into place. Positions file tells us that particle is in the 4132 cell and if we look at this particle in paraview it will be in 4132 cell alright, but actually it is in 246747 cell if we look at its coordinates in positions file. I uploaded "ParticlePositionErrorScaledText.jpg" screenshot to illustrate this. Then I traced this particle and found moment when its cell started to diverge. Particle moved from 4132 cell to 4107 cell but its cell stayed 4132 in "positions" file.

Most frustrating thing about this bug is that I can't replicate it. When I tried to continue calculation from time step before cell transition, particle correctly changed cell. This bug steadily happens when I start calculations, but I don't know bugged particle id beforehand which makes debugging difficult.

Bug happened on OpenFOAM v3.0.1 by The OpenFOAM Foundation.

Edit:

I didn't know that OpenFOAM random number generator seeded by constant, so I can replicate this bug. After some time of study I found out that though initial particle cell is correct, initial tet is not. OpenFOAM split every cell on several tets even if mesh was tetrahedral. I used injection model based on PatchInjection where particle position set by patchInjectionBase::setPositionAndCell function. On patchInjectionBase.C:218 line there is comment which says:

Code:

            // The position is between the face and cell centre, which could

            // be in any tet of the decomposed cell, so arbitrarily choose the

            // first face of the cell as the tetFace and the first point after

            // the base point on the face as the tetPt.  The tracking will pick

            // the cell consistent with the motion in the first tracking step

            tetFaceI = mesh.cells()[cellOwner][0];

            tetPtI = 1;

"The tracking will pick the cell consistent with the motion in the first tracking step" apparently it won't. This function set first cell tet as tet where particle is and somehow it is right most of the time. Particles spawned in different tet generate this error. OpenFOAM don't save particle tet IDs so when I continue the calculation it set tetIDs for existent particles and previously bugged particle moves correctly.

This modification of setPositionAndCell function in injection model fixes the error in serial case:

Code:

    bool continueLoop = false;

    do

    {

        patchInjectionBase::setPositionAndCell

        (

            this->owner().mesh(),

            this->owner().rndGen(),

            position,

            cellOwner,

            tetFaceI,

            tetPtI

        );

        if(position != pTraits<vector>::max)

        {

            continueLoop = !this->findCellAtPosition

            (

                cellOwner,

                tetFaceI,

                tetPtI,

                position,

                false

            );

        }

    } while (continueLoop);

Parallel run is still broken somehow but it is not necessarily related to this error.

Edit 2:

Fixed parallel issue. In fact I made my injection model based on PatchFlowRateInjection not on PatchInjection. Parallel calculation on PatchInjection working just fine and comparing this two models I find an error on PatchFlowRateInjection.C:173.

Code:

error

            nParcelsToInject > 0

         && (

               nParcels - scalar(nParcelsToInject)

             > this->owner().rndGen().position(scalar(0), scalar(1))

            )



correct

            nParcelsToInject > 0

         && (

               nParcels - scalar(nParcelsToInject)

             > this->owner().rndGen().globalPosition(scalar(0), scalar(1))

            )

Henry Weller fixed this error for PatchInjection model by https://github.com/OpenFOAM/OpenFOAM...8bdcd1e4d1fa80 commit in OpenFOAM-dev repo and it seems he forgot to fix it for PatchFlowRateInjection.

My previous modification of setPositionAndCell function breaks parallel execution. I returned setPositionAndCell to previous state and added "parcel.initCellFacePt();" to the end of setProperties function in injection model. For PatchFlowRateInjection it looks like this:

Code:

template<class CloudType>

void Foam::PatchFlowRateInjection<CloudType>::setProperties

(

    const label,

    const label,

    const scalar,

    typename CloudType::parcelType& parcel

)

{

    // set particle velocity to carrier velocity

    parcel.U() = this->owner().U()[parcel.cell()];



    // set particle diameter

    parcel.d() = sizeDistribution_->sample();

    parcel.initCellFacePt();

}

And I increased the distance in which initCellFacePt searches injected particle. Without it initCellFacePt throws error on some injected particles.

file particleI.H:711

Code:

if (!mesh_.pointInCellBB(position_, oldCellI, 0.3))

After this changes icoUncoupledKinematicParcelFoam works fine (at least I didn't catch any more errors) in parallel and serial case.

Greetings alchem,

Many thanks for the detailed report and updates!

Quote:

Originally Posted by alchem (Post 610593)

This has been fixed before OpenFOAM 4.0 was released, thanks to the following bug report and feedback: http://bugs.openfoam.org/view.php?id=2111 - But it wasn't fixed in 3.0.x; the commit where this was fixed is this: https://github.com/OpenFOAM/OpenFOAM...06c91dea90f572

As for the fix you've used in "setProperties", namely by using "initCellFacePt" and further increasing the search distance in that method, is something I'm unable to assess if this is really necessary or not and in under which conditions. If you can provide more details, preferably a test case with which we can reproduce this error, it would make it a lot easier to fully assess if that is the best solution or not.

Nonetheless, I'll try to keep this in mind, specially due to another recent bug report: http://bugs.openfoam.org/view.php?id=2199

Best regards,
Bruno

Greetings Bruno,

Thanks for your reply!

Error occurred in tube which deliver particles to coaxial nozzle so I made simple tube mesh to reproduce error and to my surprise on new mesh error never occurred. Original case is too large to post it here so I split mesh, made test case only with its small part and finally caught the error.
Here is the test case Attachment 50001. To reproduce the error execute following commands:

Code:

icoUncoupledKinematicParcelFoam

decomposePar

Error in decomposePar output:

Code:

--> FOAM FATAL ERROR: 

position (-0.0322173 0.0342116 0.0584572)

    for requested cell 24

    If this is a restart or reconstruction/decomposition etc. it is likely that the write precision is not sufficient.

    Either increase 'writePrecision' or set 'writeFormat' to 'binary'



    From function void Foam::particle::initCellFacePt()

    in file /home/alchem/OpenFOAMRep/OpenFOAM-dev/src/lagrangian/basic/lnInclude/particleI.H at line 718.

I made new mesh with gmsh and old mesh was made with Gambit meshing software and converted to OpenFOAM with fluent3DMeshToFoam utility. It seems that Gambit mesh is broken but checkmesh output looks fine. Maybe Gambit meshes somehow incompatible with OpenFOAM.

Code:

Checking geometry...

    Overall domain bounding box (-0.0357063 0.0180697 0.0549235) (-0.0146907 0.0357063 0.0641)

    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)

    Mesh has 3 solution (non-empty) directions (1 1 1)

    Boundary openness (-8.12545e-18 2.25481e-17 3.83927e-17) OK.

    Max cell openness = 1.80104e-16 OK.

    Max aspect ratio = 4.5683 OK.

    Minimum face area = 2.78812e-07. Maximum face area = 2.4357e-06.  Face area magnitudes OK.

    Min volume = 7.39938e-11. Max volume = 1.01678e-09.  Total volume = 4.26142e-07.  Cell volumes OK.

    Mesh non-orthogonality Max: 51.0966 average: 17.7988

    Non-orthogonality check OK.

    Face pyramids OK.

    Max skewness = 0.617267 OK.

    Coupled point location match (average 0) OK.

    Face tets OK.

    Min/max edge length = 0.000632301 0.00313668 OK.

    All angles in faces OK.

    All face flatness OK.

    Cell determinant (wellposedness) : minimum: 0 average: 1.12875

 ***Cells with small determinant (< 0.001) found, number of cells: 62

  <<Writing 62 under-determined cells to set underdeterminedCells

    Concave cell check OK.

    Face interpolation weight : minimum: 0.233926 average: 0.45583

    Face interpolation weight check OK.

    Face volume ratio : minimum: 0.305358 average: 0.846504

    Face volume ratio check OK.

Best regards,
Alexey

Hello guys,

What is the usage for the following two lines:

Code:

               nParcels - scalar(nParcelsToInject)

             > this->owner().rndGen().globalPosition(scalar(0), scalar(1))

I can guess this condition is used for determining if one more parcel is needed to be injected. However, I am not sure why we need to write this condition as the above way. I find the function of globalPosition:

https://github.com/OpenFOAM/OpenFOAM...cachedRandom.C

Code:

(>template<>

  Foam::scalar Foam::cachedRandom::globalPosition

  (

      const scalar& start,

      const scalar& end

  )

  {

      scalar value = -GREAT;



      if (Pstream::master())

      {

          value = scalar01()*(end - start);

      }



      Pstream::scatter(value);



      return start + value;

  }

So in that case, it always returns 0+1?

Thank you if you could give some hints.

Quote:

Originally Posted by alchem (Post 614840)

Code:

icoUncoupledKinematicParcelFoam

decomposePar

Error in decomposePar output:

Code:

--> FOAM FATAL ERROR: 

position (-0.0322173 0.0342116 0.0584572)

    for requested cell 24

    If this is a restart or reconstruction/decomposition etc. it is likely that the write precision is not sufficient.

    Either increase 'writePrecision' or set 'writeFormat' to 'binary'



    From function void Foam::particle::initCellFacePt()

    in file /home/alchem/OpenFOAMRep/OpenFOAM-dev/src/lagrangian/basic/lnInclude/particleI.H at line 718.

Code:

Checking geometry...

    Overall domain bounding box (-0.0357063 0.0180697 0.0549235) (-0.0146907 0.0357063 0.0641)

    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)

    Mesh has 3 solution (non-empty) directions (1 1 1)

    Boundary openness (-8.12545e-18 2.25481e-17 3.83927e-17) OK.

    Max cell openness = 1.80104e-16 OK.

    Max aspect ratio = 4.5683 OK.

    Minimum face area = 2.78812e-07. Maximum face area = 2.4357e-06.  Face area magnitudes OK.

    Min volume = 7.39938e-11. Max volume = 1.01678e-09.  Total volume = 4.26142e-07.  Cell volumes OK.

    Mesh non-orthogonality Max: 51.0966 average: 17.7988

    Non-orthogonality check OK.

    Face pyramids OK.

    Max skewness = 0.617267 OK.

    Coupled point location match (average 0) OK.

    Face tets OK.

    Min/max edge length = 0.000632301 0.00313668 OK.

    All angles in faces OK.

    All face flatness OK.

    Cell determinant (wellposedness) : minimum: 0 average: 1.12875

 ***Cells with small determinant (< 0.001) found, number of cells: 62

  <<Writing 62 under-determined cells to set underdeterminedCells

    Concave cell check OK.

    Face interpolation weight : minimum: 0.233926 average: 0.45583

    Face interpolation weight check OK.

    Face volume ratio : minimum: 0.305358 average: 0.846504

    Face volume ratio check OK.

Best regards,
Alexey

Hello,

globalPosition(start,end) returns random number between start and end. Unlike position(start, end) it returns the same number on all threads in a parallel case.

Pstream::master() - returns true if we are in a master thread
scalar01() - returns random number between 0 and 1
Pstream::scatter(value) - send value to slave threads
Finally return start + value returns same random value between start and end on all parallel threads

It seems we need to use globalPosition instead of position because otherwise different random numbers on different threads will break global number of particles consistency.

Greetings to all!

My apologies, but I never managed to look into this since my last post here. Fortunately a massive overhaul was done recently on OpenFOAM-dev that will hopefully solve this issue. See the following commit for more details: https://github.com/OpenFOAM/OpenFOAM...a9bc8c1aff0b85

And as explained in the commit, don't expect to be able to continue a simulation from a previous version with this new version. However, you should be able to start a new simulation without problems.

I don't know if and when OpenFOAM 5 will be released, but if and when it is, it will also have this feature.

Either way, you can already test this by installing and using OpenFOAM-dev, by following the:

official instructions at https://openfoam.org/version/dev/
or the unofficial instructions at https://openfoamwiki.net/index.php/I...x/OpenFOAM-dev

Best regards,
Bruno