Programming the Advection Scheme

From Gerris

(Difference between revisions)
Jump to: navigation, search
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.)
← Previous diff
Revision as of 03:27, 20 November 2007
Zaleski (Talk | contribs)
(A total rewrite. More in-depth analysis of the BCG scheme and its implementation.)
Next diff →
Line 1: Line 1:
-== How the equations are solved ==+== A large-scale view of the code ==
-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 ?+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 ? For this we first look at the large scale view, encompassing the basic operations of the code, initialization, velocity advection and tracer advection. We then turn to the detailed analysis of tracer advection.
-=== simulation_run() ===+=== 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. +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. One way to do it would be to see what is executed line by line using a tool such as [http://en.wikipedia.org/wiki/Ddd 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: We start with some more bookkeeping:
<source lang="c"> <source lang="c">
Line 48: Line 48:
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> </source>
 +
=== Inside the time loop === === Inside the time loop ===
Now we are fully inside the time loop. Now we are fully inside the time loop.
Line 105: Line 106:
</source> </source>
-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. +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 coding of tracer advection ==
 + 
 +=== Basic scheme, structures and modules for tracer advection ===
 + 
 +==== The Bell-Collela-Glaz scheme and its adaptation ====
 + 
 +Gerris implements a variant of [http://en.wikipedia.org/wiki/Godunov's_scheme Godunov's scheme] for advection that is second order in time and space. The relevant reference is :
 + 
 +(BCG) J. B. Bell, P. Colella, H. M. Glaz - A second-order projection method for the incompressible Navier-Stokes equations, ''J. Comput. Phys.'' '''85''', 257-283, (1989).
 + 
 +The BCG scheme is further simplified following :
 + 
 +D. F. Martin, An adaptive cell-centered pro jection method for the incompressible Euler
 +equations, Ph.D. thesis, University of California, Berkeley (1998).
 + 
 +D. F. Martin, P. Colella, A cell-centered adaptive projection method for the
 +incompressible euler equations, ''J. Comput. Phys.'' '''163''', 271-312 (2000).
 + 
 +The advection scheme is also discussed in Stephane's [http://gfs.sf.net/gerris.pdf Journal of Computational Physics article], Section 5, page 21.
 + 
 +To fully understand the code, you need to be familiar with the relevant parts of the above references. In these course notes we will not reproduce the corresponding developments, as this is a course about programming and not numerical schemes, but the live course sessions should include a short presentation of the version of the BCG scheme used by Gerris.
 + 
 +==== advance_tracers() module ====
The tracer is updated though the function <code>advance_tracers</code>. Click on it, do M-. to navigate and observe the function: The tracer is updated though the function <code>advance_tracers</code>. Click on it, do M-. to navigate and observe the function:
Line 136: Line 159:
</source> </source>
We see that we have two sub-functions <code>gfs_tracer_vof_advection</code> and 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.+<code>gfs_tracer_advection_diffusion</code>. One uses the VOF scheme to advect non-diffusive interfaces. The second advects traditional tracers such as temperature. To understand the advection scheme we need to look at several important structures.
-=== gfs_tracer_advection_diffusion ===+==== GfsVariableTracer object ====
 + 
 +This an object class similar to those we have seen before. The full header for that
 +class is:
 +<source lang="c">
 +/* GfsVariableTracer: header */
 + 
 +typedef struct _GfsVariableTracer GfsVariableTracer;
 + 
 +struct _GfsVariableTracer {
 + /*< private >*/
 + GfsVariable parent;
 + 
 + /*< public >*/
 + GfsAdvectionParams advection;
 +};
 + 
 +#define GFS_VARIABLE_TRACER(obj) GTS_OBJECT_CAST (obj,\
 + GfsVariableTracer,\
 + gfs_variable_tracer_class ())
 +#define GFS_IS_VARIABLE_TRACER(obj) (gts_object_is_from_class (obj,\
 + gfs_variable_tracer_class ()))
 + 
 +GfsVariableClass * gfs_variable_tracer_class (void);
 +</source>
 +As we shall see below, however, we do not need the full header to understand the functionality of this object. All we need to know is
 +<source lang="c">
 +struct _GfsVariableTracer {
 + /*< private >*/
 + GfsVariable parent;
 + 
 + /*< public >*/
 + GfsAdvectionParams advection;
 +};
 +</source>
 +and even more economically, ''all we need'' is to know about the member ''advection''.
 + GfsAdvectionParams advection;
 +So we turn to the defintion of its type:
 + 
 +==== GfsAdvectionParams struct ====
 +<source lang="c">
 +struct _GfsAdvectionParams {
 + gdouble cfl, dt; /* self-explanatory */
 + GfsVariable * v, * fv, ** u, ** g; /* v: typically, the GfsVariable that we are advecting.
 + fv: typically the flux of the GfsVariable.
 + u[c] is the component of velocity in direction c,
 + c goes form 0 to 2 in 3D, from 0 to 1 in 2D .
 + g[c] is the component of a gradient variable.
 + */
 + GfsCenterGradient gradient; /* A function pointer, so that the function
 + par->gradient(cell, c, v->i) should return
 + the value of the c component of the gradient
 + of variable v at the center of the cell. */
 +
 + gboolean use_centered_velocity; /* self-explanatory */
 + GfsUpwinding upwinding; /* self-explanatory */
 + GfsFaceAdvectionFluxFunc flux; /* A function pointer, that is used to compute fluxes */
 + GfsAdvectionScheme scheme; /* The scheme used */
 + gboolean average; /* Whether to use averaging in gfs_advection_update() */
 +};
 +</source>
 +(The comments are not in the source code, they were added for this course.) First we notice that this structure
 +is not a Gerris Object: it does not have a parent, or a related class, and is just a C structure. However it has an initialisation function, and of course it carries a lot of information. Vectors are of the form <code> ** u</code>, so that <code> u[c]</code> is a pointer to the ''c'' component of the vector ''u''.
 + 
 +==== gfs_tracer_advection_diffusion() module ====
