I'm learning about Stable Fluids from the document Stable Fluids and I have a question regarding the implementation of the project function.
In the code, the divergence is calculated as follows:
div[IX(i, j)] = -0.5 * h * (u[IX(i+1, j)] - u[IX(i-1, j)] + v[IX(i, j+1)] - v[IX(i, j-1)]);
My understanding of the definition of divergence is:
$$
\nabla \cdot \mathbf{v} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}
$$
As the distance between the two points is 2 * h, I would expect the first term to be approximated by:
(u[IX(i+1, j)] - u[IX(i-1, j)]) / (2 * h)
rather than:
(u[IX(i+1, j)] - u[IX(i-1, j)]) / 2 * h
However, when I implement my version using /(2*h), the code does not work correctly and variables such as u, v, and p become NaN after a few time steps. In contrast, using *h as in the document works perfectly.
Can someone explain why multiplying by h is correct in this context?
