(advection)=
# Advection Equation

```{note} 

**The learning objectives of this section are:**
- describe the role of advection and identify the advection equation
- explain the use of the upwind method

```

In this section we discuss the numerical solution of the 1D advection equation as given by

$$
  \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial x} = 0
$$ (advection1)

This first order hyperbolic PDE describes
the propagation of a disturbance or pulse $c(x,t)$ with a constant speed $u$ through a fluid in an infinite long 1D domain.
Examples are a single wave ($c$ being the wave amplitude), or a pollutant ($c$ as concentration).

Given any initial condition $c(x,0) = c^0(x)$ the solution at any *later* time is given
by $c(x,t) = c^0(x-ut)$. This can be verified by computing the $t-$ and $x-$derivatives of the solution and substituting
into the equation above. The function $c^0(x-ut)$ is a solution for
**right-travelling waves**, since the graph of $c^0(x-ut)$
is simply the graph of $c^0(x)$ shifted to the right by $ut$ spatial units. As time
increases, the profile $c^0(x)$ moves to the right at speed $u$.
Thus solutions are propagated, or **advected**, with constant
speed $u$ without *deforming*, that is, no change in shape. If $u > 0$, propagation is from left to right. If $u < 0$, propagation is from right to left, that is, **left-travelling waves**.

Advection is an example of transport by fluid motion. Another example is **convection** which is the transport of a fluid (water or air) in response to heat, for instance, boiling water.
Although, these two terms are interchangeable in mathematical sense, <u>advection</u> is the usual term.

The above described problem is an initial value problem (IVP). Now, let us consider a 1D domain with a finite length $L$, that is, $0 \leq x \leq L$.
In this case we need to impose some boundary conditions. A boundary condition is needed at a certain boundary location where the information *enters* the domain.
This depends on the propagation velocity. If $u > 0$, the ingoing transport is at the left boundary, $x=0$, while if $u < 0$, the ingoing transport is at the right boundary, $x=L$.

As an example, we consider Eq. {eq}`advection1` with $u = 1$ m/s. A well posed problem is then formulated as follows:

$$
  \begin{align}
    &\frac{\partial c}{\partial t} + u \frac{\partial c}{\partial x} = 0\, , \quad x \in (0,L)\, ,\quad t > 0 \\
    \\
    &c(x,0) = c^0(x)\, ,\quad x \in [0,L] \\
    \\
    &c(0,t) = f(t)\, ,\quad t > 0 \\
  \end{align}
$$ (advection2)

where $f(t)$ described the incoming wave-like disturbance; this function is usually expressed as a Fourier time series.
This IBVP is called the **wave equation** and is commonly applied in wave modelling but also in sediment transport.

```{note}

Note that the above IBVP {eq}`advection2` would become ill posed if $u<0$.

```

In contrast to the diffusion equation {eq}`diffusion1`, the solution to the advection equation is not necessarily smooth. In other words, $c(x,t)$ is allowed to be discontinuities or
non-smooth. Examples are the propagation of shock waves, hydraulic jumps and tidal bores.
Although the wave equation {eq}`advection1` appears deceptively simple, the numerical solution to this equation is quite a challenge, as we will discover below.

## Discretization

We reconsider the computational grid as discussed in Section {ref}`diffusion`.
Thus, we have our uniform mesh consisting of nodes $x_0,\cdots,x_{M}$ and $\Delta x = x_m - x_{m-1}$.
We apply the method of lines, so that we first discretize in space and then integrate over time.

A possible candidate for the spatial approximation would be central differences
(cf. Eq. {eq}`cdf1`), as follows

$$
  \frac{dc_m(t)}{dt} + u \frac{c_{m+1}(t) - c_{m-1}(t)}{2\Delta x} = 0\, , \quad m=1,\ldots,M-1\, , \quad t > 0
$$

This is a semi discretization of Eq. {eq}`advection1` for the *inner* grid points of the computational domain.

:::{card} Exercise

Verify yourself that this spatial discretization is second order accurate.

:::

The resulted system of first order ODEs with the unknowns $c_m(t)$ can be written as

$$
  \frac{d\mathbf{c}}{dt} = \mathbf{A} \mathbf{c}
$$ (vecadv)

with $\mathbf{c} = (c_1,\cdots,c_{M-1})^{\text{T}}$ being the vector representing the unknowns and
matrix $\mathbf{A}$ is the $(M-1) \times (M-1)$ discretization matrix as given by

