---
title: 'HemeLB: High-Performance Vascular Flow Simulator'
url: https://www.emergentmind.com/topics/hemelb
type: topic
---

# HemeLB: High-Performance Vascular Flow Simulator

Searching arXiv for HemeLB papers to ground the article in the cited literature.
HemeLB is a specialized, massively-parallel lattice-Boltzmann environment for high-resolution simulation of blood flow in patient-specific cerebrovascular geometries, developed to handle extremely sparse vascular domains without wasted memory or compute effort, to scale to tens of thousands of cores, and to provide real-time, in situ visualization and interactive steering for clinical use [1209.3972]. Across later work, the same code base is described as an MPI-parallel, distributed-memory C++ solver for three-dimensional, patient-specific vascular geometries, with extensions for GPU acceleration, arterial–venous self-coupling, and uncertainty quantification and visualization on exascale systems [1211.0205].

## 1. Design objectives and software structure

HemeLB was conceived for blood-flow simulation in realistic, patient-derived vascular trees, including cerebral, retinal, and peripheral vasculature, with particular emphasis on sparse geometries in which only a small fraction of the enclosing Cartesian grid is fluid. The original aims reported for the code include full-3D simulation of realistic vascular trees at high resolution, efficient strong scaling to very large CPU counts, interactive steering and visualization, and support for coupling either to complementary codes or to itself in arterial–venous configurations [2010.04144].

The software is described as an MPI-parallel, distributed-memory C++ code released under LGPL. A notable implementation choice is static polymorphism at compile time for velocity set, collision operator, and boundary-condition class, so that there is no virtual-call overhead in the lattice-Boltzmann inner loop. The supported compile-time variants explicitly mentioned are D3Q15, D3Q19, and D3Q27 velocity sets; LBGK and MRT collision operators; and SBB, GZS, BFL, and JY boundary-condition classes [1211.0205].

Domain construction is handled by preprocessing tools rather than by runtime mesh generation. One described workflow starts from a triangulated vascular surface extracted from imaging, applies Marching Cubes to obtain a closed lumen surface, skeletonizes that surface with the Mean-Curvature Skeleton method, clips more than 1,500 small vessel tips by planes orthogonal to the local skeleton direction, and finally converts the clipped surface to a block-structured lattice in which each block is $16\times16\times16$ sites and each site is tagged as fluid, wall, or outside. A compressed table then stores the mapping from block index to fluid-site and boundary-site lists, so that runtime iteration touches only active sites [2010.04144].

A recurrent distinction in the literature is that HemeLB is not presented as a general-purpose lattice-Boltzmann framework. Rather, it is specialized for the sparsity, boundary complexity, and turnaround-time constraints of patient-specific hemodynamics. This suggests that many of its algorithmic choices—indirect addressing, graph-based partitioning, coalesced halo exchanges, and in situ visualization—are best understood as co-optimized for sparse vascular anatomy rather than for dense, regular domains.

## 2. Numerical formulation and boundary treatment

In the implementations emphasized in the CPU, GPU, and exascale papers, HemeLB uses a three-dimensional single-relaxation-time lattice-Boltzmann method on a D3Q19 lattice to solve the incompressible Navier–Stokes equations for a Newtonian fluid. The continuum equations are written as mass and momentum conservation,
$$
\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,
$$
$$
\partial_t(\rho u_\alpha) + \partial_\beta(\rho u_\alpha u_\beta)
= -\partial_\alpha p + \partial_\beta[\eta(\partial_\beta u_\alpha + \partial_\alpha u_\beta)],
$$
with $p=C_s^2\rho$, $\eta=\rho C_s^2(\tau-\Delta t/2)$, and $C_s=1/\sqrt{3}$ [2202.11770].

The discrete update consists of collision followed by streaming. In one form used in the code description,
$$
f_i(\mathbf{x}+\mathbf{c}_i\Delta t,t+\Delta t)
=
f_i(\mathbf{x},t)-\Omega\bigl[f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t)\bigr],
$$
where $\Omega=\Delta t/\tau$ in BGK form, $c_s^2=1/3$, and the macroscopic fields are recovered from
$$
\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.
$$
The D3Q19 equilibrium distribution is given in the cited work by the standard second-order expansion in $\mathbf{u}$ with weights $w_0=1/3$, $w_{1-6}=1/18$, and $w_{7-18}=1/36$ [2508.15420].

Boundary treatment is a central issue because the target geometries are curved and are rarely aligned with the lattice. HemeLB implements simple bounce-back (SBB), Bouzidi-Firdaouss-Lallemand interpolation (BFL), Guo-Zheng-Shi correction (GZS), and Junk-Yang (JY). In the comparative study performed within HemeLB, SBB is exact second-order only on aligned planar walls and degrades to first order on curved or non-aligned surfaces; BFL and GZS both give second-order convergence in space; and JY is reported as poorly stable at larger Reynolds number in complex geometries [1211.0205].

