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/)
-   -   non-exact gradient values (http://www.cfd-online.com/Forums/openfoam-programming-development/104863-non-exact-gradient-values.html)

elisabet July 17, 2012 10:39

non-exact gradient values
 
1 Attachment(s)
Hi OpenFOAMers!
I've recently found a really strange behaviour on the calculation of the gradients on version 1.6-ext:

Starting from a uniform field, called P (partial pressure), which is uniform and equals 1, even at the boundaries, I want to evaluate the gradient of P, which is expected to be zero. However, it is not!!

Details on the case and code:

1.- the mesh is uniform, hence non-orthogonal corrections are not required. CheckMesh is o.k.
2.- boundary conditions of P are own made based on fixedValue type but, for this simple verification case, the value is set to 1 (and can be checked when looking at P values on the boundaries)
3.- two different strategies are used for the evaluation of the gradient:
3.a.
Code:

sndPdx=(fvc::snGrad(P) * mesh.magSf());
    dPdx= fvc::reconstruct(sndPdx);

which I guess is better for non-orthogonal meshes
3.b.
Code:

dPdx= fvc::grad(P);
In the attached figure you can find P map, dPdx_sn (from option 3.a.) and dPdx (from option 3.b.). As the mesh is orthogonal, both options yield to the same results. But the gradient is NOT zero!!!!

What's wrong? Am I going crazy?

any help is appreciated....

elisabet

akidess July 17, 2012 11:26

The gradient is (Phi_P - Phi_f)/mag(d). I'll assume d=1e-4m, and you get a small error in the computation of the difference of 1e-16 (that's about the precision you can get using double if I'm not mistaken). Dividing that by d leads to 1e-12, about the order of errors you are seeing.

Further reading: http://en.wikipedia.org/wiki/Machine_epsilon

elisabet July 18, 2012 07:38

Thanks akidess!

Obviously that was the point! sometimes, when one doesn't get the right solution begins to be suspicious of everything... sorry!

now the bad new: I must have a mistake somewhere else in my code... let's see...

elisabet

akidess July 18, 2012 07:45

No worries and good hunting!


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