4. The advection term

Tracer advection

In this chapter we treat the computation of the so-called advection term. It is the under-lined section of the equations,

$$\frac{\partial \mathbf{u}}{\partial t} = \underline{-(\mathbf{u} \cdot \nabla)\mathbf{u}} - \nabla p + \nu \nabla^2 \mathbf{u}.$$

You may also hear this term to be called “the non-linear term” or the “convective term”. In essence, this term describes the transport of the velocity field by itself. It will prove useful and instructive to consider a helper equation: The advection equation for a scalar field $s$,

$$\frac{\partial s}{\partial t} = -(\mathbf{u} \cdot \nabla) s,$$

describes the evolution of a tracer field ($s$) advected by a velocity field ($\mathbf{u}$). Clearly, if we could substitute $s \leftarrow \mathbf{u}$, we retrieve our desired non-linear term. For now, we focus on implementing a function called tracer_advection, which computes the tendency for any tracer (s) based on a velocity field.

Computing the advection term

We open a file called ns.h, indicating that this is our Navier-Stokes-equations problem-solver file, and start by writing,

#include "common.h"

// A function for the advection term of a field `s`
void tracer_advection (double s[N][N], double ux[N][N], 
                       double uy[N][N], double ds_dt[N][N]) {
    ...
}

Here we have declared a C-language function. In this case, it does not follow a standard input-to-output structure. Rather, the input is implied to be the fields s, ux, uy and ds_dt. The values of the latter will be computed and overwritten by this function (output) whereas it will use, but leave the former three fields unchanged (input). The dimension of these field are hard-coded to be 2 with the double [N][N]-array syntax, conforming with our earlier design in common.h (from the previous chapter).

but what will be the body of this function? It is good to write out the components of the previous equation as the C-compiler does not understand vector notation.

$$\frac{\partial s}{\partial t} = -(u_x \frac{\partial s}{\partial x} + u_y\frac{\partial s}{\partial y})$$

Although the velocity component values are readily available from the input to this function, the derivatives of $s$ need to be approximated,

...
void tracer_advection (double s[N][N], double ux[N][N], 
                       double uy[N][N], double ds_dt[N][N]) {
        foreach() {
            double ds_dx = ...
            double ds_dy = ...
            val(ds_dt, 0, 0) = -(val(ux, 0, 0)*ds_dx + val(uy,0,0)*ds_dy);
        }
}

It turns out that the standard 3-point 2-nd order accurate stencil for the derivative, as tested in the previous chapter, i.e.;

...
  // Not to use
  double dsdx = (val(s, -1, 0) - val(s, 1, 0))/(2*Delta);
...

is not suitable for explicit time integration of the advective tendency. This is due to an unfortunate accumulation of the errors introduced by this approximation over time. Instead, we will adopt a so-called upwinding strategy, where the stencil is biased based on the local flow direction. More specifically, we use a 2-point stencil where the gradient is estimated with upstream values of the field $s$. Conceptually, this is the gradient that is being advected towards the current cell.

...
void tracer_advection (double s[N][N], double ux[N][N], 
                       double uy[N][N], double ds_dt[N][N]) {
        foreach() {
            double ds_dx = val(ux, 0, 0) > 0 ?
                          (val(s, 0, 0) - val(s, -1, 0))/Delta :
                          (val(s, 1, 0) - val(s, 0, 0))/Delta;
            double ds_dy = val(uy, 0, 0) > 0 ?
                          (val(s, 0, 0) - val(s, 0, -1))/Delta :
                          (val(s, 0, 1) - val(s, 0, 0))/Delta;
            val(ds_dt, 0, 0) = -(val(ux, 0, 0)*ds_dx + val(uy,0,0)*ds_dy);
        }
}

It is not likely that adding so many lines of code goes without typos or other errors. We should at least test the syntax and crate a file, (advection_test.c) which #includes our new code,

#include "ns.h"

int main() {
    ;
}

Now we can test if gcc finds any issues,

$ gcc advection_test.c -lm
$ ./a.out

No errors? Then, are we done now with this term? No! We should test our tracer_advection function qualitatively and quantitatively.

Before we do that, its time for an antr’acte…