The reported benchmarks span Poiseuille flow in a tilted cylinder, Womersley flow, and Dean flow in a torus at physiologically relevant Reynolds, Womersley, and Dean numbers. For Poiseuille flow, SBB exhibits $\epsilon_u\propto D^{-1}$ while BFL and GZS exhibit $\epsilon_u\propto D^{-2}$. For Womersley flow at $\alpha=12$, the reported errors are approximately $3$–$4\%$ for SBB and approximately $1$–$1.5\%$ for BFL and GZS, with GZS slightly worse than BFL in unsteady flow. In curved-pipe Dean-flow tests, SBB shows prominent stair-casing artifacts near walls, whereas BFL and GZS recover smooth vortex centers [1211.0205].

The same study also addresses a common assumption that higher-order collision models or larger velocity sets necessarily dominate boundary treatment in importance. Within the tested moderate-Reynolds, single-component Newtonian regime, the choice between LBGK and MRT and between D3Q15, D3Q19, and D3Q27 does not significantly affect accuracy or convergence order, whereas the boundary condition does. The explicit recommendation given there is BFL as the optimal compromise for unsteady, complex-geometry hemodynamics: second-order accurate, stable at moderate Reynolds number, and significantly cheaper than GZS [1211.0205].

## 3. Sparse-geometry data structures and communication model

HemeLB’s core architectural optimization is directed at the fact that vascular domains are extremely sparse. Rather than storing the entire Cartesian lattice, the code discards empty regions and stores only active fluid sites in a one-dimensional array with precomputed neighbor pointers, an indirect-addressing scheme in which non-fluid nodes consume no memory or CPU time [1209.3972]. Later descriptions restate the same design imperative in terms of contiguous arrays for non-empty blocks and compressed tables for active-site traversal, with the explicit observation that realistic arterial or venous networks may contain only $5$–$10\%$ fluid cells in the regular grid that encloses them [2010.04144].

Parallel decomposition is graph-based. At start-up, HemeLB constructs a graph of active sites and uses ParMETIS, specifically `V3_PartKway`, to partition that graph dynamically over $p$ cores while balancing both computation and communication cost [1209.3972]. The communication scheme is correspondingly tailored to sparse nearest-neighbor exchange: all halo-exchange, visualization broadcast/reduce, and steering messages for a time step are bundled into a single `MPI_Isend`/`MPI_Irecv` per neighbor, followed by a single `MPI_Wait` synchronization. This coalesced, asynchronous pattern is reported to minimize latency overhead and improve scaling [1209.3972].

I/O and visualization are integrated into the runtime structure rather than treated as a separate workflow. A small subset of cores, the read-group, with a default of 32, performs parallel MPI-IO for checkpointing. Optionally, one core may be reserved at each time step to generate a ray-traced image and broadcast it to remote clients, at up to approximately $10$–$15$ frames per second; the steering client can send parameter updates back to the simulation [1209.3972].

Memory-management work for very large simulations extends these mechanisms to file-system constraints and node-local sharing. To circumvent the `INT_MAX` limit in 32-bit-based MPI I/O routines, HemeLB implements chunked large-count reads through an `MPI_File_read_x` wrapper or by patching OMPIO to accept `MPI_Count` directly. The same human-scale work also uses MPI-3 shared-memory windows so that one copy of the block-offset table is held per node rather than per rank [2010.04144].

## 4. Performance characteristics and predictive runtime model

The 2012 performance study reports linear or near-linear scaling up to $8\,192$ cores for medium-sized domains of approximately $20$ million sites and up to $16\,384$ cores for large domains of approximately $81$ million sites. At $32\,768$ cores, HemeLB reaches a peak of $29.5$ billion site updates per second. For a highly sparse network with fluid fraction approximately $5\%$, the reported slowdown relative to a denser cylinder with $65\%$ fluid fraction is only $11\%$, which is attributed to the indirect-addressing design [1209.3972].

The same paper gives explicit measurements for steering and in situ rendering. Rendering one image per $100$ LB steps on $2\,048$ cores adds at most approximately $10\%$ overhead. Streaming to a remote client at approximately $4.6$ frames per second, corresponding to about $32$ steps per image, while allowing bidirectional steering adds approximately $28\%$ overhead [1209.3972].

