Programming the Advection Scheme

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 04:07, 6 November 2007
Popinet (Talk | contribs)
(Fixed broken link)
← Previous diff
Revision as of 19:08, 8 November 2007
Zaleski (Talk | contribs)
(A total rewrite of this section, SourceClass removed, concentrates on how advection is coded.)
Next diff →
Line 1: Line 1:
-== The core of a fluid simulation ==+== 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 ?
-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 <code>simulation.c</code> specifically the following lines. We start with some bookkeeping: +=== simulation_run() ===
 + 
 +Let us assume that you are familiar with Stephane's [http://gfs.sf.net/gerris.pdf 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 <code>simulation.c</code>, specifically in <code>simulation_run()</code>. Let us look at this function.
 +We start with some more bookkeeping:
<source lang="c"> <source lang="c">
static void simulation_run (GfsSimulation * sim) static void simulation_run (GfsSimulation * sim)
Line 27: Line 31:
} }
</source> </source>
- +Then we start the various steps of the integration of a time step.:
<source lang="c"> <source lang="c">
gfs_simulation_set_timestep (sim); gfs_simulation_set_timestep (sim);
Line 43: Line 47:
gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim); gts_container_foreach (GTS_CONTAINER (sim->events), (GtsFunc) gfs_event_do, sim);
 +</source>
 +=== Inside the time loop ===
 +Now we are fully inside the time loop.
 +<source lang="c">
gfs_simulation_set_timestep (sim); gfs_simulation_set_timestep (sim);
Line 97: Line 105:
</source> </source>
-Now that we are more advanced, we are ready for a tip.+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.
-== Miscellaneous == +=== Advance tracers ===
- +
-=== 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":+
 +The tracer is updated though the function <code>advance_tracers</code>. Click on it, do M-. to navigate and observe the function:
<source lang="c"> <source lang="c">
-GfsEvent * parentdata = GFS_EVENT (event);+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;
 + }
 +}
</source> </source>
 +We see that we have two sub-functions <code>gfs_tracer_vof_advection</code> and
 +<code>gfs_tracer_advection_diffusion</code>. One uses the VOF scheme to advect non-diffusive interfaces. The second advects traditional tracers such as temperature.
-In this way, we load the information of the "event" object coming from+=== gfs_tracer_advection_diffusion ===
-its parent GfsEvent into this variable and we can either use it or+
-manipulate if we want. For example+
-<code> parentdata->step </code> will give us the temporal step which have been+
-read in the input.+
-=== <code>GFS_VARIABLE</code> macro ===+We concentrate for now on the second type of advection and, clicking on the corresponding function we see:
- +
-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: +
<source lang="c"> <source lang="c">
-GFS_VARIABLE (cell, v->i) = sin (2.*(*m)*M_PI*pos.x)*cos (2.*(*m)*M_PI*pos.y);+/**
-</source>+ * gfs_tracer_advection_diffusion:
-This expression replaces+ * @domain: a #GfsDomain.
-<source lang="c">+ * @par: the advection parameters.
-GFS_STATE (cell)->v = sin (2.*init->m*M_PI*pos.x)*cos (2.*init->m*M_PI*pos.y);+ *
-</source>+ * Advects the @v field of @par using the current face-centered (MAC)
-which was used in the tutorial. + * velocity field.
 + */
 +void gfs_tracer_advection_diffusion (GfsDomain * domain,
 + GfsAdvectionParams * par)
 +{
 + GfsSourceDiffusion * d;
-To use this macro, there are two points of view. + g_return_if_fail (domain != NULL);
 + g_return_if_fail (par != NULL);
-* In the "magical" approach we just now that <code>GFS_VARIABLE(somecell,somevariable->i)</code>allows us to write to the variable <code>somevariable</code>in the cell <code>somecell</code>. 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. + gfs_domain_timer_start (domain, "tracer_advection_diffusion");
-* In the rational approach, we would like to read the code and understand the definition of <code>GFS_VARIABLE</code>. However the definition of the <code>GFS_VARIABLE</code> macro is a bit baffling. It can be found in <code>fluid.h</code>. (It is a good idea to try the following: find any place in the code where the <code>GFS_VARIABLE</code> 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.)+ if ((d = source_diffusion (par->v))) {
-<source lang="c">+ GfsVariable * rhs;
-#define GFS_VARIABLE(cell, index) ((&GFS_STATE (cell)->place_holder)[index])+
-</source>+
-Notice that we have a special member of the cell state that tells us where the variable is stored. The macro <code>GFS_STATE</code> is defined by+
-<source lang="c">+
-#define GFS_STATE(cell) ((GfsStateVector *) (cell)->data)+
-</source>+
-=== The <code> FttCell</code> class ===+ 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);
 + }
