Basilisk source code (http://basilisk.fr/src/)

root / src / README.page

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
# Solvers

## [Saint-Venant](saint-venant.h)
$$
  \partial_t \int_{\Omega} \mathbf{q} d \Omega =
  \int_{\partial \Omega} \mathbf{f} (
  \mathbf{q}) \cdot \mathbf{n}d \partial
  \Omega - \int_{\Omega} hg \nabla z_b
$$
$$
  \mathbf{q} = \left(\begin{array}{c}
    h\\
    hu_x\\
    hu_y
  \end{array}\right), 
  \;\;\;\;\;\;
  \mathbf{f} (\mathbf{q}) = \left(\begin{array}{cc}
    hu_x & hu_y\\
    hu_x^2 + \frac{1}{2} gh^2 & hu_xu_y\\
    hu_xu_y & hu_y^2 + \frac{1}{2} gh^2
  \end{array}\right)
$$

* [Semi-implicit scheme](saint-venant-implicit.h)
* [Multiple layers](multilayer.h)
$$
\partial_th + \partial_x\sum_{l=0}^{nl-1}h_lu_l = 0
$$

$$
\partial_t(h\mathbf{u}_l) + \nabla\cdot\left(h\mathbf{u}_l\otimes\mathbf{u}_l + 
\frac{gh^2}{2}\mathbf{I}\right) = 
- gh\nabla z_b - \partial_z(h\mathbf{u}w) + \nu h\partial_{z^2}\mathbf{u}
$$

* [Green-Naghdi](green-naghdi.h)
$$
  \partial_t \int_{\Omega} \mathbf{q} d \Omega =
  \int_{\partial \Omega} \mathbf{f} (
  \mathbf{q}) \cdot \mathbf{n}d \partial
  \Omega - \int_{\Omega} hg \nabla z_b +
  h \left( \frac{g}{\alpha}\nabla \eta - D \right)
$$
$$
\alpha h\mathcal{T} \left( D \right) + hD = b
$$
$$
b = \left[ \frac{g}{\alpha} \nabla \eta +\mathcal{Q}_1 \left( u \right)
\right]
$$

## [Systems of conservation laws](conservation.h)
$$
\partial_t\left(\begin{array}{c}
    s_i\\
    \mathbf{v}_j\\
 \end{array}\right) + \nabla\cdot\left(\begin{array}{c}
    \mathbf{F}_i\\
    \mathbf{T}_j\\
 \end{array}\right) = 0
$$

* [Compressible gas dynamics](compressible.h)
* [All Mach compressible flows](all-mach.h)

$$
\partial_t\mathbf{q} + \nabla\cdot(\mathbf{q}\mathbf{u}) = 
- \nabla p + \nabla\cdot(\mu\nabla\mathbf{u}) + \rho\mathbf{a}
$$
$$
\partial_t p + \mathbf{u}\cdot\nabla p = -\rho c^2\nabla\cdot\mathbf{u}
$$

## Navier--Stokes
* [Streamfunction--Vorticity formulation](navier-stokes/stream.h)
$$
\partial_t\omega + \mathbf{u}\cdot\nabla\omega = \nu\nabla^2\omega
$$
$$
\nabla^2\psi = \omega
$$

* ["Markers-And-Cells" (MAC or "C-grid") formulation](navier-stokes/mac.h)
* [Centered formulation](navier-stokes/centered.h)
$$
\partial_t\mathbf{u}+\nabla\cdot(\mathbf{u}\otimes\mathbf{u}) = 
\frac{1}{\rho}\left[-\nabla p + \nabla\cdot(2\mu\mathbf{D})\right] + \mathbf{a}
$$
$$
\nabla\cdot\mathbf{u} = 0
$$

* [Two-phase interfacial flows](two-phase.h)
    + [Momentum conservation option](navier-stokes/conserving.h)
* [Momentum-conserving two-phase interfacial flows](momentum.h)

## Electrohydrodynamics

* [Ohmic conduction](ehd/implicit.h)
$$
\partial_t\rho_e = \nabla \cdot(K \nabla \phi)
$$
$$
\nabla \cdot (\epsilon \nabla \phi) = - \rho_e
$$
* [Ohmic conduction of charged species](ehd/pnp.h)
$$
\partial_tc_i = \nabla \cdot( K_i c_i \nabla \phi)
$$
* [Electrohydrodynamic stresses](ehd/stress.h)
$$
M_{ij} = \varepsilon (E_i E_j - \frac{E^2}{2}\delta_{ij})
$$

## Viscoelasticity

* [The log-conformation method for some viscoelastic constitutive 
models](log-conform.h)

$$
\rho\left[\partial_t\mathbf{u}+\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\right] = 
- \nabla p + \nabla\cdot(2\mu_s\mathbf{D}) + \nabla\cdot\mathbf{\tau}_p
+ \rho\mathbf{a}
$$
$$
\mathbf{\tau}_p = \frac{\mu_p \mathbf{f_s}(\mathbf{A})}{\lambda}
$$
$$
\Psi = \log\mathbf{A}
$$
$$
D_t \Psi = (\Omega \cdot \Psi -\Psi \cdot \Omega) + 2 \mathbf{B} +
\frac{e^{-\Psi} \mathbf{f}_r (e^{\Psi})}{\lambda}
$$

## Other equations
* [Hele-Shaw/Darcy flows](hele-shaw.h)
$$
\mathbf{u} = \beta\nabla p
$$
$$
\nabla\cdot(\beta\nabla p) = \zeta
$$
* [Advection](advection.h)
$$
\partial_tf_i+\mathbf{u}\cdot\nabla f_i=0
$$
    + [Volume-Of-Fluid advection](vof.h)
* [Interfacial forces](iforce.h)
$$
\phi\mathbf{n}\delta_s
$$
    + [Surface tension](tension.h)
    + [Reduced gravity](reduced.h)
* [Reaction--Diffusion](diffusion.h)
$$ 
\theta\partial_tf = \nabla\cdot(D\nabla f) + \beta f + r
$$ 
* [Poisson--Helmholtz](poisson.h)
$$
\nabla\cdot (\alpha\nabla a) + \lambda a = b
$$
* [Runge--Kutta time integrators](runge-kutta.h)
$$
\frac{\partial\mathbf{u}}{\partial t} = L(\mathbf{u}, t)
$$
* [Signed distance field](distance.h)

# General orthogonal coordinates

When not written in vector form, some of the equations above will
change depending on the choice of coordinate system (e.g. polar rather
than Cartesian coordinates). In addition, extra terms can appear due
to the geometric curvature of space (e.g. equations on the sphere). An
important simplification is to consider only [orthogonal
coordinates](http://en.wikipedia.org/wiki/Orthogonal_coordinates). In
this case, consistent finite-volume discretisations of standard
operators (divergence etc...) can be obtained, for any orthogonal
curvilinear coordinate system, using only a few additional geometric
parameters.

![Metric scale factors](figures/metric.svg)

The [face vector](/Basilisk C#face-and-vertex-fields) *fm* is the
*scale factor* for the length of a face i.e. the physical length is
$fm\Delta$ and the scalar field *cm* is the scale factor for the area
of the cell i.e. the physical area is $cm\Delta^2$. By default, these
fields are constant and unity (i.e. the Cartesian metric).

Several metric spaces/coordinate systems are predefined:

* [Axisymmetric](axi.h)
* [Spherical](spherical.h)
* [Radial/cylindrical](radial.h)

# Data processing

* [Various utility functions](utils.h): timing, field statistics,
  slope limiters, etc.
* [Tagging connected neighborhoods](tag.h)
* [Counting droplets](examples/atomisation.c#counting-droplets)

# Output functions

* [Multiple fields interpolated on a regular grid (text format)](output.h#output_field-multiple-fields-interpolated-on-a-regular-grid-text-format)
* [Single field interpolated on a regular grid (binary format)](output.h#output_matrix-single-field-interpolated-on-a-regular-grid-binary-format)
* [Portable PixMap (PPM) image output](output.h#output_ppm-portable-pixmap-ppm-image-output)
* [Volume-Of-Fluid facets](fractions.h#interface-output)
* [Basilisk snapshots](output.h#dump-basilisk-snapshots)
* [Basilisk View](view.h)
* [Gerris simulation format](output.h#output_gfs-gerris-simulation-format)
* [ESRI ASCII Grid format](output.h#output_grd-esri-ascii-grid-format)
* [VTK format](vtk.h)

# Input functions

* [Basilisk snapshots](output.h#dump-basilisk-snapshots)
* [Gerris simulation format](input.h#input_gfs-gerris-simulation-format)
* [ESRI ASCII Grid format](input.h#input_grd-raster-format-esri-grid)
* [Portable Gray Map (PGM) images](input.h#input_pgm-importing-portable-gray-map-pgm-images)

# Interactive Basilisk View

* [bview](bview): a script to start the client/server visualisation pipeline.
* [bview-server.c](bview-server.c): the server.
* [bview-client.py](bview-client.py): the client.

# Miscellaneous functions/modules

* [Controlling the maximum runtime](maxruntime.h)

# Tracking floating-point exceptions

On systems which support [signaling
NaNs](http://en.wikipedia.org/wiki/NaN#Signaling_NaN) (such as
GNU/Linux), Basilisk is set up so that trying to use an unitialised
value will cause a floating-point exception to be triggered and the
program to abort. This is particularly useful when developing adaptive
algorithms and/or debugging boundary conditions.

To maximise the "debugging potential" of this approach it is also
recommended to use the `trash()` function to reset any field prior to
updates. This will guarantee that older values are not mistakenly
reused. Note that this call is quite expensive and needs to be turned on
by adding `-DTRASH=1` to the compilation flags (otherwise it is just
ignored).

Doing

~~~bash
ulimit -c unlimited
~~~

before running the code will allow generation of `core` files which
can be used for post-mortem debugging (e.g. with gdb).

## Visualising stencils

It is often useful to visualise the values of fields in the stencil
which triggered the exception. This can be done using the `-catch`
option of `qcc`.

We will take this code as an example:

~~~c
#include "utils.h"

int main()
{
  init_grid (16);
  scalar a[];
  trash ({a});
  foreach()
    a[] = x;
  vector ga[];
  gradients ({a}, {ga});
}
~~~

Copy and paste this into `test.c`, then do

~~~bash
ulimit -c unlimited
qcc -DTRASH=1 -g -Wall test.c -o test -lm
./test
~~~

you should get

~~~
Floating point exception (core dumped)
~~~

Then do

~~~bash
gdb test core
~~~

you should get

~~~
...
Core was generated by `./test'.
Program terminated with signal 8, Arithmetic exception.
#0  0x0000000000419dbe in gradients (f=0x7fff5f412430, g=0x7fff5f412420)
    at /home/popinet/basilisk/wiki/src/utils.h:203
203		  v.x[] = (s[1,0] - s[-1,0])/(2.*Delta);
~~~

i.e. the exception occured in the [gradients()](utils.h#gradients)
function of [utils.h]().

To visualise the stencil/fields which lead to the exception do

~~~bash
qcc -catch -g -Wall test.c -o test -lm
./test 
~~~

you should now get

~~~
Caught signal 8 (Floating Point Exception)
Caught signal 6 (Aborted)
Last point stencils can be displayed using (in gnuplot)
  set size ratio -1
  set key outside
  v=0
  plot 'cells' w l lc 0, 'stencil' u 1+3*v:2+3*v:3+3*v w labels tc lt 1 title columnhead(3+3*v), 'coarse' u 1+3*v:2+3*v:3+3*v w labels tc lt 3 t ''
Aborted (core dumped)
~~~

Follow the instructions i.e.

~~~bash
gnuplot
gnuplot> set size ratio -1
gnuplot> set key outside
gnuplot> v=0
gnuplot> plot 'cells' w l lc 0, 'stencil' u 1+3*v:2+3*v:3+3*v w labels tc lt 1 title columnhead(3+3*v), 'coarse' u 1+3*v:2+3*v:3+3*v w labels tc lt 3 t ''
~~~

With some zooming and panning, you should get this picture

![Example of stencil/field causing an exception](figures/catch.png)

The red numbers represent the stencil the code was working on when the
exception occured. It is centered on the top-left corner of the
domain. Cells both inside the domain and outside (i.e. ghost cells)
are represented. While the field inside the domain has been
initialised, ghost cell values have not. This causes the `gradients()`
function to generate the exception when it tries to access ghost cell
values.

To initialise the ghost-cell values, we need to apply the boundary
conditions i.e. add

~~~c
  boundary ({a});
~~~

after initialisation. Recompiling and re-running confirms that this
fixes the problem.

Note that the blue numbers are the field values for the parent cells
(in the quadtree hierarchy). We can see that these are also
un-initialised but this is not a problem since we don't use them in
this example.

The `v` value in the gnuplot script is important. It controls which
field is displayed. `v=0` indicates the first field allocated by the
program (i.e. `a[]` in this example), accordingly `ga.x[]` and
`ga.y[]` have indices 1 and 2 respectively.

# See also

* [Tips]()
* [Built-in profiling](README.trace)
* [Built-in memory profiling](README.mtrace)
* [Performance profiling with Paraver](README.paraver)
* [Profiling round-off errors with CADNA](README.cadna)