src/README

    Solvers

    Saint-Venant

    \displaystyle \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 \displaystyle \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)

    \displaystyle \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 \displaystyle \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) \displaystyle \alpha h\mathcal{T} \left( D \right) + hD = b \displaystyle b = \left[ \frac{g}{\alpha} \nabla \eta +\mathcal{Q}_1 \left( u \right) \right]

    • Multilayer hydrostatic

    \displaystyle \begin{aligned} \partial_t h_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k & = 0,\\ \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) \end{aligned}

    Systems of conservation laws

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

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

    \displaystyle \frac{\partial (f \rho_i)}{\partial t } + \nabla \cdot (f \rho_i \mathbf{u}) = 0 \displaystyle \frac{\partial (\rho_i \mathbf{u})}{\partial t } + \nabla \cdot ( \rho_i \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \tau_i' \displaystyle \frac{\partial [\rho_i (e_i + \mathbf{u}^2/2)]}{\partial t } + \nabla \cdot [ \rho_i \mathbf{u} (e_i + \mathbf{u}^2/2)] = -\nabla \cdot (\mathbf{u} p_i) + \nabla \cdot \left( \tau'_i \mathbf{u} \right)

    \displaystyle \rho_ic_{p,i}\frac{DT_i}{Dt} = \beta_iT_i\frac{Dp_i}{Dt} - \nabla\cdot\mathbf{q}_i \displaystyle \left(\frac{\gamma_i}{\rho_ic_i^2} - \frac{\beta_i^2T_i}{\rho_ic_{p,i}}\right)\frac{Dp_i}{Dt} = -\frac{\beta_i}{\rho_ic_{p,i}}\nabla\cdot\mathbf{q}_i - \nabla\cdot\mathbf{u}_i

    \displaystyle \begin{aligned} \partial_t h_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k & = 0,\\ \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) - \mathbf{{\nabla}} (h \phi)_k + \left[ \phi \mathbf{{\nabla}} z \right]_k,\\ \partial_t (hw)_k + \mathbf{{\nabla}} \cdot \left( hw \mathbf{u} \right)_k & = - [\phi]_k,\\ \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k + \left[ w - \mathbf{u} \cdot \mathbf{{\nabla}} z \right]_k & = 0 \end{aligned}

    Electrohydrodynamics

    Viscoelasticity

    \displaystyle \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} \displaystyle \mathbf{\tau}_p = \frac{\mu_p \mathbf{f_s}(\mathbf{A})}{\lambda} \displaystyle \Psi = \log\mathbf{A} \displaystyle D_t \Psi = (\Omega \cdot \Psi -\Psi \cdot \Omega) + 2 \mathbf{B} + \frac{e^{-\Psi} \mathbf{f}_r (e^{\Psi})}{\lambda}

    Other equations

    The Basilisk C preprocessor

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

    Metric scale factors

    The face vector 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:

    Data processing

    Output functions

    Input functions

    Basilisk View

    Miscellaneous functions/modules

    Tracking floating-point exceptions

    On systems which support signaling NaNs (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

    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:

    #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

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

    you should get

    Floating point exception (core dumped)

    Then do

    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() function of utils.h.

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

    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.

    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

    Example of stencil/field causing an exception

    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

      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.

    Tracing permissions

    Some recent systems disallow tracing of processes for security reasons. The symptom will be an error message from gdb looking like:

    Attaching to process 9351
    Could not attach to process.  If your uid matches the uid of the target
    process, check the setting of /proc/sys/kernel/yama/ptrace_scope, or try
    again as the root user.  For more details, see /etc/sysctl.d/10-ptrace.conf
    ptrace: Operation not permitted.

    To enable tracing (which weakens your system’s security), you need to do:

    sudo sh -c 'echo 0 > /proc/sys/kernel/yama/ptrace_scope'

    See also