Basilisk C

“Basilisk C” is the extension of the C programming language used to write code in Basilisk. The goal of this extension is to provide a simple set of extra C-like constructs useful for writing discretisation schemes on Cartesian grids.

Part I covers the basic concepts which are necessary to setup simple simulations using pre-defined solvers.

Part II gives further details and is necessary if you want to setup more complex simulations and/or write your own discretisation schemes.

The description of the language below assumes a basic knowledge of C (e.g. the first four chapters of The C programming language).

Part I: Basic concepts

Fields and stencils

Field constructs are used to store quantities discretised spatially. They can be seen as a generalisation of C arrays. Three types of fields can be defined: scalar, vector and tensor.

Declaration, allocation and deallocation

A new field is declared as an automatic variable like this:

scalar a[];

Alternatively, explicit allocation and deallocation can be done like this:

scalar a = new scalar;
...
delete ({a});

where {a} is a list (with a single element in this example).

Iterators

Field values are manipulated using iterators, for example:

foreach()
  a[] = 1.;

where foreach() iterates over all the stencils defining the discretised field.

Iterators also define several variables (implicitly) such as:

double x, y, z; // coordinates of the center of the stencil
double Δ;   // size of the stencil cell

Stencils

Stencils are used to access field values and their local neighbours. By default Basilisk guarantees consistent field values in a 3x3 neighbourhood (in 2D). This can be represented like this

3x3 stencil

3x3 stencil

Neighbouring values are accessed using the indexing scheme pictured on the figure. For example, computing an approximation of the Laplacian

b=2a

could be written

scalar a[], b[];
...
foreach()
  b[] = (a[1,0] + a[0,1] + a[-1,0] + a[0,-1] - 4.*a[])/sq(Δ);

where a[] is a shortcut for a[0,0].

Vectors and tensors

Vector and tensor fields are used in a similar way. Vector fields are a collection of D scalar fields (where D is the dimension of the spatial discretisation) and tensor fields are a collection of D vector fields. Each of the components of the vector or tensor fields are accessed using the x, y or z field of the corresponding structure. For example:

vector v[];
tensor t = new tensor;
...
foreach() {
  v.x[] = 1.;
  t.x.x[] = (v.x[1,0] - v.x[-1,0])/Δ;
  t.y.x[] = (v.y[1,0] - v.y[-1,0])/Δ;
}
...
delete ({t});

Boundary conditions

Stencils located close enough to the boundaries of the domain will extend beyond it. The stencil values outside the domain (often called ghost values) need to be initialised. These values can be set in order to provide the discrete equivalents of various boundary conditions.

The default boundary condition is symmetry and takes into account whether the fields are scalars, vectors or tensors. If we consider a boundary, a scalar field a and a vector field v, symmetry can be expressed on the boundary as:

na=0 nvt=0 vn=0

with n and t the normal and tangential directions to the boundary, respectively.

Ghost values usually depend on the values inside the domain (for example when using symmetry conditions). It is thus necessary to update them when values inside the domain are modified and before any operation on stencils. This can be done by calling

boundary ({a,v});

where {a,v} is the list of fields for which boundary conditions need to be updated. If boundary conditions need to be applied to all fields, one can use the predefined list all.

A more realistic version of the tensor initialisation example above would thus be:

vector v[];
tensor t[];

foreach() {
  v.x[] = f(x,y);
  v.y[] = g(x,y);
}
boundary ({v});

foreach() {
  t.x.x[] = (v.x[1,0] - v.x[-1,0])/Δ;
  t.y.x[] = (v.y[1,0] - v.y[-1,0])/Δ;
}
boundary ({t});

Boundary conditions can be changed for each scalar field using the following syntax:

a[top] = a[];

where a[top] is the ghost value of the scalar field a immediately outside the top (respectively bottom, right, left) boundary. In this example we have set a symmetry condition on scalar field a on the top boundary (i.e. the default boundary condition). This corresponds to a Neumann condition (i.e. a condition on the normal derivative of field a). To set a Dirichlet condition instead (i.e. a condition on the value of the function on the boundary), one could use for example:

a[left] = 2.*ab - a[];

where ab is the Dirichlet value on the (left) boundary.

The library also provides functions for common boundary conditions. For example the expression for the Dirichlet condition above could also be written

a[left] = dirichlet(ab);

and the Neumann condition

na=an

would be

a[left] = neumann(an);

(the equivalent discrete expression is left as an exercise).

Note also that using these pre-defined expressions is necessary to obtain automatic homogeneous boundary conditions, which is important if scalar a is the solution of a Poisson problem.

For vector fields, boundary conditions are defined in a coordinate system local to the boundary where the x and y components are replaced by the normal n and tangential t components i.e. imposing no-flux of a vector v through the top and left boundary, together with a no-slip boundary condition would be written

v.n[top] = dirichlet(0);
v.t[top] = dirichlet(0);
v.n[left] = dirichlet(0);
v.t[left] = dirichlet(0);

In three dimensions, there is one normal direction n and two tangential directions t and r, the +z and -z boundaries are called front and back. Taking the left (-x) or right (+x) boundary as an example, we have the correspondence

u.n -> u.x
u.t -> u.y
u.r -> u.z

The relations for the other boundaries (top/bottom, front/back) are obtained by rotation of the indices.

Periodic boundaries

Periodic boundary conditions can be imposed on the right/left and top/bottom boundaries using for example

int main()
{
  ...
  periodic (right); 
  ...
}

All existing fields and all the fields allocated after the call will be periodic in the right/left direction.

Boundary conditions on specific fields can still be set to something else. For example, one could use

int main()
{
  ...
  periodic (right);
  p[left]  = dirichlet(0);
  p[right] = dirichlet(1);
  ...
}

to impose a pressure gradient onto an otherwise periodic domain.

Conversely it is possible to impose periodic conditions only on specific fields using e.g

p[top] = periodic();

Events

Numerical simulations often need to perform actions (outputs for example) at given time intervals. Because the timestep used to integrate the numerical scheme can vary, for example due to stability requirements, it is generally not trivial to ensure that specific time intervals will be respected. To solve this problem Basilisk C provides events.

The overall syntax of events is similar to that of for loops in C. For example:

event name (t = 1; t <= 5; t += 1) {
  ...
}

where name is the user-defined name of the event, t = 1 specifies the starting time, t <= 5 is the condition which must be verified for the event to carry on and t += 1 is the iteration operator. The event can happen either at specified times t or at a specified number of timesteps i.

Beyond this similarity with for loops, the syntax for events is more flexible. The order of the initialisation, condition and iteration statements can change and/or some of the statements can be omitted. For example, these are all valid events:

event name (t = 1; t += 1; t <= 5) {
  // same as previously
  ...
}

event example1 (t = 2) {
  // executes once at t = 2
  ...
}

event example2 (i++) {
  // executes at every timestep
  ...
}

event example3 (t = 1; t *= 2) {
  // executes at t = 1,2,4,8,16,...
  ...
}

In addition, lists can be used to set specific times/timesteps:

event example4 (t = {1.2, 3, 7.6}) {
  ...
}

and the special end keyword can be used for events which should happen after completion of the simulation:

event example5 (t = end) {
  ...
}

It is also possible to stop a simulation based on the return value from an event. For example:

event stop (i += 10) {
  if (my_stopping_condition())
    return 1; // any value different from zero will stop the simulation
  // the default is to return zero
}

event finalize (t = end) {
  // this will happen after the simulation has completed due to 
  // the "stop" event above
}

Part II: More advanced concepts

Lists

It is often useful to perform the same operation on several fields. Basilisk C provides a simple extension which allows iterations over lists of fields.

Declaration and allocation

An automatic list (of scalars) can be declared and allocated like this:

scalar * list = {a,b,c,d};

Automatic type casting

