9-point stencil “equivalent” for advection equation So I inherited from some people a code that solves the advection-diffusion-reaction equation for a particular system. The original code was first implemented in 1D which worked fine in cartesian coordinates. However, now that I'm trying to implement it in 2d I'm having a few problems.
The equation is:
$$\frac{\partial u}{\partial t} + \nabla \cdot \left( \boldsymbol{v} u - D\nabla u \right) = f$$
where $u$ are the primitive variables being evolved, $D$ is a diffusion coefficient, and $f$ is source+sink. 
By separating the spatial discretization from the temporal through a method of lines and using a cartesian coordinates it becomes relatively straightfrward to extend a model from 1D -> 2D because the spatial components are additive. However, the problem is that the genuine solutions to the equations I'm solving tends be somewhat spherically symmetric (but not completely), but the grid is cartesian so the solution tends to show grid orientation errors in the cartesian grid, i.e. the solution looks kinda "square" instead of looking like a "circle". A CPU costly solution to this grid orientation error is to make the grid finer, but it's getting to the point the computations become very slow.
I was able to solve the 2D grid orientation errors for the diffusion part of the advection-difussion-reaction equations by discretizing the laplacian operator into a 9-point stencil instead of the original 5-point stencil, the latter 5-point stencil is acquired by just naively adding the 1D laplacian discretations of the x and y components. The 9 point stencil adds spatially "diagonal" components to the discretization. However, the advection part still has some pretty bad grid orientation problems.
I was wondering if there is a "9-point stencil" analog to the 2D advection problem. Maybe a stencil that takes into account the "diagonal" cross-term spatial components somehow.
 A: The method you are looking for is called corner transport upwind (CTU) (Google search).
Typically, we get the finite difference schemes via Taylor expansion, keeping 1st order terms and ignoring the rest (as either small, or naturally via subtraction a la central difference). The CTU scheme keeps these 2nd order terms, so the expansion is
\begin{align}
Q\left(t+\Delta t,x,y\right)&=Q-\Delta t\left[uQ_{,x}+vQ_{,y}\right] \\
&+\frac12\Delta t^2\left[u^2Q_{,xx}+vuQ_{,xy}+uvQ_{,yx}+v^2Q_{,yy}\right]+\mathcal{O}\left(\Delta t^3\right)
\end{align}
where the comma denotes partial differentiation. To first order, you keep the $uv$ cross terms in addition to the top line (i.e., ignoring the $Q_{,xx}$ and $Q_{,yy}$ terms). This would result in something like,
\begin{align}
Q_{ij}^{n+1}&=Q_{ij}^n-\frac{u\Delta t}{\Delta x}\left[Q_{ij}^n-Q_{i-1j}^n\right] -\frac{v\Delta t}{\Delta y}\left[Q_{ij}^n-Q_{ij-1}^n\right] \\
&\qquad+uv\frac{2\Delta t^2}{\Delta x\Delta y}\left[\left(Q_{ij}^n-Q_{ij-1}^n\right)+\left(Q_{i-1j}^n-Q_{i-1j-1}^n\right)\right.\\
&\qquad\qquad\qquad\qquad\left.+\left(Q_{ij}^n-Q_{i-1j}^n\right)+\left(Q_{ij-1}^n-Q_{i-1j-1}^n\right)\right]
\end{align}
To get second order accuracy, the whole $\Delta t^2$ term needs to be kept. This results in something like (moving to central differences, rather than backwards)
\begin{align}
Q_{ij}^{n+1}&=Q_{ij}^n-\frac{u\Delta t}{2\Delta x}\left[Q_{i+1j}^n-Q_{i-1j}^n\right] -\frac{v\Delta t}{2\Delta y}\left[Q_{ij+1}^n-Q_{ij-1}^n\right] \\
&\qquad +\frac{u^2\Delta t^2}{2\Delta x^2}\left(Q_{i+1j}^n-2Q_{ij}^n+Q_{i-1j}^n\right)+\frac{v^2\Delta t^2}{2\Delta y^2}\left(Q_{ij+1}^n-2Q_{ij}^n+Q_{ij-1}^n\right)\\
&\qquad+uv\frac{2\Delta t^2}{\Delta x\Delta y}\left[\left(Q_{ij}^n-Q_{ij-1}^n\right)+\left(Q_{i-1j}^n-Q_{i-1j-1}^n\right)\right.\\
&\qquad\qquad\qquad\qquad\left.+\left(Q_{ij}^n-Q_{i-1j}^n\right)+\left(Q_{ij-1}^n-Q_{i-1j-1}^n\right)\right]
\end{align}
(you'll want to check indices, I've likely screwed them up here). This is stable so long as the following CFL condition is satisfied:
$$
\max\left(u\frac{\Delta t}{\Delta x},\,v\frac{\Delta t}{\Delta y}\right)\leq1
$$
This is also fairly common in the hydrodynamic codes I've seen & used. Randy LeVeque's Finite Volume Methods for Hyperbolic Problems (Amazon link) is probably an excellent resource for you, Chapter 20 (I think) is where he discusses this method. 
The stencils for each of these are shown below, adapted from these lecture notes (NB: 45-page PDF) on conservation laws from James Rossmanith at Iowa State University,

So you can see how the 1st order is actually a 6-point stencil and the 2nd order is a 9-point stencil.
