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.

• Can you write the particular equation you are trying to solve? – Muphrid Jul 24 '15 at 0:19
• Edited my question and added the explicit equation. – mathdummy Jul 24 '15 at 0:24

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.

• I hope you already had all of those equations typed up somewhere and could just copy them over. Otherwise, that's way too much work MathJax'ing an answer! – tpg2114 Jul 24 '15 at 1:39
• Actually, I did type them out, but they're so similar that I was able to copy & paste most of the (hence the suggestion to check indices). – Kyle Kanos Jul 24 '15 at 1:40
• Damn, this is exactly what I was looking for! My google-fu wasn't very good. Thanks! – mathdummy Jul 24 '15 at 2:02
• Hi again. I realized that this CTU scheme might be problematic to numerical scheme I have at the moment. I'm using method of lines and runge kutta fourth order, which means that I need a scheme where space and time can be decoupled so that they're separately computed (method of lines). So still this question is open. Thanks. – mathdummy Jul 28 '15 at 1:57
• @mathdummy: As I understand MOL, you are discretizing space & then transforming PDEs into ODEs to use RK integration scheme. Perhaps there is a reason, but I do not see why you couldn't use Taylor expansion of CTU in MOL. – Kyle Kanos Jul 28 '15 at 3:06