A semi-analytic performance model is provided for the time per LB step on $p$ cores. The calculation term is
$$
T_{calc}=\frac{N/p}{\tau},
$$
where $N$ is the total number of lattice sites, $p$ is the number of compute ranks excluding one steering core, and $\tau$ is the single-core update rate in pure compute-only mode. The communication term is
$$
T_{comm}=\log_2(p)\,\lambda+\frac{S_x}{\sigma},
$$
where $\lambda$ is the per-message latency, $\sigma$ is the network bandwidth per core, and $S_x$ is the geometry-dependent halo-exchange volume per core per step. The reported load-imbalance factors are $\zeta_{calc}=1.04$ and $\zeta_{comm}=1.5$, and the runtime monitoring overhead is approximately $O_{monitoring}\simeq 0.06$ [1209.3972].

For the halo volume, the cited geometry-dependent fits are
$$
S_{cyl}=1898\,(N/p)^{0.4827},
$$
$$
S_{bif}=942\,(N/p)^{0.5955},
$$
$$
S_{net}=1176\,(N/p)^{0.6134}.
$$
When visualization is enabled, the per-step cost is written as $T_{step\_vis}=T_{step}+T_{images}$, with an empirical fit on $2\,048$ cores
$$
T_{images}\approx 21.6\,k^{-0.76}\;[\mathrm{s}],
$$
where $k$ is the number of LB steps per rendered image [1209.3972].

This model is used to define an operating regime rather than merely to fit benchmark curves. The reported efficiency sweet spot is $5\,000\le N/p\le 500\,000$ sites per core, where compute and communication are balanced and near-peak SUPS per core are obtained. The associated operational recommendations are to keep $N/p$ in that range, reserve one core for steering when interactive control is required, render images in situ at approximately 100-step intervals to keep overhead below $10\%$, write full-domain snapshots sparingly, map processes to a single network island where possible, and use advance reservation or urgent-compute policies when clinical decision timescales matter [1209.3972].

## 5. Load balancing, decomposition refinements, and self-coupling

As HemeLB moved toward larger and sparser problems, decomposition quality became a first-order concern. One refinement is weighted decomposition, in which lattice sites are classified as bulk, wall, in/outlet, or wall-plus-in/outlet, with integer weights passed to ParMETIS so that the partitioner accounts for differing per-site costs. The rounded weights reported for production use are $w_{bulk}=4$, $w_{wall}=8$, $w_{inlet}=16$, and $w_{wi}=16$, obtained from least-squares calibration of pure-compute timings on Intel SandyBridge and AMD Interlagos systems [1410.4713].

A second refinement is Morton-order site sorting. HemeLB optionally sorts sites by a three-dimensional Morton, or binary Z-order, curve and uses `ParMETIS_V3_PartGeomKway` with the pre-sorted vertex list and weights. The intended effect is improved locality and reduced communication surface area of the resulting partitions [1410.4713].

The benefits and trade-offs are explicitly quantified. On a bifurcation geometry with $650\,492$ sites and $10\%$ fill at $1\,024$ cores, weighted decomposition reduces the compute-load imbalance from $I_{calc}=1.15$ to $1.02$, a reduction of approximately $85\%$, and combining weighting with the space-filling-curve ordering gives $I_{calc}\approx1.01$ with reduced communication time. On a much sparser aneurysm geometry with $5\,667\,778$ sites and $1.5\%$ fill at $8\,192$ cores, weighting similarly reduces compute imbalance but can increase communication time on HECToR by approximately $20$–$25\%$, whereas on ARCHER the communication penalty disappears and becomes a gain. The paper therefore frames weighted decomposition as a machine-dependent trade-off between compute balance and communication balance rather than as a universally monotone optimization [1410.4713].

Human-scale simulation work adds a separate axis of extension: self-coupling. HemeLB can instantiate multiple worlds in the same MPI job, such as an arterial tree and a venous tree, each with its own domain decomposition, LB updates, and internal MPI communication. Designated master ranks exchange boundary data at coupled inlets and outlets every $N$ coupling steps. The coupling strategy uses a $1\!:\!N$ mapping between arterial outlets and venous inlets, assigns a random capillary pressure drop $\Delta P_i/P_0\in[0.3,0.7]$, enforces mass conservation $\sum_i q_i=q_0$, and exchanges only averaged $\rho$ and $\mathbf{u}$ values, so the communication volume per synchronization is on the order of $10$–$100$ doubles [2010.04144].

The same paper reports strong-scaling results at human scale. On Blue Waters, for a circle-of-Willis geometry with $10.15\times10^9$ lattice sites, scaling from $18\,000$ to $170\,000$ cores gives a speed-up of $128.7\times$ for a $170.7\times$ increase in core count, corresponding to $75\%$ parallel efficiency; at $288\,000$ cores, the speed-up remains $189.8\times$, about $66\%$ efficiency. On SuperMUC-NG, scaling from $864$ to $309\,696$ ranks on the same approximately $10^{10}$-site model gives near-linear speed-up up to $190\times$, with computational and communication efficiency both above $97\%$, overall efficiency at or above $87\%$, and $5\,000$ time steps completed in $83$ seconds on $309\,696$ cores [2010.04144].

