3. Finite differences

Like with time, the spatial variability of the solution also needs to be represented discretely. The purpose of the spatial discretization is two fold: it defines the data structure for the variables that are being time integrated and it enables to compute the right-hand-side for the time-integrator problem,

$$\frac{\partial \phi}{\partial t} = \mathcal{L}(\phi).$$

This entails to evaluate the various spatial derivatives that appear in the Navier-Stokes equations. We choose to adopt the so-called finite difference representation of our fields: Any smoothly varying scalar field (e.g. $\phi$) is then stored as a collection of its point values.

$$\phi_n = \{\phi(t_n, \mathbf{x}_1), \phi(t_n, \mathbf{x}_2), ..., \phi(t_n, \mathbf{x}_N)\}.$$

We will only consider a two-dimensional domain, ($\mathbf{x} = \{x, y\}$) and the points are distributed on a regular and square Cartesian grid. In principle, each field (i.e. $u_x, u_y$ and $p$) can be stored for its own set of point locations that differs from the others. However, it is most easy (and minimalist) to have a so-called co-located grid, where the point locations for all variables coincide.

Next, the concept of indexing is introduced, A domain of size $L_0 \times L_0$ can be discretized with $N \times N$ points, giving a grid spacing of $\Delta = L_0 N^{-1}$. Each entry corresponds to a location in physical space, but also in the computer memory. We can use integer Cartesian indexing $i, j$ to refer to some location in space,

$$\phi_{i,j} = \phi\left(x = X_0 + \Delta \left(i + \frac{1}{2}\right), y = Y_0 + \Delta \left(j + \frac{1}{2}\right)\right)$$

where $\{X_0,Y_0\}$ is the origin of the axis.

We will start coding our solver with declaring the solution fields in some file name solver.c. As the various solver components are added, this file will become a coding mess, so beware that it will need to be restructured later.

#include <stdio.h>

#define N 100

double ux[N][N], uy[N][N], p[N][N]

int main() {

}

Here we have declared the fields ux, uy and p on a N x N grid. The value of the macro N is chosen small enough for quick testing, and large enough to do simple tests. The field values can be set and read in a loop over the grid points.

...
int main() {
    for (int i = 0; i < N; i++) 
        for (int j = 0; j < N; j++) {
            ux[i][j] = 0; // example
            printf ("%g\n", p[i][j]); 
        }
}

The double loop (over i, and j) will be used many times in our solver. Furthermore, the indexing with i and j may become tedious. As such we can provide ourselves with some macros that makes coding easier and less prone to errors.

#include <stdio.h> 
#define foreach() for (int _i = 0; _i < N; _i++) \ 
                      for (int _j = 0; _j < N; _j++) 
#define val(s, x_ind, y_ind) s[_i + x_ind][_j + y_ind]

The purpose of these underscore indexing (e.g. _i) is to reduce the chance that the this variable conflicts with future variables. Using these definitions, the previous sample code becomes more clear,

...
int main() {
    foreach() {
        val(ux, 0, 0) = 0;
        printf ("%g\n", val(p, 0, 0));
    }
}   

