# Programming the Advection Scheme

### From Gerris

Revision as of 03:27, 20 November 2007Zaleski (Talk | contribs) (A total rewrite. More in-depth analysis of the BCG scheme and its implementation.) ← Previous diff |
Revision as of 20:19, 15 October 2012Dustinl (Talk | contribs) (Added a description of tag navigating for vim users) Next diff → |
||

Line 1: |
Line 1: | ||

== A large-scale view of the code == | == 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. | + | 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 to 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-diffusion. We then turn to the detailed analysis of tracer advection and diffusion. |

+ | |||

+ | === Equation and numerical Scheme === | ||

+ | |||

+ | The present section should is brief guide of how Gerris treats the advection-diffusion equation of a tracer, <math>\alpha</math>, | ||

+ | |||

+ | <math> | ||

+ | \frac{\partial \alpha}{\partial t} + \nabla \cdot ( \mathbf{u} \alpha) = \nabla \cdot ( D \nabla \alpha)+S_\alpha + \nabla \cdot \mathbf{s}_\alpha | ||

+ | </math> | ||

+ | |||

+ | where <math>\mathbf{u}</math> is the fluid velocity and <math>D<\math> is the diffusion coefficient. We have added for generality two types of sources: <math>S_\alpha</math> stands for typical sources included, while <math>\mathbf{s}_\alpha</math> would be for generic source vector flux. | ||

+ | |||

+ | Gerris adopts a scheme of the form, | ||

+ | |||

+ | <math> | ||

+ | \frac{\alpha^{n+1} -\alpha^{n}}{\Delta t} + \nabla \cdot ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) = (1-\beta) \nabla \cdot ( D^n \nabla \alpha^n)+\beta \nabla \cdot ( D^{n+1} \nabla \alpha^{n+1})+S_\alpha^n + \nabla \cdot \mathbf{s}^n_\alpha | ||

+ | </math> | ||

+ | |||

+ | where <math>\beta</math> is a factor that has to be between 0.5 and 1.0: <math>\beta=</math>0.5 (the default value) corresponds to the Crank-Nicholson scheme which is second order accurate in time; <math>\beta=</math>1.0 corresponds to the Euler Forward scheme which although is first order accurate in time is more robust against oscillations. The equation is latterly reordered in a form close to a Poisson's equation, | ||

+ | |||

+ | <math> | ||

+ | \alpha^{n+1} - \nabla \cdot ( \beta \: \Delta t \: D^{n+1} \nabla \alpha^{n+1}) = \alpha^{n} - \nabla \cdot (\Delta t \: \mathbf{u}^{n+1/2} \alpha^{n+1/2}) + \nabla \cdot ( (1-\beta) \: \Delta t \: D^n \nabla \alpha^n)+S_\alpha^n | ||

+ | </math> | ||

+ | |||

+ | Adopting the cell as the control volume, <math>d \Omega</math>, assuming constant value of the variable <math>\alpha</math> in the cell and integrating the above discrete equation over that control volume, it is obtained the following, | ||

+ | |||

+ | * For the temporal term: | ||

+ | <math> | ||

+ | \iiint_\Omega \alpha \: d \Omega = a_v \: \alpha \: \: d \ell_1 d \ell_2 d \ell_3 \, = a_v \: \alpha \: g_c \, h^3 . | ||

+ | </math> where <math>h</math> is the size of the cell. For generality we have used a general system of orthogonal curvilinear coordinates given by <math>(x_1,x_2,x_3)</math>. If <math>gm_i</math> are the metric factor of each coordinate <math>x_i</math> (<math>d \ell_i=gm_i d x_i</math>), we collect them in a single metric factor <math>g= gm_1 gm_2 gm_3</math>. The subscript <math>c</math> denotes that the metric factor is to be computed at the center of mass of the cell. <math>a_v</math> is the volume fraction (it is always 1. except in the case that the cell is cut by a solid segment, i.e a mixed cell). | ||

+ | |||

+ | * For the advection term: | ||

+ | To this term it is applied additionally the divergence theorem in order to calculate it by means of the advection fluxes through the cell's contour, | ||

+ | |||

+ | <math> | ||

+ | \begin{split} | ||

+ | \iiint_\Omega \nabla \cdot ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) \: \Delta t \: d \Omega &= \iint_\Sigma \Delta t ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) \cdot \mathbf{n} \: d \Sigma \\ &= \sum_{Contour} {u_f}^{n+1/2} \alpha_f^{n+1/2} \Delta t \, s_f g_f h^2 | ||

+ | \end{split} | ||

+ | </math> | ||

+ | |||

+ | where <math>{u_f}</math> and <math>{\alpha_f}</math> are the normal velocity to the face (or solid segment) and the value of the tracer at that face (or solid segment), respectively; <math>s_f</math> is the face fraction (it is always 1. except in the case that the face were cut by a solid segment). <math>g_f</math> is, again, a metric factor denoting the subscript <math>f</math> that is evaluated at the face. | ||

+ | |||

+ | *For the diffusion terms: | ||

+ | These terms are treated similarly to the advection term resulting in, | ||

+ | <math> | ||

+ | \begin{split} | ||