$$
   \mathbf{A}=   \frac{u}{2\Delta x} \left (
   \begin{array}{rrrrrrr}    0       & -1 &  0      &          & \cdots &        &   0       \\
                             1       &  0 & -1      &          &        &        &           \\
                             0       &  1 &  0      & -1       &        &        &           \\
                             \vdots  &    & \ddots  & \ddots   & \ddots &        &  \vdots   \\
                                     &    &         &  1       &  0     & -1     &   0       \\
                                     &    &         &          &  1     &  0     &  -1       \\
                             0       &    &         & \cdots   &        &  1     &   0
   \end{array} \right )
$$

```{admonition} Skew symmetry
:class: dropdown
This matrix is a **skew symmetric** matrix ($\mathbf{A} = -\mathbf{A}^T$) implying that the nonzero eigenvalues are all purely imaginary.
Mathematically, this system is regarded as *undamped* while its solution consists of harmonic (oscillatory) components.
```

The next step is to integrate the above system of ODEs with respect to time. Although the method of lines is flexible in the sense that
the choice of time integration is independent from the choice of space discretization, we must be careful since the resulting scheme
may not be stable. For example, let us choose the forward Euler scheme, we get

$$
  \frac{c^{n+1}_m - c^{n}_m}{\Delta t} + u \frac{c^n_{m+1} - c^n_{m-1}}{2\Delta x} = 0\, , \quad m=1,\ldots,M-1\, , \quad n = 0, 1, 2, \ldots
$$ (ftcs3)
 
This explicit scheme is again the FTCS scheme but applied to the wave equation. In contrast to the diffusion equation (cf. Eq. {eq}`ftcs`),
this FTCS scheme {eq}`ftcs3` is **unconditionally unstable**, that is, it is unstable for any time step size. Therefore, the FTCS scheme is useless
for any advection-type equation.

This instability is related to the fact that the influence of the disturbance only comes from upstream and not from downstream.
This is a specific feature for hyperbolic problems. (By contrast, parabolic and elliptic problems do not have a preference of direction.)
More specifically, assuming that the propagation velocity
is positive, $u>0$, the solution $c_m$ at $x_m$ should only be affected by $c_{m-1}$ at $x_{m-1}$ (upstream) and not by the solution downstream, that is, $c_{m+1}$.
However, since this does occur with scheme {eq}`ftcs3`, it inevitably leads to unlimited growth.

```{note}

A possible solution to this instability problem is to apply the backward Euler scheme. This **implicit** scheme is always stable (see Chapter {ref}`numerical_modelling`). We then obtain
the so-called BTCS (**B**ackward in **T**ime, **C**entral in **S**pace) scheme, as follows

$$
  \frac{c^{n+1}_m - c^{n}_m}{\Delta t} + u \frac{c^{n+1}_{m+1} - c^{n+1}_{m-1}}{2\Delta x} = 0
$$

However, this is very unusual since the numerical solution is less straightforward to obtain. This can be seen as follows.

Let us consider the system of ODEs {eq}`vecadv` and subsequently apply backward Euler. This yields

$$
  \frac{\mathbf{c}^{n+1}-\mathbf{c}^n}{\Delta t} = \mathbf{A} \, \mathbf{c}^{n+1} \quad \Rightarrow \quad \left ( \mathbf{I} - \Delta t \mathbf{A} \right ) \, \mathbf{c}^{n+1} = \mathbf{c}^n
$$

with $\mathbf{c}^n = (c^n_1,\cdots,c^n_{M-1})^{\text{T}}$ being the vector representing the numerical solution at time step $n$ and
$\mathbf{A}$ is the discretization matrix as given above. Matrix $\mathbf{I}$ is the identity matrix.
We see that the solution $\mathbf{c}^{n+1}$ cannot be computed immediately. Instead, we need to solve a linear system of $M-1$ equations for the $M-1$ unknowns at the new time step.
In addition, the rank of the matrix $\mathbf{A}$ is expressed in terms of $M = L/\Delta x$ which can be quite large, especially if the mesh size $\Delta x$ is chosen very small.

This observation generally holds for any implicit scheme.

```

## The upwind method

Since we are dealing with information being propagated along a certain direction, it would be better to propose a scheme that gathers information to be propagated
in the same direction. For instance, the backward finite difference approximation, Eq. {eq}`bdfs`, applied to $\partial c/\partial x$ will gather information from the
left to propagate, while the forward finite difference operator, Eq. {eq}`fdfs`, gathers from the right. Hence,
our knowledge of which direction the wave should propagate indicates the **upwind** direction.
If $u > 0$, then propagation of the disturbance is from left to right. So, the upwind direction is *to the left*, that is, the information is obtained from the **backward**
direction. For $u < 0$, the upwind direction is *to the right* and so, the information is from the **forward** direction.

The **upwind difference scheme** is given by

$$
  \frac{\partial c}{\partial x} =
  \begin{cases} \frac{c_{m}-c_{m-1}}{\Delta x}\, , \quad \text{if} \, \, u > 0 \\
                 \\
                 \frac{c_{m+1}-c_{m}}{\Delta x}\, , \quad \text{if} \, \, u < 0
  \end{cases}
