|
[Sponsors] |
July 21, 2016, 01:40 |
Particle tracking error
|
#1 |
New Member
Join Date: May 2012
Posts: 4
Rep Power: 14 |
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. 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; 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); 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)) ) 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(); } file particleI.H:711 Code:
if (!mesh_.pointInCellBB(position_, oldCellI, 0.3)) Last edited by alchem; July 27, 2016 at 05:06. Reason: Don't want to reply to my own post |
|
August 20, 2016, 19:07 |
|
#2 | |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,980
Blog Entries: 45
Rep Power: 128 |
Greetings alchem,
Many thanks for the detailed report and updates! Quote:
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
__________________
|
||
August 21, 2016, 10:27 |
|
#3 |
New Member
Join Date: May 2012
Posts: 4
Rep Power: 14 |
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 tubeCalculation.zip. To reproduce the error execute following commands: Code:
icoUncoupledKinematicParcelFoam decomposePar 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. Alexey |
|
October 18, 2016, 19:31 |
|
#4 | |
Senior Member
Join Date: Jan 2013
Posts: 372
Rep Power: 14 |
Hello guys,
What is the usage for the following two lines: Code:
nParcels - scalar(nParcelsToInject) > this->owner().rndGen().globalPosition(scalar(0), scalar(1)) 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; } Thank you if you could give some hints. Quote:
|
||
October 19, 2016, 08:27 |
|
#5 |
New Member
Join Date: May 2012
Posts: 4
Rep Power: 14 |
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. |
|
May 6, 2017, 16:30 |
|
#6 |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,980
Blog Entries: 45
Rep Power: 128 |
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:
Best regards, Bruno
__________________
|
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
[OpenFOAM] Native ParaView Reader Bugs | tj22 | ParaView | 270 | January 4, 2016 11:39 |
OpenFOAM without MPI | kokizzu | OpenFOAM Installation | 4 | May 26, 2014 09:17 |
Compile problem | ivanyao | OpenFOAM Running, Solving & CFD | 1 | October 12, 2012 09:31 |
Ansys Fluent 13.0 UDF compilation problem in Window XP (32 bit) | Yogini | Fluent UDF and Scheme Programming | 7 | October 3, 2012 07:24 |
Version 15 on Mac OS X | gschaider | OpenFOAM Installation | 113 | December 2, 2009 10:23 |