# Gerris Flow Solver Programming Course for Dummies

(Difference between revisions)
 Revision as of 23:46, 29 October 2007Popinet (Talk | contribs) (Use tag for formatting C code)← Previous diff Revision as of 23:49, 29 October 2007Zaleski (Talk | contribs) (Wrote the material for Session 1)Next diff → Line 1: Line 1: - This wiki is an unfinished project ! + == Preamble == - == Preface == + This course material is about Gerris, a general-purpose fluid mechanics code developped by Stephane Popinet at NIWA, Wellington, New Zealand. The course is actually being taught as of writing (october 2007) in Paris by Stephane Zaleski. + The intended audience is typical first year science or engineering graduate students with either very little experience of C or with some Fortran knowledge, but willing to work hard and learn. The student should know simple C data types, pointers and functions but not structures. - This course is about Gerris, a general-purpose fluid mechanics code developped by Stephane Popinet at NIWA, Wellington, New Zealand . + Gerris is a free, GPL-licensed, open source code available at [[http://gfs.sf.net http://gfs.sf.net]] . - The intended audience is typical first year graduate students with very little experience of C or some Fortran knowledge, but willing to work hard and learn. + == Session 1 == - + - Gerris is a free, open source code available at [[http://gfs.sf.net]] . + - + - == Course 1 == + Line 80: Line 77: - +
struct Point {                                                                                                                                                                                                                                                                                                                                                                                                 struct Point {
char  name;                                                                                                                                                                                                                                                                                                                                                                                                    char  name;
double  x, y;                                                                                                                                                                                                                                                                                                                                                                                                  double  x, y;
};                                                                                                                                                                                                                                                                                                                                                                                                             };
-                                                                                                                                                                                                                                                                                                                                                                                                              +
An example of usage An example of usage - +
﻿main()                                                                                                                                                                                                                                                                                                                                                                                                        ﻿main()
{                                                                                                                                                                                                                                                                                                                                                                                                              {
-                                                                                                                        struct Point my_point;    /* declaration */

-                                                                                                                        my_point.x = 0. ;                                                                                                                                                                                                                                                                              +                                                                                                               struct Point my_point;    /* declaration */
-                                                                                                                        my_point.y = 1.;                                                                                                                                                                                                                                                                               +
-                                                                                                                        my_point.name = ‘A’;                                                                                                                                                                                                                                                                           +                                                                                                               my_point.x = 0. ;
+                                                                                                               my_point.y = 1.;
+                                                                                                               my_point.name = ‘A’;
}                                                                                                                                                                                                                                                                                                                                                                                                              }
-                                                                                                                                                                                                                                                                                                                                                                                                              +
name, x and y are members of the structure of type “Point” called my_point. We also give it a name that can be exported (printed, passed to other functions) as a character. This example shows why it is useful to use structures to store several relevant informations or data together. name, x and y are members of the structure of type “Point” called my_point. We also give it a name that can be exported (printed, passed to other functions) as a character. This example shows why it is useful to use structures to store several relevant informations or data together. Line 104: Line 102: ==== An example: the structure GfsNorm in Gerris ==== ==== An example: the structure GfsNorm in Gerris ==== - +

struct _GfsNorm {                                                                                                                                                                                                                                                                                                                                                                                              struct _GfsNorm {
Line 110:                                                                                                                                                                                                                                                                                                                                                                                                               Line 108:
};                                                                                                                                                                                                                                                                                                                                                                                                             };
typedef struct _GfsNorm GfsNorm                                                                                                                                                                                                                                                                                                                                                                                typedef struct _GfsNorm GfsNorm
-                                                                                                                                                                                                                                                                                                                                                                                                              +
From ''domain.h''. From ''domain.h''. - +
﻿GfsNorm gfs_domain_norm_residual (GfsDomain * domain,                                                                                                                                                                                                                                                                                                                                                         ﻿GfsNorm gfs_domain_norm_residual (GfsDomain * domain,
FttTraverseFlags flags,                                                                                                                                                                                                                                                                                                                                                                                        FttTraverseFlags flags,
Line 144:                                                                                                                                                                                                                                                                                                                                                                                                               Line 142:
return n;                                                                                                                                                                                                                                                                                                                                                                                                      return n;
}                                                                                                                                                                                                                                                                                                                                                                                                              }
-                                                                                                                                                                                                                                                                                                                                                                                                              +
From ''domain.c'' From ''domain.c'' Line 153: Line 151: * glib functions ''g_return_val_if_fail'' * glib functions ''g_return_val_if_fail'' + + + ===== Note ===== + + This is what the [[http://library.gnome.org/devel/glib/unstable/index.html glib documentation]] says about gpointer + +
typedef void* gpointer;
+
+                                                                                                               An untyped pointer. gpointer looks better and is easier to use than void*.
+
==== Structures that are related to other structures ==== ==== Structures that are related to other structures ==== + + More complex relationships between structures create the need for ''inheritance''. + + Below one defines two structures which have some similarities. + The ''Square'' is a kind of ''Quadrangle''. + +
+                                                                                                               Point A, B, C, D;
+                                                                                                               };
+
+
+ + Notice that we use typedef: this creates an "alias" for the structure name and will + make it able to refer to itself in its list of member names. + +
+                                                                                                               ﻿/* Related structure */
+
+                                                                                                               struct _Square {
+                                                                                                               Point A, B, C, D;
+                                                                                                               };
+
+                                                                                                               typedef struct _Square Square
+
+ + In the following ''pseudo-code'' we use the structures, for instance we define + a function that computes the area. + +
+                                                                                                               {
+                                                                                                               Triangle ABC;
+                                                                                                               ABC = create_triangle(X.A,X.B,X.C);
+                                                                                                               …
+                                                                                                               return area_Triangle (ABC) + area_Triangle (CDA)  ;
+                                                                                                               }
+
+ + If we use squares, we may need a very similar function to compute the area of + squares. +
+                                                                                                               area_Square ( Square Y )
+                                                                                                               {
+                                                                                                               …
+                                                                                                               return area_Triangle (ABC) + area_Triangle (CDA)  ;
+                                                                                                               }
+
+ + + === Elementary classes === + + * We want to avoid such repetitions. + + * Notion of parent structure: Quadrangle is the parent of Square. + One shall from now on use the “class” terminology. A class is a structure with some + functions attached to it, the functions are called methods. This should become clearer in + what follows. + + ==== Simple example : structures with inheritance ==== + +
+                                                                                                               struct _Square {
+                                                                                                               double longueur_du_cote;
+                                                                                                               };
+
+                                                                                                               typedef struct   _Square    Square
+
+                                                                                                               main()
+                                                                                                               {
+                                                                                                               Square MySquare ;   /* Déclaration of Square */
+                                                                                                               double area;
+
+                                                                                                               Initialize_Square(&MySquare) ;  /*Initialisation */
+                                                                                                               area = area_Square ( MySquare ) ;
+                                                                                                               }
+
+ + Problem: for each class we need a new function giving the area + area_Square, area_Quadrangle, etc… + So we will define the function "area_square" quite simply from "area_quadrangle". + Functions that apply to a class are called the methods of the class. + + In C++ the methods are listed when the class is defined. + In the Gerris/Glib Object system the situation is somewhat different. + A structure like the above "Square" structure is created. + Then a separated C structure is created, the ''class structure'' + that contains the methods and other information about the class. + + ---- + + This is the skeleton of such a class in C: + +
+                                                                                                               struct _Square {
+                                                                                                               double side_length;
+                                                                                                               };
+
+                                                                                                               typedef struct  _Square Square
+                                                                                                               typedef struct  _SquareClass  SquareClass
+
+                                                                                                               struct _SquareClass {
+                                                                                                               /*< private >*/
+
+                                                                                                               /*< public >*/
+                                                                                                               /* add extra methods here */
+                                                                                                               };
+
+ + ---- + + ===== The gtstemplate tool ===== + + * Automatic class creation: Usage of the Gts library utility function gtstemplate to create a Gts Object : + +
+                                                                                                               % gtstemplate
+                                                                                                               Usage: gtstemplate [OPTIONS] Class ParentClass
+                                                                                                               Options:
+                                                                                                               [--no-extra-data]
+                                                                                                               [--no-extra-method]
+
+ + ---- + Here is the result: + +
+
+                                                                                                               typedef struct _Square         Square;
+
+                                                                                                               struct _Square {
+                                                                                                               /*< private >*/
+
+                                                                                                               /*< public >*/
+                                                                                                               /* add extra data here (if public) */
+                                                                                                               };
+
+                                                                                                               typedef struct _SquareClass    SquareClass;
+
+                                                                                                               struct _SquareClass {
+                                                                                                               /*< private >*/
+
+                                                                                                               /*< public >*/
+                                                                                                               /* add extra methods here */
+                                                                                                               };
+
+                                                                                                               #define SQUARE(obj)            GTS_OBJECT_CAST (obj,\
+                                                                                                               Square,\
+                                                                                                               Square_class ())
+                                                                                                               #define SQUARE_CLASS(klass)    GTS_OBJECT_CLASS_CAST (klass,\
+                                                                                                               SquareClass,\
+                                                                                                               Square_class())
+                                                                                                               #define IS_SQUARE(obj)         (gts_object_is_from_class (obj,\
+                                                                                                               Square_class ()))
+
+                                                                                                               SquareClass * Square_class  (void);
+                                                                                                               Square * Square_new    (SquareClass * klass);
+
+                                                                                                               /* Square: Object */
+
+                                                                                                               static void Square_class_init (SquareClass * klass)
+                                                                                                               {
+                                                                                                               /* define new methods and overload inherited methods here */
+
+                                                                                                               }
+                                                                                                               static void Square_init (Square * object)
+                                                                                                               {
+                                                                                                               /* initialize object here */
+                                                                                                               }
+
+                                                                                                               SquareClass * Square_class (void)
+                                                                                                               {
+                                                                                                               static SquareClass * klass = NULL;
+
+                                                                                                               if (klass == NULL) {
+                                                                                                               GtsObjectClassInfo Square_info = {
+                                                                                                               "Square",
+                                                                                                               sizeof (Square),
+                                                                                                               sizeof (SquareClass),
+                                                                                                               (GtsObjectClassInitFunc) Square_class_init,
+                                                                                                               (GtsObjectInitFunc) Square_init,
+                                                                                                               (GtsArgSetFunc) NULL,
+                                                                                                               (GtsArgGetFunc) NULL
+                                                                                                               };
+                                                                                                               klass = gts_object_class_new (GTS_OBJECT_CLASS (Quadrangle_class ()),
+                                                                                                               &Square_info);
+                                                                                                               }
+
+                                                                                                               return klass;
+                                                                                                               }
+
+                                                                                                               Square * Square_new (SquareClass * klass)
+                                                                                                               {
+                                                                                                               Square * object; object = Square (gts_object_new (GTS_OBJECT_CLASS (klass)));
+
+
+
+                                                                                                               return object;
+                                                                                                               }
+
+
+ + ===== An simpler, actual example ===== + + Typical class usage in Gerris is less complex. An actual example from Gerris: the SourceScalar class from source.h : +
+
+                                                                                                               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 . + + ===== Object manipulation: the read method ===== + + Many Gfs objects have a read method. + From source.c +
+                                                                                                               static void source_scalar_read (GtsObject ** o, GtsFile * fp)
+                                                                                                               {
+                                                                                                               GfsSourceScalar * source;
+                                                                                                               GfsDomain * domain;
+
+                                                                                                               (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 ()));
+
+ 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);
+                                                                                                               }
+
+ + == Session 2 == + + To be written.

