CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Documentation of trackToFace subroutine (http://www.cfd-online.com/Forums/openfoam-programming-development/123442-documentation-tracktoface-subroutine.html)

boeleman September 12, 2013 19:56

Documentation of trackToFace subroutine
 
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 September 23, 2013 19:05

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 September 28, 2013 12:42

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 September 28, 2013 17:35

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 October 1, 2013 17:31

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 November 6, 2013 04:13

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 November 13, 2013 20:35

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 February 27, 2014 06:20

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 May 19, 2014 09:40

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 May 19, 2014 09:55

Hi,

what kind of reconstructPar errors are you getting?

cbpf May 19, 2014 10:45

Thank you for the fast answer.

The error is:

--> FOAM Warning :
From function void Foam::particle::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::printStack(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::particle::initCellFacePt() in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/liblagrangian.so"
#3 Foam::Cloud<Foam::passiveParticle>::initCloud(bool ) in "/home/celio/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/lib/liblagrangian.so"
#4 Foam::reconstructLagrangianPositions(Foam::polyMes 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 May 28, 2014 11:35

Quote:

Originally Posted by boeleman (Post 461877)
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:

Code:

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.


All times are GMT -4. The time now is 05:30.