Back to basics with PDEPE

Recently I got introduced to a PDE solver in MATLAB and thought it was a wonderful tool. The function in question is pdepe and on MATLAB’s webpage it claims to be a one-dimensional elliptic-parabolic PDE solver. By wrapping our governing equation, initial conditions and boundary conditions into functions with predefined function signatures, we were able to solve nonlinear system of PDEs (!) quickly. The total amount of code we wrote was less than 100 lines.

Using PDEPE
Despite the parabolic-elliptic caveat given on the MATLAB website, one may be tempted to try solving a hyperbolic equation. Trying it out is easy. Let’s solve a linear advection equation for a cosine wave:

    \begin{align*} \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} &= 0 \quad , \quad 0 \leq x \leq L \\ u(x, 0) &= \begin{cases} \cos(x) \quad , \quad 0 \leq x \leq 2\pi \\ 1 \quad , \quad otherwise \end{cases} ,\\ u(0, t) &= 1 \quad , \quad u(L, t) = 1 . \end{align*}

One first casts the governing equation into the form given by MATLAB:

    \begin{align*} c \frac{\partial u}{\partial t} = x^{-m}\frac{\partial}{\partial x}\left(x^m f\right) + s . \end{align*}

Where c, f and s are functions of \left(x, t, u, \partial u/ \partial x\right). We then write a function to compute c, f and s:

function [c, f, s] = advectionEquation(x, t, u, dudx)
    c = 1;
    f = -u;
    s = 0;
end

For the boundary conditions, we need to identify p and q such that at the boundaries,

    \begin{align*} p + qf = 0 . \end{align*}

It’s a bit annoying that MATLAB uses the flux f instead of \partial u/\partial x. If we had a Neumann boundary condition (i.e. \partial u /\partial x = 0 at L), we would have a bit of a trouble implementing it for our case. Writing p and q on the boundaries into the required function:

function [pl, ql, pr, qr] = boundaryConditions(xl, ul, xr, ur, t)
    pl = ul - 1;
    ql = 0;
    pr = ur - 1; 
    qr = 0;
end

The initial conditions is straightforward:

function u0 = initialCondition(x)
    if x <= 2*pi
        u0 = cos(x);
    else
        u0 = 1;
    end
end

Putting everything into a script:

% Parameters
m = 0;
x = linspace(0, 7*pi, 500);
t = linspace(0, 12, 5);

% Solve
sol = pdepe(...
    m,...
    @advectionEquation,....
    @initialCondition,...
    @boundaryConditions,...
    x,...
    t);
u = sol(:,:,1);

% Plot results
figure
for iT = 1:length(t)
    plot(x, sol(iT, :, 1), 'linewidth', 2)
    hold on
end
grid on
title('Advection of cosine wave')
xlabel('x')
ylabel('u')

Perfect! At first glance, one sees the wave moving smoothly to the right with no problems. However, look closer at the trough of the wave. I’ll help you with some markers:

We see that the shape is not preserved, the peak is getting smaller and smaller. this is likely due to the numerical viscosity diffusing the solution. However, the results are still encouraging.

Dealing with discontinuities
With the favorable results, let’s move on and try a more challenging case: Advection of a step function.

    \begin{align*} \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} &= 0 \quad , \quad -L \leq x \leq L \\ u(x, 0) &= \begin{cases} 1 \quad , \ x \leq 0 \\ 0 \quad , \quad x > 0 \end{cases} ,\\ u(-L, t) &= 1 \quad , \quad u(L, t) = 0 . \end{align*}

The only difference here is the initial condition and boundary conditions. The required MATLAB functions are similar to the one given earlier so I will leave them out.

Argh! Our friend spurious oscillation rears its ugly head! The theoretical discontinuity locations are marked with black markers and we can see that the drop occurs at the correct location, but there are some wild oscillation occurring before that. You can also see that the discontinuity gets a little smeared out. The quickest and dirtiest way to remedy such a situation would be to add some damping (artificial viscosity). Our equation now becomes:

    \begin{align*} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} + b \frac{\partial^2 u}{\partial x^2} = 0. \end{align*}

We need to tweak the function in MATLAB:

function [c, f, s] = advectionEquation(x, t, u, dudx)
    a = 1;
    b = -0.003; % Found through trial and error
    c = 1;
    f = -b*dudx;
    s = -a*dudx;
end

Things are starting to get a little fudgy. One has to tweak the coefficient b around a little. We want it to be big enough to remove the oscillations but small enough not to affect the results too great. Also, note that b should be <0 otherwise we will be adding “negative diffusion” (its “negative” because diffusion terms are typically written on the RHS) which destabilizes the solution. I eventually settled on a value of -0.003 such that it is small enough to have just a little overshoot at the oscillation.

Great! Now the solution looks really good. However, we have to remember what we did to arrive at this solution: the addition of a nonphysical diffusion term. This term is highly dependent on the solution and mesh size. How much to add is really a dark art and left to the coder’s discretion. Also, in this case it is relatively easy to determine how much diffusion to add because it is a pure advection equation – there should be as little smearing of the solution as possible and we already know what the solution looks like. If we had an advection-diffusion equation, and especially if the terms are nonlinear, it becomes increasingly difficult to gauge how much artificial diffusion to add.

These oscillations were a big problem in the computational fluid dynamics community in the 70s. The mathematical character of compressible flows are highly hyperbolic, discontinuities can form in the domain even from smooth initial and boundary conditions. Shockwaves are a prime example. Highly sophisticated theories and mathematical treatments gave rise to the total variation diminishing (TVD) schemes which make use of flux limiters. Despite the complicated sounding name, a flux limiter is just a way of automatically switching the discretisation from second order to first order in the solution close to discontinuities. This helps to remove the wiggles close to discontinuities while keeping the solution accurate elsewhere.

Another highly popular CFD scheme is the JST scheme, named after the people who created it: Jameson-Schmidt-Turkel. It’s basically what we did with the b coefficient, except that they devised a set of rules to compute how much diffusion to add. Artificial viscosity is added to the second- and fourth-derivatives. And of course the rigor that went into it is at a whole other level.

One thing I really want to find out is how much smearing our “fudged” solution creates. How would the second order advection equation with artifical viscosity compare to a pure advection equation solved with upwind discretization? To answer this, I wrote a backward Euler solver and will look at the results in the next post.

Leave a reply

Your email address will not be published.