+ | \iiint_\Omega \nabla \cdot (\Delta t\: D \nabla \alpha) \: d \Omega &= \iint_\Sigma \Delta t\: D \nabla \alpha \cdot \mathbf{n} \: d \Sigma \\ &= \sum_{Contour} \Delta t\: D_f g_f s_f h (h \nabla^f /alpha) | ||

+ | \end{split} | ||

+ | </math> | ||

+ | |||

+ | Reordering the equation, it results finally in the following expression, | ||

+ | |||

+ | <math> | ||

+ | \begin{split} | ||

+ | \alpha^{n+1} &- \frac{1}{a_v g_c h^2}\sum_{contour} \frac{\Delta t \: D_f \,g_f \, s_f \, \beta}{gm_f} \left(h \nabla^f \alpha^{n+1} \right) =\\ \alpha^{n} &- \frac{1}{a_v g_c} \sum_{contour} \frac{\Delta t \, s_f g_f}{h} {u_f}^{n+1/2} \alpha_f^{n+1/2} + S_\alpha^n \, \Delta t + \frac{1}{a_v g_c} \sum_{Contour} s_f^n \Delta t \, s_f g_f/h \right \\ &+ \frac{1}{a_v g_c h^2} \sum_{contour} \frac{\Delta t \: D_f g_f \, s_f \, (1-\beta)}{gm_f} \left(h \nabla^f \alpha^{n} \right)_f | ||

+ | \end{split} | ||

+ | </math> | ||

+ | |||

+ | where we have take into account, by means of <math>gm_f</math>, that a orthogonal curvilinear coordinates system could affect to the expression of the normal gradient. | ||

=== 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. 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. | + | 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 by [[running Gerris through a debugger]]. 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 31: |
Line 89: | ||

} | } | ||

</source> | </source> | ||

- | Then we start the various steps of the integration of a time step.: | + | 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 106: |
Line 164: | ||

</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. For this we shall use an interesting feature of the text editors. |

+ | |||

+ | == Tip: Navigate using the TAGS file == | ||

+ | In the above example, we have seen that we can gain some information about the macros, classes, structures and functions of the code by navigating in the code from one definition to the next. This navigation is not always immediately enlightening: we must resort to it only when we have enough general understanding of the code. To perform this navigation in a more efficient way than <code>grep</code>, we can use the following trick: | ||

+ | * Create the TAGS file : | ||

+ | <pre> | ||

+ | % cd src | ||

+ | % make tags | ||

+ | </pre> | ||

+ | |||

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

+ | # In emacs: | ||

+ | ## 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: | ||

+ | ## Open any *.c or *.h file with vim. | ||

+ | ## To inform vim of the TAGS file, execute the following command: <code>set tags+=<PATH_TO_GERRIS_SRC>/TAGS</code>, where <code><PATH_TO_GERRIS_SRC></code> is the path to the directory containing the Gerris source code (and the TAGS file). Better yet, add the command to your .vimrc file to perform this automatically. | ||

+ | ## Position the cursor on a function or variable. | ||

+ | ## Ctrl-] (in normal/command mode) will open the appropriate source file in a new buffer and navigate to the function's or variable's definition. Your location prior to executing this command is pushed into the "tag stack." | ||

+ | ## Ctrl-T "pops" the tag stack; i.e., Ctrl-T returns you to the location prior to the last executed Ctrl-]. | ||

+ | You can now try to use this trick on your own editor. It is a good idea to try the following: find any place in the code where the <code>GFS_VALUE</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, then repeat the operation to display the declarations and definitions above. | ||

+ | |||

+ | From now on, the reader or listener to the course should open his text editor and navigate freely in the code, trying to find the implementation of the numerical methods. We guide the reader below for one particular scheme: the tracer advection. | ||

+ | |||

+ | == The coding of tracer advection-diffusion == | ||

- | == The coding of tracer advection == | ||

=== Basic scheme, structures and modules for tracer advection === | === Basic scheme, structures and modules for tracer advection === | ||

Line 163: |
Line 246: | ||

==== GfsVariableTracer object ==== | ==== GfsVariableTracer object ==== | ||

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

class is: | class is: | ||

<source lang="c"> | <source lang="c"> | ||

Line 242: |
Line 325: | ||

GfsSourceDiffusion * d; | 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))) { | if ((d = source_diffusion (par->v))) { | ||

GfsVariable * rhs; | GfsVariable * rhs; | ||

- | rhs = gfs_temporary_variable (domain); | + | par->fv = rhs = gfs_temporary_variable (domain); |

gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, | gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, | ||

- | (FttCellTraverseFunc) gfs_cell_reset, rhs); | + | (FttCellTraverseFunc) copy_v_rhs, par); |

- | variable_sources (domain, par, rhs, NULL); | + | variable_sources (domain, par, rhs, NULL, NULL); |

variable_diffusion (domain, d, par, rhs, NULL); | variable_diffusion (domain, d, par, rhs, NULL); | ||

- | gts_object_destroy (GTS_OBJECT (rhs)); | + | gts_object_destroy (GTS\_OBJECT (rhs)); |

} | } | ||

else { | 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"); | ||

} | } | ||

</source> | </source> | ||

