Large-scale lattice Boltzmann simulations of complex fluids: advances through the advent of computational grids

0.21 MB·English

by  J. Harting

Large-scale lattice Boltzmann simulations of complex
fluids: advances through the advent of computational
grids
BY JENS HARTING1, JONATHANCHIN2,MADDALENA VENTUROLI2,3, AND PETERV. COVENEY2
1InstituteforComputationalPhysics,UniversityofStuttgart,Pfaffenwaldring27,D-70569Stuttgart,
5 Germany
0
2CentreforComputationalScience,ChristopherIngoldLaboratories,UniversityCollegeLondon,20
0
GordonStreet,LondonWC1H0AJ,UK
2
3SchlumbergerCambridgeResearch,HighCross,MadingleyRoad,CambridgeCB30EL,UK
n
a
J

1
DuringthelasttwoyearstheRealityGridprojecthasallowedustobeoneofthefewscientificgroups
1
involvedinthedevelopmentofcomputationalgrids.Sincesmoothlyworkingproductiongridsarenotyet

] available,wehavebeenabletosubstantiallyinfluencethedirectionofsoftwaredevelopmentandgrid
C
deploymentwithintheproject.Inthispaperwereviewourresultsfromlargescalethree-dimensional
D
latticeBoltzmannsimulationsperformedoverthelasttwoyears.Wedescribehowtheproactiveuseof
. computationalsteeringandadvancedjobmigrationandvisualizationtechniquesenabledustodoour
s
c scientificworkmoreefficiently.Theprojectsreportedoninthispaperarestudiesofcomplexfluidflows
[
undershearorinporousmedia,aswellaslarge-scaleparametersearches,andstudiesofthe

1 self-organisationofliquidcubicmesophases.
v
Keywords:Lattice-Boltzmann,complexfluids,gridcomputing,computationalsteering
1
2
0
1 1. Introduction
0
5 In recent years there has emerged a class of fluid dynamical problems, called “complex fluids”, which
0 involvebothhydrodynamicfloweffectsandcomplexinteractionsbetweenfluidparticles.Computationally,
/
s suchproblemsare toolargeandexpensiveto tacklewith atomisticmethodssuch asmoleculardynamics,
c yettheyrequiretoomuchmoleculardetailforcontinuumNavier-Stokesapproaches.
:
v Algorithmswhichwork atan intermediateor “mesoscale” levelof descriptionin orderto solve these
i problemshavebeendevelopedinresponse,includingDissipativeParticleDynamics(Espa˜nolandWarren,
X
1995;Juryet al., 1999;Flekkøyetal., 2000),Lattice Gas Cellular Automata(RivetandBoon,2001),the
r
a StochasticRotationDynamicsofMalevanetsandKapral(MalevanetsandKapral,1998;Hashimotoetal.,
2000;Sakaietal.,2000),andtheLatticeBoltzmannEquation(Succi,2001;Benzietal.,1992;Loveetal.,
2003).Inparticular,theLatticeBoltzmannmethodhasbeenfoundhighlyusefulforsimulationofcomplex
fluidflowsinawidevarietyofsystems.Thisalgorithm,describedinmoredetailbelow,isextremelywell
suitedtoimplementationonparallelcomputers,whichpermitsverylargesystemstobesimulated,reaching
hitherto inaccessible physicalregimes. We describe some of these calculations, and also attempts to take
parallelcomputingtoanewscale,bycouplingseveralsupercomputerstogetherintoacomputationalgrid,
whichinturnpermitseasyuseoftechniquessuchascomputationalsteering,codemigration,andreal-time
visualization.
Article submitted to Royal Society TEXPaper
2 J.Harting,J.Chin,M.Venturoli,andP.V.Coveney
Theterm“simplefluid”usuallyreferstoafluidwhichcanbedescribedtoagooddegreeofapproxima-
tion bymacroscopicquantitiesonly,such asthe densityfield ρ(x), velocityfield v(x), andperhapstem-
peratureT(x).Suchfluidsaregovernedbythewell-knownNavier-Stokesequations(Faber,1995),which,
beingnonlinear,aredifficulttosolveinthemostgeneralcase,withtheresultthatnumericalsolutionofthe
equationshasbecomeacommontoolforunderstandingthebehaviourof“simple”fluids,suchaswateror
air.Conversely,a“complexfluid”isonewhosemacroscopicflowisaffectedbyitsmicroscopicproperties.
A good exampleof such a fluid is blood:as it flows throughvessels (of ordermillimetres wide and cen-
timetreslong),itissubjectedto shearforces,whichcauseredbloodcells (ofordermicrometreswide)to
alignwiththeflowsothattheycanslideoveroneanothermoreeasily,causingthefluidtobecomelessvis-
cous;thischangeinviscosityinturnaffectstheflowprofile.Hence,themacroscopicbloodflowisaffected
bythemicroscopicalignmentofitsconstituentcells.Otherexamplesofcomplexfluidsincludebiological
fluidssuchas milk,cellorganellesandcytoplasm,aswellas polymersandliquidcrystals. Inallof these
cases,thedensityandvelocityfieldsareinsufficienttodescribethefluidbehaviour,andinordertounder-
standthisbehaviour,itisnecessarytotreateffectswhichoccuroveraverywiderangeoflengthandtime
scales. This length and time scale gap makes complex fluids even more difficult to model than “simple”
fluids.Whilenumericalsolutionsofthemacroscopicequationsarepossibleformanysimplefluids,sucha
levelof descriptionmay notexistfor complexfluids, yetsimulationof everysingle moleculeinvolvedis
computationallyinfeasible.
Inamixturecontainingmanydifferentfluidcomponents,anamphiphileisakindofmoleculewhichis
composedof two parts, each partbeing attracted towardsa differentfluid component.For example,soap
moleculesareamphiphiles,containingaheadgroupwhichisattractedtowardswater,andatailwhichisat-
tractedtowardsoilandgrease;analogousmoleculescanalsobeformedfrompolymers.Ifmanyamphiphile
moleculesarecollectedtogetherinsolution,theycanexhibithighlyvariedandcomplicatedbehaviour,of-
ten assembling to form amphiphile mesophases, which are complex fluids of significant theoretical and
industrialimportance.Someofthesephaseshavelong-rangeorder,yetremainabletoflow,andarecalled
liquidcrystalmesophases.Ofparticularinteresttousarethosewithcubicsymmetry,whosepropertieshave
beenstudiedexperimentally(SeddonandTempler,1995;SeddonandTempler,1993;CzeslikandWinter,
2002)inlipid-watermixtures(SeddonandTempler,1995),diblockcopolymers(Shefelbineetal.,1999),and
inmanybiologicalsystems(Landh,1995).
Overthelastdecade,significantefforthasbeeninvestedinunderstandingcomplexfluidsthroughcom-
putationalmesoscale modellingtechniques.These techniquesdo notattemptto keeptrack of the state of
everysingleconstituentelementofasystem,nordotheyuseanentirelymacroscopicdescription;instead,an
intermediate,mesoscalemodelofthefluidisdeveloped,coarse-grainingmicroscopicinteractionsenough
thattheyarerenderedamenabletosimulationandanalysis,butnotsomuchthattheimportantdetailsare
lost. SuchapproachesincludeLattice GasAutomata(RivetandBoon,2001;Frisch etal.,1986;Rothman
andKeller,1988;Love,2002),theLatticeBoltzmannequation(Succi,2001;Benzietal.,1992;McNamara
andZanetti,1988;HigueraandJime´nez,1989;Higueraetal.,1989;ShanandChen,1993;Lamuraetal.,
1999;Chenetal.,2000;ChinandCoveney,2002),DissipativeParticleDynamics(HoogerbruggeandKoel-
man, 1992; Espa˜nol and Warren, 1995; Jury et al., 1999), or the Malevanets-Kapral Real-coded Lattice
Gas(Malevanetsand Kapral, 1998;Malevanetsand Yeomans,2000;Hashimoto et al., 2000;Sakai et al.,
2000).Recently-developedtechniques(Garciaetal.,1999;Delgado-BuscalioniandCoveney,2003)which
usehybridalgorithmshaveshownmuchpromise.
AllsimulationsdescribedinthispaperusethelatticeBoltzmannalgorithm,whichisapowerfulmethod
for simulating fluid dynamics. This is due to the ease with which boundary conditions can be imposed,
andwithwhichthemodelmaybeextendedtodescribemixturesofinteractingcomplexfluids.Ratherthan
tracking the state of individual atoms and molecules, the method describes the dynamics of the single-
particledistributionfunctionofmesoscopicfluidpackets.
Article submitted to Royal Society
Large-scalelatticeBoltzmannsimulationsofcomplexfluids:advancesthroughtheadventofcomputationalgrids3
Inacontinuumdescription,thesingle-particledistributionfunctionf (r,v,t)representsthedensityof
1
fluidparticleswithpositionrandvelocityvattimet,suchthatthedensityandvelocityofthemacroscop-
icallyobservablefluidaregivenbyρ(r,t) = f (r,v,t)dv andu(r,t) = f (r,v,t)vdv respectively.
1 1
In the non-interacting,long mean free path limit, with no externally applied forces, the evolution of this
R R
functionisdescribedbyBoltzmann’sequation,
(∂ +v·∇)f =Ω[f ]. (1.1)
t 1 1
While the left hand side describes changes in the distribution function due to free particle motion, the
righthandsidemodelspairwisecollisions.ThiscollisionoperatorΩisanintegralexpressionthatisoften
simplified(Bhatnagaretal.,1954)tothelinearBhatnagar-Gross-Krook(BGK)form
1
Ω[f]≃− f −f(eq) . (1.2)
τ
h i
TheBGKcollisionoperatordescribestherelaxation,ataratecontrolledbyacharacteristictimeτ,towards
aMaxwell-Boltzmannequilibriumdistributionf(eq).Whilethisisadrasticsimplification,itcanbeshown
thatdistributionsgovernedbytheBoltzmann-BGKequationconservemass,momentum,andenergy(Succi,
2001),andobeyanon-equilibriumformoftheSecondLawofThermodynamics(Liboff,1990).Moreover,it
canbeshown(ChapmanandCowling,1952;Liboff,1990)thatthewell-knownNavier-Stokesequationsfor
macroscopicfluidflowareobeyedoncoarselengthandtimescales(ChapmanandCowling,1952;Liboff,
1990).InalatticeBoltzmannformulation,thesingle-particledistributionfunctionisdiscretizedintimeand
space.Thepositionsronwhichf (r,v,t)isdefinedarerestrictedtopointsr onalattice,andthevelocities
1 i
v are restricted to a set c joiningpointson the lattice. The density of particlesat lattice site r travelling
i
with velocityc , at timestep t is givenby f (r,t) = f(r,c ,t), while the fluid’sdensity andvelocityare
i i i
givenbyρ(r)= f (r)andu(r)= f (r)c .Thediscretizeddescriptioncanbeevolvedintwosteps:
i i i i i
the collisionstep, whereparticlesat each lattice site are redistributedacrossthe velocityvectors,andthe
P P
advection,wherevaluesofthepost-collisionaldistributionfunctionarepropagatedtoadjacentlatticesites.
Bycombiningthesesteps,oneobtainsthelattice-Boltzmannequation(LBE)
1
f (r,t+1)−f (r,t)=Ω[f]=− [f (r,t)−N (ρ,u)], (1.3)
i i i i
τ
whereN =N (ρ(r),u(r))isapolynomialfunctionofthelocaldensityandvelocity,whichmaybefound
i i
bydiscretizingthewell-knownMaxwell-Boltzmannequilibriumdistribution.Ourimplementationusesthe
Shan-Chen approach (Shan and Chen, 1993), by incorporating an explicit forcing term in the collision
operatorinordertomodelmulticomponentinteractingfluids.ShanandChenextendedf totheformfσ,
i i
whereeachcomponentisdenotedbyadifferentvalueofthesuperscriptσ,sothatdensityandmomentum
ofacomponentσ aregivenbyρσ = fσ andρσuσ = fσc .Thefluidviscosityνσ isproportional
i i i i i
to(τσ −1/2)andtheparticlemassismσ.ThisresultsinalatticeBGKequation(1.3)oftheform
P P
1
fσ(r,t+1)−fσ(r,t)=− [fσ−N (ρσ,vσ)] (1.4)
i i τσ i i
Thevelocityvσ isfoundbycalculatingaweightedaveragevelocityu′ andthenaddingatermtoaccount
forexternalforces:
ρσ ρσ τσ
u′ = uσ / , vσ =u′+ Fσ. (1.5)
τσ τσ ρσ
! !
σ σ
X X
Inordertoproducenearest-neighbourinteractionsbetweencomponents,theforcetermassumestheform
Fσ =−ψσ(x) g ψσ¯(x+c )c , (1.6)
σσ¯ i i
σ¯ i
X X
Article submitted to Royal Society
4 J.Harting,J.Chin,M.Venturoli,andP.V.Coveney
whereψσ(x) = ψσ(ρσ(x))isaneffectivechargeforcomponentσ;g isacouplingconstantcontrolling
σσ¯
the strength of the interaction between two components σ and σ¯. If g is set to zero for σ = σ¯, and a
σσ¯
positivevalueforσ 6= σ¯ then,in theinterfacebetweenbulkregionsofeachcomponent,particlesexperi-
enceaforceinthedirectionawayfromtheinterface,producingimmiscibility.Intwo-componentsystems,
it is usuallythe case thatg = g = g . Amongstotherthings,this modelhasbeenused to simulate
σσ¯ σ¯σ br
spinodaldecomposition(ChinandCoveney,2002;Gonza´lez-Segredoetal.,2003),polymerblends(Martys
andDouglas,2001),liquid-gasphasetransitions(ShanandChen,1994),andflowinporousmedia(Mar-
tys and Chen, 1996).Amphiphilicfluids may be treated by introducinga new speciesof particle with an
orientationaldegreeoffreedom,whichismodelledbyavectordipolemomentd(Chenetal.,2000)with
magnituded .Thedipolefieldd(x,t)representstheaverageorientationofanyamphiphilepresentatsite
0
x.Duringadvection,valuesofd(x,t)arepropagatedaccordingto(tildesdenotepost-collisionvalues)
ρs(x,t+1)d(x,t+1)= f˜s(x−c ,t)d˜(x−c ,t), (1.7)
i i i
i
X
Duringcollision,thedipolemomentsevolveinaBGKprocesscontrolledbyadipolerelaxationtimeτ :
d
1
d˜(x,t)=d(x,t)− d(x,t)−d(eq)(x,t) . (1.8)
τ
d
h i
Theequilibriumdipolemomentd(eq) ≃ βd h/3isalignedwiththecolourfieldhwhichcontainsacom-
0
ponenthc duetocolouredparticles,andaparths dueto dipoles.Withqσ beinga colourcharge,suchas
+1forredparticles,−1forblueparticles,and0foramphiphileparticles,onegets
hc = qσ ρσ(x+c )c , (1.9)
i i
σ i
X X
hs(x,t)= fs(x+c ,t)θ ·d (x+c ,t)+fs(x,t)d (x,t) . (1.10)
 i i j i j i i 
i j6=0
X X
Thesecond-ranktensorθ isdefinedintermsoftheunittensorIandlatticevectorc asθ =I−Dc c /c2.
j j j j j
Inthepresenceofanamphiphilicspecies,theforceonacolouredparticleincludesanadditionaltermFσ,s
to accountforthe colourfield dueto the amphiphiles.By treatingan amphiphilicparticle asa pairof oil
andwaterparticleswithaverysmallseparationd,introducingaconstantg tocontrolthestrengthofthe
σs
interactionbetweenamphiphilesandnon-amphiphilesandTaylor-expandingind,itcanbeshownthatthis
termisgivenby
Fσ,s(x,t)=−2ψσ(x,t)g d˜(x+c ,t)·θ ψs(x+c ,t). (1.11)
σs i i i
i6=0
X
Whileamphiphilesdonotpossessanetcolourcharge,theyalsoexperienceaforceduetothecolourfield,
consistingofapartFs,c duetoordinaryspecies,andapartFs,s duetootheramphiphiles:
Fs,c =2ψs(x,t)d˜(x,t)· g θ ψσ(x+c ,t), (1.12)
σs i i
σ i6=0
X X
4D
Fs,s = − g ψs(x) d˜(x+c )·θ ·d˜(x)c (1.13)
c2 ss i i i
Xi n
+ d˜(x+c )d˜(x)+d˜(x)d˜(x+c ) ·c ψs(x+c ).
i i i i
h i o
Article submitted to Royal Society
Large-scalelatticeBoltzmannsimulationsofcomplexfluids:advancesthroughtheadventofcomputationalgrids5
Whiletheformoftheinteractionsseemsstraightforwardatamesoscopiclevel,itisessentiallyphenomeno-
logical, and it is not necessarily easy to relate the interaction scheme or its coupling constants to either
microscopic molecular characteristics, or to macroscopic phase behaviour. The phase behaviour can be
verydifficulttopredictbeforehandfromthesimulationparameters,andbrute-forceparametersearchesare
oftenresortedto(Boghosianetal.,2000).
2. Technical projects
Our three-dimensional lattice Boltzmann code, LB3D, is written in Fortran 90 and designed to run on
distributed-memoryparallelcomputers,usingMPIforcommunication.Ineachsimulation,thefluidisdis-
cretizedontoacuboidallattice,eachlatticepointcontaininginformationaboutthefluidinthecorresponding
regionofspace.Eachlatticesiterequiresaboutakilobyteofmemoryperlatticesitesothat,forexample,a
simulationona1283 latticewouldrequirearound2.2GBmemory.Thehigh-performancecomputingma-
chinesonwhichmostofthesimulationworkisperformedaretypicallyratherheavilyusedThesituation
frequentlyarisesthatwhileasimulationisrunningononemachine,CPUtimebecomesavailableonanother
machinewhichmaybeabletorunthejobfasterorcheaper.TheLB3Dprogramhastheabilityto“check-
point”itsentirestatetoafile.Thisfilecanthenbemovedtoanothermachine,andthesimulationrestarted
there,evenifthenewmachinehasadifferentnumberofCPUsorevenacompletelydifferentarchitecture.
Ithasbeenverifiedthatthesimulationresultsareindependentofthemachineonwhichthecalculationruns,
sothatasinglesimulationmaybemigratedbetweendifferentmachinesasnecessarywithoutaffectingits
output.Asaconservativeruleofthumb,thecoderunsatover104 latticesiteupdatespersecondperCPU
onafairlyrecentmachine,andhasbeenobservedtohaveroughlylinearscalinguptoorder103 compute
nodes. A 1283 simulation contains around 2.1×106 lattice sites; running it for 1000 timesteps requires
aboutanhourofrealtime,splitacross64CPUs.Thelargestsimulationweperformeduseda10243lattice.
Theoutputfromasimulationusuallytakestheformofasinglefloating-pointnumberforeachlatticesite,
representing,forexample,thedensityofaparticularfluidcomponentatthatsite.Therefore,adensityfield
snapshotfroma1283systemwouldproduceoutputfilesofaround8MB.Writingdatatodiskisoneofthe
bottlenecksinlargescalesimulations.Ifonesimulatesa10243system,eachdatafileis4GBinsize.LB3D
isabletobenefitfromtheparallelfilesystemsavailableonmanylargemachinestoday,byusingtheMPI-IO
based parallelHDF5 data format(HDF5, 2003).Our code is veryrobustregardingdifferentplatformsor
clusterinterconnects:evenwithmoderateinter-nodebandwidthsitachievesalmostlinearscalingforlarge
processorcountswiththeonlylimitationbeingtheavailablememorypernode.Theplatformsourcodehas
beensuccessfullyusedonincludevarioussupercomputersliketheIBMpSeries,SGIAltixandOrigin,Cray
T3E,CompaqAlphaclusters,NECSX6,aswellaslowcost32-and64-bitLinuxclusters.However,due
tocompilerormachinepeculiaritiesitisatimeconsumingtasktoachieveoptimumperformanceonmany
differentplatforms.PortingacomplexFortrancodelikeLB3Dtonewplatformsisoftenverydifficultand
time-consumingwithoutthe assistance of welltrainedstaff at the correspondingcomputercentres. Some
of theseproblemsare duetoportabilityissues withthe Fortranlanguage.Also,tuninga codeto takefull
advantageofthemachineonwhichitrunsrequiresconsiderableknowledgeofthelocalsystem’squirks.It
is hopedthatsome ofthe portabilityissues couldbe solvedin futurebywell-designedmiddleware.Such
issues include the fact that location, size, and duration of temporary filespace change from machine to
machine,asdothemethodsforinvokingcompilersandbatchqueues.
LB3D has successfully been used to study various problems like spinodal decomposition with and
withoutshear(Gonza´lez-Segredoetal.,2003;Hartingetal.,2004b),flowinporousmedia(Hartingetal.,
2004b),theself-assemblyofcubicmesophasessuchasthe’P’-phase(NekoveeandCoveney,2001a)inbi-
narywater-surfactantsystems,orthecubicgyroidphaseinternaryamphiphilicsystems(Gonza´lez-Segredo
and Coveney,2004b;Gonza´lez-Segredoand Coveney,2004a).Before we were able to take advantageof
Article submitted to Royal Society
6 J.Harting,J.Chin,M.Venturoli,andP.V.Coveney
computationalsteeringtechniques,ourworkusuallyinvolvedlargescale parametersearchesorganisedas
taskfarming jobs in order to find the areas of interest of the available parameter space. The technique of
computationalsteering(Chinetal.,2003;Brookeetal.,2003;Loveetal.,2003)hasbeenusedsuccessfully
insmaller-scalesimulationstooptimizeresourceusage.Typically,theprocedureforrunningasimulation
oftheself-assemblyofamesophasewouldbetosetuptheinitialconditions,andthensubmitabatchjob
to run for a certain, fixed number of timesteps. If the timescale for structural assembly is unknown then
the initial number of timesteps for which the simulation runs is, at best, an educated guess. It is not un-
common to examine the results of such a simulation once they return from the batch queue, only to find
thata simulationhas notbeenrunfor sufficienttime (in whichcase it mustbe tediouslyresubmitted),or
thatitranfortoolong,andthemajorityofthecomputertimewaswastedonsimulationofanuninteresting
equilibriumsystemshowingnodynamicalbehaviour.Anotherunfortunatescenariooftenoccurswhenthe
phasediagramofasimulatedsystemisnotwellknown,inwhichcaseasimulationmayevolveawayfrom
asituationofinterest,wastingfurtherCPU time.Computationalsteering,theabilitytowatchandcontrol
acalculationasitruns,canbeusedtoavoidthesedifficulties:asimulationwhichhasequilibratedmaybe
spotted and terminated, preventing wastage of CPU time. More powerfully,a simulation may be steered
throughparameterspaceuntilitisunambiguouslyseentobeproducinginterestingresults:thistechniqueis
verypowerfulwhensearchingforemergentphenomena,suchastheformationofsurfactantmicelles,which
arenotclearlyrelatedtotheunderlyingsimulationparameters.SteeringisperformedusingtheRealityGrid
steeringlibrarywhichhasbeendevelopedbycollaboratorsattheUniversityofManchester.Thelibrarywas
builtwiththeintentionofmakingitpossibletoaddsteeringcapabilitiestoexistingsimulationcodeswith
asfewchangesaspossible,andinasgeneralamanneraspossible.Oncetheapplicationhasinitializedthe
steering library and informedit which parametersare to be steered, then after everytimestep of the sim-
ulation,it is possibleto performtaskssuch as checkpointingthe simulation,saving outputdata,stopping
thesimulation,orrestartingfromanexistingcheckpoint.Whena steeredsimulationisstarted,a Steering
GridService(SGS)isalsocreated,torepresentthesteerablesimulationontheGrid.TheSGSpublishesits
locationtoaRegistryservice,sothatsteeringclientsmayfindit.Thisdesignmeansthatitispossiblefor
clientstodynamicallyattachtoanddetachfromrunningsimulations.
Successfulcomputationalsteeringrequiresthatthesimulationoperatorshaveagoodunderstandingof
what the simulation is doing, in real time: this in turn requiresgoodvisualization capabilities. Each run-
ningsimulationemitsoutputfilesaftercertainperiodsofsimulationtimehaveelapsed.Theperiodbetween
outputemission is initially determinedby guessing a timescale over which the simulation will change in
asubstantialway;however,thisperiodisasteerableparameter,sothattheoutputratecanbeadjustedfor
optimumvisualizationwithoutproducinganexcessiveamountofdata.TheLB3Dcodeitselfwillonlyemit
volumetric datasets as described above; these must then be rendered into a human-comprehensibleform
throughtechniquesincludingvolume-rendering,isosurfacing,ray-tracing,slice planes, andFouriertrans-
forms.Theprocessofproducingsuchcomprehensibledatafromtherawdatasetsisitselfcomputationally
intensive,particularlyifitistobeperformedinrealtime,asrequiredforcomputationalsteering.Forthis
reason, we use separate visualizationclustersto renderthe data. Outputvolumesare sent fromthe simu-
lation machine to the remote visualization machine, so that the simulation can proceed independentlyof
thevisualization;thesearethenrenderedusingtheopensourceVTK(Schroederetal.,2003)visualization
library into bitmap images, which can in turn be multicast over the AccessGrid, so that the state of the
simulationcanbeviewedbyscientistsaroundtheglobe.Inparticular,thiswasdemonstratedbyperforming
and interactingwith a simulation in frontof a live worldwide audience, as part of the SCGlobal track of
theSuperComputing2004conference.TheRealityGridsteeringarchitecturewasdesignedinasufficiently
generalmanner that visualization services can also be represented by Steering Grid Services: in order to
establishaconnectionbetweenthevisualizationprocessandthecorrespondingsimulation,thesimulation
SGScanbefoundthroughtheRegistry,andtheninterrogatedfortheinformationrequiredtoopenthelink.
Article submitted to Royal Society
Large-scalelatticeBoltzmannsimulationsofcomplexfluids:advancesthroughtheadventofcomputationalgrids7
In order to be able to deploy the above described components as part of a usable simulation Grid, a
substantialamountofcoordinationisnecessary,sothattheenduserisabletolaunchanentiresimulation
pipeline,containingmigratablesimulation,visualization,andsteeringcomponents,fromaunifiedinterface.
This requiresa system for keepingtrack of which services are available, which componentsare running,
taking care of the checkpoints and data which are generated, and to harmonize communication between
thedifferentcomponents.ThiswasachievedthroughthedevelopmentofaRegistryservice,implemented
using the OGSI::Lite (McKeown, 2003) toolkit. The RealityGrid steering library(Chin et al., 2003)
communicateswiththerestoftheGridbyexposingitselfasa“GridService”.ThroughtheRegistryservice,
steeringclientsareabletofind,dynamicallyattachto,communicatewith,anddetachfromsteeringservices
tocontrolasimulationorvisualizationprocess.
Largelatticesrequireahighlyscalablecode,accesstohighperformancecomputing,terascalestorage
facilitiesandhighperformancevisualisation.LB3D providesthefirstofthese, whiletheothersarebeing
delivered by the major computing centres. We expect to be able to run our simulations in an even more
efficientwayduetothesignificantworldwideeffortbeinginvestedinthedevelopmentofreliablecompu-
tationalgrids.Theseareacollectionofgeographicallydistributedanddynamicallyvaryingresources,each
providingservicessuchascomputecycles,visualization,storage,orevenexperimentalfacilities.Themajor
differencebetweencomputationalgridsandtraditionaldistributedcomputingisthetransparentsharingand
collective use of resources, which would otherwise be individualand isolated facilities. Perhaps at some
pointcomputationalgridswillofferinformationtechnologywhatelectricitygridsofferforotheraspectsof
ourdailylife:atransparentandreliableresourcethatiseasytouseandconformstocommonlyagreedstan-
dards(FosterandKesselman,1999;Bermanetal.,2003).Robustandsmartmiddlewarewillfindthebest
availableresourcesinatransparentwaywithouttheuserhavingtocareabouttheirlocation.Unfortunately,
reliableandrobustcomputationalgridsarenotavailableyet.Weusedvariousdifferentdemonstrationgrids
whichwereassembledespeciallyforagiveneventorwereintendedforuseasprototypingplatformsrather
thanusableproductiongrids.ThesemainlyincludedgridscouplingmajorcomputeresourcesintheUKand
thebiggestefforttookplacewithintheTeraGyroidproject(Picklesetal.,2004;Blakeetal.,2004)where
the mainmachinesofthe UK’snationalHPCcentreswerecoupledwith theTeraGridfacilitiesinthe US
throughacustomhigh-performancenetwork.Intotal,about5000CPUswerepartofthisgrid.Collaborative
steeringsessions withactiveparticipantsontwo continentsandobserversworldwidewere madepossible
throughthisapproach.
3. Scientific projects
(a) Complexfluidsundershear
Inmanyindustrialapplications,complexfluidsaresubjecttoshearforces.Forexample,axialbearings
areoftenfilledwithfluidtoreducefrictionandtransportheatawayfromthemostvulnerablepartsofthe
device.Itisveryimportanttounderstandhowthesefluidsbehaveunderhighshearforces,inordertobeable
tobuildreliablemachinesandchoosetheproperfluidfordifferentapplications.Inoursimulationsweuse
Lees-Edwardsboundaryconditions,whichwereoriginallydevelopedformoleculardynamicssimulations
in1972(LeesandEdwards,1972)andhavebeenusedinlatticeBoltzmannsimulationsbydifferentauthors
before(WagnerandYeomans,1999;WagnerandPagonabarraga,2002;Hartingetal.,2004b).Weapplied
ourmodeltostudy thebehaviourofbinaryimmiscibleandternaryamphiphilicfluidsunderconstantand
oscillatory shear. In the case of spinodal decomposition under constant shear, the first results have been
published in (Harting et al., 2004b). The phase separation of binary immiscible fluids without shear has
been studied in detail by different authors, and LB3D has been shown to model the underlying physics
successfully(Gonza´lez-Segredoetal.,2003).Inthenon-shearedstudiesofspinodaldecompositionithas
beenshownthatlatticesizesneedtobelargeinordertoovercomefinitesizeeffects:1283wastheminimum
Article submitted to Royal Society
8 J.Harting,J.Chin,M.Venturoli,andP.V.Coveney
acceptablenumberoflatticesites(Gonza´lez-Segredoetal.,2003).Forhighshearrates,systemsalsohave
tobeverylongbecause,ifthesystemistoosmall,thedomainsinterconnectacrossthez = 0andz = nz
boundariestoforminterconnectedlamellaeinthedirectionofshear.Suchartefactsneedtobeeliminated
from our simulations. Figure 1 shows an example from a simulation with lattice size 128x128x512.The
volumerenderedblueandredareasdepictthedifferentfluidspeciesandthearrowsdenotethedirectionof
shear.Inthecaseofternaryamphiphilicfluidmixturesundershearweareinterestedintheinfluenceofthe
presence of surfactantmoleculeson the phase separation.We also study the stress response and stability
of cubic mesophasessuch as the gyroidphase(Gonza´lez-SegredoandCoveney,2004b)or the “P”-phase
(Nekoveeand Coveney,2001b)undershear. Such complex fluidsare expectedto exhibitnon-Newtonian
properties(seebelow).Computationalsteeringhasturnedouttobeveryusefulforcheckingonfinitesize
effects during a sheared fluid simulation, since the human eye is extremely good at spotting the sort of
structures indicative of such effects. Implementingan algorithm to automatically recognize“unphysical”
behaviourisahighlynontrivialtaskincomparison.
Figure1.Spinodaldecompositionundershear.Differentlycolouredregionsdenotethemajorityofthecorresponding
fluid.Thearrowsdepictthemovementoftheshearedboundaries(movieavailableinonlineversion).
(b) Flowinporousmedia
Studyingtransportphenomenain porousmediais ofgreatinterestinfieldsrangingfromoilrecovery
and water purification to industrial processes like catalysis. In particular, the oilfield industry uses com-
plex, non-Newtonian,multicomponentfluids (containingpolymers, surfactants and/or colloids, brine, oil
and/or gas), for processes like fracturing, well stimulation and enhanced oil recovery. The rheology and
flow behaviourof these complex fluids in a rock is differentfrom their bulk properties. It is thereforeof
considerableinteresttobeableto characteriseandpredicttheflowofthese fluidsinporousmedia.From
thepointofviewofa modellingapproach,the treatmentofcomplexfluidsinthree-dimensionalcomplex
geometries is an ambitious goal since the lattice has to be large enough to resolve individual structures.
TheadvantageoflatticeBoltzmann(orlatticegas)techniquesisthatcomplexgeometriescanbemodelled
withease.SynchrotronbasedX-raymicrotomography(XMT)imagingtechniquesprovidehighresolution,
three-dimensionaldigitisedimagesofrocksamples.ByusingthelatticeBoltzmannapproachincombina-
tion with these high resolutionimagesof rocks, notonly is it possible to computemacroscopictransport
coefficients, such as the permeability of the medium, but information on local fields, such as velocity or
fluiddensities,canalsobeobtainedattheporescale,providingadetailedinsightintolocalflowcharacter-
isationandsupportingtheinterpretationofexperimentalmeasurements(Auzeraisetal.,1996).TheXMT
techniquemeasuresthelinear attenuationcoefficientfromwhichthe mineralconcentrationandcomposi-
tionoftherockcanbecomputed.Morphologicalpropertiesofthevoidspace,suchasporesizedistribution
Article submitted to Royal Society
Large-scalelatticeBoltzmannsimulationsofcomplexfluids:advancesthroughtheadventofcomputationalgrids9
Figure2.Renderingof4.9µmresolutionX-raymicrotomographicdataofa5123sampleofBentheimersandstone.
Theporespaceisshowninred,whiletherockisrepresentedinblue.
and tortuosity,can be derivedfromthe tomographicimage of the rock volume,and the permeabilityand
conductivityoftherockcanbecomputed(Spanneetal.,1994).Thetomographicdataarerepresentedbya
reflectivitygreyscalevalue,wherethelinearsizeofeachvoxelisdefinedbytheimagingresolution,which
isusuallyontheorderofmicrons.Byintroducingathresholdtodiscriminatebetweenporesitesandrock
sites, theseimagescanbereducedtoa binary(0’sand1’s)representationoftherockgeometry.Utilizing
the latticeBoltzmannmethod,single phaseormultiphaseflowcanthenbedescribedin theserealporous
media.
LatticeBoltzmannandlatticegastechniqueshavealreadybeenappliedtostudysingleandmultiphase
flow through three-dimensionalmicrotomographicreconstruction of porous media. For example, Martys
andChen(MartysandChen,1996)andFerre´olandRothman(Ferre´olandRothman,1995)studiedrelative
permeabilities of binary mixtures in Fontainebleausandstone. These studies validated the model and the
simulationtechniques,butwerelimitedtosmalllatticesizes,oftheorderof643.Simulatingfluidflowin
real rock samples allows us to compare simulation data with experimentalresults obtained on the same,
or similar, pieces of rock. For a reasonable comparison, the size of the rock used in lattice Boltzmann
simulations should be of the same order of magnitude as the system used in the experiments, or at least
largeenoughto capturetherock’stopologicalfeatures.Themoreinhomogeneousthe rock,thelargerthe
samplesizeneedstobeinordertodescribethecorrectporedistributionandconnectivity.Anotherreason
for needingto use large lattice sizes is the influence of boundaryconditionsand lattice resolution on the
accuracy of the lattice Boltzmann method. It has been shown (see for example (He et al., 1997), (Chen
andDoolen,1998)andreferencestherein)thattheBhatnagar-Gross-Krook(BGK)(Bhatnagaretal.,1954)
approximationof the lattice Boltzmann equation which is commonly used causes so-called bounce-back
boundariesto becomeinaccurate,resultingineffectssuchasthe computedpermeabilitybeinga function
of the viscosity.Thiseffectcan be limited byloweringthe viscosityand increasingthelattice resolution.
ToaccuratelydescribehydrodynamicbehaviourusinglatticeBoltzmannsimulations,theKnudsennumber,
whichrepresentstheratioofthemeanfreepathofthefluidparticlesandthecharacteristiclengthscaleof
thesystem(suchastheporediameter),hastobesmall.Iftheporesareresolvedwithaninsufficientnumber
oflatticepoints,finitesizeeffectsarise,leadingtoaninaccuratedescriptionoftheflowfield.Inpractice,
atleastfivetotenlatticesitesareneededtoresolveasinglepore.Therefore,inordertobeabletosimulate
realisticsamplesizes,weneedlargelatticesoftheorderof5123.
UsingLB3D,weare abletosimulatedrainageandimbibitionprocessesina 5123 subsampleofBen-
theimer sandstone X-ray tomographicdata. The whole set of XMT data representedthe image of a Ben-
theimersampleofcylindricalshapewithdiameter4mmandlength3mm.TheXMTdatawereobtainedat
theEuropeanSynchrotronResearchFacility(Grenoble)ataresolutionof4.9µm,resultinginadatasetof
approximately816x816x612voxels.Figure2showsasnapshotofthe5123 subsystem.We comparesim-
Article submitted to Royal Society
10 J.Harting,J.Chin,M.Venturoli,andP.V.Coveney
Figure3.Anoriginallyfullyfluidsaturatedrockisbeinginvadedbyanotherimmisciblefluidusingabodyforcegaccn
=0.003.Theoilslowlypushestheotherfluidcomponentoutoftherockporesuntiltherockisfullysaturatedbyoilat
t=30000.Forbettervisabilityonlytheinvadingfluidisshown(movieavailableinonlineversion).
ulated velocity distributionswith experimentallyobtained magneticresonanceimaging (MRI) data of oil
andbrineinfiltrationintosaturatedBentheimerrockcore(Sheppardetal.,2003).Therocksampleusedin
theseMRIexperimentshadadiameterof38mmandwas70mmlongandwasimagedwitharesolutionof
280microns.Thesystemsimulatedwassmaller,butstillofasimilarorderofmagnitudeandlargeenough
torepresenttherockgeometry.Ontheotherhand,thehigherspaceresolutionprovidedbythesimulations
allowsadetailedcharacterisationoftheflowfieldintheporespace,henceprovidingausefultooltointer-
prettheMRIexperiments,forexampleinidentifyingregionsofstagnantfluid.Figure3showsanexample
fromabinaryinvasionstudy.Arockwhichisinitiallyfullysaturatedwith“water”(blue),isbeinginvaded
by “oil”(red)fromtherightside. Thelattice size is5123 andthe forcinglevelis setto g = 0.003.In
accn
figure3,onlytheinvadingfluidcomponentisshown,i.e.onlyareaswhereoilisthemajoritycomponent
arerendered.Periodicboundaryconditionsareapplied,andfluidleavingthesystemontheleftsideiscon-
verted to oil before re-enteringon the opposite side. After 5000timesteps, the oil has invadedaboutone
quarterofthesystemalreadyandafter25000timestepsonlysmallregionsoftherockporespacearestill
filledwithwater.After30000timesteps,thewatercomponenthasbeenfullypushedoutoftherock.This
exampleonlycoversbinary(oil/water)mixturesofNewtonianfluids,sincethisisafirstandnecessarystep
intheunderstandingofmultiphasefluidflowinporousmedia(Hartingetal.,2004b).However,weareable
to study the flow of binary immiscible fluids with an additionalamphiphiliccomponentin porousmedia
andexpectresultstobepresentedelsewhereinthenearfuture.
(c) Thecubicgyroidmesophase
ItwasrecentlyshownbyGonza´lezandCoveney(Gonza´lez-SegredoandCoveney,2004b)thatthedy-
namical self-assembly of a particular amphiphile mesophase, the gyroid, can be modelled using the lat-
tice Boltzmann method. This mesophase was observed to form from a homogeneous mixture, without
any external constraints imposed to bring about the gyroid geometry,which is an emergenteffect of the
mesoscopic fluid parameters. It is important to note that this method allows examination of the dynam-
ics of mesophase formation, since most treatments to date have focussed on properties or mathematical
Article submitted to Royal Society