HemeLB: High-Performance Vascular Flow Simulator

• HemeLB is a specialized lattice-Boltzmann environment that simulates patient-specific blood flow in sparse, complex vascular geometries using optimized data structures.
• It employs MPI-parallelism with graph-based partitioning and static polymorphism to achieve strong scaling on tens of thousands of cores with minimal memory overhead.
• The tool integrates GPU acceleration, in situ visualization, and uncertainty quantification to support clinical and exascale hemodynamic simulations.

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 (Groen et al., 2012). 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 (Nash et al., 2012).

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 (McCullough et al., 2020).

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 (Nash et al., 2012).

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 (McCullough et al., 2020).

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}$ (Zacharoudiou et al., 2022).

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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$0 with weights $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$1, $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$2, and $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$3 (Xue et al., 21 Aug 2025).

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 (Nash et al., 2012).

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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$4 while BFL and GZS exhibit $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$5. For Womersley flow at $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$6, the reported errors are approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$7–$$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$8 for SBB and approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$9–$$\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)],$$0 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 (Nash et al., 2012).

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 (Nash et al., 2012).

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 (Groen et al., 2012). 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 $$\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)],$$1–$$\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)],$$2 fluid cells in the regular grid that encloses them (McCullough et al., 2020).

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 $$\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)],$$3 cores while balancing both computation and communication cost (Groen et al., 2012). 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 (Groen et al., 2012).

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 $$\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)],$$4–$$\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)],$$5 frames per second; the steering client can send parameter updates back to the simulation (Groen et al., 2012).

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 (McCullough et al., 2020).

4. Performance characteristics and predictive runtime model

The 2012 performance study reports linear or near-linear scaling up to $$\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)],$$6 cores for medium-sized domains of approximately $$\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)],$$7 million sites and up to $$\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)],$$8 cores for large domains of approximately $$\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)],$$9 million sites. At $p=C_s^2\rho$0 cores, HemeLB reaches a peak of $p=C_s^2\rho$1 billion site updates per second. For a highly sparse network with fluid fraction approximately $p=C_s^2\rho$2, the reported slowdown relative to a denser cylinder with $p=C_s^2\rho$3 fluid fraction is only $p=C_s^2\rho$4, which is attributed to the indirect-addressing design (Groen et al., 2012).

The same paper gives explicit measurements for steering and in situ rendering. Rendering one image per $p=C_s^2\rho$5 LB steps on $p=C_s^2\rho$6 cores adds at most approximately $p=C_s^2\rho$7 overhead. Streaming to a remote client at approximately $p=C_s^2\rho$8 frames per second, corresponding to about $p=C_s^2\rho$9 steps per image, while allowing bidirectional steering adds approximately $\eta=\rho C_s^2(\tau-\Delta t/2)$0 overhead (Groen et al., 2012).

A semi-analytic performance model is provided for the time per LB step on $\eta=\rho C_s^2(\tau-\Delta t/2)$1 cores. The calculation term is

$\eta=\rho C_s^2(\tau-\Delta t/2)$2

where $\eta=\rho C_s^2(\tau-\Delta t/2)$3 is the total number of lattice sites, $\eta=\rho C_s^2(\tau-\Delta t/2)$4 is the number of compute ranks excluding one steering core, and $\eta=\rho C_s^2(\tau-\Delta t/2)$5 is the single-core update rate in pure compute-only mode. The communication term is

$\eta=\rho C_s^2(\tau-\Delta t/2)$6

where $\eta=\rho C_s^2(\tau-\Delta t/2)$7 is the per-message latency, $\eta=\rho C_s^2(\tau-\Delta t/2)$8 is the network bandwidth per core, and $\eta=\rho C_s^2(\tau-\Delta t/2)$9 is the geometry-dependent halo-exchange volume per core per step. The reported load-imbalance factors are $C_s=1/\sqrt{3}$0 and $C_s=1/\sqrt{3}$1, and the runtime monitoring overhead is approximately $C_s=1/\sqrt{3}$2 (Groen et al., 2012).

For the halo volume, the cited geometry-dependent fits are

$C_s=1/\sqrt{3}$3

$C_s=1/\sqrt{3}$4

$C_s=1/\sqrt{3}$5

When visualization is enabled, the per-step cost is written as $C_s=1/\sqrt{3}$6, with an empirical fit on $C_s=1/\sqrt{3}$7 cores

$C_s=1/\sqrt{3}$8

where $C_s=1/\sqrt{3}$9 is the number of LB steps per rendered image (Groen et al., 2012).

