pisoFoam floating point error - GAMG
 

Hi Foamers,

I'm on simulating a pipe (a quite big mesh, ~4mio cells) flow with pisoFoam and I have some troubles there. The simulation starts fine, but after a certain interations, it just crashes with an floating point error. The Co number seems to be fine. I first tried it with the PCiCG solver but the calculation took an eternity, so I switched to the GAMG solver.

This is the protocol:
Any ideas?

Hi Erik

First of all, the problem is within your turbulence model and your have bounding problems on k, and the largest value exceed 100. Have you tried upwinding your turbulence convection term? (I have changed from linear to Minmod and that has removed all problems with boundedness).
Secondly your Courant number is larger than 1, hence you would suspect to have an unstable solution, i.e. it might help lowering your time step or start using adaptive time stepping controlled by a given Courant number, e.g. 0.8.
Third, the final residual on the pressure seems to be rather large compared to the residual on the velocity, however due to the limited number of iterations, I assume that you have reached convergence to the tolerance you have specified?
I have read somewhere that 3 PISO-loops ensure that the error in the pressure correction scheme is as small as other errors in the discretization, so you might increase the number to 3. Unfortunately I have forgotten the reference.

I hope this helps,

Niels

Hello Niels,

thank you very much for your fast response!
I have to admit, that I'm still quite new to OpenFoam, so that I'm not sure about all of your answers and don't know exactly what to do.
1. What do you mean by "upwinding the turbulence convection term"? You think about changing the divSchemes in the fvSchemes file, right? My entries are:So I'll change the linear entries to upwind or minmod? Are there some more information on minmod? I haven't found anything about it, neither in the programmers nor in the guide.

2. I thought, that it is enough, if the mean Co number is <1.
I've only posted the final three steps out of 17. The first ten were all <1, but then it began tp climb slightly.
I've searches the forum for the adjustable runtime and will try it right on.

3. Am I right, that the residual for p and U is defined in the fvSolution? My entries are Do I have to increase the tolerance here?

Where can I adjust the PISO loops? Is it in the fvSolution at PISOChanging nCorrectors to 3?


Please excuse me for all the questions.

I realy appreciate your help!
Erik

Hi Erik

No problem, I am glad to help.

@1: You could e.g. change to div(phi,k) Gauss Minmod;. The clever thing is that are you misspelling the name, then the entire list of possible entries are written to the terminal. Regarding information on Minmod do a search on Google/Wikipedia or look into books on e.g. the solution of hyperbolic problems.

@2: No, unfortunately not, the stability criterion is strictly for the largest Courant number. This off course makes the meshing process a vital point in the modeling process, as poor meshes might lead to far larger computational times than strictly needed to obtain a given accuracy.

@3: Yes, you are right, however please note that you also specify a relative tolerance in p, pFinal and U, so if it is satisfied prior to tolerance, then the iteration terminates, this happens obviously in your case. Typically I have relTol = 0 for U and pFinal and sometimes also on p. Hence you could start trying that and I would recommend it (at least the first two).

@PISO correctors: Yes

Good luck:)

Niels

Hi,

I would add to Niels' posts some other hints.
(1) the mesh. as Niels mentioned, mesh is very important. So you should do a checkMesh on your case to see, if the mesh is good. If, e.g. your non-orthogonality is Max: 36.468446 average: 5.3279988, as in my case, you are OK. If the Max flies as high as 70 or more, you should put nNonOrthogonalCorrectors 2; into your PISO{} section in fvSolution. Better is to do the mesh again and better.

(2) the solution should converge in every time step, so if your U velocities reads:
GAMG: Solving for Ux, Initial residual = 0.132126, Final residual = 8.87388e-07, No Iterations 1 it is not OK. The final residual tells you only the linear solver residual, not the one you should be worried about. But that should get better if you follow what Niels suggests.

(3) To my experience, it is enough to use multigrid solver for pressure (as it is in most commercial codes anyway). If you have 4 mil. case you should speed-up by going to parallel (if you can).

(4) your k problem could be also indicating incorrect values for the inlet/ or other BCs
also this is usually the epsilon that blows up.


