Lattice Documentation

Quick Start

Build from source

# Dependencies (Ubuntu/Debian)
sudo apt install build-essential cmake pkg-config \
  libsdl2-dev libsdl2-ttf-dev mesa-common-dev \
  libgl1-mesa-dev libegl1-mesa-dev

# Build
cmake -B build -S simulation -DCMAKE_BUILD_TYPE=Release
cmake --build build -j$(nproc)

# Run interactively
cd simulation
../build/3d_fluid_simulation_car --wind=1.0

Headless render

mkdir -p frames
../build/3d_fluid_simulation_car \
  --wind=2.0 --duration=10 \
  --model=assets/3d-files/ahmed_25deg_m.obj \
  --output=frames --grid=256x128x128

Web frontend

cd website
npm ci && npm run dev
# Open http://localhost:3000

How It Works

Each frame runs five GPU compute passes in sequence:

1. Collision -- relax distribution functions toward equilibrium (MRT by default; BGK, regularized, and entropic operators available, with optional Smagorinsky SGS turbulence model)
2. Streaming -- propagate distributions to neighboring cells
3. Force computation -- accumulate drag/lift via momentum exchange at solid boundaries
4. Particle advection -- move tracer particles through the velocity field
5. Rendering -- draw particles as points or streamline trails

The LBM solver uses a D3Q19 lattice (19 velocity directions in 3D). Solid geometry is voxelized from OBJ meshes using ray casting. Boundary conditions are Zou-He velocity inlet, Zou-He pressure outlet, and Bouzidi interpolated bounce-back on solid walls. The Bouzidi scheme stores the fractional distance from each fluid cell to the actual mesh surface along every lattice link, then applies quadratic or linear interpolation during streaming. This gives second-order wall accuracy without refining the grid.

Two refinements keep the inflow and outflow physical. The inlet can inject synthetic turbulence (random Fourier modes at a configurable intensity) to mimic the 1-2% residual turbulence of a real wind tunnel. Before the outlet, a sponge layer over the last eighth of the domain ramps viscosity up and relaxes densities toward the reference state, so outgoing pressure waves are absorbed instead of reflecting -- without it, the inlet and outlet form a streamwise acoustic resonator whose standing wave contaminates the drag signal on coarse grids.

Architecture

lattice/
  simulation/
    src/main.c           # Render loop, CLI, GL setup
    src/lbm.c            # LBM grid creation, force readback
    src/ml_predict.c     # ML inference forward pass (pure C)
    shaders/
      lbm_collide.comp   # BGK + Smagorinsky + BCs
      lbm_stream.comp    # Pull-based streaming
      lbm_force.comp     # Momentum exchange drag
      particle.comp      # Particle advection + collision
      trail_update.comp  # Streamline trail history
    lib/                 # Headers
    assets/              # OBJ models, fonts
  website/
    src/app/page.tsx     # Main simulator UI
    src/app/api/         # Next.js API routes
    src/components/      # React components
  ml/
    train.cpp            # Training driver
    data_gen.py          # Modal sweep for training data
    framework/           # C++ autodiff, layers, optimizer

CLI Reference

API Reference

The web frontend talks to a Modal GPU worker via a REST endpoint.

POST /api/render

{
  "wind_speed": 1.0,
  "duration": 10,
  "model": "ahmed25",
  "viz_mode": 1,
  "collision_mode": 2,
  "reynolds": 0
}

{
  "status": "complete",
  "video_url": "https://...",
  "cd_value": 0.295,
  "cl_value": -0.012,
  "cd_series": [0.31, 0.30, 0.295],
  "effective_re": 480,
  "grid_size": "256x128x128"
}

The worker runs on an NVIDIA T4 GPU via Modal. Cold starts take 30-60 seconds; warm starts respond in ~15 seconds for a 10-second simulation.

Visualization Modes

Press V to cycle through modes in interactive mode. Streamline mode (7) stores a 32-frame position trail per particle and renders as GL_LINE_STRIP with velocity-based coloring and alpha fade.

Physics Details

Governing Equations

Fluid flow is governed by the Navier-Stokes equations for incompressible flow:

where uis the velocity field, p is pressure, ρ is density, and ν is kinematic viscosity. The Reynolds number Re = UL/ν characterizes the flow regime. For automotive aerodynamics at highway speeds, Re reaches 10⁶ or higher, indicating fully turbulent flow.

Lattice Boltzmann Method

Rather than solving Navier-Stokes directly, the LBM evolves particle distribution functions on a D3Q19 lattice (3 dimensions, 19 velocities): 1 rest, 6 face-connected, and 12 edge-connected neighbors. The BGK collision operator relaxes distributions toward equilibrium:

