How to properly implement a hyperbolic TVD flux limiter

hyperbolic-equationsnumerical methodstransport-equation

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:

We start with the transport equation:

$$
\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$:
my result

Not what I expected…

EDIT:

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

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

enter image description here

Best Answer

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

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

Related Question