Especially when we need to access neighboring field values, which is important to quantify the variability of the solution at some point. However, at the edges of our domain neighbor points do not always exists. The most minimalist way of dealing with this is by choosing periodic boundary conditions for our solver. This can be achieved by wrapping the indices if they are smaller than 0 (e.g. $-1 \rightarrow N - 1$ or larger than N - 1 (e.g. $N \rightarrow 0$). Then, the grid point at the edges “see” the data at the opposite boundary.

...
#define WRAP(i) (i < 0 ? N + i : (i >= N ? i - N : i))  
#define val(s, x_ind, y_ind) s[WRAP(_i + x_ind)][WRAP(_j + y_ind)]
...

Where C’s ternary operator is used to define the WRAP() macro. We can test if it works,

...
foreach()
    val(p, 0, 0) = 1.;
    
foreach() 
    if (val(p, 1, 1) != 1 || p( -1, -1) != 1) {
        puts ("Error");
        return 1;
    }
puts ("Succes!");
...

In the second loop, we check if the top-right or bottom-left neighbor value is indeed the value we had set it to be in the first loop. If it is not the case, it prints an error message and quits the program (using return in the main() function).

$ gcc solver.c
$ ./a.out
Succes!

A test for finite differencing

It is important to test our code a bit more quantitatively in the context of finite differences. Consider a dummy scalar field on a domain of size $[-5, 5]\times [-5, 5]$,

$$\phi(x,y) = e^{-x^2 - y^2}.$$

We can analytically differentiate this function,

$$\frac{\partial \phi}{\partial x} = -2x\ e^{-x^2 - y^2},$$

$$\frac{\partial \phi}{\partial y} = -2y\ e^{-x^2 - y^2}.$$

the test considers estimating these derivatives numerically using our code. From high school mathematics we may remember that the spatial derivative could be defined like so:

$$\frac{\partial \phi}{\partial x} = \mathrm{lim}_{h \rightarrow 0} \frac{\phi(x + h) - \phi(x)}{h}.$$

However, for this moment, there is no reason to bias in any specific direction (i.e using $x+h$, not $x - h$). We could also write a more accurate centered approximation,

$$\frac{\partial \phi}{\partial x} = \mathrm{lim}_{h \rightarrow 0} \frac{\phi(x + h) - \phi(x - h)}{2h}.$$

Instead of taking the limit of $h$ to 0, our solver will have to do with $h = \Delta$. Further, to initialize the fields, we should also define some more useful macros and variables.

    ...
    #include <math.h>
    
    #ifndef N
    #define N 100
    #endif
    ...
    #define Delta (L0/N)
    #define x (X0 + Delta*(_i + 0.5))
    #define y (Y0 + Delta*(_j + 0.5))
    #define sq(x) ((x)*(x))
    double L0 = 10;
    double X0 = -5, Y0 = -5;
    ...
    int main() {
        double phi[N][N]; // A dummy field
        foreach()
            val (phi, 0, 0) = exp(-sq(x) - sq(y));
                
        doubble L2 = 0;
        foreach() {
            double dphi_dx = (val(phi, 1, 0) - val (phi, -1, 0))/(2*Delta);
            double dphi_dy = (val(phi, 0, 1) - val (phi, 0, -1))/(2*Delta);
            double errorx = -2*x*exp(-sq(x) - sq(y)) - dphi_dx; 
            double errory = -2*y*exp(-sq(x) - sq(y)) - dphi_dy;
            L2 += sq(Delta)*sqrt(sq(errorx) + sq(errory));
        }
        printf ("%d %g\n", N, L2);
    }

Because we have chosen to set the number of cells in each dimension (N) as a macro (using the #define directive), it is not easily possible to make convergence-test loop as we did in the previous chapter. However, the additional #ifndef conditional statement allows to define N at compilation time. Say for N = 32 we can do

$ gcc -DN=32 solver.c -lm
$ ./a.out
32 0.334759

If we want to do a convergence study to check our coding, we should write a command-line bash script. Open a new file called check.sh,

#!/bin/bash
FILE=L2_error
rm $FILE
for ((R = 16 ; R <= 512 ; R *= 2))
do
    gcc -DN=$R solver.c -lm
    ./a.out >> $FILE
done

Then in the terminal, make it an executable and run this script,

$ chmod +x check.sh
$ ./check.sh
rm: cannot remove 'L2_error': No such file or directory

The rm error message can be ignored: the rm command tried to remove a file it could not find, as it was not produced yet. The bash script uses the >> operator, which redirects the output from the a.out program that would normally be seen in the terminal to a plain text file. You may check the contents of the newly created file called L2_error, which should contain the convergence-study data. Go ahead and plot it on a log–log scale, and try to find a fit. It may look something like:

If all went well, you should see that the so-called 3-point central finite difference is second-order accurate.

The code in solver.c is getting messy…

The solver.c code is burdened with a lot of moderately useful #define statements that effectively defeat the purpose of the somewhat cleaner code that follows. We can offload the code that is not directly related to the main function of the program to elsewhere. I therefore propose to create a file called common.h where we keep commonly used utilities. I fill it with:

#include <stdio.h>
#include <math.h>

#ifndef N
#define N 100
#endif

#define Delta (L0/N) 
#define x (X0 + Delta*(_i + 0.5))
#define y (Y0 + Delta*(_j + 0.5))
#define sq(x) ((x)*(x))

#define foreach() for (int _i = 0; _i < N; _i++) \
                      for (int _j = 0; _j < N; _j++) 
#define WRAP(i) (i < 0 ? N + i : (i >= N ? i - N : i))  
#define val(s, x_ind, y_ind) s[WRAP(_i + x_ind)][WRAP(_j + y_ind)]

//Default domain size
double L0 = 1.;
double X0 = 0, Y0 = 0;

The solver.c code can become more accessible and focussed now:

#include "common.h"

int main() {
  X0 = Y0 = -5.;
  L0 = 10.;
  double phi[N][N];
  foreach()
    val (phi, 0, 0) = exp(-sq(x) - sq(y));
  double L2 = 0;
  foreach() {
    double dphi_dx = (val(phi, 1, 0) - val(phi, -1, 0))/(2*Delta);
    double dphi_dy = (val(phi, 0, 1) - val(phi, 0, -1))/(2*Delta);
    double errorx = -2*x*exp(-sq(x) - sq(y)) - dphi_dx; 
    double errory = -2*y*exp(-sq(x) - sq(y)) - dphi_dy;
    L2 += sq(Delta)*sqrt(sq(errorx) + sq(errory));
  }
  printf ("%d %g\n", N, L2);
}

Mind that the common.h file in placed between quotation marks ("...") instead of the <...> markers that were used for the system header files (stdio.h and math.h). This hints the compiler where to look for these files. This alteration should not change any of the functionality of the code, and the previous script should still work, as long as solver.c is in the same folder as common.h. We are effectively creating a rudimentary interface for coding with our flow solver. It is not really user friendly (yet), because it has many quirks. For example, try to declare a variable named x. One could say this is an old-fashioned style of “academic” programming, where the programmer is typically the only user of the code.

Continue to chapter 4, where we start approximating the $\mathcal{L}(\mathbf{u})$ vector field.