# 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++) {
[i][j] = 0; // example
ux("%g\n", p[i][j]);
printf }
}
```

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(ux, 0, 0) = 0;
val("%g\n", val(p, 0, 0));
printf }
}
```

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(p, 0, 0) = 1.;
val
()
foreachif (val(p, 1, 1) != 1 || p( -1, -1) != 1) {
("Error");
puts return 1;
}
("Succes!");
puts ...
```

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(phi, 0, 0) = exp(-sq(x) - sq(y));
val
= 0;
doubble L2 () {
foreachdouble 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;
+= sq(Delta)*sqrt(sq(errorx) + sq(errory));
L2 }
("%d %g\n", N, L2);
printf }
```

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 $FILEfor ((R = 16 ; R <= 512 ; R *= 2))
do
-DN=$R solver.c -lm
gcc ./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() {
= Y0 = -5.;
X0 = 10.;
L0 double phi[N][N];
()
foreach(phi, 0, 0) = exp(-sq(x) - sq(y));
val double L2 = 0;
() {
foreachdouble 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;
+= sq(Delta)*sqrt(sq(errorx) + sq(errory));
L2 }
("%d %g\n", N, L2);
printf }
```

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.