Lists can combine elements of different types (e.g. scalar fields and vector fields). In this case the type of the resulting list is always that of the element with the lowest dimension. For example, these are valid list declarations:

scalar a[];
vector v[];
tensor t[];
...
scalar * list1 = {a,v,t};
vector * list2 = {v,t};

however trying

vector * list3 = {a,v,t};
tensor * list4 = {v,t};

will give a compilation error.

Explicit type casting

In a similar way it is sometimes useful to be able to explictly cast lists of different types. For example a list of vector fields can be cast into a list of scalar fields using:

vector * l1 = {a,b,c};
scalar * l2 = (scalar *) l1;

Explicit type casting is only allowed from a type of higher dimension down to a type of lower dimension (i.e. tensors to vectors and vectors to scalars).

List iterators

To iterate over all the elements of a list use

scalar * list = {a,b,c,d};
...
for (scalar s in list)
  dosomething (s);

It is also possible to iterate over two lists (of the same length) simultaneously using e.g.:

scalar * l1 = ...;
vector * l2 = ...;
...
scalar s; vector v;
for (s,v in l1,l2)
  dosomething (s, v);

Finally note that it is much more efficient to do:

foreach()
  for (scalar s in list)
    s[] = ...;

rather than

for (scalar s in list)
  foreach()
    s[] = ...;

because fields values for each scalar are usually stored close to one another (this is similar to ordering loops properly when accessing multidimensional arrays in C or Fortran).

Face and vertex fields

Discretisation schemes often rely on special arrangements of discretisation variables relative to the underlying grid (this is sometimes called variable staggering). Basilisk provides support for the three most common types of staggering: centered (the default), face and vertex staggering.

Example of centered, face and vertex staggering

Example of centered, face and vertex staggering

How to implement such staggering is largely a matter of convention. For example, one could just allocate standard fields in Basilisk and take as convention that the corresponding variables are staggered. While this would work in simple cases, some operations performed by Basilisk (such as interpolation and boundary conditions) need to know that these fields are staggered. Basilisk C introduces two special keywords for this purpose: face and vertex. For example, the fields in the example above would be declared

scalar p[];
face vector u[];
vertex scalar ω[];

Boundary conditions and stencils

Boundary conditions and stencils are necessarily different when considering cell-centered and staggered variables.

For vertex fields, Basilisk does not apply any boundary conditions. All the vertex values need to be defined using foreach_vertex().

For face fields, a standard call to boundary() will ensure that the stencil values in the picture below are consistent.

Stencils for face and vertex fields

Stencils for face and vertex fields

In some cases, stencil values are not required in the direction tangential to the face vector field (i.e. all values other than the normal components u.x[], u.y[], u.x[1,0] and u.y[0,1]). It is then only necessary to call

boundary_flux ({u});

to ensure consistency of the normal components.

Complex domains

By default the simulation domain is a square box with right, left, top and bottom associated boundary conditions. It is possible to define domains of arbitrary shape, with an arbitrary number of associated boundary conditions, using the mask() function. This function associates a given boundary condition to each cell of the grid.

For example, if we start with the initial square domain with 0x1 and 0y1, and want to turn the domain into a rectangle with 0y0.5, we could use

mask (y > 0.5 ? top : none);

The argument of the function is the value of the boundary condition to assign to each cell. In this example, all grid points for which y>0.5 will be assigned the (pre-defined) top boundary condition, the boundary condition of all other grid points will be unchanged (the none value is just ignored).

New types of boundary conditions can be declared using e.g.

bid circle;

where circle is a user-defined identifier. For example, a circular boundary with a no-slip boundary condition for a vector field u could be defined using

bid circle;
u.t[circle] = dirichlet(0);
...
mask (sq(x - 0.5) + sq(y - 0.5) < sq(0.5) ? circle : none);

All other boundaries would retain there (default) symmetry (i.e. slip) boundary conditions for u.

Field attributes

