Polynomial in cellLimited<cubic> seems dodgy
 

Hello :)

After having some convergence problems in a Godunov-type scheme I programmed in OpenFOAM, I have been investigating the source code for the gradScheme cellLimited<cubic> in src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/cellLimitedGrad/gradientLimiters/cubicGradientLimiter.H

The paper it is based on (Michalak and Ollivier-Gooch, 2008) says that the limiter function should be a cubic polynomial P(r) with P(r = 0) = 0, P(r = t) = 1, P'(r = 0) = 1, P'(r = t) = 0. By my calculations, this should be:

P(r) = ar^3 + br^2 + r
a = (t - 2)/(t^3)
b = - (3at^2 + 1)/(2t)
https://www.wolframalpha.com/input/?...+and+t+%3D+1.5

But the source code says:

P(r) = ar^3 + br^2 + r
a = 2/t^2 - 2/t^3
b = - (3/2)at
https://www.wolframalpha.com/input/?...+and+t+%3D+1.5

The Wolfram|Alpha graphs (shown for t = 1.5) suggest that P'(r = t) is definitely not zero for the polynomial that is in the OpenFOAM code. Coincidentally, I initially ended up with the same coefficients as in the OpenFOAM code when I accidentally said the derivative of r with respect to r was 0 and not 1 when incorporating the constraint P'(r = t), then I plotted the graphs and realised something was wrong.

Is this a mistake in the code or do I just not understand it properly? :o

Thanks,
Shannon

Perhaps worth dropping a bug report?

I submitted a bug report https://develop.openfoam.com/Develop.../-/issues/2113 and it is getting fixed :)

Well done - the system works! Thanks for putting the effort in to get this fixed for us all.