Discretization method for advection equation without numerical diffusion

Question

Given the advection equation for an incompressible flow field $$\frac{\partial c}{\partial t} + \mathrm{Pe} \frac{\partial c}{\partial x} = 0$$

what would the best method be for discretizing this without introducing any numerical diffusion or oscillations? Specifically when we have step changes in boundary conditions, and time (and space) dependent velocity $v(x,t)$.

The book Numerical Methods for Problems with Moving Fronts by Bruce A. Finlayson goes into great detail on this problem when you have step changes as boundary conditions.

He recommends filtered leapfrog as the best finite difference method: $$c^{n+1}_i = \frac{\alpha}{2}\left(c_i^n + c_i^{n+2}\right) + (1-\alpha)c_i^{n-1} - \frac{\mathrm{Pe}\Delta t}{\Delta x}\left(c_{i+1}^n - c_{i-1}^n\right)$$

(filtered means $\alpha = 1$) and Taylor-Galerkin as the best finite element method:

$$ \frac{1}{6}\left(c_{i+1}^{n+1} - c_{i+1}^{n}\right) + \frac{2}{3}\left(c_{i}^{n+1} - c_{i}^{n}\right) + \frac{1}{6}\left(c_{i-1}^{n+1} - c_{i-1}^{n}\right) = -\frac{\mathrm{Pe}\Delta t}{\Delta x} \left(c^{n}_{i+1} - c^{n}_{i-1}\right) + \frac{\mathrm{Pe}^2\Delta t^2}{\Delta x^2}\left(c^{n}_{i+1} - 2c_i^n + c^{n}_{i-1}\right)$$

but I was wondering what the consensus is today?

I have tested the mentioned methods, and they both suffer from numerical and oscillation. But perhaps this can be managed by adjusting $\Delta t$ and $\Delta x$ in relation to the velocity?

This answer to a related question regarding Crank-Nicholson and the advection equation states that

Crank-Nicolson is not necessarily the best method for the advection equation. It is second order accurate and unconditionally stable, which is fantastic. However it will generate (as with all centered difference stencils) spurious oscillation if you have very sharp peaked solutions or initial conditions.

but gives no alternative methods.

Answer 1

If your aim is to solve one-dimensional advection equation with variable velocity having minimal numerical diffusion and no unphysical oscillations and you have no other requirements (mass conservation?!) then in my opinion the so called semi-Lagrangian schemes may be your choice. By the way, using the notation Pe in your equation is a little bit confusing to me as it reminds a "Peclet" number that relates advection and diffusion.

The semi-Lagrangian methods are based on the fact that the solution is constant along characteristic curves (the "characteristics"), therefore to determine $c_i^{n+1}$ you have to find a characteristic $X(t)$ that ends in time $t^{n+1}$ in the grid point $x_i$. You can obtain it numerically by solving ordinary differential equation $$ \dot X(t) = v(X(t),t) , \,\, X(t^{n+1})=x_i \,, \,\, t \in (t^n,t^{n+1}) $$ back in time. Computing then the value of numerical solution at time $t^n$ in the position $X(t^n)$ will give you the value of $c_i^{n+1}$. This last task shall be done by some proper numerical interpolation of nodal values $c_j^n$ that will give you no unphysical values (e.g. a linear interpolation).

If the characteristics will cross the boundary of your domain (interval) for some $\tilde t \in (t^n,t^{n+1})$ then you must exploit the boundary condition, typically a given value at time $\tilde t$.

This class of methods has a great advantage of no formal restriction on the time step of method. Of course, you must work with velocity $v(x,t)$ that is smooth enough. If you have sinks or sources in your velocity field, e.g. points where two characteristics meet or start, then some careful treatment in the method shall be realized.

If you for some reasons prefer purely Eulerian type of methods (e.g. finite difference or finite volume or finite element method) then you should be aware that if you allow no unphysical oscillations then you can have at most second order accurate method, so to obtain accurate numerical solution ("less numerical diffusion"), fine grids shall be used.

Addendum

The method called "Taylor-Galerkin finite element method" (FEM) in the question is closely related to the so called Lax-Wendroff method for finite volume (or difference) method. In fact, if one uses a mass lumping in the left hand side of FEM, so it takes the form $c_i^{n+1}-c_i^n$, then it is identical to Lax-Wendroff method for constant Pe. This method is second order accurate, but it results in numerical oscillations. To avoid them, a so called limiter must be used that reduces the scheme in some grid points to the direction of first order accurate scheme. This topic is very well described in e.g. book of R. LeVeque on Finite Volume Methods for Hyperbolic Problems.

A better choice than Lax-Wendroff concerning the accuracy (a so called phase error) is the so called Fromm method that includes in the right hand side also the numerical values $c_{i-2}^n$ or $c_{i+2}^n$ depending on the sign of Pe. To have no oscillations, again a limiter must be used, but this is the case for any second order accurate scheme. This method would be my choice if its restriction on the time step is not an issue (the CFL condition).

More advanced option than the Fromm scheme with a limiter is the approach e.g. of Chi-Wang Shu and collaborators On maximum-principle-satisfying high order schemes for scalar conservation laws or Maximum-principle-satisfying high order finite volume weighted essentially nonoscillatory schemes for convection-diffusion equations. This would be my recommendation if more time to the implementation can be invested.