Transcript Simulering och visualisering av en Eulerisk fluid p˚a grafikkortet

ITN, Norrk¨oping
20 mars 2011
Simulering och visualisering av en
Eulerisk fluid p˚a grafikkortet
M ODELLERINGS - OCH S IMULERINGSPROJEKT
TNM085
Medlemmar:
Christopher Birger
Henrik B¨acklund
Erik Englesson
Anders Hedblom
Per Karlsson
Kontakt:
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
Sammanfattning
Fluidsimuleringar blir allt vanligare n¨ar avancerad mjukvara driver h˚ardvaruutvecklingen fram˚at.
Att kunna modellera och simulera fluider som r¨ok, vatten och eld har en stor betydelse inom
m˚anga omr˚aden, s˚a som film, spel och forskning, n¨ar det blir sv˚art eller om¨ojligt att anv¨anda
en riktig fluid.
Det finns givetvis olika s¨att att genomf¨ora detta p˚a, beroende p˚a vilket resultat som ska
uppn˚as – en grafiskt snygg men snabb simulering eller en mer realistisk s˚adan. Alla fluidsimuleringar bygger p˚a Navier-Stokes ekvation och det finns i princip tv˚a syns¨att f¨or att l¨osa
dessa. Antingen anv¨ands partiklar (Langrageisk) eller s˚a diskretiseras fluiden p˚a ett rutn¨at
(Eulerisk). P˚a senare tid har man a¨ ven anv¨ant en kombination av dessa. I den Euleriska metoden
tittar man p˚a fixa punkter i rymden och hur fluidkvantiteter som hastighet och tryck a¨ ndrar sig
i tiden, ist¨allet f¨or att f¨olja partiklar. I den h¨ar rapporten implemeteras en Eulerisk l¨osning.
Eftersom vi har valt att simulera vatten beh¨ovs n˚agot s¨att f¨or att representera dess yta. Till detta
anv¨ands level set-metoden.
N¨amnda metoder a¨ r mycket ber¨akningstunga och har varit n¨ast intill om¨ojlig att utf¨ora i
realtid. Tack vare dagens kraftfulla grafikkort har det dock b¨orjat bli fullt m¨ojligt att simulera
och rendera en 3D-fluid i realtid.
Den h¨ar rapporten tar upp teori och tillv¨agag˚angss¨att som anv¨ants f¨or att skapa och visualisera
en vattenliknande fluid med den Eulariaska metoden p˚a datorns GPU. All kod a¨ r skriven av oss
sj¨alva och inga externa bibliotek har anv¨ants. Projektets grundutformning var att, i grupp, skapa
en modell av ett fysikaliskt fenomen och visualisera en simulering av denna. Arbetet h¨olls i
kursen TNM085 - Modellbygge och Simuleringsprojekt p˚a Link¨opings universitet, v˚arterminen
2011.
F¨orord
Att fullf¨olja detta projekt har varit v˚ar st¨orsta utmaning hittills, och motg˚angarna har varit
m˚anga. Eftersom id´en med att simulera en Eulerisk fluid p˚a grafikkortet i realtid har blivit
bepr¨ovad av f˚a, har det varit sv˚art att f˚a tag p˚a svar till v˚ara fr˚agor.
Vi vill ge ett stort tack till PhD Andreas S¨oderstr¨om, som tog sig tid att ses och ge f¨orklaringar
p˚a v˚ara fr˚agest¨allningar och g¨ora projektet genomf¨orbart p˚a s˚a kort tid.
Inneh˚allsf¨orteckning
1
Inledning
1.1 Syfte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1
2
Metod
2.1 Beteckningar . . . . . . . .
2.2 Begr¨ansad rymd . . . . . . .
2.3 Interpolation . . . . . . . . .
2.4 Tidssteg . . . . . . . . . . .
2.5 Implementering . . . . . . .
2.5.1 GPU-programmering
2.5.2 Offline-rendering . .
2.5.3 Realtidssimulering .
2.6 Mjukvara . . . . . . . . . .
.
.
.
.
.
.
.
.
.
2
2
2
3
4
4
5
5
6
6
.
.
.
.
.
.
.
7
7
7
8
9
9
12
12
.
.
.
.
13
13
13
14
15
3
4
.
.
.
.
.
.
.
.
.
Grunderna hos en fluid
3.1 Navier-Stokes ekvation . . . .
3.2 Yttre krafter . . . . . . . . . .
3.3 Advektion . . . . . . . . . . .
3.4 Extrapolering av hastigheterna
3.5 Projektion . . . . . . . . . . .
3.6 Divergensfrihet . . . . . . . .
3.6.1 Dirichlets randvillkor .
Level set-metoden
4.1 Implicita ytor
4.2 β- och γ-band
4.3 Advektion . .
4.4 Reinitiering .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
Resultat
16
5.1 Realtid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
5.2 Utrenderat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
6
Avslutning
18
6.1 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
6.2 Vidare f¨orb¨attringar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Litteraturf¨orteckning
21
Figurer
2.1
2.2
2.3
2.4
Rymdens uppdelning i mindre voxlar. . . . . . . . . .
Definitionen av en MAC-ruta (Marker and Cell) i 2D.
Interpolation i 2D. . . . . . . . . . . . . . . . . . .
Transformation fr˚an 3D-textur till 2D. . . . . . . . .
3.1
En punkt sp˚aras ett tidsteg bak˚at med den bl˚a hastighetsvektorn. I den punkten
interpoleras en hastighet fram med hj¨alp av omkringliggande hastigheter (r¨oda
vektorer). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Bild till exempel (3.17). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Bild p˚a innan och efter Dirichlets randvillkor uppfyllts . . . . . . . . . . . . . 12
3.2
3.3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
3
4
5
4.1
Allt inom den svarta gr¨ansen a¨ r fluid. Det gr¨ona f¨altet a¨ r β-band och det
streckade a¨ r γ-band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
5.1
5.2
Resultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Resultat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Kapitel 1
Inledning
Det finns m˚anga anledningar till varf¨or man skulle vilja visualisera simuleringsbara modeller av
fluider. Specialeffekter till film, milj¨oer i spel eller avancerade tester inom andra industrier, bara
f¨or att n¨amna n˚agra. Exempel p˚a effekter som g˚ar att a˚ terskapa v¨aldigt bra med fluidsimuleringar
a¨ r r¨ok och eld, men i det h¨ar projektet handlar det om att skapa vatten, med ett realistisk
utseende och beteende.
Det finns olika s¨att p˚a hur man kan skapa en datorgenererad fluid. En metod som har blivit
omtyckt inom branschen f¨or specialeffekter, a¨ r den partikelbaserade SPH-metoden (Smoothed
Particle Hydrodynamics). Den bygger p˚a att man har en stor m¨angd partiklar, var och en med
fysikaliska egenskaper som massa och kollisionshantering. En f¨ordel med SPH a¨ r att metoden
ger ett visuellt tilltalande resultat i scener d¨ar det a¨ r mycket st¨ank. Ett problem a¨ r dock att en
SPH-fluid aldrig kan ligga riktigt still. Dess yta kommer att ha en konstant vibration. D¨arf¨or
l¨ampar sig exemplevis inte SPH i en scen d¨ar vattnet ska ligga still.
Ett annat tillv¨agag˚angss¨att a¨ r med den s˚a kallade Euleriska metoden. Den baserar sig p˚a
att, ist¨allet f¨or att anv¨anda partiklar, beskriva fluiden med hj¨alp av en fix yttre dom¨an uppdelad
i celler/voxlar. F¨or att beskriva ytan anv¨ands level set-metoden [1]. Det a¨ r den h¨ar metoden
som har anv¨ants i detta projekt, och dess olika steg och problem, som kan uppst˚a, beskrivs i
rapporten.
1.1
Syfte
Syftet med denna rapport a¨ r bland annat att:
¨ f¨orst˚aelsen f¨or hur en fluid kan vara uppbyggd
• Oka
• F¨orklara den Euleriska metoden och v¨ardefulla anv¨andningsomr˚aden f¨or den
• Dokumentera f¨or- och nackdelar med metoden
• Ge insikt i vilka problem som kan uppst˚a
1
Kapitel 2
Metod
2.1
Beteckningar
Beteckningar som anv¨ands genom rapporten a¨ r:
Skal¨arer – gemener (t), Matriser – versaler (A), Vektorer/vektorf¨alt – gemener med str¨ack (~u)
2.2
Begr¨ansad rymd
F¨or att kunna utf¨ora alla numeriska ber¨akningar p˚a datorn m˚aste en begr¨ansad rymd definieras
samt diskretiseras. Diskretisering utf¨ors genom att volymen delas in i s˚a kallade voxlar (figur
2.1). Om en volym med dimensionerna 1x1x1m skall delas in i 64x64x64 voxlar f˚as en voxelstorlek,
∆h, p˚a 1/64 m. Finare diskretisering ger h¨ogre uppl¨osning men flera ber¨akningar.
Varje voxel har en mittpunkt, sex sidor och 12 kanter med l¨angden ∆h. I denna implementation
anv¨ands ett s˚a kallat MAC-grid (Marker and Cell). Detta inneb¨ar att hastigheter sparas i mitten
av varje sida medan tryck sparas i mitten av varje voxel. Denna metod valdes d˚a precisionen i
utr¨akning av tryckgradienter samt divergens bli mer exakta [2].
2
Figur 2.1: Rymdens uppdelning i mindre voxlar.
Figur 2.2: Definitionen av en MAC-ruta (Marker and Cell) i 2D.
2.3
Interpolation
Ett problem som uppst˚ar vid en diskretisering som denna a¨ r att det enbart finns diskreta v¨arden
med ett visst mellanrum att jobba med. Oftast kommer man bli tvungen att h¨amta v¨arden som
inte finns i positionerna mellan de diskreta gr¨anserna. En interpolation blir tvungen att utf¨oras.
Det vill s¨aga, v¨arden runt omkring positionen i fr˚aga viktas och sl˚as ihop. Ett exempel p˚a enkel
linj¨ar interpolation i tv˚a dimensioner kan ses i figur 2.3.
I tre dimensioner kan enkel trilinj¨ar interpolation anv¨andas, men f¨or ett mer korrekt v¨arde
b¨or interpolation av h¨ogre ordning anv¨andas.
3
Figur 2.3: Interpolation i 2D.
I projektet anv¨ands en Monoton Catmull-Rom-interpolation [3] som a¨ r av ordning tv˚a.
Detta a¨ r en kubisk interpolation som utnyttjar derivator vid viktningen. Dessa derivator tas
fram genom central differensapproximation.
2.4
Tidssteg
Tidssteget ∆t v¨aljs s˚a att det garanterar ett CFL-tillst˚and (Courant-Friedrichs-Levy’s tillst˚and)
[1]:
∆h
(2.1)
∆t = k ·
max(ux + uy + uz )
Detta medf¨or att information maximalt kan transporteras en voxel per tidssteg. Konstanten k
anv¨ands f¨or att skala tidsteget i (2.1). F¨or att vara s¨aker p˚a att CFL-tillst˚andet uppfylls b¨or k
ligga mellan noll och ett.
Ett problem som kan uppst˚a med detta stabilitetsvillkor a¨ r att n¨ar hastigheten hos fluiden
o¨ kar, sjunker tidssteget proportionellt mot denna. En mycket h¨og hastighet ger ett mycket litet
tidssteg. Det tar dock fortfarande lika l˚ang tid f¨or datorn att utf¨ora ber¨akningarna under ett
tidssteg.
2.5
Implementering
Hela arbetet a¨ r kodat med C++. CUDA C [4] har i sin tur anv¨ants f¨or att kunna implementera
ber¨akningarna p˚a grafikkortets processor. Fluiden kan visuliseras antingen i realtid med hj¨alp
av OpenGL eller renderas ut offline. En iteration i fluidsimuleringen ser ut som f¨oljande:
1. R¨akna ut ett tidssteg
2. Extrapolera hastighetsf¨altet i β-bandet
3. L¨agg till externa krafter till hastighetsf¨altet
4. Advektera hastighetsf¨altet och level setet
4
5. Reinitisering av β- och γ-band
6. L¨os tryckekvation f¨or den nya fluiddom¨anen
7. Uppdatera hastighetsf¨altet genom att subtrahera tryckgradienten.
Stegen beskrivs i kapitel 3 och 4.
2.5.1
GPU-programmering
F¨or att programmera instruktioner till grafikkortet har Nvidia CUDA anv¨ants som till˚ater mer
parallellism a¨ n vad en CPU g¨or. F¨or att anropa en funktion (en s˚a kallad kernel) p˚a GPU’n,
k¨ors koden parallellt o¨ ver block. Totalt kan ca 65 000 block k¨oras, d¨ar varje block i sin tur kan
delas upp i tr˚adar. Antal tr˚adar beror p˚a grafikkortet och f¨or att g¨ora fluidl¨osaren kompatibel
med en st¨orre m¨angd grafikkort, anv¨ands 256 tr˚adar per block. N¨astan alla algoritmer k¨ors
o¨ ver en specifik kvantitet (t.ex. hastigheter, tryck eller level set). Varje kernel anropas med
tillr¨ackligt m˚anga block s˚a att varje tr˚ad f˚ar en egen cell/yta att arbeta med. CUDA till˚ater att
binda minnesplatser p˚a grafikkortet till en textur f¨or snabb inl¨asning av data, det vill s¨aga, om
data p˚a en minnesplats a¨ ndras s˚a kommer texturen dynamiskt att uppdateras. Anv¨andning av
texturer a¨ r ett kraftfullt verktyg d˚a lokalitet kan utnyttjas. Om ett texturv¨arde har l¨asts i en punkt,
g˚ar det snabbt att l¨asa i samma textur i n¨arliggande punkter. Tyv¨arr g˚ar det inte i dagsl¨aget att
dynamiskt skriva till 3D-texturer i CUDA, n˚agot som hade varit optimalt f¨or en simulering i tre
dimensioner. Ist¨allet a¨ r den p˚at¨ankta 3D-texturen uppdelad i skivor som l¨aggs ut som en enda
stor 2D-textur. Lokalitet bevaras d˚a fortfarande i tv˚a dimensioner men f¨orloras i den dimension
som skivorna representerar.
Figur 2.4: Transformation fr˚an 3D-textur till 2D.
2.5.2
Offline-rendering
All offline-rendering a¨ r utf¨ord i 3Ds Max, som f˚ar in en 3D-modell att jobba med. 3D-modellen
r¨aknas ut f¨or varje bildruta med metoden Marching cubes [5], d¨ar utvalda punkter i voxelrymden
anv¨ands f¨or att skapa polygoner.
5
2.5.3
Realtidssimulering
Med hj¨alp av str˚alar kan level setet renderas ut till sk¨armen. F¨or att uppn˚a realtidssimulering
m˚aste denna process vara snabb, n˚agot som g¨or det sv˚art att sp˚ara str˚alar (raytracing) och f˚a en
fotorealistisk rendering. Ist¨allet implementeras en enkel raycaster som marscherar i en str˚ales
riktning genom level setet tills den tr¨affat en punkt som ligger p˚a ytan [6]. D¨ar anv¨ands enkel
trilinj¨ar interpolation av normaler f¨or att f˚a mjuka o¨ verg˚angar mellan voxlar. Scenen renderas
med valfri shadingmetod (f¨orslagsvis Phong-shading [7]) i tv˚a pass f¨or att a˚ terge transparens
och refraktioner. I det f¨orsta passet renderas omgivningen och den lagras i en textur. I det andra
passet renderas fluiden och beroende p˚a givet transparens-v¨arde, interpoleras de tv˚a passens
renderingar. F¨or att simulera refraktioner anv¨ands en fuskl¨osning som utnyttjar normalens
riktning mot kameran [8].
2.6
Mjukvara
• Microsoft Visual Studio 2010
• Autodesk 3Ds Max
• MATLAB
6
Kapitel 3
Grunderna hos en fluid
3.1
Navier-Stokes ekvation
F¨or att kunna simulera en fluid beh¨ovs en modell som beskriver dess beteende matematiskt.
Modellen som anv¨ands a¨ r den s˚a kallade Navier-Stokes-ekvationen (3.1) som inneh˚aller termerna
advektion, tryck, yttre krafter och viskositet.
1
∂~u
+ ~u · ∇~u + ∇~p = ~g + ν∇ · ∇~u
∂t
ρ
(3.1)
Ekvationen (3.1) m˚aste dock f¨orenklas vid simulering. Den ska ocks˚a se till att fluidens volym
bevaras. D¨arf¨or l¨aggs ekvation (3.2) till som ett villkor p˚a att fluiden ska vara okomprimerbar
det vill s¨aga att hastighetsf¨altet ska vara divergensfritt.
∇ · ~u = 0
(3.2)
En ytterligare f¨orenkling p˚a Navier-Stokes-ekvationen a¨ r att ta bort viskositetstermen
(ν∇ · ∇~u) i ber¨akningarna. Detta p˚a grund av att den inte a¨ r viktig vid vattensimulering i och
med att de numeriska felen ger nog med viskositet a¨ nd˚a. De termer som a˚ terst˚ar kommer att
f¨orklaras senare i rapporter.
Den sista f¨orenklingen a¨ r att dela upp termerna i (3.1) och l¨osa de var f¨or sig vilket kommer
underl¨atta ber¨akningarna. Detta ben¨amns splitting.
3.2
Yttre krafter
F¨or att fluiden ska kunna p˚averkas utifr˚an beh¨ovs delen i Navier-Stokes-ekvationen som behandlar
yttre krafter. Utan denna skulle fluiden aldrig r¨ora p˚a sig. Yttre krafter kan innefatta t.ex.
gravitation eller vind. I det h¨ar projektet a¨ r den enda yttre kraften gravitation. Eftersom fluiden
inte har n˚agon virtuell massa kan vi se gravitationskraften som ren acceleration, som i sin
tur kan omvandlas till hastighet genom att integreras o¨ ver tiden. D˚a hela simuleringen sker i
diskreta steg f˚as den nya hastigheten genom ekvation (3.3), d¨ar ~g a¨ r gravitationen:
~un+1 = ~un + ~g · 4t
7
(3.3)
3.3
Advektion
Advektionsdelen i Navier-Stokes ekvation handlar om att propagera hastigheterna i hastighetsf¨altet.
N¨ar fluiden r¨or p˚a sig m˚aste hastighetsvektorerna, i de voxlar som fluiden sveper o¨ ver, f˚a nya
v¨arden, s˚adana att fluiden beh˚aller sin hastighet. Rent matematiskt skulle det kunna st¨allas upp
som (3.4):
∂~qn
)
(3.4)
~qn+1 = advektera(~qn , 4t,
∂t
S˚a hur g˚ar det till? Plockas advektionsdelen ut ur (3.1) kan den st¨allas upp p˚a f¨oljande vis:
∂~un+1
= −(~un · ∇) · ~un
∂t
(3.5)
Eftersom hur sm˚a tidssteg som helst inte kan tas a¨ r det, i detta fall, inte l¨ampligt att anv¨anda sig
av vanlig Euler-integration d˚a metoden a¨ r instabil. Till slut skulle v¨arden bli extremt felaktiga,
n¨ar tidssteget blir stort, och fluiden skulle inte l¨angre bete sig normalt.
Vad som ist¨allet utf¨ors a¨ r att en fiktiv partikel f˚ar sp˚aras ett tidssteg bak˚at i rymden. Den h¨ar
metoden a¨ r stabil f¨or alla tidssteg och kallas Semi-Lagrangeisk metod [9] – semi, just f¨or att
n˚agon verklig partikel aldrig skapas. D˚a metoden anv¨ander sig av tv˚a tidpunkter f¨or att fungera,
kr¨avs ocks˚a tv˚a kopior av hastighetsf¨altet; ett vars v¨arden uppdateras, och anv¨ands som sikte
bak˚at i tiden samt ett d¨ar gamla v¨arden h¨amtas ifr˚an.
Figur 3.1: En punkt sp˚aras ett tidsteg bak˚at med den bl˚a hastighetsvektorn. I den punkten
interpoleras en hastighet fram med hj¨alp av omkringliggande hastigheter (r¨oda vektorer).
Att bara g˚a ett enstaka steg bak˚at i tiden och plocka ett v¨arde d¨ar a¨ r oftast inte s¨arskilt
precist. F¨or att o¨ ka precision har en integrationsmetod av h¨ogre ordning implementerats, Runge-Kutta
ordning 4. Men Eulers stegmetod ger vid bak˚atsteget ett acceptabelt resultat.
N¨ar en punkt bak˚at i tiden har identifierats g¨aller det att f˚a ut en hastighet fr˚an just den
punkten. H¨ogst troligt a¨ r att en punkt ligger mellan voxelytor d¨ar hastigheter finns sparade,
vilket inneb¨ar att interpolation (se avsnitt 2.4) kr¨avs.
F¨or att f˚a ett a¨ nnu b¨attre resultat i advektionssteget anv¨ands a¨ ven en metod som kallas
Back-and-Forth Error Compensation and Correction (BFECC) [10]. Vad detta inneb¨ar a¨ r att
hastigheten som anv¨ands f¨or att advektera bak˚at i tiden korrigeras med ett antal steg, enligt
8
ekvation (3.6).
φˆn+1 = Advektera(φn )
φˆn = Advektera−1 (φˆn+1 )
φn = (3φn − φˆn )/2
2
n+1
φ
(3.6)
= Advektera(φn2 )
D¨ar φ a¨ r hastigheten som ska advekteras. Advektera−1 betyder att advektionen skall utf¨oras
fram˚at i tiden ist¨allet f¨or bak˚at.
3.4
Extrapolering av hastigheterna
En oklarhet kvarst˚ar; vad h¨ander om advektionsdelen i (3.1) blir tillbedd att h¨amta ett v¨arde
utanf¨or fluiden, d¨ar inga hastigheten a¨ r definierade? Det inneb¨ar problem. Man vill allts˚a se
till att det alltid finns hastigheter en liten bit utanf¨or fluidens yta. Hastigheter extrapoleras hela
tiden ut fr˚an fluiden s˚a att det alltid finns v¨arden att plocka s˚a l˚angt ut som till gr¨ansen p˚a
level setets β-band (se avsnitt 4.2). Vilka hastigheter som a¨ n befann sig utanf¨or fluiden blir
o¨ verskrivna under extrapoleringen. Ekvation (3.7) anv¨ands f¨or att utf¨ora detta repetitivt.
∂~u
= −ˆ
n · ∇~u
∂t
(3.7)
d¨ar ~u a¨ r hastighetsf¨altet och n
ˆ a¨ r normalerna till ytan. Detta kan kombineras med level setets
reinitiering (se avsnitt 4.4), som vet var de n¨armsta punkterna p˚a ytan finns och en hastighet
kan h¨amtas.
3.5
Projektion
Projektionen handlar om att l¨osa tryckdelen (3.8) av Navier-Stokes-ekvationen som g¨or hastighetsf¨altet
divergensfritt. D¨ar h¨ogerledet i (3.8) best˚ar av densiteten ρ och tryckf¨altet p~. I det h¨ar tillst˚andet
har advektion samt externa krafter applicerats p˚a ~u. Projektionen a¨ r den ber¨akningstyngsta
delen av simuleringen i och med att trycket i vektorf¨altet a¨ ndras samtidigt som divergensvillkoret
(3.2) och randvillkoret ska uppfyllas.
1
~un+1 = ~u − ∆t ∇~p
ρ
(3.8)
Eftersom p~ i ekvation (3.8) a¨ r ok¨and m˚aste detta tryckf¨alt best¨ammas. Helmholtz-Hodge
dekomposition m˚aste utnyttjas vilken s¨ager att ett vektorf¨alt kan delas upp i ett divergensfritt
samt rotationsfritt f¨alt.
~u = ~ud + k∇~p
(3.9)
D¨ar ~ud i a¨ r den divergensfria delen och ∇~p a¨ r den rotationsfria delen i ekvation (3.9). B˚ade
~ud och ∇~p a¨ r ok¨anda. Om vi l˚ater divergensoperatorn operera p˚a (3.9) f˚as:
∇ · ~u = ∇ · ~ud + k∇2 p~ = k∇2 p~
(3.10)
Eftersom den ok¨anda ~ud termen a¨ r divergensfri f¨orkortas den bort fr˚an ekvation (3.10) och
nu kan p~ l¨osas ut.
9
∆t 2
∇ p~ = ∇ · ~u
(3.11)
ρ
N¨asta steg a¨ r att l¨osa Poissons ekvationen (3.11). F¨or att kunna g¨ora det numeriskt m˚aste
ekvationen diskretiseras och vi b¨orjar med h¨ogerledet som visas i ekvation (3.12) d¨ar u, v, w a¨ r
x-, y- och z-komponenterna f¨or vektorerna i vektorf¨altet ~u. Divergensen r¨aknas bara ut i celler
som inneh˚aller fluid.
ui+1/2,j,k − ui−1/2,j,k vi,j+1/2,k − vi,j−1/2,k wi,j,k+1/2 − wi,j,k−1/2
+
+
(3.12)
(∇ · ~u)i,j,k ≈
∆h
∆h
∆h
N¨asta steg a¨ r att diskritisera v¨ansterledet som blir f¨oljande
∆t 2
∆t −6pi,j,k + pi+1,j,k + pi,j+1,k + pi,j,k+1 + pi−1,j,k + pi,j−1,k + pi,j,k−1
− ∇ p~ =
ρ
ρ
∆h2
(3.13)
Den riktiga betydelsen av v¨ansterledet a¨ r att laplaceoperatorn beskriver v¨axlingen mellan de
olika elementen i voxeln, som bearbetas, och dess grannar. Den kan a¨ ven skrivas i vektorform
vilket kommer beh¨ovas n¨ar Poissons ekvation ska l¨osas.
H¨ogerledet uttryckt i vektorer:


pi,j,k
pi+1,j,k 


pi,j+1,k 


∆t
∆t
pi,j,k+1 
−6
1
1
1
1
1
1
(3.14)
− ∇2 p~ =


ρ
ρ∆h2
pi−1,j,k 


pi,j−1,k 
pi,j,k−1
Med h¨ogerledet och v¨ansterledet diskretiserade a˚ terst˚ar det att l¨osa Poissons ekvation (3.11)
och f¨or det m˚aste a¨ ven den skrivas i vektorform/matrisform:
A~x = ~b
(3.15)
I det linj¨ara systemet (3.15) a¨ r A = ∇2 , ~x = p~ och ~b = ∇ · ~u. Detta ger att ~x ska l¨osas ut, det
vill s¨aga inversen till A ska hittas.
F¨or att f˚a en inblick i vad som egentligen h¨ander i ekvation 3.14 s˚a f¨orklaras
detta med ett
exempel. Det intressanta i ekvationen a¨ r radvektorn −6 1 1 1 1 1 1 i ekvation (3.14),
som refererar till om en voxel och dess grannar inneh˚aller fluid, luft eller solid. Ettor betyder att
voxlarna inneh˚aller fluid eller solid och nollor om voxlarna inneh˚aller luft. Om en voxel utg¨ors
av en solid beh¨ovs ett uttryck (3.16) f¨or att ber¨akna trycket i voxeln och samtidigt uppr¨atth˚alla
Neumans villkor [11]. Om en voxel till h¨oger a¨ r solid, f˚as trycket i denna voxel genom:
pi+1,j,k = pi,j,k +
ρ∆h
(ui+1/2,j,k − usolid )
∆t
(3.16)
Ekvation (3.16) ers¨atter d˚a pi+1,j,k i ekvation (3.14).
Ett 2D-exempel ser ut som f¨oljande:


pi,j
∆t −4 1 0 1 1 

ρ∆h2

10

pi,j
+ ρ∆h
(ui+1/2,j − usolid )

∆t

pi,j+1


pi−1,j
pi,j−1
(3.17)
D¨ar pi+1,j a¨ r en solid och pi,j+1 a¨ r fylld med luft vilket representeras av en nolla i radvektorn.
Figur 3.2: Bild till exempel (3.17).
11
3.6
Divergensfrihet
Ett villkor, som a¨ r mycket viktigt n¨ar fluider som vatten simuleras, a¨ r att den ska vara inkompressibel.
Volymen ska med andra ord vara konstant hela tiden. Ett villkor som kanske ses som sj¨alvklart
men som faktiskt inte a¨ r helt trivialt.
Divergens i ett vektorf¨alt kan ses som fl¨odet i en viss punkt. Om divergensen a¨ r positiv
¨ fl¨odet mot punkten st¨orre a¨ n fl¨odet
str¨ommar det ut mer fr˚an punkten a¨ n vad det kommer in. Ar
ut a¨ r divergensen negativ. Under simuleringen m˚aste det g¨alla att hela fluiden a¨ r divergensfri,
det vill s¨aga, differensen mellan fl¨odet in och ut a¨ r lika med noll.
P˚a grund av att det blir en hel del approximationer under simuleringen a¨ r det sv˚art att beh˚alla
en konstant volym p˚a fluiden. Ofta komprimeras den och det ser ut som att v¨atska f¨orsvinner
sp˚arl¨ost. Ett trick f¨or att l¨osa detta problem att g¨ora en direkt modifiering p˚a Poisson-ekvationen
under tryckber¨akningarna [12]. Om fluidens volym approximativt kan m¨atas under varje tidpunkt,
genom att r¨akna hur m˚anga voxlar som har ett negativt v¨arde, kan denna j¨amf¨oras med volymen
vid initiering. Differensen kommer att ge hur mycket volym som saknas och denna differens
skalas med en konstant och anv¨ands i divergensdelen i Poisson-ekvationen.
3.6.1
Dirichlets randvillkor
Anledningen till att en fluid inte kan str¨omma in i en solid a¨ r att soliden alltid kan mots¨atta
ett lika h¨ogt tryck mot fluiden som fluiden utg¨or p˚a den. Detta m˚aste s˚aklart a¨ ven g¨alla i
simuleringen om vattnet ska stanna kvar i akvariet. F¨or att till¨ampa detta anv¨ands Dirichlets
randvillkor (3.18) [11], som s¨ager att hastigheter i n¨arliggande soliders motsatta normalriktning
m˚aste vara noll. Med andra ord, det f˚ar inte existera n˚agra vektorkomposanter, i hastighetsf¨altet,
som pekar in mot en solid.
Figur 3.3: Bild p˚a innan och efter Dirichlets randvillkor uppfyllts
N¨ar soliden i fr˚aga har en hastighet s˚a g¨aller:
~u · n
ˆ = ~usolid · n
ˆ
(3.18)
d¨ar n
ˆ a¨ r normalen till soliden, ~u a¨ r fluidens hastigheter och ~usolid a¨ r solidens hastighet. Eftersom
det h¨ar projektet endast innefattar simulering med stillast˚aende solider kan f¨oljande ekvation
allts˚a anv¨andas:
~u · n
ˆ=0
(3.19)
Det a¨ r o¨ nskv¨art att Dirichlets randvillkor a¨ r uppfyllt under hela simuleringen. D¨arf¨or k¨ors
operationen flera g˚anger f¨or varje bildruta, mellan de olika utr¨akningarna.
12
Kapitel 4
Level set-metoden
Level set-metoden handlar om att f¨orklara fluidens form och position med implicita ytor. Det
h¨ar kapitlet f¨orklarar endast de allra grundligaste stegen f¨or att anv¨anda metoden.
4.1
Implicita ytor
Med hj¨alp av level set-metoden kan fluidens form hela tiden representeras som en implicit yta.
En implicit yta har bland annat den f¨ordel att den inte kan skapa tvetydighet om var ytan a¨ r
genom att korsa sig sj¨alv. N¨ar tv˚a ytor o¨ verlappas sl˚as de ihop och bildar en ny yta. Avst˚andet
till fluidens yta ges genom ekvationerna (4.1) och (4.2):
Distansinne = {x(t) ∈ <2 : φ(x(t)) ≤ 0}
(4.1)
Distansute = {x(t) ∈ <2 : φ(x(t)) > 0}
(4.2)
Detta ger att n¨armsta avst˚andet φ till fluidens yta kan hittas oavsett position i rymden. Ett
negativt φ-v¨arde betyder att punkten a¨ r inuti fluiden och ett positivt v¨arde utanf¨or. N¨ar φ a¨ r
exakt noll ligger punkten p˚a ytan.
4.2
β- och γ-band
Vid simulering a¨ r det endast intressant att veta var fluidens yta befinner sig. D¨arf¨or beh¨ovs det
bara ett definierat level set i omgivningen av ytan. Ett β-band begr¨ansar level setets utstr¨ackning
utanf¨or fluidens yta, d¨ar konstanten β avg¨or hur l˚ang denna utstr¨ackning a¨ r. Beroende p˚a
integrationsmetod i advektionen (se avsnitt 4.3) kr¨avs olika v¨arden p˚a β. F¨or ett vanligt Euler-steg
r¨acker det med att β har ett v¨arde motsvarande 2 · ∆h. F¨or att f¨ors¨akra sig om att korrekta
derivator finns i β-bandet utvidgas level setet med ytterligare ett band, γ-bandet, som i regel a¨ r
smalt.
13
Figur 4.1: Allt inom den svarta gr¨ansen a¨ r fluid. Det gr¨ona f¨altet a¨ r β-band och det streckade
a¨ r γ-band.
4.3
Advektion
Enligt avsnitt 2.5 uppfylls CFL-tillst˚andet med valt tidssteg. Detta a¨ r viktigt f¨or att level setet
inte ska kollapsa. F¨or att propagera ytan i tiden utf¨ors ett advektionssteg utifr˚an det divergensfria
hastighetsf¨altet med ett enkelt Euler-steg. Advektion utf¨ors bara p˚a celler som ligger innanf¨or
β-bandet.
φ(n + 1) − φ(n)
∂φ
≈
(4.3)
∂t
∆t
∂φ
= −~u · ∇φ
(4.4)
∂t
Gradienten till φ beh¨ovs allts˚a r¨aknas ut och f¨or det anv¨ands ett ”upwind scheme”, som i en
dimension ser ut som:
(
φ+
if ux < 0
∂φ
x = (φi+1,j,k − φi,j,k )/∆x,
≈
(4.5)
−
∂x
φx = (φi1,j,k − φi−1,j,k )/∆x, if ux > 0
Detta s¨ager att information endast f¨ardas i hastighetsf¨altets riktning och derivator tas fram
genom att bara kolla p˚a skillnader d¨ar fluiden har befunnit sig.
14
4.4
Reinitiering
B˚ade Eulersteget och approximeringen av gradienten introducerar fel i φ(n + 1). Eftersom
advektionen dessutom endast utf¨ors i punkter inom β-bandet (se figur 4.1) initieras ett nytt
β-band i varje iteration i fluidsimuleringen. Detta kan g¨oras med en geometrisk metod [13],
som tar fram vilka punkter som ligger i n¨arheten av ytan. Utifr˚an dessa punkter utvidgas level
setet iterativt tills hela β-bandet a¨ r fyllt. D¨arefter skapas a¨ ven ett nytt γ-band.
15
Kapitel 5
Resultat
5.1
Realtid
Resultatet blev en vattenliknande fluid i tre dimensioner. Fluiden liknar vatten, inte bara visuellt
utan ocks˚a rent beteendem¨assigt, vilket fortfarande a¨ r knepigt att uppn˚a med denna metod. Vad
som vanligtvis intr¨affar a¨ r att fluiden f˚ar en alldeles f¨or h¨og viskositet och liknar mer glycerol
eller sirap. Eftersom metoden inte a¨ r partikelbaserad ges dock inga riktiga st¨ank-effekter, som
riktiga vattenr¨orelser brukar ge upphov till. Bild p˚a fluiden i realtid:
Figur 5.1: Resultat
16
5.2
Utrenderat
Bild p˚a hur fluiden ser ut n¨ar den f˚att renderas ut i 3Ds Max, med fysikaliskt korrekt hdri-belysning
och ray-tracing-metoden Mental ray:
Figur 5.2: Resultat
17
Kapitel 6
Avslutning
6.1
Diskussion
Fluidsimulering a¨ r minst sagt en komplicerad process. Det finns en hel del information att
hitta om level set-metoden, b˚ade v¨alskrivna rapporter och hemsidor. Det knepiga och mest
motstr¨aviga under resans g˚ang har dock varit att implementera allt p˚a datorns grafikkort. Det
blir ett helt annorlunda t¨ank n¨ar ber¨akningar ska parallelliseras och optimeras f¨or att samma
funktion ska kunna ge olika resultat p˚a olika k¨arnor, v¨arden ska sparas i texturer ist¨allet f¨or
arrayer och s˚a vidare.
Ett kvarst˚aende problem a¨ r simuleringens alla approximationer. Numeriska fel, i f¨oljd av
diskretiseringarna, g¨or att det kommer att bli vissa ej helt korrekta beteenden hos fluiden,
bland annat volymf¨orlusten. Det som begr¨ansar hur fina ber¨akningar som kan g¨oras a¨ r datorns
prestanda.
6.2
Vidare f¨orb¨attringar
Det finns ett stort antal f¨orb¨attringar, som kan g¨oras med denna fluidsimulering, f¨or att ge den
ett a¨ n mer vattenliknande beteende. Det finns a¨ ven a˚ tg¨arden som kan g¨ora den mer fysikaliskt
korrekt.
Det f¨orsta som kommer p˚a tanke a¨ r problemet med volymf¨orluster. S¨attet detta a¨ r l¨ost p˚a
a¨ r ett mycket overkligt s˚adant. Att skapa volym ur ingenting genom att f¨orst¨ora en fysikaliskt
korrekt ekvation a¨ r inte optimalt. Det a¨ r dock en simpel l¨osning som fungerar bra rent visuellt,
och i m˚anga fall handlar datagrafik bara om att n˚agonting ska se bra ut. Den h¨ar metoden skulle
inte vara acceptabel i forskningssammanhang.
N¨ar det handlar om att f˚a till effekter som st¨ank och separationer av fluiden blir man i princip
tvungen att bygga vidare p˚a arbetet med en metod som kallas FLIP (fluid-implicit partikel). Det
a¨ r en Eulerisk fluidl¨osare som anv¨ander partiklar f¨or advektion. [14]
F¨or varje kernel anropar CUDA tillr¨ackligt m˚anga block s˚a att varje tr˚ad f˚ar enbart en cell
att arbeta p˚a. Det hade varit intressant att unders¨oka om det hade g˚att snabbare i realtid om varje
tr˚ad f˚att jobba med flera celler a˚ t g˚angen [15]. Styrkan i denna metod hade varit m¨ojligheten
till att lagra v¨arden i ”Shared Memory” vilket a¨ r ett minne som varje block delar mellan sig och
dessutom a¨ r betydligt snabbare a¨ n att l¨asa fr˚an textur och globalt minne [16].
Realtidsrenderingen lider tyv¨arr lite av artefakter som orsakats av cellstorleken, trots att
trilinj¨ar interpolation anv¨ants. Det borde vara m¨ojligt att innan varje rendering skapa en statisk
3D-textur av level setet, raycasta den och anv¨anda trikubisk interpolation av normaler vid ytan.
18
D˚a g˚ar det att anv¨anda grafikkortets egna process av filtrering (som a¨ r lika snabb som vanlig
texturl¨asning) och anv¨anda endast 8 texturl¨asningar ist¨allet f¨or 64 som trikubisk interpolation
vanligtvis kr¨aver [17].
19
Litteraturf¨orteckning
[1] Osher S, Fedkiw R. Level Set methods and Dynamic Implicit Surfaces. Los Angeles:
Springer; 2003
[2] Bridson R. Fluid Simulation for Computer Graphics. s. 21-25. Wellesly: A K Peters; 2008
[3] Fedkiw R, Stam J, Jensen H. W. Visual Simulation of Smoke. Standford University. http:
//physbam.stanford.edu/˜fedkiw/papers/stanford2001-01.pdf;
Senast h¨amtad 2011-03-15
[4] Sanders J, Kandrot E. CUDA by Example - An Introduction to General-Purpose GPU
Programming. Boston: Addison-Wesley; 2010
[5] Bourke P. Polygonising a Scalar Field. 1994. http://paulbourke.net/
geometry/polygonise/; Senast h¨amtad 2011-03-08
[6] Hadwiger M, Sigg C, Scharsach H, B¨uhler K, Gross M. Real-TIme Ray-Casting and
Advanced Shading of Discrete Isosurfaces. Eurographics: M. Alexa, J. Marks; 2005
[7] Watt A. 3D Computer graphics. s. 171-183. Harlow: Pearson, Addison-Wesley; Third
edition 2010
[8] Crane K, Tariq S. Real-Time Simulation and Rendering of 3D-fluids. Nvidia. http://
http.developer.nvidia.com/GPUGems3/gpugems3_ch30.html; Senast
h¨amtad 2011-03-08
[9] Stam J. Stable Fluids. Alias Wavefront. http://www.dgp.toronto.edu/
people/stam/reality/Research/pdf/ns.pdf; Senast h¨amtad 2011-03-08
[10] Liu Y, Dupont T. F. Back and Forth Error Compensation and Correction Methods for
Removing Errors Induced by Uneven Gradients of the Level Set Function. s. 311–324. J.
Comput. Phys; 2003
[11] L¨ath´en G, Nilsson O, S¨oderstr¨om A. Fluid Simulation. LiU, TNM079, Modeling and
Animation, Lab6; 2010
[12] S¨oderstr¨om A. Memory Efficient Methods for Eulerian Free Surface Fluid Animation.
Link¨opings universitet, ITN; 2010
[13] Tsai Y.-H. R. Rapid and Accurate Computation of the Distance Function Using Grids. s.
175-195. J. Computing Physics; 2002
[14] Zhu Y, Bridson R. Animating Sand as Fluid. ACM Trans. Graph; 2005
20
[15] Brandvik T, Pullan G. Acceleration of a 3D Euler Solver Using Commodity Graphics
Hardware. Cambridge, University of Cambridge; 2008. http://www.eng.cam.ac.
uk/˜gp10006/research/Brandvik_Pullan_2008a_DRAFT.pdf;
Senast
h¨amtad 2011-03-13
[16] Nvidia. CUDA Programming Guide. 2008. http://developer.download.
nvidia.com/compute/cuda/2_0/docs/NVIDIA_CUDA_Programming_
Guide_2.0.pdf; Senast h¨amtad 2011-03-13
[17] Sigg C, Hadwiger C. Fast Third-Order Texture Filtering. Nvidia. http://http.
developer.nvidia.com/GPUGems2/gpugems2_chapter20.html; Senast
h¨amtad 2011-03-08