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!
