# How do I properly implement a hyperbolic TVD flux limiter?

I'm trying to implement a TVD flux limiting scheme for the hyperbolic transport equation. However, I am currently stuck at implementing the limiter. How would I properly do this?

Here's how far I got:

$$\partial_t u(x,t) = -\partial_xF\left(u,x,t\right)$$

Using a flux conservative discretization, we get:

$$\frac{u^{k+1}_i-u^k_i}{\Delta t} = -\frac{1}{\Delta x_i}\left(F^k_{i+1/2}-F^k_{i-1/2}\right)$$

The right and left side flux are calculated by mixing of a high-order scheme $$f^H$$ and a lower-order scheme $$f^L$$:

$$F^k_{i\pm1/2}=\left(1-\phi_{i\pm1/2}\right)f^L_{i\pm1/2}+\phi_{i\pm1/2}f^H_{i\pm1/2}$$

The functions $$\phi_{i\pm1/2}$$ are called the flux limiters, and that's where my problem is! The articles I find, all assume $$F$$ to be positive and use the following approach:

$$\phi_{i+1/2}=\phi\left(r_i\right),\quad r_i=\frac{F^k_i-F^k_{i-1}}{F^k_{i+1}-F^k_i}\\ \phi_{i-1/2}=\phi\left(r_{i-1}\right),\quad r_{i-1}=\frac{F^k_{i-1}-F^k_{i-2}}{F^k_i-F^k_{i-1}}$$

For $$\phi\left(r\right)$$ they use one of the following functions: https://en.wikipedia.org/wiki/Flux_limiter

I tried to implement TVD flux limiting with the low-order Osher scheme and the high-order Lax-Wendroff scheme, however my solution is somehow broken. How can I fix this?

My schemes: $$f^L_{i+1/2}=\min{\left(0,F^k_{i+1}\right)}-\min{\left(0,F^k_i\right)}\\ f^L_{i-1/2}=\max{\left(0,F^k_{i-1}\right)}-\max{\left(0,F^k_i\right)}\\ f^H_{i+1/2}=\frac{F^k_{i+1}+F^k_i}{2}-\frac{\Delta t}{2\Delta x_i}A^k_{i+1/2}\left(F^k_{i+1}-F^k_i\right)\\ f^H_{i-1/2}=\frac{F^k_i+F^k_{i-1}}{2}-\frac{\Delta t}{2\Delta x_i}A^k_{i-1/2}\left(F^k_i-F^k_{i-1}\right)$$

(The $$A^k_{i\pm1/2}$$ are the Jacobian $$\partial_u F\left(u,x,t\right)$$ at $$x_{i\pm1/2},u^k_{i\pm1/2}$$.)

My function: $$F\left(u,x,t\right)=u$$

My initial condition: $$u\left(x,t=0\right)=\begin{cases} 1,-\frac{1}{2}\le x\le\frac{1}{2}\\ 0,\text{everywhere else} \end{cases}$$

My result after 60 iterations with $$\Delta t=0.02$$:

Not what I expected...

EDIT:

If I use $$\phi_{i-1/2}=\phi\left(r_i\right)$$, the solution is way better:

But breaks down for negative function values: $$F\left(u,x,t\right)=-u$$

Let us describe the flux-limiter scheme for the linear advection equation $$u_t = -u_x$$, as done in the book (1) p. 176, §16.2 (see also the article (2)). The numerical flux is defined as a combination of the upwing and Lax-Wendroff fluxes: $$F_{i+1/2}^k = u_i^k + \phi(r_i) \left(1 - \frac{\Delta t}{\Delta x_i}\right) \frac{u_{i+1}^k - u_{i}^k}{2} \, , \qquad r_i =\frac{u_i^k - u_{i-1}^k}{u_{i+1}^k - u_{i}^k}\, ,$$ where $$\phi$$ is deduced from a limiter function --e.g. the minmod limiter $$\phi(r) = \max(0, \min(1,r))$$-- and the smoothness indicator $$r_i$$ is the ratio of consecutive gradients.