- | At this stage the reader may be a little lost: where is the function that advects the tracer ? We do see the lines | + | It can be seen that a temporary variable, rhs (or equivalently, par-> fv), is created. It is in this variable, rhs, where the right hand side of the advection-diffusion equation (1) will be stored. |

+ | At this stage the reader may be a little lost: where is the function that advects and diffuses the tracer? We do see the lines | ||

<source lang="c"> | <source lang="c"> | ||

gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, | gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, | ||

- | (FttCellTraverseFunc) gfs_cell_reset, rhs); | + | (FttCellTraverseFunc) copy_v_rhs, par); |

</source> | </source> | ||

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. | 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. | ||

Line 274: |
Line 353: | ||

=== gfs_domain_cell_traverse() module === | === 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. | + | This is an important function which is very frequently used in the code. It allows to perform a sequence of calls of a given function, say ''func'', one call for each cell of the tree. It is used to perform operations local to each cell. |

- | 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. | + | In a sense it is the equivalent of a '''do loop''' over an NxN array in a classical fortran or C code. |

- | Here is an excerpt of its documentation as it appears in the C source code: | + | The arguments of the function are (''cell, data'') where ''cell'' is the current cell in the tree, and ''data'' is the data passed as last argument of gfs_domain_cell_traverse(). |

- | <source lang="c"> | + | In the example above, the function ''func'' is <code> gfs_cell_reset </code> and the ''data'' is <code>rhs</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. | + | |

- | */ | + | |

- | </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. | + | |

+ | We can understand this function better if have a look at the | ||

+ | [http://gerris.dalembert.upmc.fr/gerris/reference Reference Manual] . Click on the [http://gerris.dalembert.upmc.fr/gerris/reference/group__GfsDomain.html#gac8604dda373e0d5f9ad8ada611bf2c39 link to the description of gfs_domain_cell_traverse]. | ||

- | 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. | + | Following the links there we can understand a lot about the function. Let us summarize the most important facts. |

+ | 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. | ||

+ | 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>copy_v_rhs</code> initialize <code> rhs </code> to the value of the current value of <math>\alpha<math> (par-> v). We must find the actual advection in another function at a deeper level. | ||

---- | ---- | ||

Line 317: |
Line 383: | ||

=== variable_sources() module === | === 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 | + | ==== Calculating the advection term ==== |

+ | |||

+ | 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_reset | ||

Line 323: |
Line 393: | ||

* par->flux (a fonction pointer, which will be elucidated below to be gfs_face_advection_flux) | * 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 ==== | + | First we assume that par-> scheme is set to '''godunov''' otherwise the advection term would be absent of the equation. It can be a little messy how Gerris uses the variables rhs and par-> fv. Before <code> variable_source </code> is called par -> fv and rhs are the same and it has been initialized with the values of <math>\alpha</math> at instant ''n''. When we jump into <code> variable_source </code>, par -> fv is created again, it has therefore, let say a new memory address, being its values reseted to zero in the following line. So, the coincidence between rhs (named sv inside the function) and par -> fv has disappeared at this point: in this point, rhs store <math>\alpha<math> at instant n and par -> fv, set to 0 at this point by means of <code> gfs_cell_reset </code>, will be serve to store the advection term. |

- | Again the call to gfs_cell_reset does not do much: it justs resets ''par->fv'' to zero. | + | The first relevant step we observe is: |

+ | <source lang ="c"> | ||

+ | gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1, | ||

+ | (FttCellTraverseFunc) gfs_cell_advected_face_values, par); | ||

+ | </source> | ||

+ | with which are set the values of the tracer <math>\alpha<math> at instant n+1/2 at ALL the cell faces (Remember that <code> gfs_domain_cell_traverse </code> works as do loop to visit the cells applying to each cell the task specified by the selected FttCellTraverseFunc). In effect, let inspect <code> gfs_cell_advected_face_values </code>. | ||

- | ==== gfs_cell_advected_face_values ==== | + | ===== gfs_cell_advected_face_values ===== |

The documentation section above this function says what it does: | The documentation section above this function says what it does: | ||

Line 365: |
Line 439: | ||

gfs_variable_class ()) | gfs_variable_class ()) | ||

</source> | </source> | ||

- | So finally the par->v variable is just the tracer variable we are advecting, which should not be surprising. | + | So finally the par->v variable is just the tracer variable we are advecting as was assumed above. |

+ | |||

Now, ''look at the source of the function''. | Now, ''look at the source of the function''. | ||

This function appears to be rather complex. | This function appears to be rather complex. | ||

Line 384: |
Line 459: | ||

