hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

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