We concentrate for now on the second type of advection and, clicking on the corresponding function we see: We concentrate for now on the second type of advection and, clicking on the corresponding function we see:
Line 183: Line 270:
(FttCellTraverseFunc) gfs_cell_reset, rhs); (FttCellTraverseFunc) gfs_cell_reset, rhs);
</source> </source>
-It involves the important function ''gfs_domain_cell_traverse''.+It involves the important function ''gfs_domain_cell_traverse''. We explain it in the next subsection. The function variable_sources is the one that advects the tracer, we explain it in the next section.
-=== gfs_domain_cell_traverse ===+=== gfs_domain_cell_traverse() module ===
-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. +This is an important function which is very frequently used in the code. In a sense it is the equivalent of a '''do loop''' in a classical fortran or C code.
-We must find the actual advection in another function at a deeper level.+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.
 +Here is an excerpt of its documentation as it appears in the C source code:
 +<source lang="c">
 +/**
 + * gfs_domain_cell_traverse:
 + * @domain: a #GfsDomain.
 + * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
 + * %FTT_POST_ORDER.
 + * @flags: which types of children are to be visited.
 + * @max_depth: the maximum depth of the traversal. Cells below this
 + * depth will not be traversed. If @max_depth is -1 all cells in the
 + * tree are visited.
 + * @func: the function to call for each visited #FttCell.
 + * @data: user data to pass to @func.
 + *
 + * Traverses the cell trees of @domain. Calls the given function for
 + * each cell visited.
 + */
 +</source>
 +When <code>FTT_PRE_ORDER</code> is used the function is applied first to the root and then to the children. It progressively descends until max_depth, or the leaves, are visited. When <code>FTT_POST_ORDER</code> is used the function is applied first to the leaves and then to the children. The flag <code>FTT_TRAVERSE_LEAFS </code> is pretty explicit. Another possible value for the flag would be <code>FTT_TRAVERSE_NON_LEAFS</code>. Notice that when the flag is set to <code>FTT_TRAVERSE_LEAFS </code> the argument <code>@order</code> has no effect. Notice also that we do not give explicit values to the various order and flag arguments, instead we use values defined in macros, which makes them more explicit.
