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...


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


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;
        num = upp-up;
        den = up-u;
    for i=1:Nx
        if abs(den(i))>=eps
            r(i) = num(i)./den(i);
    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;

% graphics
hold on
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

  • $\begingroup$ 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... $\endgroup$ – jan.sende Aug 13 at 10:08
  • $\begingroup$ 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? $\endgroup$ – jan.sende Aug 13 at 10:49
  • $\begingroup$ @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). $\endgroup$ – Harry49 Aug 13 at 11:59
  • $\begingroup$ 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. $\endgroup$ – jan.sende Aug 13 at 12:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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