Parallel multiscale simulations of a brain aneurysm

17 Pages·2013·3.59 MB·English

by  Leopold Grinberg

JournalofComputationalPhysics244(2013)131–147
ContentslistsavailableatSciVerseScienceDirect
Journal of Computational Physics
journal homepage: www.elsevier.com/locate/jcp
Parallel multiscale simulations of a brain aneurysm
Leopold Grinberga, Dmitry A. Fedosovb, George Em Karniadakisa,⇑
aDivisionofAppliedMathematics,BrownUniversity,Providence,RI02912,USA
bInstituteofComplexSystemsandInstituteforAdvancedSimulation,ForschungszentrumJülich,Jülich52425,Germany
a r t i c l e i n f o a b s t r a c t
Articlehistory: Cardiovascularpathologies,suchasabrainaneurysm,areaffectedbytheglobalbloodcir-
Availableonline1September2012 culationaswellasbythelocalmicrorheology.Hence,developingcomputationalmodels
forsuchcasesrequiresthecouplingofdisparatespatialandtemporalscalesoftengoverned
Keywords: bydiversemathematicaldescriptions,e.g., bypartialdifferential equations (continuum)
Thrombosis andordinarydifferentialequationsfordiscreteparticles(atomistic).However,interfacing
Bloodmicrorheology atomistic-based with continuum-based domain discretizations is a challenging problem
Parallelcomputing
thatrequiresbothmathematicalandcomputationaladvances.Wepresenthereahybrid
Atomistic-continuumcoupling
methodologythatenabledustoperformthefirstmultiscalesimulationsofplateletdepo-
Dissipativeparticledynamics
sitionsonthewallofabrainaneurysm.Thelargescaleflowfeaturesintheintracranialnet-
Spectralelements
Domaindecomposition work are accurately resolved by using the high-order spectral element Navier–Stokes
solverNejTar.Thebloodrheologyinsidetheaneurysmismodeledusingacoarse-grained
stochastic molecular dynamics approach (the dissipative particle dynamics method)
implementedintheparallelcodeLAMMPS.Thecontinuumandatomisticdomainsoverlap
withinterfaceconditionsprovidedbyeffectiveforcescomputedadaptivelytoensurecon-
tinuityofstatesacrosstheinterfaceboundary.Atwo-wayinteractionisallowedwiththe
time-evolving boundary of the (deposited) platelet clusters tracked by an immersed
boundarymethod.Thecorrespondingheterogeneoussolvers(NejTarandLAMMPS)are
linkedtogetherbyacomputationalmultilevel message passing interface thatfacilitates
modularityandhighparallelefficiency.Resultsofmultiscalesimulationsofclotformation
insidetheaneurysminapatient-specificarterialtreearepresented.Wealsodiscussthe
computational challenges involved and present scalability results of our coupled solver
onupto300Kcomputerprocessors.Validationofsuchcoupledatomistic-continuummod-
elsisamainopenissuethathastobeaddressedinfuturework.
(cid:2)2012ElsevierInc.Allrightsreserved.
1.Introduction
Cerebralaneurysmsoccurinupto5%ofthegeneralpopulation,witharelativelyhighpotentialforruptureleadingto
strokesinabout30,000Americanseachyear[1].Therearenoquantitativetoolstopredicttheruptureofaneurysmsand
noconsensusexistsamongmedicaldoctorswhenexactlytooperateonpatientswithcerebralaneurysms.
Thebiologicalprocessesprecedingtheaneurysmrupturearenotwellunderstood.Thereareseveraltheoriesrelatingthe
wallshearstresspatternsandpressuredistributionwithinaneurysmandpathologicalchangesatcellularleveloccurringat
thearterialwalllayers.Realisticsimulationofsuchprocessesmustbebasedonresolvingconcurrentlythemacro-(centime-
ter)aswellasthemicro-(sub-micron)scaleflowfeatures,andalsotheinteractionofbloodcellswiththeendothelialcells
⇑
Correspondingauthor.
E-mailaddresses:[email protected],[email protected](G.E.Karniadakis).
0021-9991/$-seefrontmatter(cid:2)2012ElsevierInc.Allrightsreserved.
http://dx.doi.org/10.1016/j.jcp.2012.08.023
132 L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147
formingtheinnerlayerofthearterialwall.Suchcomplexprocessesareclearlymultiscaleinnaturenecessitatingtheuseof
differentmathematicalmodelstoresolveeachscale.
Inthecurrentwork,weconsidertheinitialformationofaplateletclotinbrainaneurysmasarepresentativeexample.
Clotsformnotonlyinbrainaneurysmsbutalsoinaorticaneuryms[2],coronaryarteries[3],veins[4],etc.Themainfocus
ofthispaperisthegeneralframeworkweproposeformultiscalemodelingofarterialbloodflow,includingmodelingofthe
initialthrombusformation.Thismultiscalemethodologyisdevelopedbasedonatomisticandcontinuumdescriptions,and
canbeusedtosimulate,forinstance,clotformationandgrowthinaorticaneurysmsorincarotidarteries.Somealgorithmic
andcomputationalaspectsofthecoupledatomistic-continuummodelcanalsobeusedinotherresearchareas,e.g.inmate-
rialsmodelingandsolidmechanics.Inthefollowingwereview(a)physiologicalaspectsrelatedtocerebralaneurysms,and
(b)recentadvancesandopenissuesinmultiscalemodeling.
1.1.Cerebralaneurysm:overview
Cerebralaneurysmsarepathological,blood-filledpermanentdilationsofintracranialbloodvesselsusuallylocatednear
bifurcationsintheCircleofWillis(CoW)[5–8],seeFig.1.Currently,prospectiveruptureratesareculledfromvariousstud-
ies,buttheseratesdonottakeintoaccountpatient-specificanatomicandphysiologicfactorswhichmayaffectrupturerisk.
Theselimitationsaremoremagnifiedwhenconsideringoneparadoxthatprevailswhenanalyzingruptureratesfromthe
literature:The prospective rate of rupture ofanterior circulationaneurysmsless than 7mm in size is extremelylow [7],
butthemajorityofrupturedaneurysmsinmostseriesareanteriorcirculationaneurysmslessthan7mminsize[10].Anal-
ysisofpatient-specificanatomicandphysiologicdatamaybeusedtoidentifyspecificpatientswhomightbeathigherrisk
forrupture.
Theprocessesofinitiation,growthandruptureofCAsarenotwellunderstood,andseveral–oftencontrasting–hypoth-
esesonthemodifiedhemodynamicshaveemergedinrecentyears,e.g.seetherecentreviewsin[5,6,8,11]andreferences
therein.Whileonewouldexpectthataneurysmruptureisassociatedwithhighpressureorhighwallshearstress(WSS)
magnitude,thisisnotthecasewithCAsasnoelevatedpeakpressureandtypicallylowWSSvaluesareobservedwithin
theaneurysm.
Histologicalobservationshaverevealedadegenerationofendothelialcells(EC)anddegradationoftheintracellularmatrix
ofthearterialwallsduetodecreasingdensityofsmoothmusclecells(SMCs)[12,13].WSSisrelatedtotheendothelialgene
expression,andinlaminarflowproducesaquiescentphenotypeprotectingfrominflammationorcellapoptosis.However,
flowinstabilitiesandoscillatoryWSScantriggercertaingenetictraitsthataffecttheelasticpropertiesofthearterialwall
[14,15].Forexample,recentstudiesreportedthatlowWSSlevelsinoscillatoryflowcancauseirregularECpatterns,poten-
tiallyswitchingfromanatheroprotectivetoanatherogenicphenotype[5,15–17].Ontheotherhand,lowWSScanalsobe
protectiveasitleadstothickeningofarterialwalls,whichthenbecomesmoretoleranttomechanicalloads.Hence,bysimply
examiningthemagnitudeofWSSwecannotarriveataconsistenttheoryofaneurysmgrowthorrupture.Tothisend,exper-
imentalstudieswithECssubjectedtoimpingingflowin[15,18]haveshowntheimportanceofWSSgradients(WSSG)(rather
thanthemagnitude)onthemigrationofECsdownstreamofthestagnationregion.
PerhapsthemostdramaticevidenceontheroleofECsandtheirinteractionwithbloodcellshasbeendocumentedina
Japaneseclinicalstudycomparingrupturedandunrupturedaneurysms[9].Rupturedaneurysmsexhibitedsignificantendo-
thelialdamageandinflammatorycellinvasioncomparedwithunrupturedaneurysms;seeFig.1right)foratypicalimage.
WeobservethattheEClayerisdrasticallyalteredandcoveredwithbloodcellsandafibrinnetwork.Similarly,itwasre-
portedin[9]thatleukocytes(WhiteCells,WCs)inthewallcouldbeassociatedwithsubarachnoidhemorrhage(SAH),and
endothelialerosionenhancesleukocyteinvasionofthewallbeforerupture.
Fig.1. Bloodflowinthebrainisamultiscaleproblem.Left:themacrodomainwherethelarge-scaleflowdynamicsismodeledbytheNavier–Stokes
equations;differentcolorscorrespondtodifferentcomputationalpatches(CourtesyofProf.J.R.Madsen,HarvardMedicalSchool).Shownintheinsetarethe
microdomainswheredissipativeparticledynamicsisappliedtomodelthemicro-scalefeatures.Right:Scanningelectronmicroscopeimageoftheinside
layerofarupturedMCAaneurysmshowingthedisruptedpatternofendothelialcellsandbloodcellsadheredtotheinter-endothelialcellgaps[9].(For
interpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147 133
DamagedECtypicallyinitiatesathrombus(bloodclot)formationprocess.Oneofthemostimportantbuildingblocksin
theclotisa2–3lminsizebloodcellcalledplatelet.Theroleofplateletshasnotbeenfullyexploredinexperimentalwork,
but clinical tests have documented the existence of spontaneous thrombosis in giant aneurysms [19–21] due to platelet
deposition.Thrombusformationwithintheaneurysmisnon-uniformduetothecomplexflowpatternsandtheinteraction
ofplateletswiththedamagedEClayer[22].Similarly,itwasfoundinendovascularstudiesthatcreationofspontaneous
thrombosisoccursafterstentplacement[23].Inotherstudies,enhancedplateletaggregabilityhasbeenobservedincerebral
vasospasmfollowinganeurysmalSAH[24].B.FurieandB.C.Furiereviewedtwoindependentpathwaystothrombusforma-
tion[25].Plateletactivationmaystartduetoexposureofbloodcellstoeithersubendothelialcollagenorthrombingenerated
bytissuefactorderivedfromthevesselwall.Theauthorsalsoemphasizedthatlargescaleflowfeaturessuchasshearand
turbulenceaffectthethrombusformationprocess.
Akeyquestionregardingtheprocessofaneurysmprogressionishowtocorrelatethebloodcelldynamicsandinterac-
tions to the complex flow patterns observed within the aneurysms [26]. To this end, there are currently two schools of
thought:high-floweffectsandlow-floweffects[5].TheformersuggeststhatlocalizedelevatedWSScauseendothelialinjury
thatcausesanoverexpressionofNOproduction,whichinturncanleadtoapoptosisofSMCsandsubsequentwallweaken-
ing.Thelatterpointstothestagnationtypeflowinthedomeoftheaneurysm,whichalsocausesirregularproductionofNO.
Thisdysfunctionofflow-inducedNOleadstotheaggregationofredbloodcellsandtheaccumulationandadhesionofplate-
letsandleukocytesalongtheintimalsurface(thefirstlayerofthearterialwall).This,inturn,maydamagetheintima,allow-
ingforinfiltrationofwhitebloodcells(WBC)andfibrin[27].Weshouldalsonotethattheresponseofdamagedaneurysm
endotheliumtostimulisuchasWSSmaybedifferentfromtheresponseofhealthyendothelialtissue.Tothisend,toadvance
ourunderstanding,weneedtoinvestigatescenariosbeyondthepuremechanicalpointofview.Inparticular,weneedto
understandtheroleofendothelialcellsandtheirinteractionswiththebloodcells,i.e.,platelets,RBCsandWBCs.
1.2.Multi-scalemodelingofbloodflow
Modelingthebloodflowasamultiscalephenomenon,usingcoupledcon-tinuum-atomisticmodels,isessentialtobetter
understandthethrombusformationprocess.Seamlessintegrationofheterogeneouscomputercodesbasedoncontinuum
models with codes that implement atomistic-level descriptions is key to the successful realization of parallel multiscale
modelingofrealisticphysicalandbiologicalsystems.
Multiscalehybridapproaches[28–35]canpotentiallyprovideanelegantsolutionfornon-feasiblemicro-scalesimula-
tionsperformedinverylargecomputationaldomains.Themainchallengeindevelopingefficienthybridmethods[36]isfor-
mulatingrobustinterfaceconditionscouplingmultipledescriptions(e.g.,continuum,atomistic,mesoscopic)suchthatthe
mainphysicalquantities(mass,momentum,andenergy)andtheirfluxesareconservedatthedomaininterface.Ingeneral,
mostexistinghybridmethodscanbeclassifiedintotwomaincategoriesofcouplingapproaches[36]:(i)thestate-exchange
method[28,32–34,37,35]and(ii)theflux-exchangemethod[29–31]ortheircombination.
Thestate-exchangemethodreliesonthestateinformation(e.g.,boundaryvelocity)whichistransferredbetweenvarious
descriptionsacrosstheinterfacesof(inmostcases)overlappingdomains.Intheoverlappingregionallparticipatingdescrip-
tionsmustbevalidandtheconstraineddynamicsisoftenimposedviaadynamicrelaxationtechnique[34].Thestate-ex-
changemethodiscloselyrelatedtothealternatingSchwarzmethod[38],wheretheintegrationincontinuumandatomistic
domainsisperformedinasequentialorsimultaneousfashionwithaconsecutivecouplingstrategybetweendifferentdo-
mains. The main drawback of the state-exchange method is that the information propagation (e.g., flow development)
throughtheinterfacemaybeunphysical,resultinginaslowconvergencetosteady-stateforsteady-flowproblemsorinsig-
nificantrestrictionsofthemethod’sapplicabilitytofastunsteadyflows.Ingeneral,thestate-exchangemethodcanproperly
handleunsteadyflowswhichcanbeconsideredasquasi-steady.
Theflux-exchangemethodisbasedontheexchangeofrelevantfluxes(e.g.,massflux,momentumflux)betweendifferent
descriptions.Eventhoughtheflux-exchangemethodnaturallyfollowstheconservationlaws,itdoesnotguaranteethecon-
tinuityofstatevariables,whichmayrequireanadditionaltreatment[39].Moreover,theflux-exchangemethodappearsto
bemorerestrictivewithrespecttoanefficienttimedecoupling[36]thanthestate-exchangemethod,whichisabsolutely
necessaryforalargedisparityoftimescalesinatomisticandcontinuumrepresentations.Thechoiceofaparticularhybrid
methodoroftheircombinationisdictatedbytheflowproblem,i.e.,consideringthemainflowcharacteristics,thedisparity
oflengthandtimescales,andtheefficiencyofthechosenalgorithm.
Imposinginterfaceconditionsatartificialboundariesofatomisticandcontinuumdomainscreatesanadditionalchallenge
incoupledsimulations.Eventhoughboundaryconditions(BC)incontinuumapproachescanoftenbeproperlyimposed,the
extractionofrequiredstateinformationfromatomisticparticle-basedsolutionisadifficulttask.Theextractionofthemean
fieldproperties(ensembleaverages)typicallyrequiresspatio-temporalsamplingofcharacteristicscomputedbyatomistic
solver.Hadjiconstantinouet.al.[40]obtainedanaprioriestimateforthenumberofsamplesrequiredtomeasuretheaverage
offlowpropertiesinacellofaselectedvolumeforafixederror.Thenumberofsamplesdependsonthesignificanceofdeter-
ministicflowproperties withrespectto thethermalfluctuationsamplitudeand canbeverylarge ifthermalfluctuations
overwhelmdeterministicflowcharacteristics.Therearealsoapproacheswhichattempttominimizethedrawbackofthe
presenceofthermalfluctuations[39]ortofilterthemoutaswewillpresentinthispaper.
Non-periodicBCsinparticlemethodsareassociatedwithparticleinsertions/deletions[41]andtheapplicationofeffective
boundaryforces[33,35,42].TheexistingBCalgorithmsforNewtonianfluidsinparticlemethodsworkquitewell;however,
134 L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147
thecorrespondingalgorithmsforcomplexfluids(e.g.,polymericandbiologicalsuspensions)arestillunderactivedevelop-
ment.Insertionofbondedparticles(e.g.polymers,redbloodcells,etc.)demandsmoresophisticatedmethodsthaninsertion
ofnon-bondedparticles(e.g.,representingbloodplasma).Themainchallengehereistominimizenon-physicaldisturbances
inlocaldensityandvelocity.
Insummary,formultiscalemodelingofbloodflowweneedtoaddressalltheaforementionedchallenges.Thefirstchal-
lenge is the derivation of mathematical models which are able to accurately predict physiological processes. The second
challengeisthedesignofrobustinterfaceconditionsrequiredtointegratethecontinuumandatomisticsimulations.The
thirdchallengeistheprocessingofnon-stationarydatafromatomisticsimulations.Duetorelativelyfrequentdataexchange
betweenthestochasticanddeterministicsolvers,statisticalaveragingshouldbeperformedondatacalculatedwithinshort
timeintervals.Thefourthchallengeisacomputationalone.Inparticular,billionsofdegreesoffreedomarerequiredtoaccu-
ratelyresolvethemicro-scaleandalsotheinteractionbetweenthebloodcellsandendothelialcellsinavolumeassmallas
1mm3.Simulatingeffectsofsmallscalesonthelargerscalesrequiresincreasedlocalresolutioninthediscretizationofcom-
putationaldomainwithinthecontinuumdescription.Forexample,clotformationstartsfromaggregationofindividualblood
cellsandtheclotmustgrowtoacertainsizebeforeaffectingthelargescaleflowdynamics.However,toaccuratelycapture
theclotgeometryanditsshading-offsegments,veryhighlocalresolutionisrequired.Moreover,tofollowthedetachedclot
segmentsanadaptivelocalhp-typemeshrefinementmustbeused.
Inthispaperweaddresssomeoftheaforementioneddifficulties.InSection2wedescribethenumericalmethodsapplied
tosimulatebloodflowandclotformationprocessatvariousscaleswithsufficientlyhighresolution.Inthesamesectionwe
describeourapproachforcouplingtheatomisticandcontinuumsolvers.Specifically,wereviewtheinterfaceconditionsand
waystoincreasethecomputationalefficiency.InSection3wepresentresultsofcoupledsimulationsofabloodflowandclot
formationinasimplifieddomainofanarteryandalsoinapatient-specificbrainvasculaturewithananeurysm.InSection4
wesummarizeourstudy,discussthecurrentlimitationsofourmethods,andprovideanoutlookforfuturedevelopment.
2.Methods
Inthissectionwedescribethenumericalmethodsusedtosimulatebloodflowatbothcontinuumandatomisticlevels.
First,anumericalmethodforsolvingtheNavier–Stokesequationsusingacontinuumdescriptionisreviewed.Specifically,
wefocusonspatio-temporaldiscretizationandonthesmoothedprofilemethod(SPM)[43,44]implementedinourresearch
codeNejTar.Second,wedescribeanumericalapproachforcoarse-grainedatomisticsimulationusingthedissipativepar-
ticledynamics(DPD)method[45,46]andaprojectionmethodofatomisticdatatocontinuumfield.Finally,couplingofthe
continuum solver NejTar with the atomistic solver DPD-LAMMPS (based on the open source LAMMPS code [47]) is
described.
2.1.Continuum-basedmodeling:spectralelementmethod
Thelarge-scaleflowdynamicsismodeledbytheNavier–Stokesequations:
@u
þu(cid:2)ru¼(cid:3)rpþmr2uþf; r(cid:2)u¼0; ð1Þ
@t
whereuisthevelocityvector,pisthepressure,misthekinematicviscosity,tisthetime,andfisaforceterm,whichwewill
discussinmoredetailfurtherbelow.TheflowproblemisdefinedinarigiddomainX ofarterialnetworks.Tomodellarge-
C
scaleflowdynamics,weassumebloodtobeanincompressibleNewtonianfluidwithconstantdensityandviscosity.
Thethree-dimensional(3D)Navier–Stokesequationsaresolvedusingtheopen-sourceparallelcodeNejTardeveloped
at BrownUniversity. NejTar employsthe spectral/hp element spatial discretization (SEM/hp) [48], which provideshigh
spatialresolutionandissuitableforsolvingunsteadyflowproblemsingeometricallycomplexdomains.Thecomputational
domainisdecomposedintopolymorphicelements.WithineachelementthesolutionofEq.(1)isapproximatedbyhierar-
chical,mixedorder,semi-orthogonalJacobipolynomialexpansions[48].Fig.2showsanillustrationofthedomaindecom-
positionandthepolynomialbasisusedinNejTar.FortimeintegrationNejTaremploysahigh-ordersemi-implicittime-
steppingscheme[49].
Tosimulatemovingobjectsortime-evolvingstructureswithinafixedcomputationaldomainweusethesmoothedpro-
filemethod[43,44].SPMbelongstothefamilyofimmersedboundarymethods,andhenceithasnorequirementsonthe
meshtobeconformalwiththeboundariesofmovingobjectsorstructures.TheforcingtermfinEq.(1)iscomputedusing
SPMsuchthatno-slipboundaryconditions(BCs)areweaklyimposedatany(virtual)surfaceCvðt;xÞ2XC.
TheEq.(1)issolvedinseveralstepsusingthenumericalschemeoutlinedbelow.Inthefollowingequationstheindexnin
theschemecorrespondstothetimesteptn¼nDt,whereDtisthesizeofatimestep;c anda arecoefficientsfortheback-
0 k
warddifferentiation,b areextrapolationcoefficients,andJeisthetimediscretizationorder.
k
ThefirststepinsolvingEq.(1)includesexplicitcomputingoftheprovisionalvelocityfieldu(cid:4):
u(cid:4)(cid:3)PJke¼(cid:3)01ðakun(cid:3)kÞ¼(cid:3)XJe(cid:3)1ðb un(cid:3)k(cid:2)run(cid:3)kÞ: ð2Þ
Dt k
k¼0
L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147 135
Fig.2. AschematicoftheunstructuredsurfacegridandthepolynomialbasisemployedinNejTar.Thecomputationaldomainisdecomposedintoon-
overlappingelements.Withineachelementthesolutionisapproximatedbymixedorder,semi-orthogonalJacobipolynomialexpansions.Theshape
functionsassociatedwiththevertex,edgeandfacemodesforafourth-orderpolynomialexpansiondefinedontriangularandquadrilateralelementsare
shownincolor.(Forinterpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
InthesecondstepwesolveimplicitlytheHelmholtzequationforeachvelocityvectorcomponent:
c u~(cid:3)u(cid:4)
0 ¼(cid:3)rpþmr2u~; ð3Þ
Dt
wherethepressurefieldpiscomputedfromthePoissonequation:
r2p¼ðDtÞ(cid:3)1r(cid:2)u(cid:4): ð4Þ
ThePoissonequationforpisderivedbyapplyingthedivergenceoperatortoEq.(3)andusingtheconditionr(cid:2)u~(cid:5)0.Inthe
caseofzeroforceterm,thesecondstepcompletesthenumericalintegrationofEq.(1)andthesolutionvectorattimestep
tþDtisdefinedtobeunþ1¼u~.Withanon-zeroforceterm,anadditionalstepisrequiredtoobtainsolutionattimestep
tþDt.
Atthethirdstep,thevelocityfieldu~ obtainedinthesecondstepiscorrectedinordertoimposeno-slipBCsatmoving
boundariesCvðt;xÞ:
unþ1(cid:3)u~ u (cid:3)u~
c ¼/nþ1 p (cid:3)rp~; ð5Þ
0 Dt Dt
wherethepressurefieldp~iscalculatedfrom:
r2p~¼ðDtÞ(cid:3)1r(cid:2)(cid:2)/nþ1ðu (cid:3)u~Þþu(cid:4)(cid:3): ð6Þ
p
ThePoissonequationforp~isderivedbyapplyingthedivergenceoperatortoEq.(5)andassumingthatr(cid:2)unþ1(cid:5)0.
Thescalarfield/ðt;xÞinEq.(6)istheso-calledindicatorfunction.Thisfunctioniszerooutsidetheregionsboundedby
Cvðt;xÞ,whileinsidetheseregions/ðt;xÞ¼1.Notethat/ðt;xÞisnotastepfunctionsuchthatitsmoothlyvariesfromzero
toonearoundCvðt;xÞwith/ðt;CvÞ¼0:5.Thevectorupðt;xÞspecifiesrequiredflowvelocitywithitsvaluesbeingimportant
onlyfor/ðt;xÞ>0.Thesmoothnessoffields/andu isveryimportant,sincehigh-orderaccuracycanonlybeachievedfor
p
smoothfields.Insection2.3wewillfurtherextendthediscussionofCv;/,andup,andshowtherelationofthesefieldsto
dataprovidedbytheatomisticsolver.
The patient-specific arterial networkconsidered in our study isvery large. Moreover,to accuratelyrepresent the wall
shearstresses(animportantcharacteristicinbiologicalflows)highspatialresolutionisrequired.Thisleadstoanextremely
largesizeofcomputationalproblem(intermsofdegreesoffreedom(DOF))withOð109ÞDOF.Toefficientlysolvesuchlarge
problemweemployedamulti-patchdomaindecompositionmethod[50].Thismethoddecomposesthefulltightly-coupled
problem into a number of small tightly-coupled problems defined in subdomains (patches), where the global continuity
(coupling)isenforcedbyprovidingproperinterfaceconditions.Thus,themulti-patchmethodsignificantlyimprovescom-
putationalefficiencyandessentiallyremovesthelimitsimposedduetolargeproblemsize.Themulti-patchapproachiswell
suitedwiththefunctionaldecompositionstrategyweemployinmultiscalesimulationssuchthatineachpatchadifferent
numericalmodelfortheflowproblemcanbeapplied.Forexample,thecomputationaldomainofbrainarteriesshownin
Fig.1isdecomposedintofourpatches,butonlyonepatchhasaninterfacewiththeatomisticsolver,wheretheSPMmethod
isemployed.ThesolutionofEq.(1)withanon-zeroforcetermimpliesthatonemorePoissonequationandthreeadditional
projectionproblemsmustbesolved,whichpracticallydoublesthecomputationaleffort.UsingtheflowsolverwiththeSPM
inasinglepatchleadstoabettermanagementofcomputationalresources.
2.2.Atomisticbasedmodeling:dissipativeparticledynamicsmethod
Tomodelbloodflowdynamicsatatomistic/mesoscopicscaleweemployacoarse-grainedstochasticmoleculardynamics
approach[51]usingthedissipativeparticledynamics(DPD)method[45,46].DPDisamesoscopicparticlemethodwitheach
136 L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147
particlerepresentingamolecularclusterratherthananindividualmolecule.TheDPDcanbeseamlesslyappliedtosimulate
bondedstructures(e.g.,polymers,bloodcells)andnon-bondedparticles(e.g.,bloodplasma).TheDPDsystemconsistsofN
pointparticlesinteractingthroughconservative,dissipative,andrandomforcesgivenrespectivelyby
(cid:4) r (cid:5) n
FCij ¼aij 1(cid:3)ricj ^rij; FDij ¼(cid:3)cxDðrijÞðvij(cid:2)^rijÞ^rij; FRij¼rxRðrijÞpffiDiffijffiffitffiffi^rij; ð7Þ
where^r ¼r =r ; v ¼v (cid:3)v; Dtisthetimestep,andr isthecutoffradiusbeyondwhichallforcesvanish.Thecoefficients
ij ij ij ij i j c
a ; c,andrdefinethestrengthofconservative,dissipative,andrandomforces,respectively.xDandxRðr Þ¼(cid:7)1(cid:3)rij(cid:8)kwith
ij ij rc
the exponent k are weight functions, and n ¼n is a normally distributed random variable. The random and dissipative
ij ji
forces form a thermostat and must satisfy the fluctuation–dissipation theorem [46] leading to the two conditions:
xDðr Þ¼(cid:2)xRðr Þ(cid:3)2 andr2¼2ck T withTbeingtheequilibriumtemperature.ThemotionofDPDparticlesisgovernedby
ij ij B
Newton’slaw:
dri¼v m dvi¼f f ¼X(cid:7)FCþFDþFR(cid:8): ð8Þ
dt i i dt i i ij ij ij
j;i–j
Theatomisticproblemisdefinedinafixednon-periodicdomainX ,wherethefluidisrepresentedbyacollectionofsolvent
A
particleswithanumberofsuspendedplatelets.ThedomainboundariesC arediscretizedintotriangularelementsT,where
A
localBCvelocitiesareimposed.Ingeneral,weimposeeffectiveboundaryforcesF ontheparticlesnearboundariesthat
eff
representsolidwallsandinflow/outflowBCs.Suchforcesimposeno-slipBCsatsolidwallsandcontroltheflowvelocities
atinflow/outflow[42].Inaddition,atinflow/outflowweinsert/deleteparticlesaccordingtolocalparticleflux[42].Fig.3
presentsatypicaldiscretizationofC .Incoupledatomistic-continuumsimulationstheBCsprescribedatC arebasedon
A A
datareceivedfromthecontinuumsolver,whichwillbedescribedinthenextsection.
CouplingofatomisticandcontinuumsolversrequiresthecalculationofaveragedpropertiesinX suchasvelocityand
A
density.Accuratecomputingoftheaveragedfieldsthroughprocessingofatomisticdatapresentsseveralchallengesinclud-
inggeometricalcomplexityoftheatomisticdomains,thermalfluctuations,andflowunsteadiness.Tocomputeaverageprop-
ertiesweadoptaspatio-temporalbinningstrategy.ThecomputationaldomainX isdiscretizedintotetrahedralelements
A
(bins) X . The time–space averaged velocity and density are computed by sampling data over predefined time intervals
b
Dt withineachbin:
s
1 nX1þNtsX
u(cid:2) ¼ un; ð9Þ
b N i
p n¼n1 i
where u(cid:2) denotes the average velocity in bin b;N is the number of particles passing through X over the time interval
b p b
Dt ¼N Dt,andun isthevelocityofaparticleiattimestepn.
s ts i
Theaveragingofatompropertiesoverbin’svolumeleadstocellcentereddata,whichisdiscontinuousatthebininter-
faces.Aswementionedabove,thespectralconvergenceinthespectralelementmethodcanbeachievedonlyforsmooth
fields.TomapthecellcentereddatatoavertexcenteredC0 continuousfieldweapplyaprojectionoperatorPcc2vc:
XNb
u(cid:2)v ¼Pcc2vcðu(cid:2)bÞ¼ ðxvbu(cid:2)bÞ; ð10Þ
b¼1
wherevrepresentsavertexofXb; Nb isthenumberofbinssharingthevertexv,andxvb arethenormalizedintegration
weights:8vPbxvb¼1.Theweightsareinverselyproportionaltothedistancebetweenthevertexvandthecenterofbin
b.Analternativemethodtocomputevertexcenteredflowpropertiesistosampletheatomisticdatawithinasphericalvol-
umecenteredattheverticesofelementsX .Bothmethodsproducesimilarresults,whilethemaindifferenceisincompu-
b
tationalcomplexity.Ingeneral,itiscomputationallymoreefficienttodetermineifaparticleislocatedinsideasphericalbin
Fig.3. AtomisticdomainXA withatriangulationofthedomainboundariesCA.Colorsrepresentthey-componentofthevelocityimposedatCA.(For
interpretationofthereferencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147 137
incomparisontothatinsideapolymorphicelement.Moreover,thenumberofverticesistypicallyafractionofthenumberof
elements.However,possibleoverlapbetweenthevolumescenteredateachvertexrequiresthatthelocationofeachparticle
istestedinanumberofsphericalbins,whichsignificantlyincreasesthecomputationalburden.Incaseofhexahedralbins,
wherebinfacesareproperlyalignedwiththecoordinatesystem,calculationofcellcentereddataandtheirprojectiononthe
vertex centered configuration is the most efficient. However, such hexahedral bins cannot be used in arbitrary complex
geometries.
TheaccuracyincomputingofaveragepropertiesiscontrolledbythesizeofabinandthetimeintervalDt.Insteadyflows
s
Dt canbechosenverylargetominimizetheeffectofthermalfluctuationsandtoaccuratelycapturethedeterministiccom-
s
ponent.Innon-stationaryflows,forinstancebloodflowconsideredinourstudy,theaveragedpropertiesaretimedependent,
anditisnotalwaysclearhowtodefinethetimeintervalDt.Thus,dataaveragedoverashortDt maybeverynoisydueto
s s
thermalfluctuations,whileaveragingoveralongDt willinadequatelyapproximatethetimevariationofthedeterministic
s
component within the atomistic fields. In order to overcome this difficulty we employ the window proper orthogonal
decomposition (WPOD) [52,53]. WPOD is a spectral analysis tool used to transform a field into M orthogonal temporal
and spatial modes. The low-order POD modes contain most of the energy and correspond to the ensemble average field,
whilethehigh-orderPODmodesreflectthermalfluctuationcomponents.ToconstructthePODmodesweapplythemethod
of snapshots [54], where each snapshot contains a spatio-temporal average of atomistic data computed over a relatively
short Dt. Fig. 4 illustrates the applicationof WPODmethod for processingof atomisticdata.The plots clearly showthat
s
WPODsubstantiallyimprovestheaccuracyinthereconstructionofthetimedependentaverageflowfieldandalsoitsderiv-
atives.Forthisproblem,theaveragedsolutionobtainedwiththeWPODtechniquewasapproximatelyoneorderofmagni-
tudemoreaccurate(inL -norm)thanthesolutionobtainedusingstandardaveraging.Oursimulationtestsindicatethatin
2
ordertoachieveacertainaccuracy,theWPODapproachallowsustobalancethenumberofPODmodesandthelengthof
timeintervalDt foraveragingofasinglesnapshot.Fig.5showsthatforMDt ¼constidenticalaccuracyinreconstructionof
s s
thedeterministiccomponentofanunsteadyflowcanbeachievedinpractice.
2.3.Coupledatomistic-continuumsolver
ToconstructasolverformultiscalesimulationswecouplethecontinuumsolverNejTarwiththeatomisticsolverDPD-
LAMMPS.AschematicrepresentationoffunctionaldecompositionimplementedinourcoupledsolverispresentedinFig.6
andanexampleofacoupledproblemisprovidedinFig.1.Eachsolverfullypreservesitsoriginalcapabilities,inparticular
themulti-patchdomaindecompositioninNejTarandmultiplereplicastocomputethestatisticalpropertiesmoreaccu-
ratelyinDPD-LAMMPS.Independenttightly-coupledproblemsaresolvedwithineachcontinuumoratomisticdomainusing
non-overlappinggroupsofprocessors.Additionalsub-groupsofprocessorsarederivedtohandletasksrequiredtoimpose
BCsatinlets/outletsandalsoatinter-domaininterfaces.
Couplingoftheatomisticandcontinuumsolversfollowstheframeworkdescribedin[35],wherethecontinuumsolverfor
Navier–StokesequationsiscoupledtoDPDandalsotomoleculardynamics(MD).Theflowdomainisdecomposedintoa
numberofoverlappingregions,whichmayemploydifferentdescriptionssuchasMD,DPD,orcontinuum.Eachsub-domain
is integrated independently. The coupling between overlapping sub-domains is performed by exchanging state variables
Fig.4. ProcessingofdatafromtheDPDsimulationofunsteadyflowinapipe:(a)–averagedsolution(streamwisevelocitycomponent,alignedwiththex
axis)computedwithastandardaveragingwithNts¼50;(b)–averagedsolutionprocessedwiththeWPOD;herethedataisreconstructedusingthefirst
twoPODmodes.(c)–exactsolution.(d–f)–gradientsofthestreamwisecomponentofthevelocityfiled.(Forinterpretationofthereferencestocolourin
thisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
138 L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147
100
>)
u
<
or (
err 10−1
−
2
L
Nts100_m40
10−2 Nts200_m20
0 50 100 150 200 250 300
t
Fig.5. ProcessingofdatafromtheDPDsimulationofunsteadyflowinapipe:accuracy.Thestreamwisevelocityiscomputedfromatomisticdatausingthe
WPODmethod.Thetimewindowoverwhichthedataisprocessediskeptconstant.Blacksolidcurve–eachsnapshotiscomputedover100timesteps,
numberofPODmodesis40.Reddot-dashcurve–eachsnapshotiscomputedover200timesteps,numberofPODmodesis20.(Forinterpretationofthe
referencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
Fig.6. Aschematicrepresentationoftaskparallelisminourcoupledsolver.ThecontinuumdomainispartitionedintoNoverlappingpatchesXCi,andthe
patchesexchangedataviainter-patchcommunicators.EachXCi hasmultipleinletsandoutlets,andeachpatchmaybeconnectedto0Dor1Darterial
networkmodelinordertomodelBCsforinlets/outlets.TwoatomisticdomainsXAiareplacedinsideXC1.EachXAicanbereplicatedseveraltimestoreduce
statisticalerror.EachXAiislinkedtoaXC1viatwointerfaces:A2DinterfacewhichconsiderstheboundariesofXAiandusesdatacomputedinXC1asBCs.A
3Dinterface,whichistailoredtotheimmersedboundarymethod;XC1 usesanexternalforcefieldcomputedinXAi.
everyDsintimeprogressionasshowninFig.7.ThetimeDsmaycorrespondtoadifferentnumberoftimestepsfordistinct
descriptions.Inthispaperwealsoextendthemethodof[35]andapplyittounsteadyflowinacomplexgeometry.
Tosetupamultiscaleproblemwithheterogeneousdescriptionswehavetodefinelengthandtimescales.Inprinciple,the
choiceofspatio-temporalscalesmaybeflexible,butitislimitedbyvariousfactorssuchasmethodapplicability(e.g.,sta-
bility,flowregime)andproblemconstraints(e.g.,temporalresolution,microscalephenomena).Forexample,aunitoflength
(L )inX correspondsto1mm,whileaunitoflength(L )inX isabout1lm.Inaddition,fluidproperties(e.g.,viscosity)in
C C A A
differentdescriptionsmaynotnecessarilybeinthesameunitsinvariousmethods.Toproperlycoupledifferentdomainswe
arerequiredtoconsistentlynon-dimensionalizethetimeandlengthscalesandtomatchnon-dimensionalnumberswhich
characterizetheflow.Forexample,tomatchtheReynoldsnumber(Re)indifferentdomainsthevelocityinX isscaledas:
A
L m
u ¼u C A; ð11Þ
A CL m
A C
L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147 139
t
communication
step K+1
Δτ
communication
step K
Δ t Δ t
NS DPD
Fig.7. Aschematicrepresentationofthetimeprogressionindifferentsub-domains.ThetimestepratioemployedinourstudyisDtNS=DtDPD¼20,and
Ds¼10DtNS¼0:0344s.
wherem andm arethekinematicfluidviscositiesintheX andX .Thetimescaleineachsub-domainisdefinedast(cid:6)L2=m
C A C A
andisgovernedbythechoiceoffluidviscosity.Inoursimulationsasingletimestepinthecontinuumsolver(Dt )corre-
NS
sponds to 20 time steps in the atomistic solver (Dt ). The data exchange between the two solvers occurs every
DPD
Ds¼10Dt ¼200Dt (cid:6)0:0344s.
NS DPD
Fig.8illustratesthesetupforcoupledatomistic-continuumsimulationofplateletaggregation.First,thecontinuumdo-
mainX iscreated.Second,theatomisticdomainX isplacedintheareaofinterest,suchthatX 2X .Third,anadditional
C A A C
sub-domainX forsamplingofatomisticdataisinsertedintoX suchthatX 2X .TheboundariesofX arediscretized
S A S A A
using2Delements(triangles),whilethevolumeofX isdiscretizedusing3Delements(bins).
S
Theinterfaceconditionsusedforatomistic-continuumcouplingarebasedontherequirementsforeachsolver.Theatom-
isticsolverrequiresalocalvelocityfluxtobeimposedateachelementT ofC .ThecoordinatesofT centersareprovidedby
A
theatomisticsolvertothecontinuumsolveratthepreprocessingstage,andinterpolationoperatorsprojectingthevelocity
fromthecontinuumfiledontothecentersofT areconstructed.Duringtimeintegrationdatacomputedbythecontinuum
solverisinterpolatedandtransferredtotheatomisticsolver.Thecontinuumsolveralsorequirestheindicatorfunction/and
thevelocityu .ThemainobjectiveofthesimulationshowninFig.8istostudytheformationofbloodclot.Theclotisformed
p
byaggregationoftheDPDparticles(platelets)onthearterialwall,i.e.,thebloodclotissimulatedbytheatomisticsolver.The
interactionbetweenthemicroandmacroscalesstronglydependsontheeffectsoftheclotgrowthontheflow.Specifically,
theclotformsanobstacleattachedtothewall,andalsosomepartsoftheclotmaybreakupandadvectdownstream.Since
theclotisformedbyanevolvingclusterofDPDparticleswhichrepresentplatelets,thelocaldensityoftheseparticlesin-
creasesincomparisontothatinbulkbloodflow.Tomimicthisobstacle(aggregate)weconstructtheindicatorfunction/in
suchawaythatithasavalueofoneintheregionwithhigh-densityofactive(sticky)plateletsandzerointheregionwith
normalplateletsdensityasinbulkflow.Tothisend,theindicatorfunctionisdefinedas:
(cid:9) (cid:4)q ðt;xÞ(cid:3)q (cid:5) (cid:10)
/ðt;xÞ¼0:5 tanh p 0 þ1 ; ð12Þ
Łq
Fig.8. Acoupledatomistic-continuumsimulation:illustrationofcomputationaldomainsandinter-domaindataexchange.AnatomisticdomainXAisfully
embeddedintothecontinuumdomainXCintheregionwhereplatelet(yellowdots)depositionissimulated.Interfacevelocityconditionsareimposedat
theboundariesCAofXA.VirtualboundaryconditionsforvelocityareimposedinsideXCusingthesmoothedprofilemethod,whichrequirestheindicator
function Uðt;xÞ and the particle velocity upðt;xÞ. The fields Uðt;xÞ and upðt;xÞ are sampled within the sub-domain XS placed inside XA. Right –
computationaldomainofananeurysmXA2XC isplacedatthewallofXC,wherethecontactsurfaceiscoloredinpurple.(Forinterpretationofthe
referencestocolourinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)
140 L.Grinbergetal./JournalofComputationalPhysics244(2013)131–147
whereq ðt;xÞistheconcentrationofactiveplatelets,q isthepredefinedthreshold(inoursimulationsq istwotothree
p 0 0
timesgreaterthanthebulkplateletdensity),andŁqisthescalingfactorwhichcontrolsthesteepnessoftheindicatorfunc-
tion.ThevaluesoftheindicatorfunctionarecomputedbytheatomisticsolverwithinX atthecoordinateswhichcorre-
S
spondtothelocationsofgridpointsinX .Duetoahigh-spatialresolutioninthecontinuumdomainthenumberofgrid
C
pointscanbeverylarge.Tocomputetheindicatorfunctionefficientlythedensityq ðt;xÞisfirstcomputedatthebincenters
p
byspatio-temporalaveragingovertimeintervalDs.Then,theprojectionoperatorPcc2vc isappliedtocomputethevertex-
centereddata.Theprojectionfromcell-centeredtovertex-centereddataisfollowedbytheinterpolationofq ðt;xÞatbin
p
verticestothegridpointsofX usinglocallinearinterpolationwithineachbin.Theindicatorfunctionisthencalculated
C
fromtheinterpolateddensityvalues.
Monitoringthesizeofthebloodclotovertimemaybeaverydifficulttask.Themajorchallengeistodefinetheclotandits
boundariesintermsoftheactiveplateletsandtheirvelocity.Thedescriptionofthebloodvolumeoccupiedbyactiveplate-
letsbyanindicatorfunctionprovidesaconvenientmethodtomeasurethevolumeoccupiedbytheclot.Forexample,we
maysimplyintegratetheindicatorfunction/ðt;xÞoverX orX insteadofcountingindividualactiveplatelets.However,
A C
suchmethodwillalsoincludetheregionsoccupiedbymovingobjectscomposedofaggregatedplatelets,whicharenota
partofthestaticbloodclotnearthewall.Thus,tomaketheestimationoftheclotgrowthmorereliable,itmightbenecessary
tointroduceweightsfortheindicatorfunctionwhichassumeavalueofoneintheregionofveryslowflowanddecayexpo-
nentiallyfastastheflowvelocityincreases.
2.4.Plateletaggregationmodel
PlateletsaremodeledbysingleDPDparticleswithalargereffectiveradius[55]thanthatofplasmaparticlesandarecou-
pledtotheplasmathroughtheDPDdissipativeinteractions.Themodelofplateletaggregationisadoptedfrom[56],where
plateletscanbeinthreedifferentstates:passive,triggered,andactivated.Inthepassivestate,plateletsarenon-adhesiveand
interactwitheachotherthroughtherepulsiveDPDforces,whichprovidetheirexcludedvolumeinteractions.Passiveplate-
letsmaybetriggerediftheyareinclosevicinityofanactivatedplateletorinjuredwall.Whenaplateletistriggered,itstill
remainsnon-adhesiveduringtheso-calledactivationdelaytime,whichischosenrandomlyfromaspecifiedtimerange.After
the selected activation delay time, a triggered platelet becomes activated and adhesive. Activated platelets interact with
other activated particles and adhesive sites, which are placed at the wall representing an injured wall section, through
theMorsepotentialasfollows
UMðrÞ¼De(cid:2)e2bðr0(cid:3)rÞ(cid:3)2ebðr0(cid:3)rÞ(cid:3); ð13Þ
whereristheseparationdistance,r isthezeroforcedistance,D isthewelldepthofthepotential,andbcharacterizesthe
0 e
interactionrange.TheMorsepotentialinteractionsareimplementedbetweeneveryactivatedplateletoradhesivesiteifthey
arewithinadefinedpotentialcutoffradiusr .TheMorseinteractionsconsistofashort-rangerepulsiveforcewhenr<r
d 0
andofalong-rangeattractiveforceforr>r suchthatr correspondsapproximatelytotheeffectiveradiusofaplatelet.
0 0
Finally,activatedplateletsmaybecomepassiveagain,iftheydidnotinteractwithanyactivatedplateletoradhesivesite
duringafiniterecoverytime.
2.5.Accuracyconsiderations
Validationofthemathematicalmodelsusedforcomplexbiologicalsystemsisextremelydifficult.Ultimately,itwould
requirethepossibilitytoperformexperimentsontheonsetandsubsequentgrowthoftheclotformationinthebrain,which
iscurrentlynotfeasible.Theaccuracyofthecoupledatomistic-continuumsolverdependsonmodelinganddiscretization
errors.Specifically,modelingerrorsinthecontinuumapproachcanberelatedtotheassumptionofrigidvesselwalls,while
realvesselwallsaredeformable,andtouncertaintiesintheinflow/outflowboundaryconditionsaswellasinmaterialprop-
erties.Attheatomisticlevel,theplateletaggregationmodelmayalsoleadtomodelingerrors.Weusedtheplateletaggre-
gationmethodpresentedin[56],whichproducedplateletaggregationratescomparablewithexperimentallyobservedones
ininvitroexperiments[57].However,theclotformationprocessmaybeverydifferentinthebrain,anditsbiomechanical
detailshavenotbeenquantifiedyet.Inaddition,interactionsbetweenredbloodcells,whitebloodcells,andfibrinogenhave
tobeincluded,andsuchinteractionsareparticularlyimportantforsimulatingadvancedstagesofclotformation.Moreover,
modeling clot formation at the initial and intermediate stages may require the use of time-dependent models such that
aggregationinteractionsaretimedependent,forexample,totreatthedepthoftheMorsepotential(D )asatime-varying
e
parameter.
Thediscretizationerrorsareattributedtotheusualtime–spacediscretizationparameters,i.e.,Dt;handpmeshrefine-
ment, errors in the DPD solver due to thermal noise, to coupling errors between continuum and atomistic descriptions,
butalsotothestatisticalerrorsincomputingensembleaverages.Thenumericalaccuracyofdifferentcomponentsinour
coupledsolverhasbeenverifiedbyperformingtimeandspace-discretizationrefinementinthecontinuumsolver.Wehave
also performed atomistic simulations using three different resolutions, with L =L ratio equal to 100, 200 and 400. Com-
A C
parableresultshavebeenobtainedwithL =L ¼200andL =L ¼400,whilethefinestresolutionhasbeenusedinthesim-
A C A C
ulationspresentedbelow.Errorsincouplingthecontinuumandatomisticdescriptionsforsimulatingstationaryflowshave