Programming the Advection Scheme

From Gerris

Contents 1 The core of a fluid simulation 2 Miscellaneous 2.1 Cast to the parent class 2.2 GFS_VARIABLE macro 2.3 The FttCell class 2.4 Tip: Navigate using the TAGS file 3 More about classes: 3.1 SourceScalar 3.1.1 SourceScalar Header 3.1.2 SourceScalar Object methods 3.1.2.1 the read method 4 The end for now

The core of a fluid simulation

Let us assume that as a computational physicists you are most interested in how the Navier-Stokes equations are read. You are familiar with Stephane's J. Comput. Phys. article and you want to see how it is implemented in code. The core of the simulation in simulation.c specifically the following lines. We start with some 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;
  }

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);
 
    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);
}

Now that we are more advanced, we are ready for a tip.

Miscellaneous

Cast to the parent class

The example shows how one can extract information from an object if one knows in which parent it is stored. For example, if we want to use the information of the event variable which came from the GfsEvent class, we can have access to it by means of the following variable "parentdata":

GfsEvent * parentdata = GFS_EVENT (event);

In this way, we load the information of the "event" object coming from its parent GfsEvent into this variable and we can either use it or manipulate if we want. For example parentdata->step will give us the temporal step which have been read in the input.

GFS_VARIABLE macro

This macro is used everywhere in the code to access the value of a GFS_VARIABLE in a cell. An example of its use is found in the init_periodic object showed in the previous session:

GFS_VARIABLE (cell, v->i) = sin (2.*(*m)*M_PI*pos.x)*cos (2.*(*m)*M_PI*pos.y);

This expression replaces

GFS_STATE (cell)->v = sin (2.*init->m*M_PI*pos.x)*cos (2.*init->m*M_PI*pos.y);

which was used in the tutorial.

To use this macro, there are two points of view.

• In the "magical" approach we just now that GFS_VARIABLE(somecell,somevariable->i)allows us to write to the variable somevariablein the cell somecell. For us dummies a question may arise: what value should the index i have ? Is it the cell number ? Well i is not an index: it is a member of each GfsVariable. We do not have to set it: it is already there in the variable.

• In the rational approach, we would like to read the code and understand the definition of GFS_VARIABLE. However the definition of the GFS_VARIABLE macro is a bit baffling. It can be found in fluid.h. (It is a good idea to try the following: find any place in the code where the GFS_VARIABLE macro is used ---that should be easy--- and then click on it. Type M-. in emacs ---or the equivalent in vim--- and get to the definition.)

#define GFS_VARIABLE(cell, index)     ((&GFS_STATE (cell)->place_holder)[index])

Notice that we have a special member of the cell state that tells us where the variable is stored. The macro GFS_STATE is defined by

#define GFS_STATE(cell)               ((GfsStateVector *) (cell)->data)

The FttCell class

This class is defined as follows:

struct _FttCell {
  /*< public >*/
  guint flags;
  gpointer data;
 
  /*< private >*/
  struct _FttOct * parent, * children;
};

The data member of this class is a gpointer. Remember from last session that it is an untyped pointer, so there is no difficulty to cast to a GfsStateVector by using

struct _GfsStateVector {
  /* temporary face variables */
  GfsFaceStateVector f[FTT_NEIGHBORS];
 
  /* solid boundaries */
  GfsSolidVector * solid;
 
  gdouble place_holder;
};

Tip: Navigate using the TAGS file

- Create the TAGS file :

% cd src
% make tags

You can then use the tags in emacs, nedit or vim.

1. In emacs:
   1. open any *.c or *.h file with emacs.
   2. Position the cursor on a function or variable.
   3. do ESC . to find its definition. (or M-. using the emacs Meta key convention )
   4. do ESC * to return to the previous location (s) .
2. In vim, nedit ... : vim, nedit users please contribute ...

More about classes:

SourceScalar

SourceScalar Header

The SourceScalar class from source.h has some interesting differences with the three previous examples :

/* GfsSourceScalar: Header */
 
typedef struct _GfsSourceScalar         GfsSourceScalar;
 
struct _GfsSourceScalar {
  /*< private >*/
  GfsSourceGeneric parent;
 
  /*< public >*/
  GfsVariable * v;
};
 
#define GFS_SOURCE_SCALAR(obj)            GTS_OBJECT_CAST (obj,\
			GfsSourceScalar,\
			gfs_source_scalar_class ())
#define GFS_IS_SOURCE_SCALAR(obj)         (gts_object_is_from_class (obj,\
			gfs_source_scalar_class ()))
 
GfsSourceGenericClass * gfs_source_scalar_class  (void);

Notice that we do not use all the possibilities in the template. The class structure is a SourceGenericClass, there is no need to create a specific SourceScalarClass .

SourceScalar Object methods

the read method

The source scalar class has a complex read method. From source.c

static void source_scalar_read (GtsObject ** o, GtsFile * fp)
{
  GfsSourceScalar * source;
  GfsDomain * domain;
 
  if (GTS_OBJECT_CLASS (gfs_source_scalar_class ())->parent_class->read)
    (* GTS_OBJECT_CLASS (gfs_source_scalar_class ())->parent_class->read) 
      (o, fp);
  if (fp->type == GTS_ERROR)
    return;

This is standard Gerris “read” code that uses the parent class read method to read first the parameters that are defined in the parent class. It will initialize the variables of the parent object.

source = GFS_SOURCE_SCALAR (*o);

This “casts” the “generic” object **o into an object of type “source_scalar”.

domain =  GFS_DOMAIN (gfs_object_simulation (source));
  if (fp->type != GTS_STRING) {
    gts_file_error (fp, "expecting a string (GfsVariable)");
    return;
  }
  source->v = gfs_variable_from_name (domain->variables, 
		fp->token->str);
  if (source->v == NULL) {
    gts_file_error (fp, "unknown variable `%s'", fp->token->str);
    return;
  }

We expect a string: the name of the variable. If the name of the variable is not found in the list of domain gfs_variables, an error results.

If it is found, the member v of the structure source becomes this gfs_variable.

if (source->v->sources == NULL)
    source->v->sources = 
      gts_container_new (GTS_CONTAINER_CLASS (gts_slist_container_class ()));
  gts_container_add (source->v->sources, GTS_CONTAINEE (source));

The variable then aquires a new source in its list of sources v->sources . gts_containers are defined in the gts library.

gts_file_next_token (fp);
}

The end for now

---

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