This model is used to define an operating regime rather than merely to fit benchmark curves. The reported efficiency sweet spot is $$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],$$0 sites per core, where compute and communication are balanced and near-peak SUPS per core are obtained. The associated operational recommendations are to keep $$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],$$1 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 $$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],$$2, 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 (Groen et al., 2012).

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 $$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],$$3, $$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],$$4, $$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],$$5, and $$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],$$6, obtained from least-squares calibration of pure-compute timings on Intel SandyBridge and AMD Interlagos systems (Groen et al., 2014).

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 (Groen et al., 2014).

The benefits and trade-offs are explicitly quantified. On a bifurcation geometry with $$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],$$7 sites and $$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],$$8 fill at $$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],$$9 cores, weighted decomposition reduces the compute-load imbalance from $\Omega=\Delta t/\tau$0 to $\Omega=\Delta t/\tau$1, a reduction of approximately $\Omega=\Delta t/\tau$2, and combining weighting with the space-filling-curve ordering gives $\Omega=\Delta t/\tau$3 with reduced communication time. On a much sparser aneurysm geometry with $\Omega=\Delta t/\tau$4 sites and $\Omega=\Delta t/\tau$5 fill at $\Omega=\Delta t/\tau$6 cores, weighting similarly reduces compute imbalance but can increase communication time on HECToR by approximately $\Omega=\Delta t/\tau$7–$\Omega=\Delta t/\tau$8, 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 (Groen et al., 2014).

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 $\Omega=\Delta t/\tau$9 coupling steps. The coupling strategy uses a $c_s^2=1/3$0 mapping between arterial outlets and venous inlets, assigns a random capillary pressure drop $c_s^2=1/3$1, enforces mass conservation $c_s^2=1/3$2, and exchanges only averaged $c_s^2=1/3$3 and $c_s^2=1/3$4 values, so the communication volume per synchronization is on the order of $c_s^2=1/3$5–$c_s^2=1/3$6 doubles (McCullough et al., 2020).

The same paper reports strong-scaling results at human scale. On Blue Waters, for a circle-of-Willis geometry with $c_s^2=1/3$7 lattice sites, scaling from $c_s^2=1/3$8 to $c_s^2=1/3$9 cores gives a speed-up of $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$0 for a $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$1 increase in core count, corresponding to $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$2 parallel efficiency; at $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$3 cores, the speed-up remains $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$4, about $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$5 efficiency. On SuperMUC-NG, scaling from $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$6 to $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$7 ranks on the same approximately $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$8-site model gives near-linear speed-up up to $$\rho=\sum_i f_i,\qquad \rho\,\mathbf{u}=\sum_i f_i\,\mathbf{c}_i,\qquad p=c_s^2\rho.$$9, with computational and communication efficiency both above $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$00, overall efficiency at or above $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$01, and $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$02 time steps completed in $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$03 seconds on $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$04 cores (McCullough et al., 2020).

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 (Zacharoudiou et al., 2022).

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 (Zacharoudiou et al., 2022).

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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$05 GPUs to $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$06 GPUs yields $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$07 efficiency, and an ExaPipe run on $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$08 GPUs, approximately two-thirds of Summit, yields $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$09 efficiency. Reported throughput is approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$10 MLUPS per node, dropping when the number of sites per GPU falls below $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$11. The paper also reports a GPU-versus-CPU speed-up of approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$12 at $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$13k 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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$14 million parallel tasks (Zacharoudiou et al., 2022).

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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$15 SUPS for a $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$16-site case and $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$17 SUPS for a $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$18-site case, with near-linear scaling up to four GPUs and efficiency at or above $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$19. The same work estimates that at approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$20 SUPS, a $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$21k-cell vascular model requires about $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$22 seconds per LB step and about $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$23 seconds for a $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$24-step cardiac cycle on four GTX 1080 Ti GPUs, while the visualization pipeline delivers frame rates above $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$25 frames per second for typical approximately $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$26k-voxel volumes (Esfahani et al., 2019).

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 $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$27 with

$$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$28

and uncertainty in quantities such as pressure and wall shear stress is visualized through Shannon entropy on marching-cubes topologies,

$$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$29

and an isosurface-crossing probability

$$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$30

The Frontier study reports near-ideal strong scaling up to more than $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$31 MPI ranks for domains of order $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$32 cells, constant MLUPS per GPU node under weak scaling, achieved per-node throughput on MI250X of order $$\partial_t \rho + \partial_\alpha(\rho u_\alpha)=0,$$33 MLUPS, and integration of sampled parameters such as inlet peak velocity and Smagorinsky constant directly into HemeLB configuration files (Xue et al., 21 Aug 2025).

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.