ENH: New gradient scheme iterativeGaussGrad

Merged Kutalmış Berçin requested to merge feature-skew-corrected-gradient into develop



Improvement of numerical fidelity only via improving gradient computations for meshes with non-zero face skewness.

(Theory) What is face skewness:

  • Level of closeness of the following two spatial points:​
    • Face intersection of the PN vector between cell centres of an owner cell and a neighbour cell​
    • Face centre of the common face of the two cells

(Theory) Why is face skewness important:

  • Prediction fidelity
    • Any deviation from face centre as the face average centre reduces the second-order accuracy (i.e. non-zero face skewness)​
    • Hence OpenFOAM has various skewness corrections for centre-to-face interpolations
  • Numerical stability
    • Reduces the diagonal dominance of the discrete Poisson operator which slows down convergence. ADI (i.e. alternating-direction implicit) type preconditioners also become less effective.

(Theory) What is the technical challenge in face skewness:

  • Omission wherever possible due to cost concerns rather than a true technical challenge​
    • Gradient computations are accounted for ~20% of a typical simulation​

(Theory) Treatments in OpenFOAM:

  • Many methods to quantify skewness
    • OpenFOAM has its own definition​
    • Skewness formulations in checkMesh and snappyHexMesh are a bit different
  • Various skewness corrections available​
    • skewCorrectedSnGrad surface-normal gradient scheme​
    • skewCorrected surface interpolation scheme​
    • No corrections on boundary skewness​
    • No corrections in Gauss-Green gradient computations​
    • Not needed for leastSquares and pointLinear based gradient schemes

Resolved bugs



Problem definition:

  • Face skewness is not accounted during any Gauss-Green gradient computations.​
  • Incorporation of face skewness into gradient computations may or may not improve the level of prediction fidelity.​

A problem solution: Iterative Gauss gradient

  • Add an outer loop around the gradient and skew-correction vector computations​
  • Add an optional relaxation factor inside the gradient computation

Test case:

  • Very difficult to design a test case where skewness is isolated from non-orthogonality​
  • Manufactured solution: A two-dimensional triangle-cell domain where theoretical gradient is square-root of 2 in each co-ordinate direction


  • Numerical stability
    • Simulation crashes or not
  • Prediction fidelity - Level of discrepancy between theoretical and numerical values:​
  • Arithmetic average of error field​
    • Coefficient of variation (CoV) of error field​
    • Relative standard deviation: std/mean​
    • Measure of amount of dispersion​
    • Low CoV: Values are close to mean​
    • High CoV: Values spread out over a wider range, potentially indicating outliers with respect to the mean​

Control variable:

  • Skewness is not the control variable since changing skewness also changes non-orthogonality in general​
  • Therefore:​
    • We fixed skewness and non-orthogonality with a fixed-geometry test case​
    • Control variable becomes the gradient scheme itself​
    • We quantified the performance of each gradient scheme with the chosen metric and compared them with each other

Mesh metrics: image










What happened?

  • New iterative Gauss scheme​
    • Reduced errors in internal fields in comparison to other Gauss schemes​
    • Reduced errors in boundary fields in comparison to least square schemes​
  • Increasing number of iterations monotonically reduced errors​
  • Application of interpolation-scheme limiters increased errors

What do the results mean in practise? (How will it change in how we do things?)​

  • The new scheme may be used to improve prediction fidelity for gradient computations on meshes with skewness and non-orthogonality​
  • More tests are needed for:​
    • Effects of extreme levels of face skewness​
    • Cost estimations​
  • The results do not suggest that the new scheme can improve numerical stability


  • No changes in existing output.
  • No changes in existing user input.


  • No finite-area
  • Needs extensive tests to assume safe for single-phase, multiphase and overset finite-volume applications​
  • No boundary treatments for Dirichlet and Neumann conditions


  • Compilation (incl. submodules):
    • linux64ClangDPInt32Opt (clang11)
    • linux64GccDPInt32Opt
    • linux64GccSPDPInt64Debug
  • Alltest: No change in output with respect to the develop HEAD + no error
Edited by Kutalmış Berçin

Merge request reports