[boeleman]

Hello,

I am trying to write a Volume of Fluid (VoF) + Lagragian Particle Tracking (LPT) solver based on OF2.1.1 interDymFoam and the basicKinematicCloud template. However, my code hangs in the trackToFace subroutine in particleTemplates.H. Using gdb I figured out that tetLambda gives back a very small value for lam, lambdaMin is set to 0 and the correction algorithm kicks in, but my code still hangs. I would like to understand the the trackToFace and tetLambda code better, but there is not a lot of documentation online. Could anyone provide references to the algorithms used?

Thanks,

Arnout

[boeleman]

Figured I'd post what I think is going on based on reading the code. Not sure this is correct:

1) A distance "dCorr" is picked based on either the parcel velocity or the maximum Courant number, whichever one is smaller.
2) findtris checks whether while going from the center of the tetrahedral, in which the parcel is located, to its current position + dCorr a face is crossed.
3) hitWallFaces checks whether a wall was hit.
4) if both findtris and hitWallFaces do not find any faces the trackToFace will return 1.
otherwise tetLambda will calculate the actual amount traveled before the face was hit.

What I am still confused about is why first an initial guess is calculated using findTris, and why not finding any faces with findTris guarantees that no faces where hit. It calculated from the center of a tet after all and not using the actual particle trajectory.

[GPesch]

Hi Arnout,

sorry for just throwing that at you (as I didnt have time to look at it closely):

I am working with the DSMC Solver and sometimes also have the problem that one of my parcels is stuck in the trackToFace Routine. I could not figure out what exactly happenend and also didnt fully understand the code of trackToFace (as it looks rather complicated )...

Someone was giving me the hint of looking into

Graham Macpherson, Niklas Nordin and Henry G. Weller: "Particle tracking in unstructured, arbitrary polyhedral meshes for use in CFD and molecular dynamics" which was released in
Commun. Numer. Meth. Engng 2009, 55: 263-273.

I didnt look very closely into it but the paper seems to actually describe the trackToFace routine.

Please contact me if you are interested and if you need help gaining access to the publication.

Georg.

[boeleman]

Hello Georg,

Thank you so much for the reference. I just downloaded it and will take a look at it. I did find a workaround for the hanging particles, which might be useful to you. In the trackToFace routine you twice find the line:

position_ += trackingCorrectionTol*(tet.centre() - position_);

in the particle position correction algorithm. trackingCorrectionTol is in the order of 10^-5. If you increase this value to, for example, 1 findTris and tetLambda will agree which other and there will be no more hanging particles (at least not for me). It is not a very elegant hack, and I am going to look at the paper to find a more elegant solution, but it does get the job done. Hopefully it will work for you too.

Cheers,

Arnout

[boeleman]

Just read the article. Things make a lot more sense now. (The reason for having trackToFace is because it is computationally cheaper to store which cell a particle is in, than to find the cell it is in by calculating it every iteration. How lambda is part of the calculation which face was crossed). However, I feel there must be a follow up article, because this article does not mention tet decomposition as a workaround to particles hanging in concave cells. Would you happen to know which article this can be found in?

[GPesch]

Hi Arnout,

sorry for not getting back to you earlier!

Unfortunately I do not know whether a follow-up article is existing! If you stumble over it at some point I would be very happy if you could let me know!

I just wanted to thank you for the little "hack" you posted: I haven't tried it yet, but I will definitely do in order to avoid the hanging particles! Thanks!

Georg

[boeleman]

Hello Georg,

Figured I might as well share the modifications I made:

https://www.dropbox.com/s/9jc27rwjp9...ngianOf211.zip

I ended up writing an extra counter, which resets particles to the center of a grid cell whenever the number of iterations is larger than 100. As long as the counter is smaller than 100, the normal correction is used. Also, I modified the calculation at the wall. My particles were getting stuck because they were larger than grid cells at the wall. Last but not least I think I found a bug in particleTemplates.C where I replaced triI with tI in the trackToFace subroutine. You can do a diff to see the details. My code runs without hanging now, so hopefully it helps you too. The article was very helpful!

Cheers,

Arnout

[GPesch]

Finally had the time to test it out! I'm glad, that you made the modifications as I was to lazy to look into the code! Works well for me—you are my hero, thanks a bunch!!

[cbpf]

Hello boeleman,

Can you say how can i use your files to overcome the lagrangian reconstructPar errors?

I'm very new in openfoam and i don't know how to recompile only that folders.

Thank you.

[GPesch]

Hi,

what kind of reconstructPar errors are you getting?

[cbpf]

Thank you for the fast answer.

The error is:

--> FOAM Warning :
From function void Foam:article::initCellFacePt()
in file particle/particleI.H at line 759
cell, tetFace and tetPt search failure at position (-0.0125212 -0.00987232 0.00017957)
for requested cell 0

--> FOAM FATAL ERROR:


FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 Foam:article::initCellFacePt() in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/liblagrangian.so"
#3 Foam::Cloud<Foam:assiveParticle>::initCloud(bool ) in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/liblagrangian.so"
#4 Foam::reconstructLagrangianPositions(Foam:olyMes h const&, Foam::word const&, Foam::PtrList<Foam::fvMesh> const&, Foam::PtrList<Foam::IOList<int> > const&, Foam::PtrList<Foam::IOList<int> > const&) in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/libreconstruct.so"
#5
in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/reconstructPar"
#6 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7
in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/reconstructPar"
Abortado (imagem do núcleo gravada)

[DSMC123]

Originally Posted by boeleman

Hello Georg,

Figured I might as well share the modifications I made:

https://www.dropbox.com/s/9jc27rwjp9...ngianOf211.zip

I ended up writing an extra counter, which resets particles to the center of a grid cell whenever the number of iterations is larger than 100. As long as the counter is smaller than 100, the normal correction is used. Also, I modified the calculation at the wall. My particles were getting stuck because they were larger than grid cells at the wall. Last but not least I think I found a bug in particleTemplates.C where I replaced triI with tI in the trackToFace subroutine. You can do a diff to see the details. My code runs without hanging now, so hopefully it helps you too. The article was very helpful!

Cheers,

Arnout

Hi Arnout,

I have been seeing similar errors with hanging during the trackToFace function and would like to try your code to sort it, but the .zip file doesn't seem to be working? I'd be grateful to see what you have done.

For interest, I was having similar problems a few years ago with small (microscale) meshes and Graham Macpherson sent me a small loop to add to particleTemplates.C that fixed that issue, but does not solve this one. The loop is:

tetAreas[0] = tet.Sa();
        tetAreas[1] = tet.Sb();
        tetAreas[2] = tet.Sc();
        tetAreas[3] = tet.Sd();
        
            //******
            for (label i = 0; i < 4; i++)
            {
                tetAreas[i] /= (mag(tetAreas[i]) + VSMALL);
            }
            //******

        FixedList<label, 4> tetPlaneBasePtIs;

        tetPlaneBasePtIs[0] = basePtI;
        tetPlaneBasePtIs[1] = f[fPtAI];
        tetPlaneBasePtIs[2] = basePtI;
        tetPlaneBasePtIs[3] = basePtI;

Where the bold is the new code and the rest will show you the correct place in the code to place it. I suspect this stopped particles being larger than the grid cells at the wall as you discussed previously. This loop never seems to have made it into any OpenFOAM release, but Graham was still working for OpenCFD when he sent me this.

Thanks,

Craig

Edit: Downloaded the .zip file now. I hope my code snippet above is at least of interest, if not use.

[DustExplosion]

I just wanted to confirm that the tetArea normalization worked for me on OF 3.0.1. No more stuck particles on micrometer meshes!

It would be nice to get this into the main code version, but I do not know how to request. Can anyone point me in the right direction?

[akidess]

Originally Posted by DustExplosion

I just wanted to confirm that the tetArea normalization worked for me on OF 3.0.1. No more stuck particles on micrometer meshes!

It would be nice to get this into the main code version, but I do not know how to request. Can anyone point me in the right direction?

Post a bug on bugs.openfoam.org, ideally with a simple test case, and the code you suggest as a patch.

[DustExplosion]

Thanks Akid,

I just did the report (http://bugs.openfoam.org/view.php?id=2376). I put it down as "major" as you basically can't use the code without fixing the bug, but low priority as not many people would be running these types of simulations.

[DustExplosion]

Just an update on this. The suggested fix broke a bunch of other functionality so it was not accepted.

[aghsin]

Originally Posted by boeleman

Hello Georg,

Figured I might as well share the modifications I made:

https://www.dropbox.com/s/9jc27rwjp9...ngianOf211.zip

I ended up writing an extra counter, which resets particles to the center of a grid cell whenever the number of iterations is larger than 100. As long as the counter is smaller than 100, the normal correction is used. Also, I modified the calculation at the wall. My particles were getting stuck because they were larger than grid cells at the wall. Last but not least I think I found a bug in particleTemplates.C where I replaced triI with tI in the trackToFace subroutine. You can do a diff to see the details. My code runs without hanging now, so hopefully it helps you too. The article was very helpful!

Cheers,

Arnout

Dear Arnout ,
would you please re-share your solution for this issue again., the file has been deleted. I have the same problem , my particles stuck in unfinited loop close to wall.
Thanks alot

[aghsin]

Originally Posted by DustExplosion

Just an update on this. The suggested fix broke a bunch of other functionality so it was not accepted.

would you please explain what are the other functionality which were ruined by this modification? thanks inadvance

[guin]

I suggest you to try it out with the current development version from the OpenFOAM Foundation. In summary: all these ugly tracking problems shall be history thanks to the recent implementation of the barycentric tracking algorithm (see https://cfd.direct/openfoam/free-sof...tric-tracking/). I am simulating parcels with extremely high propagation speeds and used to have similar troubles with the old tracking method.

Sent from my HUAWEI RIO-L01 using CFD Online Forum mobile app