(Difference between revisions)
+                                                                                                                                                                            the function ''par->flux'' . This is not a fixed name for a function but a pointer to a function that defines the advection flux. We cannot click on it to determine how it is defined ! Doing a grep on flux reveals many lines. Doing a
-                                                                                                   % cd src                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             +
-                                                                                                   % make tags                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          +                                                                                                                                                                            % grep par\-\>flux *.c

+                                                                                                                                                                            flux = u_face T_face
+
+ and we know how the code implements the scheme ! == The end for now == == The end for now ==

## How the equations are solved

Our first two courses were about the general notion of object-oriented programming and the implementation of these notions in some typical Gerris classes. Now we turn the heart of the matter: how does Gerris solve the Navier-Stokes equations ?

### simulation_run()

Let us assume that you are familiar with Stephane's Journal of Computational Physics article and you want to see how it is implemented in code. On way to do it would be to see what is executed line by line using a tool such as ddd. In this course we lead you directly to the interesting part. We skip the bookkeeping and initialisation stuff that is executed first in gerris. The core of the simulation is in `simulation.c`, specifically in `simulation_run()`. Let us look at this function. We start with some more bookkeeping:

`static void simulation_run (GfsSimulation * sim){  GfsVariable * p, * pmac, * res = NULL;  GfsDomain * domain;  GSList * i;   domain = GFS_DOMAIN (sim);   p = gfs_variable_from_name (domain->variables, "P");  g_assert (p);  pmac = gfs_variable_from_name (domain->variables, "Pmac");  g_assert (pmac);   gfs_simulation_refine (sim);  gfs_simulation_init (sim);   i = domain->variables;  while (i) {    if (GFS_IS_VARIABLE_RESIDUAL (i->data))      res = i->data;    i = i->next;  }`

Then we start the various steps of the integration of a time step.:

`gfs_simulation_set_timestep (sim);  if (sim->time.i == 0) {    gfs_approximate_projection (domain,				&sim->approx_projection_params,				&sim->advection_params,				p, sim->physical_params.alpha, res);    advance_tracers (domain, sim->advection_params.dt/2.);  }  while (sim->time.t < sim->time.end &&	 sim->time.i < sim->time.iend) {    GfsVariable * g[FTT_DIMENSION];    gdouble tstart = gfs_clock_elapsed (domain->timer);     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);`

### Inside the time loop

Now we are fully inside the time loop.

`gfs_simulation_set_timestep (sim);     gfs_predicted_face_velocities (domain, FTT_DIMENSION, &sim->advection_params);     gfs_variables_swap (p, pmac);    gfs_mac_projection (domain,    			&sim->projection_params,     			&sim->advection_params,			p, sim->physical_params.alpha, g);    gfs_variables_swap (p, pmac);     gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_half_do, sim);     gfs_centered_velocity_advection_diffusion (domain,					       FTT_DIMENSION,					       &sim->advection_params,					       g,					       sim->physical_params.alpha);     if (gfs_has_source_coriolis (domain)) {      gfs_poisson_coefficients (domain, sim->physical_params.alpha);      gfs_correct_normal_velocities (domain, 2, p, g, 0.);      gfs_correct_centered_velocities (domain, 2, g, sim->advection_params.dt);      gfs_source_coriolis_implicit (domain, sim->advection_params.dt);      gfs_correct_normal_velocities (domain, 2, p, g, 0.);      gfs_correct_centered_velocities (domain, 2, g, -sim->advection_params.dt);    }     gfs_domain_cell_traverse (domain,			      FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,			      (FttCellTraverseFunc) gfs_cell_coarse_init, domain);    gfs_simulation_adapt (sim);     gfs_approximate_projection (domain,   				&sim->approx_projection_params,     				&sim->advection_params, p, sim->physical_params.alpha, res);     advance_tracers (domain, sim->advection_params.dt);     sim->time.t = sim->tnext;    sim->time.i++;     gts_range_add_value (&domain->timestep, gfs_clock_elapsed (domain->timer) - tstart);    gts_range_update (&domain->timestep);    gts_range_add_value (&domain->size, gfs_domain_size (domain, FTT_TRAVERSE_LEAFS, -1));    gts_range_update (&domain->size);  }  gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);    gts_container_foreach (GTS_CONTAINER (sim->events),			 (GtsFunc) gts_object_destroy, NULL);}`