Equilibrium Distribution

The Maxwell-Boltzmann equilibrium distribution is:

with weights w₀ = 1/3, w_1-6 = 1/18, w_7-18 = 1/36 and speed of sound c_s² = 1/3.

Collision Kernel

#version 430 core
layout(local_size_x = 8, local_size_y = 8, local_size_z = 8) in;

uniform float tau;

void main() {
    uvec3 pos = gl_GlobalInvocationID;
    uint idx = pos.x + pos.y*NX + pos.z*NX*NY;

    // Compute density and velocity
    float rho = 0.0;
    vec3 u = vec3(0.0);
    for (int i = 0; i < 19; i++) {
        float fi = f[i][idx];
        rho += fi;
        u += fi * e[i];
    }
    u /= rho;

    // BGK collision
    float omega = 1.0 / tau;
    for (int i = 0; i < 19; i++) {
        float feq = equilibrium(i, rho, u);
        f[i][idx] -= omega * (f[i][idx] - feq);
    }
}

Smagorinsky SGS Model

For turbulent flows, the Smagorinsky subgrid-scale model computes a local eddy viscosity from the non-equilibrium stress tensor:

where C_s is the Smagorinsky constant (default 0.1, tunable via --smagorinsky) and |Π| is the Frobenius norm of the non-equilibrium stress. This increases τ locally in turbulent regions, preventing instability without damping the entire domain.

Entropic Collision

As an alternative to MRT, the entropic LBM (--entropic) enforces the discrete H-theorem each step: the relaxation is written as f' = f + αβ(f^eq− f) and α is solved per cell (Newton-Raphson, 2-3 iterations) so that the entropy H = ∑ f_i ln(f_i/w_i) never decreases. In under-resolved regions α drops below the BGK value of 2, adding exactly enough local dissipation to keep the scheme stable -- no eddy-viscosity constant to tune.

Drag Computation

Surface forces use the Mei-Luo-Shyy momentum exchange method, which pairs with the Bouzidi boundary to account for the actual wall position rather than assuming it sits on the grid midpoint:

where f_i* is the post-collision distribution heading toward the wall and f_ī is the Bouzidi-reflected distribution coming back. The force is accumulated on the GPU every LBM step and each reported sample is the mean over its full sampling window, which suppresses the discretization noise that dominates instantaneous momentum-exchange snapshots on coarse grids. The drag and lift coefficients are then:

ML Surrogate Model

A geometry-aware neural network predicts Cd and Cl for any input mesh in under a millisecond. It gives an instant estimate before the LBM has time to converge and lets the browser return a number with no server round trip. The model reads the shape of the body directly rather than a fixed label, so it generalizes to geometries it was never trained on.

Reading geometry

Every mesh is first placed in the same frame the solver uses. The body is centered, its longest axis is rotated onto the flow direction, and it is scaled into a unit cube. From that canonical form the encoder measures nineteen scale invariant descriptors that capture how the body displaces air. These include the three aspect ratios of the bounding box, solidity and convexity, slenderness and fineness, frontal fill, wetted area ratio, sphericity, the fore aft shift of the volume centroid, and an eight station profile of cross section area along the flow axis. The result is a fixed length vector that describes any mesh, which is what lets one model serve an Ahmed body, a full car, and an SUV at the same time.

Architecture

The network takes twenty z-score normalized inputs, being the nineteen geometry descriptors plus the base ten logarithm of the Reynolds number, and outputs two values, Cd and Cl. Reynolds is the only flow input the model needs, because at a fixed Reynolds number the force coefficients no longer depend on the wind speed. Two hidden layers with ReLU activations keep the network small enough to evaluate instantly in the browser.

Linear(20, 64) -> ReLU -> Linear(64, 64) -> ReLU -> Linear(64, 2)

inputs   19 shape descriptors + log10(Reynolds)
outputs  Cd, Cl
about 5,600 parameters total

Training data

The model is trained on 1,421 high fidelity CFD samples drawn from four public reference datasets. These are AhmedML for the Ahmed body, WindsorML for the Windsor body, DrivAerML for the DrivAer road car, and a sample of SHIFT-SUV for a detailed SUV. Each source reports its force coefficients against the frontal area of that specific geometry, so the Cd convention stays consistent as the shapes vary. Meshes are fetched and encoded on Modal, which keeps the terabyte scale geometry off the local machine.

AhmedML     499 Ahmed body variants     CC-BY-SA
WindsorML   349 Windsor body variants   CC-BY-SA
DrivAerML   484 DrivAer car variants    CC-BY-SA
SHIFT-SUV    89 SUV variants (sample)   CC-BY-NC
            ----
            1,421 geometries

