[card]

Hi all,

I posted this in the Running/Solving forum, but it appears that no one is interested. I was reluctant to post this as a bug, but I really don't understand how a uniform freestream with appropriate BC's does not return absolute zero residuals in icoFoam. Also, when making a small change to the grid the behaviour is different. I have a very detailed post here:

Strange behaviour

In the post, I have a simple test case with everything you need. An expert would take no more than five minutes to run it and see what I am saying. I've set this up in Fluent and the behaviour is as expected: zero residuals forever. The interesting part is if you run foamCalc for the divergence and the magnitude of the gradient, you'll find that it is absolute zero except in a few locations along the boundaries. Also, the gradient is absolute zero except along two boundaries. Along those two boundaries, the behaviour is symmetric. Thanks to anyone who atleast takes a look. Happy new year.

[luca_g]

Dear David,

I looked at your post and downloaded the case, so I offer you my opinion: there's nothing wrong in icoFoam (nor in the OpenFOAM libraries): it is just a matter of reference scale, although I agree that the resulting behaviour can be misleading and might not be ideal in this case.

Basically all the residuals you see printed by any solver are NORMALISED residuals based on a normalisation factor which is computed in :
OpenFOAM-1.4.1/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

If you edit the file:
OpenFOAM-1.4.1/.OpenFOAM-1.4.1/controlDict
and set the debug switch lduMatrix to 2 and run all your cases again you will get a printout of the normalisation factor as well.
Then, you will notice that the normalisation factors are always either 1e-20 (which a thresold set in the code) or of the order of machine epsilon (1e-12 to 1e-14).
Thus, despite the "normalised residuals" appear to be of order 1, the true residuals are actually of order epsilon and thus no further convergence is really possible.
However, because the linear solvers decide whether or not to try solving based on the "normalised residual", they will try and solve the linear systems. I agree this behavior might be misleading since you maybe thought that the "tolerance" defined in the fvSolution dictionary was a truly absolute residual, while it actually specify a normalised residual (relTol beign a relative normalised residual).

See (e.g.) code in:
OpenFOAM-1.4.1/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTests.C
OpenFOAM-1.4.1/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C

Thus, although the behavior for a very special case might look odd, you shouldn't expect any problem in a real case, when the starting solution is not an exact one.

Setting nu=0 gives perfectly 0 residuals because the computation of the divergence alone results in a perfect cancellation of the fluxes, while the computation of the laplacian (in the case nu> 0) is zero only to machine epsilon.

Kind regards,

Luca

[card]

Hi Luca,

Thank you so very much for looking at the post and the useful debugging information. I do see that the residuals are normalized, but ultimately that's not the concern. Even normalized, the residuals should be 0 as all the fluxes should provide perfect cancellation as you say. With nu=0, we see this to be true. But, if you scale the grid slightly, say 1.01 as in my original post and keep nu=0, the residuals exist and you do not get perfect cancellation. This change in behaviour is what is puzzling me for the moment. So the laplacian calculation plays no role, only the divergence. In using foamCalc to calculate the divergence, the output file shows zero everywhere except in one location along the boundary where it is 1.e-14 or so. I know this is small, but if it is zero everywhere else as it should be and this is as perfect as a grid as one could hope for, I fail to see why this happens. With central difference and no NVD stuff happening I'm confused. Does this make sense? Many, many thanks for your excellent post.

All the best,

David

[listmg]

David,

It would seem that your problem is still related to precision. I reformulated your case slightly to always be scaled to a 1 meter box. I then progressed through grid sizes (1x1, 2x2, 3x3, ..., 10x10) confirming that all grids with power of 2 coordinates gave "clean" results in checkMesh vs. the other grids (3x3, 5x5, etc) gave "small" perturbations in checkMesh. Whilst these are all OK, they do introduce enough error that the fluxes do not exactly cancel and thus the problem is over-specified. This results in the non-zero residuals.

When one of the boundaries is set to zeroGradient in U and the remaining fixed U boundaries are set to have p as zeroGradient, the solution behaves as expected with some error in each of the velocity components (presumably also from the precision).

Some of the divergence behavior still is intriguing. When setting U to uniform (1 0 0) at time step 0, `foamCalc div U` will give uniform 0 for most of these cases when 1.0 is used to scale the grid, but there are some non-zero entries when 1.1 is used. This still looks to be a precision issue.

Lastly, I modified the case to run with potentialFoam. It gave unacceptable results on the USRT case with center values of velocity hovering near zero and sharp gradients at the boundaries. When one of the boundaries was relaxed as afore mentioned, potentialFoam worked as expected giving U = (1 0 0) +/- some error. Although I agree USRT should be perfect and potentialFoam should have run with the boundaries as-were, the solver clearly felt the problem was overspecified.

Regards,
Mike