For arbitrary wave speeds $$a$$ such that $$u_t = -au_x$$, the method has numerical flux $$F_{i+1/2}^k = a\frac{u_i^k + u_{i+1}^k}{2} - |a|\frac{u_{i+1}^k - u_{i}^k}{2} + \phi(r_i) \left(\operatorname{sgn} a - a\frac{\Delta t}{\Delta x_i}\right) a \frac{u_{i+1}^k - u_{i}^k}{2} \, ,$$ and the smoothness indicator is the ratio of slopes in the upwind direction $$r_i = \frac{u_{I+1}^k - u_{I}^k}{u_{i+1}^k - u_{i}^k}, \qquad I = i - \operatorname{sgn} a \, .$$ Below is a Matlab implementation of it

% parameters
N = 400;
CFL = 0.95; % Courant number
a = 1;
tend = 1;

dx = 4/N;
x = ((-2-2*dx):dx:(2+2*dx))';
Nx = length(x);
dt = CFL*dx/abs(a);

% initialization
t = 0;
u0 = @(x) (-0.5<x).*(x<0.5);
u = u0(x);

% iterations
while (t+dt)<tend
up = circshift(u,-1);
upp = circshift(up,-1);
um = circshift(u,+1);

phi = zeros(Nx,1);
r = ones(Nx,1);
if a>=0
num = u-um;
den = up-u;
else
num = upp-up;
den = up-u;
end
for i=1:Nx
if abs(den(i))>=eps
r(i) = num(i)./den(i);
end
end
phi = min([ones(Nx,1), r],[],2);
phi = max([zeros(Nx,1), phi],[],2);

Fp05 = 0.5*(a*(u+up)-abs(a)*(up-u)) + 0.5*(1-CFL)*abs(a)*(up-u).*phi;
Fm05 = circshift(Fp05,+1);

u = u - dt/dx*(Fp05 - Fm05);
u(1) = u(3);
u(2) = u(3);
u(end-1) = u(end-2);
u(end) = u(end-2);

t = t + dt;
end

% graphics
figure;
plot(x,u,'bo');
hold on
plot(x+a*t,u0(x),'k-');
xlim([x(1) x(end)]);
ylim([0 2]);


with the following output

Note that the analytical solution $$u(x,t) = u(x- at,0)$$ of $$u_t = -au_x$$ translates the initial data with speed $$a$$. For the prescribed rectangular initial data, we have $$u(x,t) = \Bbb I_{at-1/2\leq x\leq at+1/2}$$ where $$\Bbb I$$ denotes the indicator function.

(1) R.J. LeVeque, Numerical Methods for Conservation Laws. Birkhäuser, 1992. doi:10.1007/978-3-0348-8629-1

(2) P.K. Sweby: "High resolution schemes using flux limiters for hyperbolic conservation laws", SIAM J. Numer. Anal. 21(5) (1984), 995-1011. doi:10.1137/0721062

• Hi, thanks for the answer. I'm trying it currently. Would you mind writing the $F_{i-1/2}^k$ as well, I don't want get the indexes wrong... – jan.sende Aug 13 at 10:08
• It's not as clear cut as you think... For example replacing $r_i$ with $r_{i-1}$ is not the right approach... Anyway, it works now! Thanks! If I interpret it right, my problem is mixing Lax-Wendroff with Osher. I should instead use Upwind for the low-order integrator... Do you by any chance have an article/link for the derivation of the arbitrary wave speed form? – jan.sende Aug 13 at 10:49
• @jan.sende Oh yes, it's the right way to do it -- just shift the index $i$ from $i$ to $i-1$ (see my Matlab implementation in the updated post). – Harry49 Aug 13 at 11:59
• I think Mathematica is trolling me. It didn't work before with $\phi\left(r_{i-1}\right)$, but now it does! But thanks for the effort of putting it into Matlab. I will check out the Book at the library in the evening. – jan.sende Aug 13 at 12:26