good luck
matej

Hi Matej,
thanks for your suggestions.
I will set the nNonOrthogonalCorrectors to 3, thank you. My non-orthogonality Max is ~80, and unfortunately there's no chance for a fast remesh of case.
I've tried the simulation with the new settings from Nils, but I stil get a floating point error. I tried an adjustableTimeStep with maxDeltaT 1e-06 and a maxCo 0.7. At the beginning it ran fine, but with more iterations, bounding k was strongly diverging.
As you mentioned, I also think, that I have some problems with the BCs for k.
Here is my k file:
I thougt about changing the outlet condition to zeroGradient and the wall BCs to kqRWallFunction - but in the tutorials it was always fixedValue for pisoFoam.

Actually this is here just a test case, for a similar calculation for the buoyantBoussinesqSimpleFoam solver. There, I do have exactly the problem, that k and epsilon are exploding, and I get the floating point error. I think, it's about the BC's on epsilon and k, but I'm not sure how to deal with it.

By the way, the parallelisation works fine on the case, there are no problems.

Thank you,
Erik

Erik,

(1) non-orthogonality 80 is a top class number ;). what is you mesh made of, polyhedrals? You may try to use some limiting like:

gradSchemes
{
default cellLimited leastSquares 1.0;
}
as well in "divSchemes" like this example:
div(phi,U) Gauss upwind cellLimited leastSquares 1.0;
you may try to use these for others too.

(2) what turb. model do you use? solving only for k, not eps. or omega...
(3) k BCs: yep, using RANS model, you should use the kqRWallFunction
on walls. Your inlet value should be reasonable according to you size of inlet,
inlet velocity and so on.

(4) your initial conditions. Do you start from zeros, or from some previous solution?
if you can, it's always good to start from some solution. What I do is usually running potentialFlow and then pisoFoam or simpleFoam first to start from some physical solution.
to run potentialFlow, it's very convenient to use pyFoam routines, as pyFoamPotentialRunner.py, which does all necessary changes in settings, runs potentialFoam and then turns everything to your previous settings.


good luck
matej

Hi Matej,

1) I tried your suggestionsEven with Gauss upwind cellLimited leastSquares 1.0 phi, I recieve this error:
Where do I have to enter the "surfaceScalarField cellLimited from objectRegistry region0"? And what do I have to enter there?

2) I'm using the LES model: oneEqEddy with the default entries from the tutorial.

3) The kqRWallFunction seems fine now. I'll check the value right on.

4) I start from zero. I've tried to perform potentialFoam but it didn't run. It still made some problems, after entering the requirded laplacian (1,p) parameters.

Best regards,
Erik

Ok. now we are getting somwhere!!

it's LES, so few things to consider, if you want reasonable LES results:

(1) Co Number should be kept around 0.2 definitely not higher

(2) there is usually no reason to run LES simulation from zero initial cond.
try harder with potential flow, or at least laminar flow (switch off turbulence in LESproperties). After you have initialization from laminar or potential flow, you need to
have turbulence properties all around the domain, which you can make with some randomisation (search the forum) or you have to run the simulation for few (4-6) volume change in the domain. only after that you can have some reasonable statistics in the flow. To run the potentialFoam, you need to change few bits in system directory, so follow the error messages, or use pyFoam which do this for you.

(3) BCs: for very fine grid (small y+ ) you may use k=0, for larger values you may use the wall function. For good night reading I'd suggest following discussion:
http://www.cfd-online.com/Forums/ope...043-les-5.html

(4) discretisation: sorry - to fast writing. no upwinding in grad schemes of course!
since you are running LES, you should not use upwind alone, which smooth your turbulent disturbances by numerical diffusion. you should use higher order schemes e.g. linearUpwind or linearUpwindV for velocity. for more info, look here: http://www.opencfd.co.uk/openfoam/doc/fvSchemes.html

http://www.cfd-online.com/Forums/ope...earupwind.html

good source of info on schemes is the source in:
$FOAM/src/finiteVolume/interpolation/surfaceInterpolation/

hope this helps
matej