Hemodynamics withlattice Boltzmann
Rupert Nashr.nash@epcc.ed.ac.uk@rupe_nash

The lattice Boltzmann
methodA lightning summary�.�molecular�edynamics
FHP hexagonallatticefromFrisch
et ofal.points.
[��].EachShownishasonmodels [��], which
restricts
particles
to a latticegrid point3
�nalstates(right)
arechosen
betweenrandomly
withequlinks lower
to itsneighbours
along
which
theparticles
stream
without
interaction
duringone
discretetime
LB prehistory - lattice gas cellular�e update rules for the occupation numbers n i (x�, t) (where istep, each link being either occupied or vacant. �e particles then collide moving to di�erent links atautomata(1980’s)
their new site.

can be written as
pdate rules for the occupation numbers n i (x�, t) (where i inde
∆ i n i = C i (n�)be written ase ∆ i n i ≡ n i (x� + c�i ∆t, t + ∆t) − n i (∆x�i,nt)i =represents
C i (n�) the free pion operator which obeys conservation of mass and momenn i ≡ n i (x� + c�i ∆t, t + ∆t) − n i (x�, t) represents the free propa
el has extremely simple dynamics, yet when coarse grained
operator which obeys conservation of mass and momentum
Figure �.� equations.
�e FHP
hexagonal
latticefromFrischet al. [��].
Shown is one
collision,
wherethe upper andsho
er-Stokes
�atthiscandisplay
thecomplexities
One
collisionrule– chooserandomly
betweentwo
outputstateslower �nal states (right) are chosen between randomly with equal probability.
as extremely simple dynamics, yet when coarse grained it aactions at the molecular scale of a�real �uid is somewhat surpris
update rules for the occupation numbers n (x , t) (where i indexes the list of discrete velocities
okes�eequations.
�at this can display the complexities shownietrized by the local density ρ, velocity u� and tempera
rium f (0)Boltzmann, parametrized by theBGK
local density
ρ, velocity u� and temperature T, giving theEquationn [��, ��]:
11(0)
(0)
��+v
⋅∇+g⋅
∇f=−(f−f),[∂]
��txc[∂ t + v� ⋅ ∇x� + g� ⋅ ∇�c ] f =τ − ( f − f ),τD�24
mplicated, nonlinear, integral term with relaxation on a characteristic time τ towards
thef(0)mm= ρ�
�c� − u��2 � .
� exp �−2πk B TD�22k B Tmm
2f = ρ��c� − u�� � .
� exp �−hydrodynamic variables
at a pointvelocity
moments of the distrib
2πk2kTB T are the (microscopic)
B3(0)∫ f d c�,ρu = ∫ f c d c�,��ariables at a point are the (microscopic) velocity mom
S = ∫ f (c c − c δ ) d c�;

n:

ρ=ααβ

α3α β2s αβ

3ely the density, momentum density and deviatoric stress tensor. Note that the N

o the same order as the Maxwellian, we have5(u α g β + g α u β )Q i α βDiscretiseinvelocity
g α c ispaceαΦ (x�, t) ≡ [ g� ⋅ ∇ ] = −ρw �+��
c iii

c s22c s4• Reduce continuous three dimensional space to discrete

set of velocities

g all these results gives the discrete velocity Boltzmann equation• DVBE:∂ t f i + c�i ⋅ ∇ f i + [ g� ⋅ ∇�c f ] i = −(0)( fi − fi )τ.

mote the singleForcerelaxation
term too timelong operator on the right-hand side of Equafor slides!

��η s = ρc s2 τ s2
η b = ρc s2 τ b .
3(�.��a)(�.��b)Everyone uses simple velocity sets
�e LBM can be extended to incorporate thermal �uctuations. An algorithm due to Ladd [��] asextended by Adhikari et al. [�] will be used here. For the simulations in this work, the three dimensional,
��een velocity (D�Q��) velocity set, shown in Figure �.�, will be used, in the implementation described• Usually(ornoted.
D3Q15 or D3Q27)
in AppendixD3Q19
A, unless otherwise�.� �e D�Q��
set is based on a simple
cubic mesh of points. �e velocities are: one zero vector;• OK Figure
as long
asvelocity
symmetric“enough”
six with �c��2 ≡ 1 along the unit vectors, and eight with �c��2 = 3 along the [1, 1, 1] type directions.

• Introduces more errors but still second order in Δx�.� S������������� �� �������������
• Except
sometimes
for turbulence?
• White& Chong, J. Comput. Physics (2011); S. K. Kang and Y. A.
Singular distributions of forces are widely used in theoretical models of �uid immersed boundaries;
Hassan, J. Comput. Phys. 232, 100 (2012).
for example an interface between gas and liquid phases can be described as a two-dimensional surface,a polymer as a one-dimensional curve [��] or a moving particle at distances large compared to its
radius as a point [��]. In particular, at low Reynolds number the velocity �eld around a body can be

6

7DVBE ⇒ LBE• Integrate along each velocity with the trapezium rule forone timestep• Make simple substitution get explicit equationsf¯i (~x + ~ci t, t +t)

f¯i (~x, t) =t⌧+

t/2(f¯i (~x, t)f¯i0 (~x, t))

8LB vs. CFD

Boltzmann
Eqn

Navier-Stokes
Equations
Discretevelocity BE

FD/FEM/FVM

Schemerecovers NS
as dx ->0LBE

Schemerecovers NS
as dx->0

9LB is strong• Possible to introduce more physics rigorously*- Particles- Multicomponent, multiphase fluids- Liquid crystals• Meshing is much simpler than FEM (only 1 grid resolutionto set)
• Parallel implementation (P<1000) easy!
• Parallel implementation (P<=250,000, efficiency >90%)
possible!
• No global communications required (unlike e.g. Poissonsolver)

10

LB is weak

1 xpc=s• Speed of sound fixedt
3- Need to carefully chose params such that Mach number <<1- Mach number important because it is the parameter in theChapman Enskog expansion- Sound takes many time steps to cross domain- Pressure takes many time steps to converge- Is there really an advantage over just solving the Poisson equation
as in most NSE based codes?

11

LB is horrible• It struggles to get to high velocities (positivity ofdistribution functions)• It struggles to get to low viscosity as has poor stability• Hence high Reynolds number requires advances
techniques: cascaded LB

HemeLB applications

(a)

13

CFD (Colour for doctors)
sprouting front

Bvelocity (m s–1)
0.04

DDB

AE

0.01

0.001

C(b)

WSS (Pa)

sprouting front

20100Figure 8. P5B simulation results: (a) velocity magnitude plotted on a cross section along the z ¼ 0 plane. Velocity shows the expected parabolic profile acrovessel diameter. Velocity is higher in the artery located at the centre of the domain, in particular close to the optic disc. Velocity magnitude quickly decreases
artery progresses towards the sprouting front and it stops being a preferential flow path at the points where its identity stops being clearly defined. (bmagnitude plotted on the model surface. Areas of preferential flow tend to experience highest WSS magnitudes. WSS is generally low across the dexcept for the arterial segment close to the optic disc and some first-order branches. WSS values higher than 20 Pa are considered unphysiological and the rexperiencing them are coloured in black. Black circles indicate regions of interest referenced in the manuscript.

14Angiogenesis(a)Downloaded from rsif.royalsocietypublishin(b)

Figureretinal
1. Murine
retinal (a)vascularplexusand2121days
daysfrom
postnatal
Mouse
networks
6 daysand 6(b)birth (panel (a)andnext(b),
Withindays, theofprimitive
vesselpathways
network remodels(Micedidrespectively).
not survive)The
stepis pruning– removal
redundant– does into
hemodynamics
play aSamplesrole? were collected, mounted and imaged as
mature vasculature.
described in §3.

sksu
Fimop(Ctosubu
m
arord

15

Pruning in action?

Downloaded from rsif.royalsocietyp

WSS (Pa)
20100

Figure 10. Traction vectors (of constant length and coloured
according to~50 µmmagnitude) on the luminal surface of the region of interest highlighted infigure 9b. The loop branch undergoing regression (upper branch) experiences

16

Download HemeLB if you want a go
• Started at UCL• Developers now at UCL, Edinburgh, Brunel, Clemson• C++, object oriented, heavy use of templates• Source code on GitHub:- https://github.com/UCL/hemelb• Docker image available too:- https://hub.docker.com/hemelb/hemelb

17Performance – Site updates per second
1011 7

10S/s 1
S/Ncores /s1010

109
SmallOriginal

108107 110106 210NBC
Original
DMAPPNBCImmediate
102103

DMAPPImmediate
310Cores104Nsites /Ncores

MedSmallMed

104

105

LargeLarge

105

106

Parallel in TimeParallelism in Time
• Parallel-in-time integration can allow the use of even moreconcurrency to reduce wall clock time
- Note: will always increase CPU time

• Parareal is the most commonly used PinT method• Has has been used with LB [Randles 2014, J Comp Phys]

19

Parareal• Treat your simulation as an initial value problem• Have a fine operator (F) and coarse one (G).
- Cost of applying G must be << F
- G can be less accurate- G can be on a coarser grid
• Apply G several time to advance to target time
• From intermediate points start multiple F simulations inparallel• Correct and iterate

20Parareal21

Parareal for HemeLB• Collaboration with Derek Groen (Computer Science, UoBLondon) and Daniel Ruprecht (Mathematics, UoLeeds)• Project began October• Will implement the fine and coarse operators with existingsimulation• Are working on coarsening and refinement operations- Requires mapping between two grids in parallel with different

decompositions• Unclear how well this will work for lattice-Boltzmann

22Public cloud for HPC?
Experience with HemeLB on AWS and Azure

Comparison: policy
Public cloud

HPC centre• Pay provider through• Apply for time throughgrant• On demand (minutesmonths)• User defines security

grant process• Batch system (minutesdays)• Centre defines security- You have root on VM- You do NOT have root

• Many services (Web, DB,• Limited services

serverless, ML, analytics)(filesystem, batch, pre/post-processing)24Comparison: hardware
Public cloud

HPC centre• Many node types• 1-3 node types• New servers constantly• New hardware every 2-5

being added• Network:years• Network:- Ethernet- Specialised low-latency high-bandwidth- (Some Azure nodes haveInfiniband)

25

Comparison: storage
Public cloud

HPC centre• Nodes usually have some

• Nodes might have locallocal• Object storage• Various databasesstorage (NVRAM?)• Shared parallel POSIX
filesystem- Relational- Columnar- NoSQL- Graph26

Comparison: build
Cloud (Azure)

HPC (Archer)• yum install ...

• module load ...

• Configure NFS and MPI• cmake• Compile dependencies• make• cmake• qsub• make install• Save image• Start pool of nodes• Submit job• Stop pool27

HemeLB performance28

Simulation cost

29

AWS Performance

Time / S

100

10Cirrus
AWS-EC2
ARCHER

110100Cores1000

AWS Cost
87

Simulation cost / GBP

65Cirrus

4AWS-EC2
ARCHER (partner)

32

10163264128
Cores256

5121024

Acknowledgements• HemeLB was started at UCL (Peter Coveney)• Major code contributors:- RWN, Miguel Bernabeu, Derek Groen, Sebastian Schmieschek,Hywel Carver, Jens Nielsen, Mayeul d’Avezac, JamesHetherington, Marco Mazzeo, Steven Manos• AWS study by Steven Steven• Azure credit provided by MS Research• All errors are mine