C++ Implementation & Code Architecture

This page documents the solver architecture, memory model, parallelization strategy, build system, and test suite. Every design decision is explained with its performance and engineering rationale.

1. Module Reference

2. Memory Model

The 2D fluid grid is stored in a single std::vector<double> with 1D indexing. This is the single most important performance decision in the code. Compare the two approaches:

// BAD: Nested vector, non-contiguous memory, pointer indirection per row
std::vector<std::vector<std::vector<double>>> f_bad(NY, std::vector<std::vector<double>>(NX, std::vector<double>(9)));

// GOOD: Flat 1D array, contiguous memory, single allocation, cache-friendly
std::vector<double> f(NX * NY * 9);
// Access at (x, y, i):  f[(y * NX + x) * 9 + i]

The flat layout ensures that the 9 distributions at a single node are stored in consecutive memory addresses. During the collision step, the CPU prefetcher can load all 9 values into cache in a single cache line operation. With nested vectors, each row access requires a separate pointer dereference and heap allocation lookup.

3. The Core Solver Loop

The execution flow for each timestep is:

1. enforce_inflow()    Zou/He velocity inlet at x = 0
2. enforce_outflow()   Convective outlet at x = NX-1
3. collide()           BGK relaxation (OpenMP parallel)
4. stream()            Propagate + bounce-back (OpenMP parallel)
5. extract_forces()    Momentum exchange on cylinder surface

The collision step, shown below, is the computational heart of the solver. It applies the BGK operator at every fluid node:

// Collision step: relax each node toward equilibrium
#pragma omp parallel for collapse(2)
for (int y = 0; y < NY; ++y) {
  for (int x = 0; x < NX; ++x) {
    int idx = y * NX + x;
    if (sys.obstacle[idx]) continue;

    double rho = 0.0, mx = 0.0, my = 0.0;
    double* f_node = &sys.f[idx * 9];
    for (int i = 0; i < 9; ++i) {
      rho += f_node[i];
      mx  += f_node[i] * cx[i];
      my  += f_node[i] * cy[i];
    }
    double u = mx / rho, v = my / rho;

    for (int i = 0; i < 9; ++i) {
      double feq = compute_equilibrium(i, rho, u, v);
      f_node[i] -= (1.0 / tau) * (f_node[i] - feq);
    }
  }
}

Key optimizations in this loop:

• Flat pointer access: f_node points to the start of the 9-value block for the current node, no index recomputation inside the inner loops
• Macro computation fused: density and momentum are computed in a single pass over the 9 directions
• Branch elision: obstacle check at node level, not inside the direction loop

4. OpenMP Parallelization

Both the collision and streaming loops are parallelized using #pragma omp parallel for collapse(2). The collapse(2) clause fuses the nested y and x loops into a single iteration space, enabling better load balancing across threads:

#pragma omp parallel for collapse(2)
for (int y = 0; y < NY; ++y) {
  for (int x = 0; x < NX; ++x) {
    // ... collision work ...
  }
}

Without collapse(2), if NY = 150 and NX = 400, only 150 chunks would be distributed across threads, causing load imbalance when some rows contain obstacle nodes. With collapse(2), 60,000 individual node tasks are distributed, giving near-perfect load balance.

Collision and streaming are separate parallel regions; they cannot be fused into a single loop because streaming reads from the post-collision state of neighboring nodes, introducing a data dependency across grid cells.

5. Force Extraction: Momentum Exchange Method

The forces on the cylinder are computed by summing momentum transfers across all fluid-obstacle boundary links. Each boundary link connects a fluid node adjacent to an obstacle node:

// For each obstacle node, check all 9 neighbor directions
for (int y = 0; y < NY; ++y) {
  for (int x = 0; x < NX; ++x) {
    int node_idx = y * NX + x;
    if (!sys.obstacle[node_idx]) continue;

    for (int i = 0; i < 9; ++i) {
      int nx = x + cx[i], ny = y + cy[i];
      // ... boundary check + periodicity ...
      int fluid_idx = ny * NX + nx;
      if (!sys.obstacle[fluid_idx]) {
        double delta = sys.f[fluid_idx * 9 + i] - sys.f[node_idx * 9 + bounce_back[i]];
        fx += cx[i] * delta;
        fy += cy[i] * delta;
      }
    }
  }
}

This method is exact in the LBM framework (no modeling error for stationary boundaries) and requires no body-fitted mesh or interpolation; the forces emerge naturally from the distribution functions at the boundary.

6. Build System & Quick Start

# Configure and build all targets
cmake -B build && cmake --build build

# Run a case via its CMake target
cmake --build build --target sim        # cylinder (LBM_Engine)
cmake --build build --target cavity     # lid-driven cavity
cmake --build build --target flat_plate # flat plate boundary layer
cmake --build build --target amr        # AMR demonstration

# Or run the executable directly with arguments
./build/LBM_Engine 200          # cylinder at Re=200
./build/LBM_Engine 100 12000    # cylinder at Re=100 for 12000 steps

# Run the Google Test suite (12 cases)
./build/LBM_Tests

# Batch sweep across cases
bash scripts/run_all_cases.sh

# Post-process VTK to JSON for the web viewer
python3 scripts/postprocess.py output/re100 --json --every 5

# Preview the documentation site
python3 -m http.server -d docs 8765

7. Build Configuration

8. Test Suite

The test suite covers 12 test cases across 8 test groups:

9. Solver Output & the ML Data Pipeline

Every simulation run emits a structured result set that doubles as training data for the neural-network surrogate. Per frame, the solver writes frame_*.json containing the velocity components (u, v), the vorticity field (omega), and the pressure field (p) on a downsampled grid; a meta.json records the grid dimensions, lattice spacing, Reynolds number, relaxation time, and inflow velocity; and forces.jsonl appends time-resolved drag and lift coefficients. Paraview VTK is written alongside for 3D inspection and offline post-processing.

scripts/postprocess.py converts the VTK output to JSON/PNG and can slice individual frames for the web gallery (with --split, --cmap, --strouhal, and --vorticity flags). The neural-network data loader (pinn/data/loader.py) then ingests those frame*.json files into NumPy arrays, samples sensor points (importance-sampled near the geometry) and collocation points for the PDE residual, and normalizes the coordinates (x, y) to [-1, 1] and the physical parameters (e.g. Reynolds number to Re_n). This JSON-to-NumPy step is the bridge that turns a multi-hour C++ simulation into a training batch.

This result pipeline is exactly what feeds the surrogate work documented next: Physics-Informed Neural Networks.

10. Continuous Integration

GitHub Actions builds and tests on both Ubuntu and macOS on every push and pull request:

.github/workflows/ci.yml
  - matrix: [ubuntu-latest, macos-latest]
  - cmake -B build
  - cmake --build build
  - ./build/LBM_Tests