11. Some improved functionalities

In this chapter we will improve our code to increase its function, and improve the user interface. These points are inspired by my wishes when composing the Gallery chapter. We consider four main topics,

1. An option to solve for a flow tracer field
2. An option for an additional body force to the Navier-Stokes equations
3. A user-interface for the time loop

(1 and 2) A tracer field and a buoyancy field

Since we already have the solving infrastructure, it would be nice to solve for the evolution of a tracer field as well. It can be implemented by adding a global field in ns.h called tracer and then update it during the advance_ns() step. Furthermore, it will prove useful to add an additional body-force momentum-source term to the vertical velocity (e.g. a gravity acceleration force). These two additions require very little new code to be added in ns.h,

...
double ux[N][N], uy[N][N], p[N][N];
double tracer[N][N], acceleration_y[N][N];
...

void advance_ns (double dt) {
...
  advance_scalar (uy_star, nu, dt);
  advance_scalar (tracer, kappa, dt);
  foreach()
    val(uy_star, 0, 0) += dt*val(acceleration_y, 0, 0);
  // Project and update
  ...

That it! See in the Gallery for its usage. But you may not recognize the clean case setup files before reading about…

(3) The User interface

I wish to implement a function called run() which will run the time loop (i.e. for (t = 0; t < t_end; t += dt)) for our solver. In its not user-friendly form we would add to ns.h,

...

double t_end = 1;

void run() {
  for (t = 0; t < t_end; t += dt) {
    // advance
    dt = dt_next(t_event);
    advance_ns(dt);
}

were we could change t_end in our .c-case files to our desire. However, it would be nice to *optionally* break into this loop and add functions as we please. To illustrate, one may want to execute some function every iteration, but the function body is only defined in a later.c`-case file. A method to do this is to use function pointers.

double t_end = 1;
// Function pointer
void (*every_iter)() = NULL;

void run() {
  for (t = 0; t < t_end; t += dt) {
    if (every_iter != NULL) // Do not excecute an undefined function
      every_iter();
    // advance
    dt = dt_next(t_event);
    advance_ns(dt);
}

Using this, one may optionally assign a function to the void every_iter() function pointer in a .c file. It could look like,

...
void my_fun() {
  printf ("%g \n", t);
}

int main () {
  ...
  every_iter = my_fun;
  run();
}

This would result in the execution of the void my_fun() function every time step. Printing the value of the time parameter for every iteration. Similar, I would like functions that execute every iter_interval iterations (i.e. solver time steps) and every t_interval of time units. The code becomes

...
// Time loop variables
double t_interval = 1e8, t_end = 1;
int iter, iter_interval = 1e8;

// Function pointers
void (*every_t_interval)() = NULL;
void (*every_iter_interval)() = NULL;
void (*every_iter)() = NULL;

// A time loop 
void run() {
  //(re)set
  iter = 0;
  t = 0;
  // time of next t_interval
  double t_event = 0;
  for (t = 0; t < t_end; t += dt) {
    // Call events if defined
    if (every_iter != NULL)
      every_iter();
    if ((iter % iter_interval) == 0 && every_iter_interval != NULL)
      every_iter_interval();
    if (fabs(t - t_event) < 1e-8) { // Allow small binary-representation error
      if (every_t_interval != NULL)
        every_t_interval();
      // update t_interval 
      t_event = min(t_end, t_event + t_interval);
    }
    // advance
    dt = dt_next(t_event);
    advance_ns(dt);
    iter++;
  }
}

Finally,

I have added a noise() macro to common.h which computes a random value between -1 and 1.

#include <stdlib.h>
#define noise() ((double)((rand()%2000) - 1000)/1000.)

We will use this for the cases in the gallery