From now on it becomes difficult to offer a traditional written description of the course. Instead, in a class the teacher should open a text editor, projecting it on screen in front of the students, and navigate through the code using the tip explained in the previous course. Correspondingly, the reader or listener to the course should open his text editor and navigate freely in the code, trying to find the implentation of the numerical methods. We guide the reader below for one particular scheme: the tracer advection-diffusion equation.

The tracer is updated though the function `advance_tracers`. Click on it, do M-. to navigate and observe the function:

`static void advance_tracers (GfsDomain * domain, gdouble dt){  GSList * i = domain->variables;  while (i) {    if (GFS_IS_VARIABLE_TRACER_VOF (i->data)) {      GfsVariableTracer * t = i->data;       t->advection.dt = dt;      gfs_tracer_vof_advection (domain, &t->advection);      gfs_domain_variable_centered_sources (domain, i->data, i->data, t->advection.dt);    }    else if (GFS_IS_VARIABLE_TRACER (i->data)) {      GfsVariableTracer * t = i->data;       t->advection.dt = dt;      gfs_tracer_advection_diffusion (domain, &t->advection);      gfs_domain_cell_traverse (domain,				FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,				(FttCellTraverseFunc) GFS_VARIABLE1 (t)->fine_coarse, t);    }    i = i->next;  }  }`

We see that we have two sub-functions `gfs_tracer_vof_advection` and `gfs_tracer_advection_diffusion`. One uses the VOF scheme to advect non-diffusive interfaces. The second advects traditional tracers such as temperature.

We concentrate for now on the second type of advection and, clicking on the corresponding function we see:

`/** * gfs_tracer_advection_diffusion: * @domain: a #GfsDomain. * @par: the advection parameters. * * Advects the @v field of @par using the current face-centered (MAC) * velocity field. */void gfs_tracer_advection_diffusion (GfsDomain * domain,				     GfsAdvectionParams * par){  GfsSourceDiffusion * d;   g_return_if_fail (domain != NULL);  g_return_if_fail (par != NULL);   gfs_domain_timer_start (domain, "tracer_advection_diffusion");   if ((d = source_diffusion (par->v))) {    GfsVariable * rhs;     rhs = gfs_temporary_variable (domain);    gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,			      (FttCellTraverseFunc) gfs_cell_reset, rhs);    variable_sources (domain, par, rhs, NULL);    variable_diffusion (domain, d, par, rhs, NULL);    gts_object_destroy (GTS_OBJECT (rhs));  }  else {    variable_sources (domain, par, par->v, NULL);    gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, par->v);  }   gfs_domain_timer_stop (domain, "tracer_advection_diffusion");}`

At this stage the reader may be a little lost: where is the function that advects the tracer ? We do see the lines

`gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,			      (FttCellTraverseFunc) gfs_cell_reset, rhs);`

It involves the important function gfs_domain_cell_traverse.

### gfs_domain_cell_traverse

In the example above, this function applies the function ` (FttCellTraverseFunc) gfs_cell_reset ` to the data ` rhs ` in each cell. It is often used in the code to perform operations local to each cell. However in the example above the call to this function does not do much: it justs resets `rhs` to zero. We must find the actual advection in another function at a deeper level.

### variable_sources (domain, par, rhs, NULL)

Having a look at this function (we do not display it here, but you should look at it now using you text editor) reveals a large number of functions. The gfs_domain_cell_traverse function is called many times to apply various functions. One of the important functions for the implementation of the u grad T term is the function par->flux . This is not a fixed name for a function but a pointer to a function that defines the advection flux. We cannot click on it to determine how it is defined ! Doing a grep on flux reveals many lines. Doing a

```% grep par\-\>flux *.c
```

reveals the definition of the flux function in the file `advection.c` in the module `gfs_advection_params_read`. There are actually several options for this function (you can see them in the code). The option can be selected in the simulation .gfs file. We can see it in `advection.c` . Only one of them is the default one. Where is the default initialized ? It is not initialized in the object initialization function, but in the module `simulation_init() ` of `simulation.c` .

`object->advection_params.flux = gfs_face_velocity_advection_flux;`

If we look in this function we see

`flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt    /ftt_cell_size (face->cell);`

so now we know that the flux has the well-known form

```flux = u_face T_face
```

and we know how the code implements the scheme !

## The end for now

Go to course top level. Go to top of page. Go to Next Session.