A literate Lattice Boltzmann Code

Adrian Kummerländer

What?

Literate Programming

…a tool for reproducible science.

Instead of imagining that our main task is to instruct a computer what to do, let us concentrate rather on explaining to human beings what we want a computer to do. —Donald Knuth (Literate Programming, 1984)

Setup

  • Emacs text edtor with org-babel
  • Literate source file can be …
  • tangled to a compilable program
  • woven into a prose document (e.g. HTML / PDF)
  • Languages may be easily mixed

LiterateLB Features

  • Self contained
  • (Smagorinsky) BGK collision
  • Workable set of boundary conditions
  • Integrated (volumetric) visualization
  • Interesting set of examples
  • Competitive GPU performance

Visualization

  • Useful both for science and art
  • Just in Time vs. Real Time
  • Inspired by reality and freed from its restrictions

Signed Distance Functions

  • \(d : \mathbb{R}^3 \to \mathbb{R}\) is shortest distance to surface
  • Iteratively approximate directed distances
  • Surface normals
  • Model for boundary configuration and visualization

Constructive Solid Geometry

  • library of geometric primitives
  • combine with boolean operators
v = translate(v, center);

return ssub(
    rounded(box(v, v3(20, 20, 20)), 5),
    sphere(v, 25),
    5
);

Fancy example

add(
  sadd(
    sub(
      rounded(box(v, v3(5, 28, 38)), 1),
      rounded(box(v, v3(6, 26, 36)), 1)
    ),
    cylinder(translate(v, v3(0,0,-45)), 5, 12),
    1
  ),
  sintersect(
    box(v, v3(5, 28, 38)),
    add(
      add(
        box(rotate_x(v, angle), v3(10, width, 100)),
        box(rotate_x(v, -angle), v3(10, width, 100))
      ),
      add(
        add(
          add(
            box(rotate_x(translate(v, v3(0,0,25)), angle), v3(10, width, 100)),
            box(rotate_x(translate(v, v3(0,0,25)), -angle), v3(10, width, 100))
          ),
          add(
            box(rotate_x(translate(v, v3(0,0,-25)), angle), v3(10, width, 100)),
            box(rotate_x(translate(v, v3(0,0,-25)), -angle), v3(10, width, 100))
          )
        ),
        add(
          add(
            box(rotate_x(translate(v, v3(0,0,50)), angle), v3(10, width, 100)),
            box(rotate_x(translate(v, v3(0,0,50)), -angle), v3(10, width, 100))
          ),
          add(
            box(rotate_x(translate(v, v3(0,0,-50)), angle), v3(10, width, 100)),
            box(rotate_x(translate(v, v3(0,0,-50)), -angle), v3(10, width, 100))
          )
        )
      )
    ),
    2
  )
);

Sphere tracing

Volumetric Visualization

Ray Marching

  • Common practice for volume rendering
  • Periodically sample some scalar value along a ray
  • … velocity norm, density, curl magnitude, shear layer visibility, various quality criteria …
  • Discretized volume rendering equation
  • Extensible to arbitrarily complex light transport

Volume rendering equation

Color \(C\) along ray \(r\) with length \(L\) with absorption \(\mu\): \[C(r) = \int_0^L c(x) \mu(x) \exp\left(-\int_0^x \mu(t) dt\right) dx\]

Simplified for speed: \[C(r) = \sum_{i=0}^N C(i \Delta x) \alpha (i \Delta x) \prod_{j=0}^{i-1} \left(1 - \alpha(j\Delta x)\right)\]

Color palette

  • Can make or break a good visualization
  • SciVisColor.org offers a variety of useful palettes
  • LiterateLB includes a selection of them

blue_orange.png

4wave_ROTB.png

Virtual Rheoscopic Fluid

rheo_paper.png

\(n_\text{shear} = 2 D u \hat{u} - 2\left(u^T D u \hat{u} \right) \hat{u}\)

Shear Layer Normal

  • fuid stress tensor \(\sigma\), strain-rate tensor \(Du\)
  • \(R_{\hat{u}} = \sigma \hat{u}\) force on infinitesimal plane in flow
  • \(R_{\hat{u}} - (R_{\hat{u}} \cdot \hat{u}) \hat{u} = 2\mu Du\hat{u} - 2\mu(\hat{u}^TDu\hat{u})\hat{u}\)
  • Neglect viscosity scalar \[n_\text{shear} = 2 D u \hat{u} - 2\left(u^T D u \hat{u} \right) \hat{u}\]
  • Construct strain-rate tensor from \(f^\text{neq}\)

Shear Layer Extraction

Rheoscopic Fluid

Examples

Flow around a sphere

Nozzle

Flow around rotating shaft