Training setup

Training runs in JAX with optax. Inputs and targets are z-score normalized and the normalizer is saved next to the weights. The objective is mean squared error on the two coefficients, optimized with AdamW under gradient clipping and a learning rate that halves on plateau.

Framework:    JAX + optax
Optimizer:    AdamW (lr=1e-3, weight_decay=1e-4, clip_norm=1.0)
Loss:         MSE on [Cd, Cl]
Batch size:   64
Schedule:     halve lr on plateau, early stopping
Train/val:    80/20 split

Accuracy

On a held out split the model reaches a mean absolute error of 0.020 on Cd and 0.101 on Cl. The more demanding test is leave one geometry out, where every shape is scored by a model that never saw it during training. Under that protocol the network reaches a mean absolute error of 0.018 on Cd against 0.024 for a baseline that always predicts the dataset mean, and 0.104 against 0.150 on Cl. Beating the mean baseline on shapes held out of training is the evidence that the model has learned to read geometry rather than memorize the training bodies.

Metric                     Model     Mean baseline
Held out Cd MAE            0.020
Held out Cl MAE            0.101
Leave one geometry out Cd  0.018     0.024
Leave one geometry out Cl  0.104     0.150

Inference

An export step writes the weights, the normalizer, and precomputed descriptors for the bundled geometries into a single JSON file. The browser loads /models/anyobj.json and runs a pure TypeScript forward pass with no dependencies, so a prediction costs well under a millisecond and needs no backend.

# Fetch the public datasets into the Modal volume
modal run ml/anyobj/modal_train.py --fetch ahmedml,windsorml,drivaerml

# Train the surrogate
modal run ml/anyobj/modal_train.py

# Export the browser model
python ml/anyobj/export_web.py

Validation

The solver uses a sphere drag test to validate the force computation pipeline (momentum exchange, Bouzidi boundary, Cd formula). At Re=100 with 16 cells across the sphere diameter, the measured Cd falls within the expected range from Clift, Grace & Weber (1978).

Sphere reference test

Why Cd differs from wind tunnel data

Published Ahmed body data (Ahmed 1984, Lienhart 2003) is measured at Re > 500,000 where the flow is fully turbulent. The LBM solver on typical grids (128-256 cells) operates at Re = 50-200 where the flow is laminar. At laminar Re the boundary layer is thicker, separation occurs earlier, and Cd is 3-10x higher than at turbulent Re. This is physically correct behavior, not a bug.

The Smagorinsky SGS turbulence model adds local eddy viscosity that partially compensates, but matching high-Re experimental data requires grids with hundreds of cells across the body -- beyond what consumer GPUs can handle in real time.

Grid size and Reynolds number

The achievable Reynolds number is limited by the grid resolution. The relaxation parameter τ must stay above 0.52 for numerical stability. This relationship determines the maximum Re for each grid:

Try it yourself

Use the calculator below to see what Reynolds number and memory footprint your grid configuration produces:

Performance

Performance depends on grid size and GPU. The LBM kernel is memory-bandwidth-bound: each cell reads and writes 19 float32 distributions per step. Typical throughput on modern GPUs:

MLUPS = Million Lattice Updates Per Second. Memory usage per cell is approximately 200 bytes (19 distributions × 2 buffers + velocity + solid + Bouzidi q). A 256×128×128 grid needs about 960 MB of GPU buffer space.

Convergence time

The drag coefficient needs 3-5 flow-throughs to converge. A flow-through is the time for a fluid element to traverse the domain: N_x / U_lattice steps. For a 256-cell domain at U = 0.05, one flow-through takes 5,120 steps. The simulation reports an exponential moving average (EMA) of Cd that stabilizes after approximately 50 samples.

Contributing

Lattice is open source under the MIT license. Contributions are welcome -- see CONTRIBUTING.md for the full style guide, testing, and PR process.

Code style

The simulation is C11, formatted with clang-format (LLVM style, 80 char limit). The website is TypeScript/React with Tailwind. Tests run in CI via GitHub Actions.

Commit sign-off

Every commit must carry a sign-off line certifying the Developer Certificate of Origin. Use git commit -s to add it automatically:

git commit -s -m "lbm: fix bounce-back at concave corners"

# Produces:
# lbm: fix bounce-back at concave corners
#
# Signed-off-by: Your Name <your@email.com>

Running tests

# Run tests
cd simulation && xvfb-run ../build/test_lbm  # LBM unit tests
xvfb-run bash test/test_integration.sh       # Integration test
cd website && npm test                        # Frontend tests
cd ml/framework && ./build/test_ml            # ML framework tests