$$

and it correctly passes information in the direction of the advective transport.
This approximation is called **first order upwinding**.

:::{card} Exercise

Verify the first order accuracy of this scheme to $\partial c/\partial x$.

:::

Application of forward Euler yields

$$
  c^{n+1}_m = c^n_m - \sigma
  \begin{cases} c^n_m - c^n_{m-1} \, , \quad \mbox{if} \, \, u > 0 \\
                 \\
                 c^n_{m+1} - c^n_m \, , \quad \mbox{if} \, \, u < 0
  \end{cases}
$$ (ftbs)

with

$$
  \sigma = \frac{u\,\Delta t}{\Delta x}
$$

defined as the **Courant number**. This number plays a very important role for hyperbolic problems.

Scheme {eq}`ftbs` is called the **FTBS** scheme, which stands for **F**orward in **T**ime, **B**ackward in **S**pace,
and is first order accurate in both time and space.

:::{card} Exercise

For the FTBS scheme with $u>0$, show that $\tau_{\Delta t,\Delta x} = \mathcal{O} (\Delta t,\Delta x)$.

```{admonition} Solution
:class: tip, dropdown

With $u>0$, the FTBS scheme is given by

$$
  \frac{c^{n+1}_m - c^n_m}{\Delta t} + u \frac{c^n_m - c^n_{m-1}}{\Delta x} = 0
$$

To check consistency we must require the solution $c(x,t)$ to be smooth, so that we can apply the following Taylor series expansions

$$
  c(x,t+\Delta t) = c(x,t) + \Delta t\, c_t(x,t) + \frac 12 \Delta t^2\, c_{tt}(x,t)
$$

and

$$
  c(x-\Delta x,t) = c(x,t) - \Delta x\, c_x(x,t) + \frac 12 \Delta x^2\, c_{xx}(x,t)
$$

Note that both series have been expanded till second order ("*first derivative + first order accuracy*").

Next, substituting both expansions into the FTBS scheme yields

$$
\begin{align*}
\tau_{\Delta t,\Delta x} &= \frac{c(x,t+\Delta t)-c(x,t)}{\Delta t} + u \, \frac{c(x,t)-c(x-\Delta x,t)}{\Delta x} \\
                         &\\
                         &=\frac{\partial c}{\partial t}(x,t)+\frac 12 \Delta t \frac{\partial^2 c}{\partial t^2}(x,t) + u \frac{\partial c}{\partial x}(x,t) -
                           \frac{u\Delta x}{2}\frac{\partial^2 c}{\partial x^2}(x,t) + \ldots\\
                         &\\
                         &=\frac 12 \Delta t \frac{\partial^2 c}{\partial t^2}(x,t) - \frac{u\Delta x}{2}\frac{\partial^2 c}{\partial x^2}(x,t) + \ldots
\end{align*}
$$

and, thus

$$
  \tau_{\Delta t,\Delta x} = \mathcal{O} \left(\Delta t, \Delta x \right)
$$

Hence, the FTBS scheme is first order in $\Delta t$ and first order in $\Delta x$.

```
:::

The above exercise demonstrates that the FTBS scheme is consistent with the wave equation {eq}`advection1`.
Next, we check the stability of the FTBS scheme. Since it is explicit, we expect that this scheme can either be
*conditionally stable* or *unconditionally unstable* (like the FTCS scheme).

Suppose that $u>0$. Then the FTBS scheme is given by

$$
  \frac{c^{n+1}_m - c^n_m}{\Delta t} + u \frac{c^n_m - c^n_{m-1}}{\Delta x} = 0
$$

To check stability, we rewrite the scheme as follows

$$
  c^{n+1}_m = c^n_m - \frac{u\,\Delta t}{\Delta x} \, \left ( c^n_m - c^n_{m-1} \right )
$$

or 

$$
  c^{n+1}_m = c^n_m - \sigma \, \left ( c^n_m - c^n_{m-1} \right )
$$

This equation is further rewritten as

$$
  c^{n+1}_m = (1-\sigma)c^n_m + \sigma c^n_{m-1}
$$

We assume by induction that $c^n_m \geq 0$ and $c^n_{m-1} \geq 0$. Then $c^{n+1}_m \geq 0$ if $1-\sigma \geq 0$ and $\sigma \geq 0$.
Hence, the stability condition of the FTBS scheme reads

$$
  0 \leq \sigma \leq 1
$$

This condition is known as the *Courant-Friedrichs-Lewy* (**CFL**) condition.

