"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](Basilisk C#part-i-basic-concepts) covers the basic concepts
which are necessary to setup simple simulations using [pre-defined
solvers](src/README#solvers).
[Part II](Basilisk C#part-ii-more-advanced-concepts) 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*](http://en.wikipedia.org/wiki/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](http://en.wikipedia.org/wiki/Automatic_variable) like this:
~~~literatec
scalar a[];
~~~
Alternatively, explicit allocation and deallocation can be done like this:
~~~literatec
scalar a = new scalar;
...
delete ({a});
~~~
where `{a}` is a [list](#lists) (with a single element in this example).
### Iterators
Field values are manipulated using *iterators*, for example:
~~~literatec
foreach()
a[] = 1.;
~~~
where `foreach()` iterates over all the *stencils* defining the
discretised field.
Iterators also define several variables (implicitly) such as:
~~~literatec
double x, y, z; // coordinates of the center of the stencil
double Delta; // 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](/src/figures/stencil.svg)
Neighbouring values are accessed using the indexing scheme pictured on
the figure. For example, computing an approximation of the Laplacian
$$b=\nabla^2 a$$
could be written
~~~literatec
scalar a[], b[];
...
foreach()
b[] = (a[1,0] + a[0,1] + a[-1,0] + a[0,-1] - 4.*a[])/sq(Delta);
~~~
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:
~~~literatec
vector v[];
tensor t = new tensor;
...
foreach() {
v.x[] = 1.;
t.x.x[] = (v.x[1,0] - v.x[-1,0])/Delta;
t.y.x[] = (v.y[1,0] - v.y[-1,0])/Delta;
}
...
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:
$$\partial_n a=0$$
$$\partial_n v_t=0$$
$$v_n=0$$
with `n` and `t` the normal and tangential directions to the boundary,
respectively.
Boundary conditions can be changed for each scalar field using the
following syntax:
~~~literatec
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:
~~~literatec
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
~~~literatec
a[left] = dirichlet(ab);
~~~
and the Neumann condition
$$
\partial_na=a_n
$$
would be
~~~literatec
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*](src/poisson.h#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
~~~literatec
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, top/bottom
and front/back boundaries using for example
~~~literatec
int main()
{
...
periodic (right);
...
}
~~~
All existing fields and all the fields allocated after the call will
be periodic in the right/left direction.
## 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:
~~~literatec
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:
~~~literatec
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:
~~~literatec
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:
~~~literatec
event example5 (t = end) {
...
}
~~~
It is also possible to stop a simulation based on the return value
from an event. For example:
~~~literatec
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:
~~~literatec
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:
~~~literatec
scalar a[];
vector v[];
tensor t[];
...
scalar * list1 = {a,v,t};
vector * list2 = {v,t};
~~~
however trying
~~~literatec
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:
~~~literatec
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
~~~literatec
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.:
~~~literatec
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:
~~~literatec
foreach()
for (scalar s in list)
s[] = ...;
~~~
rather than
~~~literatec
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](src/figures/staggering.svg)
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
~~~literatec
scalar p[];
face vector u[];
vertex scalar omega[];
~~~
### 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()*](Basilisk C#foreach_vertex).
For *face* fields stencil values in the picture below are consistent.
![Stencils for *face* and *vertex* fields](src/figures/stencil_face.svg)
## Layers
For some problems for which the aspect ratio (vertical/horizontal) is small --
like the ocean or the atmosphere -- it might be useful to use a custom vertical
coordinate. To use this extension, one must define
~~~literatec
#define LAYERS 1
~~~
in the header. The number of layers is set by the variable nl (by default
`nl = 1` ). You can then declare layered variables via
~~~literatec
scalar h = new scalar[nl];
vector u = new vector[nl];
face vector hu = new face vector[nl];
...
~~~
Such layered variables correspond to `nl` stacked 1D or 2D fields. You need to
cleanup these fields at the end of the execution because they are not automatic
variables:
~~~literatec
event cleanup (t = end, last)
{
delete ((scalar *){h, u, hu});
}
~~~
This layered module comes with an iterator `foreach_layer()` which must be
included within a global field iterator (e.g. `foreach()` etc.)
~~~literatec
foreach()
foreach_layer()
a[] = 1.;
~~~
Within a `foreach()` iterator, you can also access a
specific layer with the last index in the bracket as shown below for layer 10 on a 2D grid.
~~~literatec
foreach()
a[0,0,10] = 1.;
~~~
## 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
~~~literatec
attribute {
int myattribute;
}
~~~
You can then access this attribute using
~~~literatec
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=\nabla a$$
Using a centered scheme, this could be written:
~~~literatec
scalar a[];
vector g[];
...
foreach() {
g.x[] = (a[1,0] - a[-1,0])/(2.*Delta);
g.y[] = (a[0,1] - a[0,-1])/(2.*Delta);
}
~~~
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:
~~~literatec
...
foreach()
foreach_dimension()
g.x[] = (a[1,0] - a[-1,0])/(2.*Delta);
~~~
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.
The `right`, `left`, `top`, `bottom`, `front` and `back` direction
indicators 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 $\nabla a$ normal to each face
can be written:
~~~literatec
scalar a[];
face vector g[];
...
foreach_face()
g.x[] = (a[] - a[-1,0])/Delta;
~~~
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](Basilisk C#face-and-vertex-fields)
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
~~~literatec
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](Basilisk C#face-and-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
~~~literatec
face vector u[];
vertex scalar omega[];
...
foreach_vertex()
omega[] = (u.y[] - u.y[-1,0] - u.x[] + u.x[0,-1])/Delta;
~~~
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:
~~~literatec
scalar a[];
...
double mu = 1., flux = 0.;
foreach_boundary (top)
flux += sq(Delta)*mu*(a[ghost] - a[])/Delta;
~~~
where *top* indicates which boundary we want to traverse (the keywords
are the same as used for [boundary conditions](#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](#parallel-programming) on `flux` in
foreach_boundary().
`foreach_boundary()` iterates over the cells which are
touching a given boundary but the coordinates `x,y,z` correspond not the center of the cells but the coordinates of a given 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](https://en.wikipedia.org/wiki/Inheritance_%28object-oriented_programming%29)
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 appearance 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 appearance 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:
~~~bash
% qcc -events ...
~~~
The resulting code will then output (on standard error) a trace looking like:
~~~bash
...
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
## Diagonalization
The `diagonalize(s)` operator replaces the central stencil value of
the scalar field `s` (i.e. `s[]`) with one and its other values with
zero, in all operations within the following C statement. This is
typically used to automatically construct Jacobi relaxation operators
for multigrid solutions of linear systems.
## Named/optional arguments in function calls
C++ allows named and optional arguments in function calls, for example:
~~~c
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:
~~~literatec
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](https://en.wikipedia.org/wiki/OpenMP) or
[MPI](https://en.wikipedia.org/wiki/Message_Passing_Interface). 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:
~~~literatec
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.
Similarly, when using MPI parallelism, each process will only hold a local
`maxi` value, whereas a global value is required.
To cover this common pattern, OpenMP provides [*reduction
operators*](https://www.openmp.org/spec-html/5.0/openmpsu107.html)
which can be used in Basilisk like this:
~~~literatec
double maximum (scalar a) {
double maxi = - 1e100;
foreach (reduction(max:maxi))
if (fabs(a[]) > maxi)
maxi = fabs(a[]);
return maxi;
}
~~~
It is also possible to reduce an array by indicating its length. For example, one could count the number of leafs cells at each level by writing:
~~~literatec
int len = grid->maxdepth + 1, cells[len] = {0};
foreach (reduction(+:cells[:len]))
cells[level]++;
~~~
Note that Basilisk generalises reduction operators so that they work
also with MPI parallelism.
Beyond reduction operations, care must be taken to avoid concurrent
write accesses when writing Basilik code. Consider the following code
which makes use of an auxilliary variable to store an intermediate
result:
~~~literatec
double norm;
foreach() {
norm = sqrt(sq(u.x[]) + sq(u.y[]));
a[] = norm;
b[] = norm/2.;
}
~~~
The code will work fine in serial (or with MPI), however it will fail
when using shared-memory parallelism (i.e. OpenMP), because `norm`
will be shared amongst threads and written to concurrently. The
solution is simple: just make `norm` a variable local to the loop body
i.e.
~~~literatec
foreach() {
double norm = sqrt(sq(u.x[]) + sq(u.y[]));
a[] = norm;
b[] = norm/2.;
}
~~~
Each thread then has its own private `norm` variable and concurrent
accesses are avoided. Note that this is also a cleaner/clearer way of
writing code, even for the serial version. Generally speaking,
variables should be declared/allocated as closely as possible to where
they are used.
By default Basilisk will parallelise all iterators. This can be a
problem, for example when writing data to a file. Serial iterations
for individual loops can be forced like this:
~~~literatec
foreach (serial)
printf ("%g %g %g\n", x, y, a[]);
~~~
## Tracing and profiling
The `trace` keyword is used to indicate that the function immediately
following should be instrumented for [built-in
profiling](src/README.trace) or [profiling with
Paraver](src/README.paraver). For example one could write
~~~literatec
trace
double maximum (scalar a) {
...
~~~