Discrete approximations of $\nabla^2{\bf v}$ 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?
 A: 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.
A: 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.