-=== variable_sources (domain, par, rhs, NULL)=== 
-Having a look at this function (we do not display it here, but you should look at it+The module <code>gfs_domain_cell_traverse</code> is a very useful and widely used function. However in the example above the call to <code>gfs_cell_reset</code> does not do much: it justs resets <code>rhs</code> to zero. We must find the actual advection in another function at a deeper level.
-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+ 
 +----
 +'''Note'''
 + 
 +As usual it is important to separate functionality from implementation. It is enough to understand how the Gerris code operate to learn about functionality. If we want to know more about the <code>gfs_domain_cell_traverse</code> function, we may feel the urge to look at its source code. This may not be a good idea, especially at this stage of the course. Anyway, should we embark on it, it requires, as usual, a lot of navigation. We explore in a sequence the functions:
 + 
 +# gfs_domain_cell_traverse
 +# gts_container_foreach
 +# static void box_traverse
 +# void ftt_cell_traverse
 +# cell_traverse_pre_order_all
 + 
 +(''Navigate on your screen through the functions.'')
 + 
 +----
 + 
 +=== variable_sources() module ===
 + 
 +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. These functions are
 + 
 +* gfs_cell_reset
 +* gfs_cell_advected_face_values
 +* par->flux (a fonction pointer, which will be elucidated below to be gfs_face_advection_flux)
 + 
 +We detail the various functions below.
 + 
 +==== gfs_cell_reset ====
 + 
 +Again the call to gfs_cell_reset does not do much: it justs resets ''par->fv'' to zero.
 + 
 +==== gfs_cell_advected_face_values ====
 + 
 +The documentation section above this function says what it does:
 +<source lang="c">
 +/**
 + * gfs_cell_advected_face_values:
 + * @cell: a #FttCell.
 + * @par: the advection parameters.
 + *
 + * Fills the face variable (@v field of #GfsFaceStateVector) of all the
 + * faces of @cell with the advected value of variable @par->v at time
 + * t + dt/2.
 + */
 +void gfs_cell_advected_face_values (FttCell * cell,
 + const GfsAdvectionParams * par)
 +</source>
 +The ''v field of GfsFaceStateVector'' is a special field defined on the faces which we will use to store the values of the tracer variable to be advected.
 + 
 +What is the variable par->v ? To find out, we need to look for where the GfsAdvectionParam variable ''par'' is initialized. On can search using grep as in the paragraph about ''par->flux'' below. Finally one discovers that initialization is done in ''variable.c'' as follows:
 +<source lang="c">
 +static void variable_tracer_init (GfsVariableTracer * v)
 +{
 + gfs_advection_params_init (&v->advection);
 + v->advection.gradient = gfs_center_van_leer_gradient;
 + v->advection.flux = gfs_face_advection_flux;
 + v->advection.v = GFS_VARIABLE1 (v);
 + v->advection.fv = NULL;
 + GFS_VARIABLE1 (v)->description = g_strdup ("Tracer");
 +}
 +</source>
 +where the ''GFS_VARIABLE1'' macro is
 +<source lang="c">
 +#define GFS_VARIABLE1(obj) GTS_OBJECT_CAST (obj,\
 + GfsVariable,\
 + gfs_variable_class ())
 +</source>
 +So finally the par->v variable is just the tracer variable we are advecting, which should not be surprising.
 +Now, ''look at the source of the function''.
 +This function appears to be rather complex.
 + 
 +<source lang="c">
 +void gfs_cell_advected_face_values (FttCell * cell,
 + const GfsAdvectionParams * par)
 +{
 + FttComponent c;
 + gdouble size;
 + GfsStateVector * s;
 + 
 + g_return_if_fail (cell != NULL);
 + g_return_if_fail (par != NULL);
 + 
 + s = GFS_STATE (cell);
 + size = ftt_cell_size (cell);
 + for (c = 0; c < FTT_DIMENSION; c++) {
 + gdouble unorm = par->use_centered_velocity ?
 + par->dt*GFS_VARIABLE (cell, par->u[c]->i)/size :
 + par->dt*(s->f[2*c].un + s->f[2*c + 1].un)/(2.*size);
 + gdouble g = (* par->gradient) (cell, c, par->v->i);
 + gdouble vl = GFS_VARIABLE (cell, par->v->i) + MIN ((1. - unorm)/2., 0.5)*g;
 + gdouble vr = GFS_VARIABLE (cell, par->v->i) + MAX ((- 1. - unorm)/2., -0.5)*g;
 + gdouble src = par->dt*gfs_variable_mac_source (par->v, cell)/2.;
 + gdouble dv;
 + 
 +#if FTT_2D
 + dv = transverse_term (cell, size, FTT_ORTHOGONAL_COMPONENT (c), par);
 +#else /* FTT_3D */
 + static FttComponent orthogonal[FTT_DIMENSION][2] = {
 + {FTT_Y, FTT_Z}, {FTT_X, FTT_Z}, {FTT_X, FTT_Y}
 + };
 + 
 + dv = transverse_term (cell, size, orthogonal[c][0], par);
 + dv += transverse_term (cell, size, orthogonal[c][1], par);
 +#endif /* FTT_3D */
 + 
 + s->f[2*c].v = vl + src - dv;
 + s->f[2*c + 1].v = vr + src - dv;
 + }
 +}
 +</source>
 +[[Image:u_LR.png|Figure 1 (from Dan Martin's PhD thesis)|300px|thumbnail]]
 +It actually implements the [[Session 3 General architecture of the code#The Bell-Collela-Glaz scheme | BCG scheme]]. In the code <code>vl</code> and <code>vr</code> are the advected "left" and ""right" value of Figure 1. Look at page 21, Eq. (19) of Stephane's [http://gfs.sf.net/gerris.pdf Journal of Computational Physics article] for the expression of the "Left State" and "Right State" of the faces (you should really be familiar with Godunov's method to follow this). The gradient ''g'' is computed using the ''gradient'' function specified in the''gradient'' member of GfsAdvectionParam par . The member is initialized in variable.c in the function
 + 
 + static void variable_tracer_init (GfsVariableTracer * v)
 + 
 +shown [[Session 3 General architecture of the code#gfs_cell_advected_face_values | above]]. The gradient is defined there as the ''gfs_center_van_leer_gradient()'' module, a self explanatory name.
 + 
 +==== par->flux ====
 + 
 +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 flux *.c
 +reveals too many lines to be useful. Doing a
<pre> <pre>
% grep par\-\>flux *.c % grep par\-\>flux *.c
</pre> </pre>
-reveals the definition of the flux function in the file <code>advection.c</code>+seems to reveal the definition of the flux function in the file <code>advection.c</code>
-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 the module <code>gfs_advection_params_read</code>. There are actually several options for this function. (''Look at the source code:'' you can see the various options.) The option can be selected in the simulation file (for instance ''simulation.gfs''). Only one of them is the default one. So where is the default initialized? (''Navigate to the initialization function of the object <code>GfsAdvectionParams</code>''). Navigation, as you see, does not help.
-of <code>simulation.c</code> .+We need to do a
 + % grep flux *.c | grep =
 +which yields two initializations, one in the module <code>simulation_init()</code>
 +of <code>simulation.c</code>:
<source lang="c"> <source lang="c">
object->advection_params.flux = gfs_face_velocity_advection_flux; object->advection_params.flux = gfs_face_velocity_advection_flux;
</source> </source>
 +and another in <code>variable.c</code> :
 +<source lang="c">
 +v->advection.flux = gfs_face_advection_flux;
 +</source>
 +The first piece of code initializes the function ''flux'' for the velocity advection "in" the simulation parameters of the ''object'' of type <code>GfsSimulation</code>. This object describes the whole Navier-Stokes resolution procedure, and its <code>AdvectionParam</code> is naturally that for velocity advection in Navier--Stokes. The other piece of code initializes it in the initialisation function of <code>GfsVariableTracer</code>.
-=== gfs_face_velocity_advection_flux ===+==== gfs_face_advection_flux ====
-If we look in this function we see +If we look in this function we see the comments
<source lang="c"> <source lang="c">
-flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt+/**
- /ftt_cell_size (face->cell);+ * gfs_face_advection_flux:
 + * @face: a #FttCellFace.
 + * @par: the advection parameters.
 + *
 + * Adds to variable @par->fv, the value of the (conservative)
 + * advection flux of the face variable through @face.
 + *
 + * This function assumes that the face variable has been previously
 + * defined using gfs_cell_advected_face_values().
 + */
</source> </source>
-so now we know+and then some source code. One part is particularly telling:
 +<source lang="c">
 +flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt*
 + gfs_face_upwinded_value (face, GFS_FACE_UPWINDING, NULL)/ftt_cell_size (face->cell);
 +</source>
 +so now we have the hint that
that the flux has the well-known form that the flux has the well-known form
<pre> <pre>
flux = u_face T_face flux = u_face T_face
</pre> </pre>
-and we know how the code implements the scheme ! 
-== The end for now ==+The face variable T_face is computed in an ''upwinded'' manner following the [[Session 3 General architecture of the code#The Bell-Collela-Glaz scheme | BCG scheme]]. It uses the following function:
 + 
 +===== gfs_face_upwinded_value =====
 + 
 +<source lang="c">
 +/**
 + * gfs_face_upwinded_value:
 + * @face: a #FttCellFace.
 + * @upwinding: type of upwinding.
 + * @u: the cell-centered velocity.
 + *
 + * This function assumes that the face variable has been previously
 + * defined using gfs_cell_advected_face_values().
 + *
 + * Returns: the upwinded value of the face variable.
 + */
 +gdouble gfs_face_upwinded_value (const FttCellFace * face,
 + GfsUpwinding upwinding,
 + GfsVariable ** u)
 +</source>
 + 
 + 
 + 
 +To see how this function works, we need to understand some of the functions defined in ''ftt.h'' and ''ftt.c''. These are part of the Fully Threaded Tree structure explained in the next session.
---- ----
-Go to [[ Gerris Flow Solver Programming Course for Dummies | course top level]]. Go to [[Session 3 General architecture of the code#The core of a fluid simulation | top of page]]. Go to [[Session 4 | Next Session]].+Go to [[ Gerris Flow Solver Programming Course for Dummies | course top level]]. Go to [[Session 3 General architecture of the code#How the equations are solved | top of page]]. Go to [[Session 4 The Fully Threaded Tree | Next Session]].

Revision as of 03:27, 20 November 2007

Contents

A large-scale view of the code

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 ? For this we first look at the large scale view, encompassing the basic operations of the code, initialization, velocity advection and tracer advection. We then turn to the detailed analysis of tracer advection.

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. One 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 coding of tracer advection

Basic scheme, structures and modules for tracer advection

The Bell-Collela-Glaz scheme and its adaptation

Gerris implements a variant of Godunov's scheme for advection that is second order in time and space. The relevant reference is :

(BCG) J. B. Bell, P. Colella, H. M. Glaz - A second-order projection method for the incompressible Navier-Stokes equations, J. Comput. Phys. 85, 257-283, (1989).

The BCG scheme is further simplified following :

D. F. Martin, An adaptive cell-centered pro jection method for the incompressible Euler equations, Ph.D. thesis, University of California, Berkeley (1998).

D. F. Martin, P. Colella, A cell-centered adaptive projection method for the incompressible euler equations, J. Comput. Phys. 163, 271-312 (2000).

The advection scheme is also discussed in Stephane's Journal of Computational Physics article, Section 5, page 21.

To fully understand the code, you need to be familiar with the relevant parts of the above references. In these course notes we will not reproduce the corresponding developments, as this is a course about programming and not numerical schemes, but the live course sessions should include a short presentation of the version of the BCG scheme used by Gerris.

advance_tracers() module

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. To understand the advection scheme we need to look at several important structures.

GfsVariableTracer object

This an object class similar to those we have seen before. The full header for that class is:

/* GfsVariableTracer: header */
 
typedef struct _GfsVariableTracer GfsVariableTracer;
 
struct _GfsVariableTracer {
/*< private >*/
GfsVariable parent;
 
/*< public >*/
GfsAdvectionParams advection;
};
 
#define GFS_VARIABLE_TRACER(obj) GTS_OBJECT_CAST (obj,\
GfsVariableTracer,\
gfs_variable_tracer_class ())
#define GFS_IS_VARIABLE_TRACER(obj) (gts_object_is_from_class (obj,\
gfs_variable_tracer_class ()))
 
GfsVariableClass * gfs_variable_tracer_class (void);

As we shall see below, however, we do not need the full header to understand the functionality of this object. All we need to know is

struct _GfsVariableTracer {
/*< private >*/
GfsVariable parent;
 
/*< public >*/
GfsAdvectionParams advection;
};

and even more economically, all we need is to know about the member advection.

 GfsAdvectionParams advection;

So we turn to the defintion of its type:

GfsAdvectionParams struct

struct _GfsAdvectionParams {
gdouble cfl, dt; /* self-explanatory */
GfsVariable * v, * fv, ** u, ** g; /* v: typically, the GfsVariable that we are advecting.
fv: typically the flux of the GfsVariable.
u[c] is the component of velocity in direction c,
c goes form 0 to 2 in 3D, from 0 to 1 in 2D .
g[c] is the component of a gradient variable.
*/

GfsCenterGradient gradient; /* A function pointer, so that the function
par->gradient(cell, c, v->i) should return
the value of the c component of the gradient
of variable v at the center of the cell. */

 
gboolean use_centered_velocity; /* self-explanatory */
GfsUpwinding upwinding; /* self-explanatory */
GfsFaceAdvectionFluxFunc flux; /* A function pointer, that is used to compute fluxes */
GfsAdvectionScheme scheme; /* The scheme used */
gboolean average; /* Whether to use averaging in gfs_advection_update() */
};

(The comments are not in the source code, they were added for this course.) First we notice that this structure is not a Gerris Object: it does not have a parent, or a related class, and is just a C structure. However it has an initialisation function, and of course it carries a lot of information. Vectors are of the form ** u, so that u[c] is a pointer to the c component of the vector u.

gfs_tracer_advection_diffusion() module

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. We explain it in the next subsection. The function variable_sources is the one that advects the tracer, we explain it in the next section.

gfs_domain_cell_traverse() module

This is an important function which is very frequently used in the code. In a sense it is the equivalent of a do loop in a classical fortran or C code. 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. Here is an excerpt of its documentation as it appears in the C source code:

/**
* gfs_domain_cell_traverse:
* @domain: a #GfsDomain.
* @order: the order in which the cells are visited - %FTT_PRE_ORDER,
* %FTT_POST_ORDER.
* @flags: which types of children are to be visited.
* @max_depth: the maximum depth of the traversal. Cells below this
* depth will not be traversed. If @max_depth is -1 all cells in the
* tree are visited.
* @func: the function to call for each visited #FttCell.
* @data: user data to pass to @func.
*
* Traverses the cell trees of @domain. Calls the given function for
* each cell visited.
*/

When FTT_PRE_ORDER is used the function is applied first to the root and then to the children. It progressively descends until max_depth, or the leaves, are visited. When FTT_POST_ORDER is used the function is applied first to the leaves and then to the children. The flag FTT_TRAVERSE_LEAFS is pretty explicit. Another possible value for the flag would be FTT_TRAVERSE_NON_LEAFS. Notice that when the flag is set to FTT_TRAVERSE_LEAFS the argument @order has no effect. Notice also that we do not give explicit values to the various order and flag arguments, instead we use values defined in macros, which makes them more explicit.


The module gfs_domain_cell_traverse is a very useful and widely used function. However in the example above the call to gfs_cell_reset does not do much: it justs resets rhs to zero. We must find the actual advection in another function at a deeper level.



Note

As usual it is important to separate functionality from implementation. It is enough to understand how the Gerris code operate to learn about functionality. If we want to know more about the gfs_domain_cell_traverse function, we may feel the urge to look at its source code. This may not be a good idea, especially at this stage of the course. Anyway, should we embark on it, it requires, as usual, a lot of navigation. We explore in a sequence the functions:

  1. gfs_domain_cell_traverse
  2. gts_container_foreach
  3. static void box_traverse
  4. void ftt_cell_traverse
  5. cell_traverse_pre_order_all

(Navigate on your screen through the functions.)


variable_sources() module

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. These functions are

  • gfs_cell_reset
  • gfs_cell_advected_face_values
  • par->flux (a fonction pointer, which will be elucidated below to be gfs_face_advection_flux)

We detail the various functions below.

gfs_cell_reset

Again the call to gfs_cell_reset does not do much: it justs resets par->fv to zero.

gfs_cell_advected_face_values

The documentation section above this function says what it does:

/**
* gfs_cell_advected_face_values:
* @cell: a #FttCell.
* @par: the advection parameters.
*
* Fills the face variable (@v field of #GfsFaceStateVector) of all the
* faces of @cell with the advected value of variable @par->v at time
* t + dt/2.
*/

void gfs_cell_advected_face_values (FttCell * cell,
const GfsAdvectionParams * par)

The v field of GfsFaceStateVector is a special field defined on the faces which we will use to store the values of the tracer variable to be advected.

What is the variable par->v ? To find out, we need to look for where the GfsAdvectionParam variable par is initialized. On can search using grep as in the paragraph about par->flux below. Finally one discovers that initialization is done in variable.c as follows:

static void variable_tracer_init (GfsVariableTracer * v)
{
gfs_advection_params_init (&v->advection);
v->advection.gradient = gfs_center_van_leer_gradient;
v->advection.flux = gfs_face_advection_flux;
v->advection.v = GFS_VARIABLE1 (v);
v->advection.fv = NULL;
GFS_VARIABLE1 (v)->description = g_strdup ("Tracer");
}

where the GFS_VARIABLE1 macro is

#define GFS_VARIABLE1(obj)            GTS_OBJECT_CAST (obj,\
GfsVariable,\
gfs_variable_class ())

So finally the par->v variable is just the tracer variable we are advecting, which should not be surprising. Now, look at the source of the function. This function appears to be rather complex.

void gfs_cell_advected_face_values (FttCell * cell,
const GfsAdvectionParams * par)
{
FttComponent c;
gdouble size;
GfsStateVector * s;
 
g_return_if_fail (cell != NULL);
g_return_if_fail (par != NULL);
 
s = GFS_STATE (cell);
size = ftt_cell_size (cell);
for (c = 0; c < FTT_DIMENSION; c++) {
gdouble unorm = par->use_centered_velocity ?
par->dt*GFS_VARIABLE (cell, par->u[c]->i)/size :
par->dt*(s->f[2*c].un + s->f[2*c + 1].un)/(2.*size);
gdouble g = (* par->gradient) (cell, c, par->v->i);
gdouble vl = GFS_VARIABLE (cell, par->v->i) + MIN ((1. - unorm)/2., 0.5)*g;
gdouble vr = GFS_VARIABLE (cell, par->v->i) + MAX ((- 1. - unorm)/2., -0.5)*g;
gdouble src = par->dt*gfs_variable_mac_source (par->v, cell)/2.;
gdouble dv;
 
#if FTT_2D
dv = transverse_term (cell, size, FTT_ORTHOGONAL_COMPONENT (c), par);
#else /* FTT_3D */
static FttComponent orthogonal[FTT_DIMENSION][2] = {
{FTT_Y, FTT_Z}, {FTT_X, FTT_Z}, {FTT_X, FTT_Y}
};
 
dv = transverse_term (cell, size, orthogonal[c][0], par);
dv += transverse_term (cell, size, orthogonal[c][1], par);
#endif /* FTT_3D */
 
s->f[2*c].v = vl + src - dv;
s->f[2*c + 1].v = vr + src - dv;
}
}
Figure 1 (from Dan Martin's PhD thesis)
Enlarge
Figure 1 (from Dan Martin's PhD thesis)

It actually implements the BCG scheme. In the code vl and vr are the advected "left" and ""right" value of Figure 1. Look at page 21, Eq. (19) of Stephane's Journal of Computational Physics article for the expression of the "Left State" and "Right State" of the faces (you should really be familiar with Godunov's method to follow this). The gradient g is computed using the gradient function specified in thegradient member of GfsAdvectionParam par . The member is initialized in variable.c in the function

static void variable_tracer_init (GfsVariableTracer * v) 

shown above. The gradient is defined there as the gfs_center_van_leer_gradient() module, a self explanatory name.

par->flux

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 flux *.c 

reveals too many lines to be useful. Doing a

% grep par\-\>flux *.c 

seems to reveal 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. (Look at the source code: you can see the various options.) The option can be selected in the simulation file (for instance simulation.gfs). Only one of them is the default one. So where is the default initialized? (Navigate to the initialization function of the object GfsAdvectionParams). Navigation, as you see, does not help. We need to do a

% grep flux *.c  | grep =

which yields two initializations, one in the module simulation_init() of simulation.c:

object->advection_params.flux = gfs_face_velocity_advection_flux;

and another in variable.c :

v->advection.flux = gfs_face_advection_flux;

The first piece of code initializes the function flux for the velocity advection "in" the simulation parameters of the object of type GfsSimulation. This object describes the whole Navier-Stokes resolution procedure, and its AdvectionParam is naturally that for velocity advection in Navier--Stokes. The other piece of code initializes it in the initialisation function of GfsVariableTracer.

gfs_face_advection_flux

If we look in this function we see the comments

/**
* gfs_face_advection_flux:
* @face: a #FttCellFace.
* @par: the advection parameters.
*
* Adds to variable @par->fv, the value of the (conservative)
* advection flux of the face variable through @face.
*
* This function assumes that the face variable has been previously
* defined using gfs_cell_advected_face_values().
*/

and then some source code. One part is particularly telling:

flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt*
gfs_face_upwinded_value (face, GFS_FACE_UPWINDING, NULL)/ftt_cell_size (face->cell);

so now we have the hint that that the flux has the well-known form

flux = u_face T_face

The face variable T_face is computed in an upwinded manner following the BCG scheme. It uses the following function:

gfs_face_upwinded_value
/**
* gfs_face_upwinded_value:
* @face: a #FttCellFace.
* @upwinding: type of upwinding.
* @u: the cell-centered velocity.
*
* This function assumes that the face variable has been previously
* defined using gfs_cell_advected_face_values().
*
* Returns: the upwinded value of the face variable.
*/

gdouble gfs_face_upwinded_value (const FttCellFace * face,
GfsUpwinding upwinding,
GfsVariable ** u)


To see how this function works, we need to understand some of the functions defined in ftt.h and ftt.c. These are part of the Fully Threaded Tree structure explained in the next session.


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

Personal tools
communication