Incorrect centered finite differences expression in colors module

[rhdtownsend]

In the compute_derivatives_at_point() subroutine of the hermite_interp.f90 module in colors, the expression used to calculate centered finite differences,

df_dx = (f(i + 1, j, k) - f(i - 1, j, k))/(2.0_dp*dx)

will give incorrect results when the grid isn't uniformly spaced -- you can appreciate this from the fact that dx = x(i+1) - x(i), and so x(i-1) is not used.

As an example, applying this formula to the function y(x) = 3x + 1 sampled at x = 0, 1, 1.5 yields a slope 3.5, rather than 3.

[Debraheem]

When I plug in your example for dx = 0.5, I get 4.5, not 3.5, so I think it's even worse. and 2.25 if using dx = 1.

Is the solution here to use a non-uniform 3 point difference formula like the one Numpy uses for their gradient function?

https://numpy.org/devdocs/reference/generated/numpy.gradient.html

[rhdtownsend]

Yeah, I must have messed up the calculation!

MSG in effect uses this formula for calculating centered finite differences, but constructs it using a different approach -- namely, by fitting a parabola through the points f(x-hs), f(x) and f(x+hd), and then evaluating the slope at f(x). This leads to

df_dx = [hs/hd f(x+hd) + (hd/hs - hs/hd) f(x) - hd/hs f(x-hs)] / (hs + hd)

...which is mathematically identical to the formula above.

[nialljmiller]

This was fixed in PR #955