## 6. GPU implementations, visualization, and exascale uncertainty analysis

HemeLB has also been reimplemented for accelerator-based systems. The large-scale GPU code is written in CUDA C++, uses a structure-of-arrays layout for the 19 distribution streams, and places the principal LBM state and auxiliary tables in GPU global or constant memory. Six CUDA kernels implement collision and streaming for distinct site classes—interior, walls, inlet, outlet, inlet-plus-wall, and outlet-plus-wall—and the host retains responsibility for MPI-driven halo exchange, buffer packing and unpacking, and pointer swaps between `f_old` and `f_new` [2202.11770].

GPU parallel execution is organized by a `StepManager` with stages `PreSend`, `Send`, `PreReceive`, `Receive`, `PostReceive`, and `EndIteration`. The implementation reorders work so that domain-edge and mid-domain kernels are launched before MPI sends, and it uses multiple CUDA streams together with asynchronous device-to-host and host-to-device transfers to overlap GPU computation with inter-node communication. Work decomposition remains ParMETIS-driven on the sparse vascular mesh, with one MPI rank driving one GPU [2202.11770].

The Summit benchmarks show that the accelerator version retains the strong-scaling behavior of the CPU code at much larger task counts. For the GPU implementation on Summit with six V100 GPUs per node, baseline strong scaling from $768$ GPUs to $6\,144$ GPUs yields $90\%$ efficiency, and an ExaPipe run on $18\,432$ GPUs, approximately two-thirds of Summit, yields $74\%$ efficiency. Reported throughput is approximately $120$ MLUPS per node, dropping when the number of sites per GPU falls below $10^6$. The paper also reports a GPU-versus-CPU speed-up of approximately $\times85$ at $12$k units when the comparison is made against raw “cores,” and notes that when GPU “cores” are equated to streaming multiprocessors, the code continues scaling beyond $30$ million parallel tasks [2202.11770].

A separate GPU acceleration and visualization effort targets hospital deployment on a workstation with four GTX 1080 Ti devices. There, HemeLB is coupled to a CUDA/OpenGL rendering pipeline and a Qt-based interface. For two patient-specific STL vascular domains, the reported performance is $15\,168\,964$ SUPS for a $296\,814$-site case and $13\,674\,962$ SUPS for a $132\,910$-site case, with near-linear scaling up to four GPUs and efficiency at or above $80\%$. The same work estimates that at approximately $1.5\times10^7$ SUPS, a $300$k-cell vascular model requires about $0.02$ seconds per LB step and about $40$ seconds for a $2\,000$-step cardiac cycle on four GTX 1080 Ti GPUs, while the visualization pipeline delivers frame rates above $30$ frames per second for typical approximately $300$k-voxel volumes [1906.11925].

More recent work places HemeLB in an exascale uncertainty-quantification and visualization workflow. In that setting, the code is combined with a large-eddy Smagorinsky model, EasyVVUQ for non-intrusive parameter sampling, and Frontier for ensemble execution. The effective viscosity is written as $\nu_{\rm eff}=\nu_0+\nu_t$ with
$$
\nu_t=(C_{\rm smag}\,\Delta)^2\,|\bar{\mathbf{S}}|,
$$
and uncertainty in quantities such as pressure and wall shear stress is visualized through Shannon entropy on marching-cubes topologies,
$$
E(q)=-\sum_{t=1}^{256} p_t(q)\log_2 p_t(q),
$$
and an isosurface-crossing probability
$$
I(q)=\sum_{t=2}^{255} p_t(q).
$$
The Frontier study reports near-ideal strong scaling up to more than $100\,000$ MPI ranks for domains of order $10^9$ cells, constant MLUPS per GPU node under weak scaling, achieved per-node throughput on MI250X of order $100$ MLUPS, and integration of sampled parameters such as inlet peak velocity and Smagorinsky constant directly into HemeLB configuration files [2508.15420].

Taken together, these strands define HemeLB as a specialized vascular LBM environment whose distinguishing features are sparse-domain execution, graph-based decomposition, interactive in situ visualization, and a progression from clinically oriented cerebrovascular simulation to human-scale, accelerator-based, and uncertainty-aware workflows. The literature does not present a single static form of the software; instead, it shows a code base repeatedly extended while preserving the same central premise that sparse vascular anatomy, communication structure, and visualization requirements must be treated as a single HPC design problem.

Source: https://www.emergentmind.com/topics/hemelb