for (c = 0; c < FTT_DIMENSION; c++) { | for (c = 0; c < FTT_DIMENSION; c++) { | ||

gdouble unorm = par->use_centered_velocity ? | gdouble unorm = par->use_centered_velocity ? | ||

- | par->dt*GFS_VARIABLE (cell, par->u[c]->i)/size : | + | par->dt*GFS_VALUE (cell, par->u[c])/size : |

par->dt*(s->f[2*c].un + s->f[2*c + 1].un)/(2.*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 g = (* par->gradient) (cell, c, par->v->i); | ||

- | gdouble vl = GFS_VARIABLE (cell, par->v->i) + MIN ((1. - unorm)/2., 0.5)*g; | + | gdouble vl = GFS_VALUE (cell, par->v) + MIN ((1. - unorm)/2., 0.5)*g; |

- | gdouble vr = GFS_VARIABLE (cell, par->v->i) + MAX ((- 1. - unorm)/2., -0.5)*g; | + | gdouble vr = GFS_VALUE (cell, par->v) + MAX ((- 1. - unorm)/2., -0.5)*g; |

gdouble src = par->dt*gfs_variable_mac_source (par->v, cell)/2.; | gdouble src = par->dt*gfs_variable_mac_source (par->v, cell)/2.; | ||

gdouble dv; | gdouble dv; | ||

Line 409: |
Line 484: | ||

</source> | </source> | ||

[[Image:u_LR.png|Figure 1 (from Dan Martin's PhD thesis)|300px|thumbnail]] | [[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 | + | 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) | static void variable_tracer_init (GfsVariableTracer * v) | ||

Line 415: |
Line 490: | ||

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. | 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 ==== | + | ===== par->flux ===== |

One of the important functions for the implementation of the ''u grad T'' term is | One of the important functions for the implementation of the ''u grad T'' term is | ||

Line 439: |
Line 514: | ||

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>. | 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_advection_flux ==== | + | ===== gfs_face_advection_flux ===== |

If we look in this function we see the comments | If we look in this function we see the comments | ||

Line 457: |
Line 532: | ||

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

<source lang="c"> | <source lang="c"> | ||

- | flux = GFS_FACE_FRACTION (face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt* | + | flux = gfs_domain_face_fraction (par->v->domain, face)*GFS_FACE_NORMAL_VELOCITY (face)*par->dt* |

gfs_face_upwinded_value (face, GFS_FACE_UPWINDING, NULL)/ftt_cell_size (face->cell); | gfs_face_upwinded_value (face, GFS_FACE_UPWINDING, NULL)/ftt_cell_size (face->cell); | ||

</source> | </source> | ||

- | so now we have the hint that | + | so now we have the hint |

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> | ||

+ | Observe that <code>gfs_domain_face_fraction </code> takes into account either the surface fraction, <math>s_f</math> and the metric factor <math>g_f</math>. | ||

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: | 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 ===== | + | ====== gfs_face_upwinded_value ====== |

<source lang="c"> | <source lang="c"> | ||

Line 487: |
Line 563: | ||

</source> | </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. | ||

+ | Let resume the description <code>source_variables</code>. At the moment we have calculated thee advection term which is stored at par->fv. | ||

+ | The next two lines we found in <code>source_variables</code> serves of to add the advection term to the variable <code>rhs</code>; first par-> v is set to sv (rhs) and then to par -> v (which is in fact in that point, sv, or in other words, rhs) is added the content of par-> fv by means of par-> update (This is again a pointer to a function, in this case <code> gfs_advection_update </code>). It is in this function where to the advection flux it is applied the weighting factor <math>1/(a_v g_c)</math> (The function <code> gfs_domain_cell_fraction () </code> returns <math>a_v g_c</math>). | ||

- | 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. | + | ==== Calculating the source term ==== |

+ | |||

+ | It is also in this module where the additional sources are computed and added to rhs. Just a few tips about sources. | ||

+ | |||

+ | If we have sources it seems sensitive that they were derived from the GfsSourceScalar class since they apply to a scalar variable, <math>\alpha<math>. We first make a grep search, grep GfsSourceScalar *, inspecting the lines showed we see that it is defined in source.h. Open it with emacs, and locating its definition we see that: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | typedef struct _GfsSourceScalar GfsSourceScalar; | ||

+ | |||

+ | struct _GfsSourceScalar { | ||

+ | /*< private >*/ | ||

+ | GfsSourceGeneric parent; | ||

+ | |||

+ | /*< public >*/ | ||

+ | GfsVariable * v; | ||

+ | }; | ||

+ | </source> | ||

+ | |||

+ | We see that its parent is GfsSourceGeneric, Navigating to it using tag we see in its definition the following lines: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | ... | ||

+ | gdouble (* mac_value) (GfsSourceGeneric *, FttCell *, GfsVariable *); | ||

+ | gdouble (* centered_value) (GfsSourceGeneric *, FttCell *, GfsVariable *); | ||

+ | gdouble (* face_value) (GfsSourceGeneric *, FttCellFace *, GfsVariable *); | ||

+ | void (* flux) (GfsSourceGeneric *, GfsDomain *, | ||

+ | GfsVariable *, GfsVariable *, | ||

+ | gdouble); | ||

+ | ... | ||

+ | </source> | ||

+ | |||

+ | These lines are the core of the source classes; they are pointers to functions. We have to provide one (or various) functions in the form shown above with which we compute the source term. | ||

+ | |||

+ | We can provide center values of the source (<math>S_\alpha<math>) as is the case of the GfsSource class: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | static void source_init (GfsSourceGeneric * s) | ||

+ | { | ||

+ | s->mac_value = s->centered_value = source_value; | ||

+ | } | ||

+ | </source> | ||

+ | |||

+ | Also, we can provide a function that compute the flux added by means of the flux vector <math>\mathbf{s_\alpha}</math> to the cell by the source in a period <math>\Delta t</math>. In other words, with the nomenclature of above the function should provide the sum | ||

+ | <math> \Delta t \sum_{contour} s_f g_f s_\alpha</math> of the discrete equation. | ||

+ | |||

+ | This is the case of the GfsSourceDiffusionExplicit class. In this case, if we take a look to the function initializing GfsSourceDiffusionExplicit: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | static void gfs_source_diffusion_explicit_init (GfsSourceGeneric * s) | ||

+ | { | ||

+ | s->mac_value = s->centered_value = NULL; | ||

+ | s->flux = source_diffusion_explicit_flux; | ||

+ | } | ||

+ | </source> | ||

+ | |||

+ | |||

+ | It is available other type of sources, GfsSourceVelocity, which suits better to the addition of sources (volume forces) to the momentum equation, the equations for the velocity components. Those interested can take a look to GfsSourceCoriolis or to GfsSourceElectric in the electric module ''electrohydro.mod''. | ||

+ | |||

+ | After this brief explanation of the sources classes, it is clear the role of the two following lines of <code>source_variable </code>: | ||

+ | |||

+ | * source -> flux are computed by means of <code> gfs_domain_variable_fluxes</code> and after weighting its value by <math>a_v g_c</math>, are added to <code>rhs</code>. | ||

+ | * source -> centered_value are computed by means of <code>gfs_domain_variable_centered_sources</code> and then directly added to <code>rhs</code>. | ||

+ | |||

+ | Summarizing, in the function <code> variable_source </code> it is where is computed the right side of the discrete equation except the diffusion term. | ||

+ | |||

+ | === variable_diffusion() module === | ||

+ | |||

+ | Tag navigating we will first take a look to GfsSourceDiffusion, | ||

+ | |||

+ | <source lang ="c"> | ||

+ | struct _GfsSourceDiffusion { | ||

+ | /*< private >*/ | ||

+ | GfsSourceScalar parent; | ||

+ | |||

+ | /*< public >*/ | ||

+ | GfsDiffusion * D; | ||

+ | }; | ||

+ | </source> | ||

+ | It is clear that in essence the GfsSourceDiffusion is a scalar source that contain only a GfsDiffusion object. Is in this object, GfsDiffusion, in where all the relevant information for the integration of the diffusion term is available: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | typedef struct _GfsDiffusion GfsDiffusion; | ||

+ | |||

+ | struct _GfsDiffusion { | ||

+ | /*< private >*/ | ||

+ | GfsEvent parent; | ||

+ | |||

+ | /*< public >*/ | ||

+ | GfsFunction * val; | ||

+ | GfsVariable * mu; | ||

+ | GfsMultilevelParams par; | ||

+ | }; | ||

+ | </source> | ||

+ | |||

+ | Once we have clear how is the structure of GfSourceDiffusion, we are in condition to analyze <code> variable_diffusion </code>. If we tag navigate to it the first thing we observe is the line: | ||

+ | |||

+ | <source lang ="c"> | ||

+ | gfs_diffusion_coefficients (domain, d, par->dt, rhoc, axi, alpha, d->D->par.beta); | ||

+ | </source> | ||

+ | |||

+ | It is self explanatory what this function do: To calculate the diffusion weighting factors <math>\Delta t \: D_f \, g_f \, s_f \, \beta /gm_f</math> at all the faces and solid segments. Additionally, since will be neccessary, the variable rhoc is filled with <math>g_c a_v</math> (this is done inside of <code> diffusion_mixed_coef </code>). | ||

+ | |||

+ | After applying boundary conditions for the tracer <math>\alpha</math> on embedded surfaces by means of <code>gfs_domain_surface_bc </code> the function <code>gfs_diffusion_rhs</code> is called. Within this function the computation of the right hand side of the discrete equation is completed since in it a do-loop over the cells (using again <code>gfs_domain_cell_traverse</code>) is performed computing the summatory of diffusion fluxes over the faces+solid of each cell and dividing the sum by the factor <math>h^2 a_v g_c</math>. Observe that the total flux (the summatory) is pondered with <math>(1-\beta)/\beta</math> to correct the fact that the diffusion coefficients has been calculated with <math>\beta</math>. | ||

+ | |||

+ | <code>gfs_diffusion</code> serves finally to solve the discrete using multilevel relaxation scheme similar to the one described in the Stephane's [http://gfs.sf.net/gerris.pdf Journal of Computational Physics article]. | ||

+ | |||

---- | ---- | ||

- | 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]]. | + | ⇐ '''previous''': [[The Fluid Domain]], ⇑ '''up''': [[Gerris Flow Solver Programming Course for Dummies]], ⇒ '''next''': [[The Gerris Object System]] |

## Revision as of 20:19, 15 October 2012

## 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 to 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-diffusion. We then turn to the detailed analysis of tracer advection and diffusion.

### Equation and numerical Scheme

The present section should is brief guide of how Gerris treats the advection-diffusion equation of a tracer, <math>\alpha</math>,

<math> \frac{\partial \alpha}{\partial t} + \nabla \cdot ( \mathbf{u} \alpha) = \nabla \cdot ( D \nabla \alpha)+S_\alpha + \nabla \cdot \mathbf{s}_\alpha </math>

where <math>\mathbf{u}</math> is the fluid velocity and <math>D<\math> is the diffusion coefficient. We have added for generality two types of sources: <math>S_\alpha</math> stands for typical sources included, while <math>\mathbf{s}_\alpha</math> would be for generic source vector flux.

Gerris adopts a scheme of the form,

<math> \frac{\alpha^{n+1} -\alpha^{n}}{\Delta t} + \nabla \cdot ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) = (1-\beta) \nabla \cdot ( D^n \nabla \alpha^n)+\beta \nabla \cdot ( D^{n+1} \nabla \alpha^{n+1})+S_\alpha^n + \nabla \cdot \mathbf{s}^n_\alpha </math>

where <math>\beta</math> is a factor that has to be between 0.5 and 1.0: <math>\beta=</math>0.5 (the default value) corresponds to the Crank-Nicholson scheme which is second order accurate in time; <math>\beta=</math>1.0 corresponds to the Euler Forward scheme which although is first order accurate in time is more robust against oscillations. The equation is latterly reordered in a form close to a Poisson's equation,

<math> \alpha^{n+1} - \nabla \cdot ( \beta \: \Delta t \: D^{n+1} \nabla \alpha^{n+1}) = \alpha^{n} - \nabla \cdot (\Delta t \: \mathbf{u}^{n+1/2} \alpha^{n+1/2}) + \nabla \cdot ( (1-\beta) \: \Delta t \: D^n \nabla \alpha^n)+S_\alpha^n </math>

Adopting the cell as the control volume, <math>d \Omega</math>, assuming constant value of the variable <math>\alpha</math> in the cell and integrating the above discrete equation over that control volume, it is obtained the following,

- For the temporal term:

<math> \iiint_\Omega \alpha \: d \Omega = a_v \: \alpha \: \: d \ell_1 d \ell_2 d \ell_3 \, = a_v \: \alpha \: g_c \, h^3 . </math> where <math>h</math> is the size of the cell. For generality we have used a general system of orthogonal curvilinear coordinates given by <math>(x_1,x_2,x_3)</math>. If <math>gm_i</math> are the metric factor of each coordinate <math>x_i</math> (<math>d \ell_i=gm_i d x_i</math>), we collect them in a single metric factor <math>g= gm_1 gm_2 gm_3</math>. The subscript <math>c</math> denotes that the metric factor is to be computed at the center of mass of the cell. <math>a_v</math> is the volume fraction (it is always 1. except in the case that the cell is cut by a solid segment, i.e a mixed cell).

- For the advection term:

To this term it is applied additionally the divergence theorem in order to calculate it by means of the advection fluxes through the cell's contour,

<math> \begin{split} \iiint_\Omega \nabla \cdot ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) \: \Delta t \: d \Omega &= \iint_\Sigma \Delta t ( \mathbf{u}^{n+1/2} \alpha^{n+1/2}) \cdot \mathbf{n} \: d \Sigma \\ &= \sum_{Contour} {u_f}^{n+1/2} \alpha_f^{n+1/2} \Delta t \, s_f g_f h^2 \end{split} </math>

where <math>{u_f}</math> and <math>{\alpha_f}</math> are the normal velocity to the face (or solid segment) and the value of the tracer at that face (or solid segment), respectively; <math>s_f</math> is the face fraction (it is always 1. except in the case that the face were cut by a solid segment). <math>g_f</math> is, again, a metric factor denoting the subscript <math>f</math> that is evaluated at the face.

- For the diffusion terms:

These terms are treated similarly to the advection term resulting in, <math> \begin{split} \iiint_\Omega \nabla \cdot (\Delta t\: D \nabla \alpha) \: d \Omega &= \iint_\Sigma \Delta t\: D \nabla \alpha \cdot \mathbf{n} \: d \Sigma \\ &= \sum_{Contour} \Delta t\: D_f g_f s_f h (h \nabla^f /alpha) \end{split} </math>

Reordering the equation, it results finally in the following expression,

<math> \begin{split}

\alpha^{n+1} &- \frac{1}{a_v g_c h^2}\sum_{contour} \frac{\Delta t \: D_f \,g_f \, s_f \, \beta}{gm_f} \left(h \nabla^f \alpha^{n+1} \right) =\\ \alpha^{n} &- \frac{1}{a_v g_c} \sum_{contour} \frac{\Delta t \, s_f g_f}{h} {u_f}^{n+1/2} \alpha_f^{n+1/2} + S_\alpha^n \, \Delta t + \frac{1}{a_v g_c} \sum_{Contour} s_f^n \Delta t \, s_f g_f/h \right \\ &+ \frac{1}{a_v g_c h^2} \sum_{contour} \frac{\Delta t \: D_f g_f \, s_f \, (1-\beta)}{gm_f} \left(h \nabla^f \alpha^{n} \right)_f \end{split}

</math>

where we have take into account, by means of <math>gm_f</math>, that a orthogonal curvilinear coordinates system could affect to the expression of the normal gradient.

### 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 by running Gerris through a debugger. 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. For this we shall use an interesting feature of the text editors.

## Tip: Navigate using the TAGS file

In the above example, we have seen that we can gain some information about the macros, classes, structures and functions of the code by navigating in the code from one definition to the next. This navigation is not always immediately enlightening: we must resort to it only when we have enough general understanding of the code. To perform this navigation in a more efficient way than `grep`

, we can use the following trick:

- Create the TAGS file :

% cd src % make tags

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

- In emacs:
- 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:
- Open any *.c or *.h file with vim.
- To inform vim of the TAGS file, execute the following command:
`set tags+=<PATH_TO_GERRIS_SRC>/TAGS`

, where`<PATH_TO_GERRIS_SRC>`

is the path to the directory containing the Gerris source code (and the TAGS file). Better yet, add the command to your .vimrc file to perform this automatically. - Position the cursor on a function or variable.
- Ctrl-] (in normal/command mode) will open the appropriate source file in a new buffer and navigate to the function's or variable's definition. Your location prior to executing this command is pushed into the "tag stack."
- Ctrl-T "pops" the tag stack; i.e., Ctrl-T returns you to the location prior to the last executed Ctrl-].

You can now try to use this trick on your own editor. It is a good idea to try the following: find any place in the code where the `GFS_VALUE`

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, then repeat the operation to display the declarations and definitions above.

From now on, the reader or listener to the course should open his text editor and navigate freely in the code, trying to find the implementation of the numerical methods. We guide the reader below for one particular scheme: the tracer advection.

## The coding of tracer advection-diffusion

### 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 is 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;

....

if ((d = source_diffusion (par->v))) {

GfsVariable * rhs;

par->fv = rhs = gfs_temporary_variable (domain);

gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,

(FttCellTraverseFunc) copy_v_rhs, par);

variable_sources (domain, par, rhs, NULL, NULL);

variable_diffusion (domain, d, par, rhs, NULL);

gts_object_destroy (GTS\_OBJECT (rhs));

}