-This class is defined as follows: + gfs_domain_timer_stop (domain, "tracer_advection_diffusion");
 +}
 +</source>
 +At this stage the reader may be a little lost: where is the function that advects the tracer ? We do see the lines
<source lang="c"> <source lang="c">
-struct _FttCell {+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- /*< public >*/+ (FttCellTraverseFunc) gfs_cell_reset, rhs);
- guint flags;+
- gpointer data;+
- +
- /*< private >*/+
- struct _FttOct * parent, * children;+
-};+
</source> </source>
 +It involves the important function ''gfs_domain_cell_traverse''.
-The data member of this class is a gpointer. Remember from last session that it is+=== gfs_domain_cell_traverse ===
-an untyped pointer, so there is no difficulty to cast to a <code>GfsStateVector</code>+
-by using+
-<source lang="c">+
-struct _GfsStateVector {+
- /* temporary face variables */+
- GfsFaceStateVector f[FTT_NEIGHBORS];+
- /* solid boundaries */+In the example above, this function applies the function <code> (FttCellTraverseFunc) gfs_cell_reset </code> to the data <code> rhs </code> 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 <code>rhs</code> to zero.
- GfsSolidVector * solid;+We must find the actual advection in another function at a deeper level.
- gdouble place_holder;+=== variable_sources (domain, par, rhs, NULL)===
-};+
-</source>+
-=== Tip: Navigate using the TAGS file ===+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
- - Create the TAGS file : <pre>+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+<pre>
-% make tags+% grep par\-\>flux *.c
</pre> </pre>
- +reveals the definition of the flux function in the file <code>advection.c</code>
-You can then use the tags in emacs, nedit or vim.+in the module <code>gfs_advection_params_read</code>. 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 <code>advection.c</code> . 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 <code>simulation_init() </code>
-# In emacs: +of <code>simulation.c</code> .
-## open any *.c or *.h file with emacs.+
-## Position the cursor on a function or variable. +
-## do ESC . to find its definition. (or M-. using the emacs Meta key convention ) +
-## do ESC * to return to the previous location (s) . +
-# In vim, nedit ... : vim, nedit users please contribute ...+
- +
-== More about classes: == +
- +
-=== SourceScalar ===+
- +
-==== SourceScalar Header ====+
- +
-The SourceScalar class from <code>source.h</code> has some interesting differences with the three previous examples : +
<source lang="c"> <source lang="c">
-/* GfsSourceScalar: Header */+object->advection_params.flux = gfs_face_velocity_advection_flux;
- +
-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);+
</source> </source>
-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 ====+=== gfs_face_velocity_advection_flux ===
-===== the read method =====+If we look in this function we see
- +
-The source scalar class has a complex read method. +
-From <code>source.c </code>+
<source lang="c"> <source lang="c">
-static void source_scalar_read (GtsObject ** o, GtsFile * fp)+flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt
-{+ /ftt_cell_size (face->cell);
- 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;+
-</source>+
-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 lang="c">+
- source = GFS_SOURCE_SCALAR (*o);+
-</source>+
-This “casts” the “generic” object **o+
-into an object of type “source_scalar”. +
-<source lang="c">+
- 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;+
- }+
-</source>+
-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. +
-<source lang="c">+
- 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));+
-</source>+
-The variable then aquires a new source in its list of sources v->sources . +
-gts_containers are defined in the gts library.+
-<source lang="c"> +
- gts_file_next_token (fp);+
-}+
</source> </source>
 +so now we know
 +that the flux has the well-known form
 +<pre>
 +flux = u_face T_face
 +</pre>
 +and we know how the code implements the scheme !
== The end for now == == The end for now ==

Revision as of 19:08, 8 November 2007

Contents

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.

Advance tracers

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.

gfs_tracer_advection_diffusion

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;

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.

Personal tools
communication