[alchem]

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:

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:

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

    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.

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:

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

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.