else {

.....

}

}

It can be seen that a temporary variable, rhs (or equivalently, par-> fv), is created. It is in this variable, rhs, where the right hand side of the advection-diffusion equation (1) will be stored. At this stage the reader may be a little lost: where is the function that advects and diffuses the tracer? We do see the lines

gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,

(FttCellTraverseFunc) copy_v_rhs, par);

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. It allows to perform a sequence of calls of a given function, say *func*, one call for each cell of the tree. It is used to perform operations local to each cell.
In a sense it is the equivalent of a **do loop** over an NxN array in a classical fortran or C code.
The arguments of the function are (*cell, data*) where *cell* is the current cell in the tree, and *data* is the data passed as last argument of gfs_domain_cell_traverse().
In the example above, the function *func* is ` gfs_cell_reset `

and the *data* is `rhs`

.

We can understand this function better if have a look at the Reference Manual . Click on the link to the description of gfs_domain_cell_traverse.

Following the links there we can understand a lot about the function. Let us summarize the most important facts.
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 `copy_v_rhs`

initialize ` rhs `

to the value of the current value of <math>\alpha<math> (par-> v). 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:

- 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

#### Calculating the advection term

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)

First we assume that par-> scheme is set to **godunov** otherwise the advection term would be absent of the equation. It can be a little messy how Gerris uses the variables rhs and par-> fv. Before ` variable_source `

