An Example of Use of the Gerris Object System
From Gerris
Contents |
A new class to initialize periodic velocity fields
In this session, we are going to attempt to create a new class that will initialize the velocity field in the form of an array of vortices. While it is possible to do this with a little C code section in the simulation file, we will do it here as an educational exercise.
A template for new object classes
GTS comes with a simple script called gtstemplate which will generate a C code template for the new class we want to create. For a summary of its syntax just type:
% gtstemplate
Usage: gtstemplate [OPTIONS] Class ParentClass
Options:
[--no-extra-data]
[--no-extra-method]
[--overload=METHOD]
As an example we are going to try and create a template for a new initial condition class called InitPeriodic. Just type:
% gtstemplate --overload=read --overload=write --overload=event \ InitPeriodic GfsGenericInit > init_periodic.c
We have just created a template for a new class called InitPeriodic derived from GfsGenericInit and where the read, write and event methods are overloaded. Fire up your favourite editor and have a look at the file generated: init_periodic.c. The first thing you see is that the file is divided in two sections starting with
/* InitPeriodic: Header */
and
/* InitPeriodic: Object */
As in the previous examples these sections correspond respectively to the declaration of the structures and functions (header part) and to the corresponding definitions (object part). If we were to use these functions and structures in a library, these two parts would be in separate files (a .h and a corresponding .c file). For what we are interested in they are fine staying in the same file. There is a function we will not need: init_periodic_new, just remove the lines declaring and defining it in the header and object sections. We are left with only one function: init_periodic_class. This function essentially registers our new class and its associated attributes, which are things like:
- the name of the class: "InitPeriodic",
- the sizes of the object and class structures: sizeof (InitPeriodic) and sizeof (InitPeriodicClass),
- function to call to initialize the class: init_periodic_class_init,
- function to call to initialize a new object: init_periodic_init,
If we now look at what the init_periodic_class function does, we see that it indeed overloads the event, read and write methods of our class with locally defined functions. The GFS_EVENT_CLASS and GTS_OBJECT_CLASS macro calls are casting operators which change the type of our new class (InitPeriodicClass) to the types of one of its parent class. GtsObjectClass is the ancestor of all object classes and we see that the read and write methods are thus defined for all objects. The parent class GfsGenericInitClass is also an alias of the GfsEvent class (as we have seen before it is a special type of event occurring once at the start of the simulation) which has an associated event method.
The goal: initial conditions
With this preliminary overall understanding of what the template does, we are ready to define our new initialization class. What we want to do is initialize the 2D velocity field using the following functions:
u(x,y) = − cos(2mπ x)sin(2mπ y) (1)
v(x,y) = sin(2mπ x)cos(2mπ y) (2)
where m is a parameter. This is an exact stationary solution of the 2D Euler equations with periodic boundary conditions (each vortex is a kind of flywheel).
Header
We need to store the value of m somewhere. A clean way to do that is to add m as a data associated with our InitPeriodic object. To do this just modify the definition of the structure in the header part of the file, like this:
struct _InitPeriodic {
/*< private >*/
GfsGenericInit parent;
/*< public >*/
gdouble m;
};
Note again that the GfsGenericInit parent; a declaration of member data is the first member in the structure. The whole object-inheritance mechanism in C relies on this data alignment.
Object
We now need to specify the value of m. The first value we want to assign is the default value of m for an object newly created. We can do that using the initialisation function init_periodic_init which is called everytime a new InitPeriodic object is created (it is the constructor function of the InitPeriodicClass).
Init method
Just add:
static void init_periodic_init (InitPeriodic * object)
{
object->m = 1.;
}
Our default value for m is now one. What we really want is being able to specify the value of m in the parameter file used by gerris. We can do that by using the read and write methods.
Write method
We will start with the write method: init_periodic_write. This function first calls the write method of the parent class of our object:
/* call write method of parent */
if (GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->write)
(* GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->write)
(o, fp);
The parent class is given by the parent_class field of any GtsObjectClass. The write method might not be defined for the parent class, so it is always a good (safe) idea to first test that it is indeed defined.
Now that we have written the data associated with the parent class, we can write our own data like this:
static void init_periodic_write (GtsObject * o, FILE * fp)
{
/* call write method of parent */
if (GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->write)
(* GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->write)
(o, fp);
/* do object specific write here */
fprintf (fp, " %g", INIT_PERIODIC (o)->m);
}
Note that we used a space to separate this new data from the data of the parent class and that we did not add anything after our own data (no space or newline). This is so that we can eventually extend (through inheritance) this class by adding more data in exactly the same way.
Read Method
We now need to define a “symmetrical” read method which will read the value of m from the parameter file. In the init_periodic_read method, we first read the parameters associated with the parent class. The line:
if (fp->type == GTS_ERROR)
return;
just checks if an error occurred while reading parameters for the parent class, in which case the method returns prematurely. To read our parameter we are going to use functions associated with the GtsFile data type (have a look at the GTS reference manual for details):
static void init_periodic_read (GtsObject ** o, GtsFile * fp)
{
/* call read method of parent */
if (GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->read)
(* GTS_OBJECT_CLASS (init_periodic_class ())->parent_class->read)
(o, fp);
if (fp->type == GTS_ERROR)
return;
if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
gts_file_error (fp, "expecting a number (m)");
return;
}
INIT_PERIODIC (*o)->m = atof (fp->token->str);
/* do not forget to prepare for next read */
gts_file_next_token (fp);
}
We first check that the current token is an integer (GTS_INT) or a floating point number (GTS_FLOAT). If it is not we set an error message describing the problem and return prematurely. If it is, we convert the string describing the current token (fp->token->str) into a floating point number (using the standard C atof function) and assign it to m. Note that we need to use the type conversion macro INIT_PERIODIC because the argument passed to init_periodic_read is a generic GtsObject. As we are dealing with a method of our InitPeriodicClass we know that this GtsObject is also a InitPeriodic object and the type conversion is then legitimate.
Event method
We can now create, read and write our new object. However we are not doing anything with it yet. Like for any other GfsEvent, the action performed is controlled by the event method. If we look at the init_periodic_event function, we see that it returns a gboolean. If TRUE this return value means that the event took place. We then call the event method of the parent class and check its return value. If TRUE we do something and return TRUE. We know that the parent class is GfsInit, an event which takes place once at the start of the simulation. What we want to do is initialise the velocity field using our formula. We can do that like this:
typedef struct {
gdouble * m;
GfsVariable * u;
GfsVariable * v;
} MyData;
static void init_velocity (FttCell * cell, MyData * data)
{
FttVector pos;
gdouble * m = data->m;
GfsVariable * v = data->v;
GfsVariable * u = data->u;
ftt_cell_pos (cell, &pos);
GFS_VALUE (cell, u) = - cos (2.*(*m)*M_PI*pos.x)*sin (2.*(*m)*M_PI*pos.y);
GFS_VALUE (cell, v) = sin (2.*(*m)*M_PI*pos.x)*cos (2.*(*m)*M_PI*pos.y);
}
static gboolean init_periodic_event (GfsEvent * event, GfsSimulation * sim)
{
InitPeriodic * inputdata = INIT_PERIODIC(event);
if ((* GFS_EVENT_CLASS (GTS_OBJECT_CLASS (init_periodic_class ())\
->parent_class)->event) (event, sim)) {
/* do object-specific event here */
MyData data;
GfsDomain * domain = GFS_DOMAIN (sim);
data.m = &(inputdata->m);
data.u = gfs_variable_from_name (domain->variables, "U");
data.v = gfs_variable_from_name (domain->variables, "V");
gfs_domain_cell_traverse (GFS_DOMAIN (sim),
FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) init_velocity,
&data);
return TRUE;
}
return FALSE;
}
The gfs_domain_cell_traverse function traverses the cell tree and calls the init_velocity function for each leaf cell (FTT_TRAVERSE_LEAFS). A GfsSimulation is an object derived from GfsDomain which justifies the GFS_DOMAIN (sim) type casting. We shall have a look at these structures and functions later in the course.
We also pass an extra parameter to init_velocity: &data. This extra argument contains some data obtained from our object through the object InitPeriodic * init .
What happens in init_velocity is straightforward. We first get the position of the center of the cell using ftt_cell_pos and then assign the values of the two components of the velocities using our formula and the m parameter defined in init.
The full source of our new object can be found here
Discussion of the example
Cast to the class
The example shows how one can extract information from an object using a cast. This applies if we want information in a parent. For example, if we want to use the information of the sim variable which came from the parent GfsDomain class, we can have access to it by means of the following variable "domain":
GfsDomain * domain = GFS_DOMAIN (sim);
In this way, we load the information of the "sim" object coming from
its parent GfsDomain into this variable and we can either use it or
manipulate if we want. For example
domain->variables will give us the variables of the simulation
sim.
GFS_VALUE macro
This macro is used everywhere in the code to access the value of a GfsVariable in a cell. An example of its use is found in the init_periodic object above:
GFS_VALUE (cell, v) = sin (2.*(*m)*M_PI*pos.x)*cos (2.*(*m)*M_PI*pos.y);
By definition, the macro call GFS_VALUE(cell,variable) allows to read or
write to the GfsVariable variable in the cell cell.
Note For somebody trying to learn naively about the code, just stopping at that may seem irrational. Perhaps we would like to read more of the code and understand the definition of GFS_VALUE. This would violate the distinction between the higher-level and the lower-level parts of the code. We can do a lot in Gerris without ever knowing about the lower-level parts.
Testing our new class
Creating a module
We have now almost all that is needed for our new object. How do we use this new piece of code with gerris? What we want to do is to create a dynamically loadable module. First of all we need to check that we can compile the code. There are a few missing headers we need to add (at the top of the file):
#include <math.h>
#include <stdlib.h>
#include <gfs.h>
The gfs.h header contains all the function declarations for the gerris library. We can now try to compile using for example:
% cc `pkg-config gerris2D --cflags` -Wall -g -O2 -c init_periodic.c
where pkg-config gerris2D –cflags (quoted using inverted quotes) defines the compilation flags needed to use the 2D version of the gerris library.
To make a proper module we also need to add the following at the end of the file:
/* Initialize module */
const gchar * g_module_check_init (void);
const gchar * g_module_check_init (void)
{
init_periodic_class ();
return NULL;
}
This just tells gerris how to initialise the module after it has been loaded. The g_module_check_init function will be called. In our case, it does only one thing which is to instantiate our new class: this is necessary so that gerris registers how to handle it.
We are now all set to create our new module. Sadly, dynamically loadable module creation is not a standardised process and the command-line arguments vary from compiler to compiler and from system to system. In the following, I will assume that you use gcc on a linux box. If you are using another system supporting dynamically loadable modules, you will need to read your local manual. On a linux box:
% cc `pkg-config gerris2D --cflags` -Wall -g -O2 -c -fPIC init_periodic.c % cc -shared init_periodic.o -lc -o libinit_periodic2D.so
should work. Note that we indicate to pkg-config that this module uses the gerris2D library. We could have built both a 2D and a 3D version of the same module. At runtime the gerris executable uses this extension to check which version to load.
Installing the module
Our module is now ready to use. We just need to install it in a directory where it will be found by gerris. If you installed gerris in /home/joe/local just type:
% cp libinit_periodic2D.so /home/joe/local/lib/gerris
Using the module
We can now use it directly in a parameter file, for example:
1 2 GfsSimulation GfsBox GfsGEdge {} {
GfsTime { end = 50 }
GfsRefine 6
GModule init_periodic
InitPeriodic {} 2
GfsOutputTime { istep = 10 } stdout
GfsOutputProjectionStats { istep = 10 } stdout
GfsOutputPPM { step = 0.1 } vorticity-%4.2f.ppm { v = Vorticity }
}
GfsBox {}
1 1 right
1 1 top
The first part of the object definition: InitPeriodic {} is the generic GfsEvent definition, the second part: 2 is the parameter m we added. What if we do not specify the right parameter? Just try to replace InitPeriodic {} 2 with InitPeriodic {} a. You should get a message like
gerris: file 'periodic.gfs' is not a valid simulation file periodic.gfs:5:18: expecting a number (m)
which tells you where the error occurred (in file periodic.gfs, line 5, character 18) together with the error message which we specified in the read method of our class.
⇐ previous: The Gerris Object System, ⇑ up: Gerris Flow Solver Programming Course for Dummies, ⇒ next: GfsView