I am writing a Navier Stokes solver. The vector field is represented as a grid with integer coordinates

I am looking at other people's computer code. I don't entirely understand the vector calculus, but if I am interpreting it correctly, I see the $\nabla^2{\bf v}$ term approximated (in the three-dimensional case) as

$\left( \begin{array}{c} \sum_{i=\pm 1}v_x(x+i,y,z)+v_x(x,y+i)+v_x(x,y,z+i) \\ \sum_{i=\pm 1}v_y(x+i,y,z)+v_y(x,y+i)+v_y(x,y,z+i) \\ \sum_{i=\pm 1}v_z(x+i,y,z)+v_z(x,y+i)+v_z(x,y,z+i) \\ \end{array} \right) - 6{\bf v}(x,y,z)$

This term needs to be scaled by the resolution of the grid but it is essentially the sum of the gradients coming into the grid cell at each edge.

First question, have I interpreted this correctly?

Second question: what is a better approximation of this term?

The above approximation derives from the simple practical fact that each cell in the grid normally has direct access only to its immediate lateral neighbours. This simplifies the code. However, constraints on the grid resolution and time step, combined with sometimes-large viscosity and diffusion parameters, means the simulation sometimes behaves badly.

I would like to try a better approximation of the diffuse term, considering diagonal neighbours and neighbours further than one grid step away. (Existing numerical tricks that I have seen seem to interact badly with my boundary conditions, so I would like to build a solution from the ground up.)

What should I use?

  • 1
    $\begingroup$ I assume $v_x$ is the $x$ component of $v$ and not a derivative? $\endgroup$ – k1next Feb 26 '13 at 17:01
  • $\begingroup$ essentially the sum of the gradients coming into the grid cell at each edge. I would say, essentially the sum of second derivatives. For example, the expression $v_x(x+1,y,z)-2v_x(x,y,z)+v_x(x-1,y,z)$ is a discrete form of the second derivative in the $x$ direction. $\endgroup$ – user53153 Feb 26 '13 at 18:54

You are right with your assumption for question 1). For simplicity, let us consider $u$ as a scalar function instead of a vector $v$.Assuming an equdistant grid in all 3 dimension, we then have \begin{align} \nabla^2u(x,y,z) = \Delta u (x,y,z)&= \partial_{xx}u+\partial_{yy}u+\partial_{zz}u\\ &\approx \frac{u(x+h,y,z)-2u(x,y,z)+u(x-h,y,z)}{h^2} \\ &+ \frac{u(x,y+h,z)-2u(x,y,z)+u(x,y-h,z)}{h^2}\\ & +\frac{u(x,y,z+h)-2u(x,y,z)+u(x,y,z-h)}{h^2}\\ &=\frac{1}{h^2}\Bigl( u(x+h)+u(x-h)+u(y+h)+u(y-h)+u(z+h)+u(z-h)\Bigr) \\ & -\frac{6}{h^2} u(x,y,z) \end{align} Which is exactly the same you have for $h=1$ and a vector $v$ instead of $u$.

derives from the simple practical fact that each cell in the grid normally has direct access only to its immediate lateral neighbours

This is correct. This are the neighboured cells of a cube.

Concerning the second question of yours, I am not sure if this is a standard approach. For finite elements, the above is quite intuitive, as it represents the flux through the faces of an element ( a cube for 3D equidistant). But one possible approach would be the following (Reference): \begin{align} f''(x) &\approx \frac{-f(x+2 h)+16 f(x+h)-30 f(x) + 16 f(x-h) - f(x-2h)}{12 h^2} \end{align} which needs to be adapted into 3D. This combines the information out of more cells than the direct neighbours. And is called 5-point stencil. Of course you would need to check whether you are too close to the boundary to use this formula (direct or indirect neighbour to boundary).

You can get some further information if you simply google for 7 point stencil or 27 point stencil. The below posted hyperlinks somehow don't work.

I also found a 9-point stencil for 1D cases here in exercise 4. As it would be a 27-point stencil in 3D Google gave me for example this.

  • $\begingroup$ Somehow my hyperlinks are not formatted correctly $\endgroup$ – k1next Feb 26 '13 at 17:35
  • $\begingroup$ Sorry I don't know why that is not working. But basically I googled 9 point stencil and 27 point stencil. I will remove the hyperlinks from my answer $\endgroup$ – k1next Feb 26 '13 at 17:38
  • $\begingroup$ The URLs were Google search URLs, apparently you copied them from Google search result page. If you instead navigated to the search result and took URL from the browser address line, this would not happen. $\endgroup$ – user53153 Feb 26 '13 at 18:52
  • $\begingroup$ What a stupid mistake. Sorry. Thank you very much for your edit. $\endgroup$ – k1next Feb 26 '13 at 19:11

How are you performing the smoothing step?

If you are integrating $\nabla^2$ explicitly, and are encountering instability, then a larger or different stencil will not help you. I would recommend instead solving the viscous step implicitly; see for instance page 3, right column, of this paper.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.