View Issue Details Jump to Notes ] Issue History ] Print ]
IDProjectCategoryView StatusDate SubmittedLast Update
0000141OpenFOAM[All Projects] Bugpublic2011-02-17 18:102011-03-21 12:47
Assigned Tohenry 
PlatformLinuxOSFedoraOS Version12
Product VersionOther 
Target VersionFixed in Version 
Summary0000141: wrong gradient calculation using leastSquares on an unstructured hexahedral mesh
DescriptionI have a solver based on interFoam in which I try to apply a volScalarField based on a radial Gaussian distribution to an area where fvc::grad(alpha1) is non-zero, i.e. the fluid interface. The code is:
    volScalarField QT = I0*mag(fvc::grad(alpha1));
I0 is a volScalarField that has the desired Gaussian distribution not limited to the interface.

When using "leastSquares" as the discretization scheme for the gradient calculation, the resulting QT field is highly erroneous. "Gauss linear" produces the expected result (no parameters other than the default GradSchemes in fvSchemes were changed, and the I0 field is the same for both cases).
Additional InformationOF-1.6.x, latest commit: 5e99c7bac54c85ae93babf825187247620968938

Grid was generated by gridPro and converted to openfoam. checkMesh reports no errors:

Mesh stats
    points: 1829427
    faces: 5433600
    internal faces: 5379840
    cells: 1802240
    boundary patches: 6
    point zones: 0
    face zones: 0
    cell zones: 0

Overall number of cells of each type:
    hexahedra: 1802240
    prisms: 0
    wedges: 0
    pyramids: 0
    tet wedges: 0
    tetrahedra: 0
    polyhedra: 0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch Faces Points Surface topology
    left 7296 7475 ok (non-closed singly connected)
    right 7296 7475 ok (non-closed singly connected)
    bottom 12288 12417 ok (non-closed singly connected)
    top 12288 12417 ok (non-closed singly connected)
    front 7296 7475 ok (non-closed singly connected)
    back 7296 7475 ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.01 -0.002 -0.01) (0.01 0.001 0.01)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (5.382788e-17 9.360708e-16 -9.564311e-17) OK.
    Max cell openness = 3.549774e-16 OK.
    Max aspect ratio = 470.5362 OK.
    Minumum face area = 4.648805e-12. Maximum face area = 2.008524e-07. Face area magnitudes OK.
    Min volume = 2.526641e-17. Max volume = 1.48396e-11. Total volume = 1.2e-06. Cell volumes OK.
    Mesh non-orthogonality Max: 32.5824 average: 6.88462
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.166914 OK.

Mesh OK.
TagsNo tags attached.
Attached Filespng file icon ls_gs_bugreport.png [^] (157,461 bytes) 2011-02-17 18:10

? file icon LSbugrep.7z.001 [^] (1,843,200 bytes) 2011-03-21 10:43
? file icon LSbugrep.7z.002 [^] (1,725,546 bytes) 2011-03-21 10:44

- Relationships

-  Notes
henry (manager)
2011-02-17 18:19

Have you tried OpenFOAM-1.7.x ?

Any work we do on this will be in 1.7.x so we will need the details of behavior of this version rather than any previous releases.
akidess (reporter)
2011-02-17 18:33

I have not. I'll test 1.7.x and report back.
akidess (reporter)
2011-02-20 21:22

Results with 1.7.x look the same as the images posted for 1.6.x.
henry (manager)
2011-02-20 21:27

Could you send a small test code and case for us to investigate this issue?
akidess (reporter)
2011-02-25 08:45

Apologies for the delayed response, I was expecting to automatically get a notification if there are new notes. I will put together a test case for you this weekend. Do you have an ftp I can upload the files to (you can email login information to A.Kidess at Even compressed the case weighs about 130 MB.
henry (manager)
2011-02-25 08:51

It would be useful if you could create a small case from what you have which reproduces the problem and post it here. If the mesh is generated automatically with a tool we supply then the pack need not include the mesh only the means to generate it.
akidess (reporter)
2011-02-25 08:55

The mesh was generated by a third party using gridPro. I'm not sure if I can get a smaller mesh , but I'll see what I can do.
akidess (reporter)
2011-03-21 10:49

Sorry I couldn't get the case and solver smaller than 3.5 MB, so I had to split the archive in two. Steps to reproduce the error:

Extract the archives by using "7za x LSbugrep.7z.001". Compile the solver (which I stripped down to the necessary), then run the solver "sp3d" in the case directory. It will run one time step and write out a solution file.

Now, in the postprocessor, set a threshold of alpha = 0.4 and display I0 / QT on the resulting domain. The field I0 will have the expected smooth distribution, while the field QT which is I0 * mag(grad(alpha1)) is not when using leastSquares to approximate gradients. Using Gauss linear the result is as expected.
henry (manager)
2011-03-21 12:19

Your case is extreme both in terms of the cell aspect ratios and the gradients being evaluated. It looks like with least-squares the round-off errors accumulate to such an extent that the gradient has a large error. This is partly due to the area weighting currently used but also the round-off errors accumulated in the tensor inversion.

There are various choices for the weighting in the least-squares; the current version being carefully constructed to work well for a range of 1D, 2D and 3D meshes but it may not be ideal for your case. A simpler approach is still used in the experimental extendedLeastSquares scheme, try
default extendedLeastSquares 0;

which behave much better for your case than leastSquares but not possibly not as well as Gauss linear. It may be for your mesh that the Gauss linear scheme is the best or for greater accuracy you could use a higher order interpolation scheme in conjunction with the Gauss gradient scheme.
akidess (reporter)
2011-03-21 12:36

The uploaded mesh was coarsened automatically and the aspect ratios became much more extreme than on the original grid (1850 vs. 470 max). As you expected extendedLeastSquares behaves MUCH better (both on the coarsend and the original grid), though not quite as good as Gauss linear.

If we can conclude that the erroneous result is due to limitations of the least squares method on this grid and not a programming bug then I'm happy with the solutions you proposed and we can close this report.
henry (manager)
2011-03-21 12:47

Yes, this is not a programming bug but a limitation of the current leastSquares implementation which accumulates round-off error on extreme meshes.

The Gauss gradient method is more robust and can be made more accurate on distorted and irregular meshes by choosing a higher-order interpolation scheme.

- Issue History
Date Modified Username Field Change
2011-02-17 18:10 akidess New Issue
2011-02-17 18:10 akidess File Added: ls_gs_bugreport.png
2011-02-17 18:19 henry Note Added: 0000253
2011-02-17 18:33 akidess Note Added: 0000254
2011-02-20 21:22 akidess Note Added: 0000256
2011-02-20 21:27 henry Note Added: 0000257
2011-02-25 08:45 akidess Note Added: 0000266
2011-02-25 08:51 henry Note Added: 0000267
2011-02-25 08:55 akidess Note Added: 0000268
2011-03-21 10:43 akidess File Added: LSbugrep.7z.001
2011-03-21 10:44 akidess File Added: LSbugrep.7z.002
2011-03-21 10:49 akidess Note Added: 0000292
2011-03-21 12:19 henry Note Added: 0000293
2011-03-21 12:36 akidess Note Added: 0000294
2011-03-21 12:47 henry Note Added: 0000295
2011-03-21 12:47 henry Status new => closed
2011-03-21 12:47 henry Assigned To => henry
2011-03-21 12:47 henry Resolution open => fixed