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, integrating time with a Runge-Kutta 4th order scheme, and using a cartesian coordinate discretation it becomes relatively straightforward 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 "transverse" 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 "transverse" cross-term spatial components somehow. I know of Corner Transport Upwind (CTU) which has a transverse component:

\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}

However I'm not sure if it can be used with MOL and Runge-Kutta 4th order because it depends explicitly on time e.g. $\Delta t$, while it seems that cell discretization schemes used for MOL only depend on the spatial components so that time can be integrated separately. However I might be wrong and CTU could be perfectly appropriate. I would be really grateful for some guidance here because I'm a noob!

  • $\begingroup$ Could you use a stencil to do a multidimensional least-squares fit? That will account for transverse components. $\endgroup$ Commented Jul 28, 2015 at 7:49
  • $\begingroup$ Can you elaborate? You mean Finite Element Methods? Can you point to a specific stencil I could use? Thanks. $\endgroup$
    – mathdummy
    Commented Jul 28, 2015 at 23:32
  • $\begingroup$ I've seen it done in finite volume before. Do you have access to dx.doi.org/10.1175/MWR-D-14-00054.1 and dx.doi.org/10.1175/2009MWR2917.1 which describe a 2D cubic upwind-biased stencil. $\endgroup$ Commented Jul 29, 2015 at 6:57
  • $\begingroup$ I saw those papers. I'm having a hard time understanding what they're doing. They have some complicated geometries and I'm basically just looking for a uniform cartesian grid stencil. $\endgroup$
    – mathdummy
    Commented Jul 30, 2015 at 8:40
  • $\begingroup$ They should be reasonably straightforward to adapt to uniform cartesian grids. I happen to be implementing the scheme from Weller & Shahrokhi 2014 on such a grid in Python, but in a somewhat simplified form. Happy to share the code once it's debugged and working. $\endgroup$ Commented Jul 30, 2015 at 11:19

1 Answer 1


You could use a stencil to do a multidimensional least-squares fit which will account for transverse components. I've seen this done in finite volume before, [Weller et al. 2009, Weller & Shahrokhi 2014] give details of using upwind-biased stencil on unstructured two-dimensional meshes, but the same techniques can be straightforwardly applied to uniform Cartesian meshes.

I have seen this technique used with flux-form advection that gives stable results, but my own tests with advective form suggest that such a scheme is not stable, likely due to Godunov's theorem.


  • [Weller et al. 2009] Voronoi, Delaunay, and block-structured mesh refinement for solution of the shallow-water equations on the sphere. Monthly Weather Review, 137, 4208-4224, dx.doi.org/10.1175/2009MWR2917.1
  • [Weller & Shahrokhi 2014] Curl-free pressure gradients over orography in a solution of the fully compressible Euler equations with implicit treatment of acoustic and gravity waves. Monthly Weather Review, 142, 4439-4457, dx.doi.org/10.1175/MWR-D-14-00054.1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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