is called par -> fv and rhs are the same and it has been initialized with the values of <math>\alpha</math> at instant *n*. When we jump into ` variable_source `

, par -> fv is created again, it has therefore, let say a new memory address, being its values reseted to zero in the following line. So, the coincidence between rhs (named sv inside the function) and par -> fv has disappeared at this point: in this point, rhs store <math>\alpha<math> at instant n and par -> fv, set to 0 at this point by means of ` gfs_cell_reset `

, will be serve to store the advection term.

The first relevant step we observe is:

gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,

(FttCellTraverseFunc) gfs_cell_advected_face_values, par);

with which are set the values of the tracer <math>\alpha<math> at instant n+1/2 at ALL the cell faces (Remember that ` gfs_domain_cell_traverse `

works as do loop to visit the cells applying to each cell the task specified by the selected FttCellTraverseFunc). In effect, let inspect ` gfs_cell_advected_face_values `

.

##### 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 as was assumed above.

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_VALUE (cell, par->u[c])/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_VALUE (cell, par->v) + MIN ((1. - unorm)/2., 0.5)*g;

gdouble vr = GFS_VALUE (cell, par->v) + 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;

}

}

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 the *gradient* 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_domain_face_fraction (par->v->domain, 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 the flux has the well-known form

flux = u_face * T_face

Observe that `gfs_domain_face_fraction `

takes into account either the surface fraction, <math>s_f</math> and the metric factor <math>g_f</math>.

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.

Let resume the description `source_variables`

. At the moment we have calculated thee advection term which is stored at par->fv.
The next two lines we found in `source_variables`

serves of to add the advection term to the variable `rhs`

; first par-> v is set to sv (rhs) and then to par -> v (which is in fact in that point, sv, or in other words, rhs) is added the content of par-> fv by means of par-> update (This is again a pointer to a function, in this case ` gfs_advection_update `

). It is in this function where to the advection flux it is applied the weighting factor <math>1/(a_v g_c)</math> (The function ` gfs_domain_cell_fraction () `

returns <math>a_v g_c</math>).

#### Calculating the source term

It is also in this module where the additional sources are computed and added to rhs. Just a few tips about sources.

If we have sources it seems sensitive that they were derived from the GfsSourceScalar class since they apply to a scalar variable, <math>\alpha<math>. We first make a grep search, grep GfsSourceScalar *, inspecting the lines showed we see that it is defined in source.h. Open it with emacs, and locating its definition we see that:

typedef struct _GfsSourceScalar GfsSourceScalar;

struct _GfsSourceScalar {

/*< private >*/

GfsSourceGeneric parent;

/*< public >*/

GfsVariable * v;

};

We see that its parent is GfsSourceGeneric, Navigating to it using tag we see in its definition the following lines:

...

gdouble (* mac_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);

gdouble (* centered_value) (GfsSourceGeneric *, FttCell *, GfsVariable *);

gdouble (* face_value) (GfsSourceGeneric *, FttCellFace *, GfsVariable *);

void (* flux) (GfsSourceGeneric *, GfsDomain *,

GfsVariable *, GfsVariable *,

gdouble);

...

These lines are the core of the source classes; they are pointers to functions. We have to provide one (or various) functions in the form shown above with which we compute the source term.

We can provide center values of the source (<math>S_\alpha<math>) as is the case of the GfsSource class:

static void source_init (GfsSourceGeneric * s)

{

s->mac_value = s->centered_value = source_value;

}

Also, we can provide a function that compute the flux added by means of the flux vector <math>\mathbf{s_\alpha}</math> to the cell by the source in a period <math>\Delta t</math>. In other words, with the nomenclature of above the function should provide the sum <math> \Delta t \sum_{contour} s_f g_f s_\alpha</math> of the discrete equation.

This is the case of the GfsSourceDiffusionExplicit class. In this case, if we take a look to the function initializing GfsSourceDiffusionExplicit:

static void gfs_source_diffusion_explicit_init (GfsSourceGeneric * s)

{

s->mac_value = s->centered_value = NULL;

s->flux = source_diffusion_explicit_flux;

}

It is available other type of sources, GfsSourceVelocity, which suits better to the addition of sources (volume forces) to the momentum equation, the equations for the velocity components. Those interested can take a look to GfsSourceCoriolis or to GfsSourceElectric in the electric module *electrohydro.mod*.

After this brief explanation of the sources classes, it is clear the role of the two following lines of `source_variable `

:

- source -> flux are computed by means of
`gfs_domain_variable_fluxes`

and after weighting its value by <math>a_v g_c</math>, are added to`rhs`

. - source -> centered_value are computed by means of
`gfs_domain_variable_centered_sources`

and then directly added to`rhs`

.

Summarizing, in the function ` variable_source `

it is where is computed the right side of the discrete equation except the diffusion term.

### variable_diffusion() module

Tag navigating we will first take a look to GfsSourceDiffusion,

struct _GfsSourceDiffusion {

/*< private >*/

GfsSourceScalar parent;

/*< public >*/

GfsDiffusion * D;

};

It is clear that in essence the GfsSourceDiffusion is a scalar source that contain only a GfsDiffusion object. Is in this object, GfsDiffusion, in where all the relevant information for the integration of the diffusion term is available:

typedef struct _GfsDiffusion GfsDiffusion;

struct _GfsDiffusion {

/*< private >*/

GfsEvent parent;

/*< public >*/

GfsFunction * val;

GfsVariable * mu;

GfsMultilevelParams par;

};

Once we have clear how is the structure of GfSourceDiffusion, we are in condition to analyze ` variable_diffusion `

. If we tag navigate to it the first thing we observe is the line:

gfs_diffusion_coefficients (domain, d, par->dt, rhoc, axi, alpha, d->D->par.beta);

It is self explanatory what this function do: To calculate the diffusion weighting factors <math>\Delta t \: D_f \, g_f \, s_f \, \beta /gm_f</math> at all the faces and solid segments. Additionally, since will be neccessary, the variable rhoc is filled with <math>g_c a_v</math> (this is done inside of ` diffusion_mixed_coef `

).

After applying boundary conditions for the tracer <math>\alpha</math> on embedded surfaces by means of `gfs_domain_surface_bc `

the function `gfs_diffusion_rhs`

is called. Within this function the computation of the right hand side of the discrete equation is completed since in it a do-loop over the cells (using again `gfs_domain_cell_traverse`

) is performed computing the summatory of diffusion fluxes over the faces+solid of each cell and dividing the sum by the factor <math>h^2 a_v g_c</math>. Observe that the total flux (the summatory) is pondered with <math>(1-\beta)/\beta</math> to correct the fact that the diffusion coefficients has been calculated with <math>\beta</math>.

`gfs_diffusion`

serves finally to solve the discrete using multilevel relaxation scheme similar to the one described in the Stephane's Journal of Computational Physics article.

⇐ **previous**: The Fluid Domain, ⇑ **up**: Gerris Flow Solver Programming Course for Dummies, ⇒ **next**: The Gerris Object System