Scalar field can have associated “attributes” which behave very much like “members” of C structures. To add a new attribute (to all scalar fields) use something like

attribute {
  int myattribute;
}

You can then access this attribute using

scalar s[];
s.myattribute = 2;

Of course you are not limited to int, you can use any declaration which would be legal in a C structure, provided the new name does not conflict with the names of existing attributes.

If the attribute is not set explicitly, it defaults to “zero”.

Pre-defined attributes

Point function calls

More field iterators

We have already seen the basic field iterator foreach(). While many numerical schemes can be implemented using only this iterator, more efficient and elegant implementations can often be designed using other iterators provided by Basilisk.

foreach_dimension()

As an example, consider the code required to compute

g=a

Using a centered scheme, this could be written:

scalar a[];
vector g[];
...
foreach() {
  g.x[] = (a[1,0] - a[-1,0])/(2.*Δ);
  g.y[] = (a[0,1] - a[0,-1])/(2.*Δ);
}

It is clear that computing g.y is almost identical to computing g.x, all that is required is a permutation of the indices. Most multi-dimensional schemes involve similar permutations. Writing these permutations manually, as we did above, is error-prone and leads to code duplication. While this may be acceptable for the simple example above, this can become very heavy for more complicated schemes. To avoid this Basilisk provides the foreach_dimension() iterator which can be used like this:

...
foreach()
  foreach_dimension()
    g.x[] = (a[1,0] - a[-1,0])/(2.*Δ);

At compilation, the code block associated with foreach_dimension() (a single line in this example) will be automatically duplicated for each dimension of the problem with the appropriate coordinate permutations: .x will become .y, .y will become .x (or .z in 3D) etc… and the field indices will be rotated correspondingly.

Any variable or function identifier ending with _x, _y and _z will also be rotated.

foreach_face()

The foreach() operator can be seen as iterating over the cells of the grid. When writing flux-based schemes it is often useful to consider iterations over the faces of each cell (because fluxes are by definition associated with cell faces).

Implementing an iteration over cell faces using only the foreach() operator is not easy, even for a one-dimensional problem, due to boundary conditions. Things also become more complex when considering multi-dimensional and adaptive grids. The foreach_face() field iterator solves these problems by providing a consistent traversal of cell faces in all cases.

For example, computing the components of a normal to each face can be written:

scalar a[];
face vector g[];
...
foreach_face()
  g.x[] = (a[] - a[-1,0])/Δ;

Note that a face is defined by convention as separating a cell from its left neighbour (hence a[-1,0] and not a[1,0], see the Face and vertex fields section above). Note also that the duplication/permutation performed by foreach_dimension() is implicitly included in foreach_face().

In some cases, it can be necessary to apply different operations to each component of a face vector field. For example let’s assume we need to initialise a face vector field with the components (y,x). This could be done using

face vector u[];
...
foreach_face(x)
  u.x[] = y;
foreach_face(y)
  u.y[] = x;

Note that the coordinates x and y correspond to the center of the face.

foreach_vertex()

When initialising vertex fields, it is necessary to traverse each vertex of the grid. For example, computing vorticity using the components of velocity defined on faces could be written

face vector u[];
vertex scalar ω[];
...
foreach_vertex()
  ω[] = (u.y[] - u.y[-1,0] - u.x[] + u.x[0,-1])/Δ;

Consistently with the other iterators, the x and y declared implicitly within foreach_vertex() loops are the coordinates of the vertex.

foreach_boundary()

It is sometimes useful to be able to traverse only the cells which are touching a given boundary, for example to compute surface fluxes and other diagnostics. While using a combination of foreach() and conditions on cell coordinates would be possible, it would also be inefficient and could be difficult when using adaptive meshes. For this reason, Basilisk provides the foreach_boundary() iterator which can be used like this for example:

scalar a[];
...
double μ = 1., flux = 0.;
foreach_boundary (top)
  flux += sq(Δ)*μ*(a[ghost] - a[])/Δ;

