|Anonymous | Login | Signup for a new account||2014-10-01 12:11 UTC|
|My View | View Issues | Change Log | Roadmap|
|View Issue Details|
|ID||Project||Category||View Status||Date Submitted||Last Update|
|0000141||OpenFOAM||[All Projects] Bug||public||2011-02-17 18:10||2011-03-21 12:47|
|Target Version||Fixed in Version|
|Summary||0000141: wrong gradient calculation using leastSquares on an unstructured hexahedral mesh|
|Description||I 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 Information||OF-1.6.x, latest commit: 5e99c7bac54c85ae93babf825187247620968938|
Grid was generated by gridPro and converted to openfoam. checkMesh reports no errors:
internal faces: 5379840
boundary patches: 6
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
tet wedges: 0
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)
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.
|Tags||No tags attached.|
|Attached Files|| ls_gs_bugreport.png [^] (157,461 bytes) 2011-02-17 18:10
LSbugrep.7z.001 [^] (1,843,200 bytes) 2011-03-21 10:43
LSbugrep.7z.002 [^] (1,725,546 bytes) 2011-03-21 10:44
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.
|I have not. I'll test 1.7.x and report back.|
|Results with 1.7.x look the same as the images posted for 1.6.x.|
|Could you send a small test code and case for us to investigate this issue?|
|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 tudelft.nl)? Even compressed the case weighs about 130 MB.|
|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.|
|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.|
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.
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.
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.
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.
|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|