The physical interpretation of this CFL condition is that the distance covered during $\Delta t$ with finite speed $u$ must be smaller than $\Delta x$, that is,
$u\,\Delta t \leq \Delta x$. This implies that the parameters $\Delta t$ and $\Delta x$ are mutually coupled: a small $\Delta x$ implies a small $\Delta t$.
The CFL condition plays an important role in the numerical solution to wave-like (or advection-like) PDEs.

:::{card} Exercise

Suppose that $u<0$. Verify that in that case the FTBS scheme is stable if

$$
 -1 \leq \frac{u\,\Delta t}{\Delta x} \leq 0
$$

:::

## The leapfrog scheme

In the previous section we have discussed a conditionally stable scheme to solve the wave equation, namely, the FTBS scheme. This explicit scheme is constructed based on the forward Euler
scheme and the upwind difference scheme. A main disadvantage of this scheme is that is only first order accurate both in time and space, which is rather inaccurate.

A widely used scheme that is second order accurate but also stable, albeit conditionally, designed for *oscillatory motions* is the **leapfrog** scheme.
This explicit scheme is devised based on central differences both in time and space and is given by

$$
  \frac{c^{n+1}_m - c^{n-1}_m}{2 \Delta t} + u \frac{c^n_{m+1} - c^n_{m-1}}{2\Delta x} = 0\, , \quad m=1,\ldots,M-1\, , \quad n = 0, 1, 2, \ldots
$$ (leapfrog)

:::{card} Exercise

Show that the leapfrog scheme {eq}`leapfrog` is second order accurate both in time and space.

```{admonition} Solution
:class: tip, dropdown

First, we replace the numerical solution $c^n_m$ by the exact one, $c(x,t)$ at point $x_m$ and time $t_n$, thus

$$
  \tau_{\Delta t,\Delta x} = \frac{c(x,t+\Delta t) - c(x,t-\Delta t)}{2 \Delta t} + u \frac{c(x+\Delta x,t) - c(x-\Delta x,t)}{2\Delta x}
$$ (leapfrog2)

Furthermore, we require the solution $c(x,t)$ to be smooth. The Taylor series for $c(x,t\pm\Delta t)$ will be expanded till the third order
since we expect second order accuracy for the approximation of the first order derivative in time (2+1=3). This also holds for the expansion
of $c(x\pm\Delta x,t)$ in space. Hence, we have

$$
  c(x,t\pm\Delta t) = c(x,t) \pm \Delta t\, c_t(x,t) + \frac 12 \Delta t^2\, c_{tt}(x,t) \pm \frac 16 \Delta t^3\, c_{ttt}(x,t)
$$

and

$$
  c(x\pm\Delta x,t) = c(x,t) \pm \Delta x\, c_x(x,t) + \frac 12 \Delta x^2\, c_{xx}(x,t) \pm \frac 16 \Delta x^3\, c_{xxx}(x,t)
$$

Next, substituting both expansions into Eq. {eq}`leapfrog2` yields

$$
\begin{align*}
\tau_{\Delta t,\Delta x} &= \frac{c(x,t+\Delta t)-c(x,t-\Delta t)}{2\Delta t} + u \, \frac{c(x+\Delta x,t)-c(x-\Delta x,t)}{2\Delta x} \\
                         &\\
                         &=\frac{\partial c}{\partial t}(x,t)+\frac 16 \Delta t^2 \frac{\partial^3 c}{\partial t^3}(x,t) + u \frac{\partial c}{\partial x}(x,t) +
                           \frac{u\Delta x^2}{6}\frac{\partial^3 c}{\partial x^3}(x,t) + \mathcal{O} \left(\Delta t^5, \Delta x^5 \right)\\
                         &\\
                         &=\frac 16 \Delta t^2 \frac{\partial^3 c}{\partial t^3}(x,t) + \frac{u}{6}\Delta x^2\frac{\partial^3 c}{\partial x^3}(x,t) +
                          \mathcal{O} \left(\Delta t^5, \Delta x^5 \right)
\end{align*}
$$

and, thus

$$
  \tau_{\Delta t,\Delta x} = \mathcal{O} \left(\Delta t^2, \Delta x^2 \right)
$$

Hence, the leapfrog scheme is second order in $\Delta t$ and second order in $\Delta x$.

```
:::

The leapfrog scheme is an example of a <u>multi-step method</u>, namely, a two-step scheme. This has a disadvantage that it requires *two* starting values while the wave equation
is provided with *one* initial condition. Hence, the leapfrog scheme is not self-starting. A solution to this problem is to apply a one-step scheme (e.g., FTCS, FTBS) at the first
time step only of the simulation.

Finally, the leapfrog scheme is stable (here shown without proof) if

$$
  |\sigma| = \frac{|u|\,\Delta t}{\Delta x} \leq 1
$$