where top indicates which boundary we want to traverse (the keywords are the same as used for boundary conditions) and ghost is the index of the corresponding ghost cell: for this example a[ghost] == a[0,1].

This example computes the (diffusive) flux of a through the top boundary. Note that for this to work in parallel you would also need to add a reduction clause on flux in foreach_boundary().

foreach_neighbor()

This iterator is local to a cell, so it must be included within a global field iterator (e.g. foreach() etc.). It will iterate over all the cells belonging to the local stencil i.e. the 5x5 stencil by default (in 2D). This can optionally be reduced to a 3x3 stencil (i.e. only the nearest neighbors) using the syntax foreach_neighbor(1).

Event inheritance

In contrast with C functions, multiple events can share the same name. This is used to implement inheritance of existing events and allows to modify and extend the functionality of existing solvers.

When two (or more) events have the same name, the order of execution of this group is modified. The following rules are applied:

  • all events in the group are executed “together”,
  • the first event in the group (in order of apparition in the source files), sets the point where the group is executed,
  • within a group, the order of execution is the inverse of the order of apparition in the source files.

Although this may seem complicated, these rules work well in practice and lead to a reasonably simple inheritance mechanism.

A trace of the order of execution of events can be obtained using the compilation option:

% qcc -events ...

The resulting code will then output (on standard error) a trace looking like:

...
events (i = 4, t = 0.00191154)
  logfile                   rising.c:131
  set_dtmax                 src/navier-stokes/centered.h:163
  stability                 src/tension.h:59
  stability                 src/vof.h:44
  stability                 src/navier-stokes/centered.h:165
  vof                       src/vof.h:174
  vof                       src/navier-stokes/centered.h:175
  tracer_advection          src/navier-stokes/centered.h:176
...

There are two event groups in this example: stability and vof. They correspond to extensions of the initial solver /src/navier-stokes/centered.h: adding the stability conditions for surface tension (at line 59 of /src/tension.h) and for VOF (at line 44 of /src/vof.h) and plugging in VOF advection (at line 174 of /src/vof.h) at the correct location (line 175 of /src/navier-stokes/centered.h).

Grid allocation and deallocation

Hierarchy of grids

Multigrids

Adaptive grids

Metric

Named/optional arguments in function calls

C++ allows named and optional arguments in function calls, for example:

void example (float f = 0, double d = 0);
...
example (f = 2); // d = 0
example (d = 4, f = 1.2);

Unfortunately this is not allowed in C99. Basilisk C includes an extension which provides this functionality (using structure initialisation in C99). Using this extension, the C++ example above would be written:

struct Example {
  float f;
  double d;
};

void example (struct Example p)
{
  float f = p.f;
  double d = p.d;
  ...
}
...
example (f = 2); // d = 0
example (d = 4, f = 1.2);

Parallel programming

Basilisk can automatically parallelise field iterations (i.e. foreach() etc…) on systems supporting OpenMP. If the write operations performed on stencils are purely local (i.e. concurrent accesses by several threads are only possible for read operations), then nothing special needs to be done to parallelise the corresponding field iteration.

A classical example where this condition is not verified are reduction operations. Consider for example the function:

double maximum (scalar a) {
  double maxi = - HUGE;
  foreach()
    if (fabs(a[]) > maxi)
      maxi = fabs(a[]);
  return maxi;
}

In this example, if the foreach() iteration is performed in parallel, several threads may try to set the maxi variable simultaneously (this is a concurrent write access). If this happens the result returned by the function will be incorrect.

To cover this common pattern, OpenMP provides reduction operators which can be used in Basilisk like this:

double maximum (scalar a) {
  double maxi = - 1e100;
  foreach (reduction(max:maxi))
    if (fabs(a[]) > maxi)
      maxi = fabs(a[]);
  return maxi;
}

Tracing and profiling

The trace keyword is used to indicate that the function immediately following should be instrumented for built-in profiling or profiling with Paraver. For example one could write

trace
double maximum (scalar a) {
...