## Preamble

This course material is about Gerris, a general-purpose fluid mechanics code developped by Stephane Popinet at NIWA, Wellington, New Zealand. The course is actually being taught as of writing (october 2007) in Paris by Stephane Zaleski.

The intended audience is typical first year science or engineering graduate students with either very little experience of C or with some Fortran knowledge, but willing to work hard and learn. The student should know simple C data types, pointers and functions but not structures.

Gerris is a free, GPL-licensed, open source code available at [http://gfs.sf.net] .

## Session 1

### Introduction

This course is about how Gerris ([[1]] ) is programmed. It is intended to help in understanding the Gerris source code and learning how to modify it usefully.

This course is not about the numerical methods in Gerris (however it would be good for students of this course to learn about them, for instance in the J. Comput. Phys. article) .

Gerris is closely linked to Gts, the GNU triangulated surface library, also written by Stephane Popinet

Gerris and Gts are programmed in a style analogous to that of Glib, Gnome and GTK . It is a style of C programming that offers several advantages:

• Most aspects of Object-Oriented Programming (OOP) , such as the existence of classes with their own methods and inheritance.
• The ability to interface to other programming languages. (As far as I know, this feature is not used in Gerris/Gts)

To implement Object-Oriented Programming, Gerris/Gts uses its own “Object system”. This system is analogous to the Glib object system (Gobject), but not identical to it. Learning more about Gobject can be very useful.

### References

• Introduction and in-depth discussion of the Gobject system (in order of decreasing ease of reading):
• Introduction to the Gerris/Gts object system :

You can find one in the [Gerris Tutorial Section on Objects].

### Useful tips

#### Navigate using the TAGS file in emacs or vim

- Create the TAGS file :
```% cd src
% make tags
```
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:
1. to be written ...

#### No order in which to read the code

There is no good order in which to read the code (I have not found it yet)

#### Keep a C precedence and associativity table nearby

A lot of macros and functions such g_assert come from the Glib. Keep a bookmark to the Glib documentation: [[2]]

### C basics

#### Introduction to Structures

```struct Point {
char  name;
double  x, y;
};
```

An example of usage

```﻿main()
{

struct Point my_point;    /* declaration */

my_point.x = 0. ;
my_point.y = 1.;
my_point.name = ‘A’;
}
```

name, x and y are members of the structure of type “Point” called my_point. We also give it a name that can be exported (printed, passed to other functions) as a character. This example shows why it is useful to use structures to store several relevant informations or data together.

#### An example: the structure GfsNorm in Gerris

```
struct _GfsNorm {
gdouble bias, first, second, infty, w;
};
typedef struct _GfsNorm GfsNorm
```

From domain.h.

```﻿GfsNorm gfs_domain_norm_residual (GfsDomain * domain,
FttTraverseFlags flags,
gint max_depth,
gdouble dt,
GfsVariable * res)
{
GfsNorm n;
gpointer data[2];

g_return_val_if_fail (domain != NULL, n);
g_return_val_if_fail (res != NULL, n);

gfs_norm_init (&n);
data[0] = res;
data[1] = &n;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, flags, max_depth,
#ifdef HAVE_MPI
domain_norm_reduce (domain, &n);
#endif /* HAVE_MPI */
gfs_norm_update (&n);

dt *= dt;
n.bias *= dt;
n.first *= dt;
n.second *= dt;
n.infty *= dt;
return n;
}
```

From domain.c

Notice the use of

• glib basic types gdouble, gpointer
• glib functions g_return_val_if_fail

##### Note

This is what the [glib documentation] says about `gpointer`

```typedef void* gpointer;

An untyped pointer. gpointer looks better and is easier to use than void*.
```

#### Structures that are related to other structures

More complex relationships between structures create the need for inheritance.

Below one defines two structures which have some similarities. The Square is a kind of Quadrangle.

```struct _Quadrangle {
Point A, B, C, D;
};

```

Notice that we use typedef: this creates an "alias" for the structure name and will make it able to refer to itself in its list of member names.

```﻿/* Related structure */

struct _Square {
Point A, B, C, D;
};

typedef struct _Square Square
```

In the following pseudo-code we use the structures, for instance we define a function that computes the area.

```area_Quadrangle (Quadrangle X )
{
Triangle ABC;
ABC = create_triangle(X.A,X.B,X.C);
…
return area_Triangle (ABC) + area_Triangle (CDA)  ;
}
```

If we use squares, we may need a very similar function to compute the area of squares.

```area_Square ( Square Y )
{
…
return area_Triangle (ABC) + area_Triangle (CDA)  ;
}
```

### Elementary classes

• We want to avoid such repetitions.
• Notion of parent structure: Quadrangle is the parent of Square.

One shall from now on use the “class” terminology. A class is a structure with some functions attached to it, the functions are called methods. This should become clearer in what follows.

#### Simple example : structures with inheritance

```struct _Square {
double longueur_du_cote;
};

typedef struct   _Square    Square

main()
{
Square MySquare ;   /* Déclaration of Square */
double area;

Initialize_Square(&MySquare) ;  /*Initialisation */
area = area_Square ( MySquare ) ;
}
```

Problem: for each class we need a new function giving the area ` area_Square, area_Quadrangle, ` etc… So we will define the function "area_square" quite simply from "area_quadrangle". Functions that apply to a class are called the methods of the class.

In C++ the methods are listed when the class is defined. In the Gerris/Glib Object system the situation is somewhat different. A structure like the above "Square" structure is created. Then a separated C structure is created, the class structure that contains the methods and other information about the class.

This is the skeleton of such a class in C:

```struct _Square {
double side_length;
};

typedef struct  _Square Square
typedef struct  _SquareClass  SquareClass

struct _SquareClass {
/*< private >*/

/*< public >*/
/* add extra methods here */
};
```

##### The gtstemplate tool
• Automatic class creation: Usage of the Gts library utility function gtstemplate to create a Gts Object :
```% gtstemplate
Usage: gtstemplate [OPTIONS] Class ParentClass
Options:
[--no-extra-data]
[--no-extra-method]
```

Here is the result:

```/* Square: Header */

typedef struct _Square         Square;

struct _Square {
/*< private >*/

/*< public >*/
/* add extra data here (if public) */
};

typedef struct _SquareClass    SquareClass;

struct _SquareClass {
/*< private >*/

/*< public >*/
/* add extra methods here */
};

#define SQUARE(obj)            GTS_OBJECT_CAST (obj,\
Square,\
Square_class ())
#define SQUARE_CLASS(klass)    GTS_OBJECT_CLASS_CAST (klass,\
SquareClass,\
Square_class())
#define IS_SQUARE(obj)         (gts_object_is_from_class (obj,\
Square_class ()))

SquareClass * Square_class  (void);
Square * Square_new    (SquareClass * klass);

/* Square: Object */

static void Square_class_init (SquareClass * klass)
{
/* define new methods and overload inherited methods here */

}
static void Square_init (Square * object)
{
/* initialize object here */
}

SquareClass * Square_class (void)
{
static SquareClass * klass = NULL;

if (klass == NULL) {
GtsObjectClassInfo Square_info = {
"Square",
sizeof (Square),
sizeof (SquareClass),
(GtsObjectClassInitFunc) Square_class_init,
(GtsObjectInitFunc) Square_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
klass = gts_object_class_new (GTS_OBJECT_CLASS (Quadrangle_class ()),
&Square_info);
}

return klass;
}

Square * Square_new (SquareClass * klass)
{
Square * object; object = Square (gts_object_new (GTS_OBJECT_CLASS (klass)));

return object;
}

```
##### An simpler, actual example

Typical class usage in Gerris is less complex. An actual example from Gerris: the SourceScalar class from `source.h` :

```/* 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 .

##### Object manipulation: the read method

Many Gfs objects have a read method. From `source.c `

```static void source_scalar_read (GtsObject ** o, GtsFile * fp)
{
GfsSourceScalar * source;
GfsDomain * domain;

(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 ()));
```

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

To be written.