Gerenderde animatie

Hier is een eerste gerenderde animatie te zien van ons hybride systeem. De beelden werden opgenomen aan 200FPS, en hier getoond aan slechts 50FPS, waardoor het water in slow motion beweegt.

Let op: Het aquarium is slechts een visualisatie. De werkelijke grens van het water water loopt gewoon verder boven het aquarium zelf.

Finale Renders

We hebben de afgelopen dagen hard gewerkt om de thesis af te krijgen. Het is ondertussen eindelijk af en verstuurd naar de drukker. We verwachten de afgedrukte versie van de acco tegen morgen. In de tussentijd dachten we nog onze finale renders online te gooien.








Je kan op de foto's klikken voor een hoog-resolutie versie.

Voor alle duidelijkheid: Het gaat hier dus om renders van onze simulatie met POVRAY, en niet interactieve beelden.

Buoyancy en renders

In de vorige blogpost werd aangehaald hoe een 2D-3D koppeling met een 2D hard boundary voor problemen kan zorgen. De voorgestelde oplossing bestaat erin om de hard boundary weg te laten. De particles moeten wel nog dienen om het oppervlakte detail te verhogen, dus deze worden met een artificiële buoyancy kracht terug naar het oppervlak geduwd.

Om deze implementatie mogelijk te maken moeten we allereerste een koppeling tussen particles en heightfield uitwerken. Wanneer een particle in het heightfield gaat, moet zijn massa worden toegevoegd. Dit houdt in dat 1 particle een soort uitmiddeling van zijn massa moet doen over verschillende heightfieldkolommen. Zowel een kernel die over alle particles itereert als een kernel die over alle heightfieldkolommen itereert kan dit op parallelle hardware niet verwezenlijken: er moet in beide gevallen op verschillende plaatsen tegelijk geschreven worden.

Om dit probleem op te lossen nemen we om te beginnen de particlediameter even groot als de breedte/lengte van een heightfieldkolom. We maken dan de abstractie dat een particle slechts in een heightfieldkolom tegelijkertijd kan zitten. Een kernel over alle particles berekent dan hoe diep elke particle in de heightfield zit. Hierna kan een kernel over alle heightfieldkolommen berekenen wat de uitgemiddelde massatoevoeging van ieder particle voor deze specifieke kolom is (elke particle verdeeld zijn massa over meerdere kolommen, afhankelijk van de SPH smoothinglength).

Particles vallen nu door de heightfield en de heightfield genereert golven. Om ze realistisch te laten bewegen moet nog een drag force geïmplementeerd worden die de particles in het watervolume afremt en een buoyancy force, die we groter veronderstellen dan ze echt is, zodat de particles terug aan het oppervlak komen.

Enkele dagen van implementatie en debuggen (onder andere een quick fix aan de heightfield en enkele uren undefined behaviour opsporen met Valgrind) leverde de gewenste resultaten:


Ook het NVIDIA Marching Cubes algoritme werd geïntegreerd in ons systeem (alsook een betere camera control, instelbare parameters, een heleboel refactors om de leesbaarheid te verhogen en extra particle emitters (in de vorm van kubussen en bollen)). De scene kan nu dus ook met een oppervlak bekeken worden:


Een samengestelde screenshot levert dan volgende concept-art op:


Wanneer we de data van Marching Cubes exporteren naar een povray formaat, zijn we ook in staat om mooie renders te maken. Een ray trace van dezelfde scene ziet er als volgt uit:


Dit type renders kunnen we gebruiken om nog een highres filmpje te maken.

Ook geavanceerdere ray traces zijn mogelijk, zoals met photon mapping:


Omwille van ruis kunnen dit soort renders niet gebruikt worden voor een filmpje. Voordat we een filmpje kunnen opnemen moeten we nog enkele artifacts uit de scene halen (zoals rechts op de render het zwarte lijntje, dit is een artifact van het heightfield). Ook werkt Marching Cubes nog niet goed wanneer particles verwijderd worden. Aan het algoritme zelf moet niet veel meer gedaan worden, alleen komt het met de huidige maximale tijdstap nog niet helemaal tot een rustsituatie (de particles blijven wat schommelen op het water). Wanneer deze zaken opgelost zijn kunnen we onze versnelling vergelijken met een particle-only simulatie en alles evalueren.

Hybride algoritmes

We ondervinden enkele problemen met de koppeling tussen het heightfield en de particles. Wat er eigenlijk gebeurt kun je zien in volgende screenshots:


De particles duwen op het heightfield waardoor er lokaal een put ontstaat. Hier komen steeds meer particles in, zodat de put groter wordt. Het resultaat zie je in het rechter prentje: niet echt de bedoeling dus.

Om dit effect te counteren hebben we enkele methodes uitgetest:

• Parameters aanpassen
  We kunnen bijvoorbeeld een grotere graviteitsconstante of kleinere timestep nemen voor het heightfield. Hierdoor zal hij stroefer bewegen maar wel tot een horizontaal evenwicht komen. Deze methode werd toegepast in het eerder geposte filmpje. We willen toch liever niet met fysische waarden prutsen dus gingen we opzoek naar andere technieken.
• Krachten van heightfield op particles
  We kunnen berekenen welke versnelling het heightfield kent en deze omzetten naar krachten. Deze krachten zorgde er wel voor dat de particles minder hard duwden, maar een horizontaal evenwicht is nog steeds niet gegarandeerd: van zodra teveel particles bij elkaar komen ontstaat er weer een put.
• Herverdeling van het water
  Tijdens de integratie fase van het heightfield moet verdwenen volume gecompenseerd worden (het volume werd eerst afgenomen op plaatsen waar particles op de heightfield duwen). Dit gebeurde eerst uniform maar we kunnen lager water meer voordeel geven. Hier zal dan sneller water toegevoegd worden waardoor het heightfield wel terug naar een horizontale toestand evolueert. De resultaten vind je hieronder.

In het eerste plaatje zie je een golf die terug naar de gemiddelde hoogte wordt getrokken (helemaal rechts is de hoogte het laagste en daar zal volumecompensatie optreden). Het heightfield zal hierdoor niet heel drastisch bewegen, zoals je ook kan zien in het rechter plaatje. De bewegingen sterven nu heel snel uit.

We denken dat we beter resultaten zullen krijgen wanneer de heightfield geen hard boundary is. We laten de particles er gewoon doorvallen en door middel van een artificiële buoyancy kracht duwen we ze terug naar het oppervlak. Hoe we dit net op de grafische kaart gaan implementeren volgt in een volgende blogpost, samen met de resultaten. Hieronder volgt nog een schets van het huidige algoritme:

Marching Cubes

Terwijl Daniel verder werkt aan de koppeling tussen het heightfield en de particles heb ik eens snel naar een render techniek gekeken. Met Marching Cubes samplen we het density veld van onze simulatie op de hoekpunten van een kubus. Ofwel valt een hoekpunt binnen het watervolume ofwel niet. Afhankelijk van de configuratie van de verkregen kubus berekenen we triangles en normalen volgens volgende tabel:



Een CUDA GPU implementatie van dit algoritme was reeds beschikbaar in de CUDA SDK. Ik heb deze geïntegreerd in ons particle systeem om het resultaat eens te bekijken. Een simulatie van 6000 particles en een 32x32x32 Marching Cubes grid draait nu aan 15 FPS. Meer particles nemen heeft niet veel effect, maar dan is er wel een groter grid nodig zodat het oppervlak niet te ruw wordt. De gridgrootte van het Marching Cubes algoritme is duidelijk de bottleneck.

Dubbele koppeling sph - heightmap

Hier is een video om de dubbele koppeling tussen heightmap en SPH particles te demonstreren.



De heightfield wordt omlaag geduwd op plaatsen met hoge druk, en compenseert dit volume verlies met de verhoging van het oppervlak op plaatsen met lage druk. Wederzijds is er een kracht die particles uit het heightfield duwt aan de hand van de normaal. Om te voorkomen dat particles in het wateroppervlak vallen is er ook een harde boundary voorzien moest de kracht van heightmap op particle niet groot genoeg zijn.

Momenteel is er nog een onopgeloste stabiliteit wanneer het heightfield te laag komt te staan. Op dit punt is het systeem namelijk numeriek instabiel wegens kleine massa per heightfield kolom. Dit probleem kan voorkomen bij grote drukvariaties in kleine volumes water en kan leiden tot ontploffingen van de heightmap.

Binnenkort gaan we een video maken van dezelfde opstelling, maar dan enkel met particles, om het verschil te demonstreren.

Verder bestaat er al een implementatie van het algoritme dat particles omzet tot heightfield om de hoeveelheid particles in het systeem te beperken. Ook hiervan zullen we binnenkort een video laten zien, nadat we hebben geëxperimenteerd met een aantal heuristieken.

Shallow Water

Onze shallow water implementatie is ondertussen af. Een voorsmaakje:


Het is de bedoeling om het wateroppervlak te samplen met particles, waardoor het zelfs niet gerendered zal worden. De interface tussen deze 2D en 3D technieken is ook bijna af, we werken deze week aan de integratie van beide technieken.

Deze shallow water implementatie gebruikt Jacobi Iteraties om het stelsel op te lossen, maar er bestaan geavanceerdere, snellere technieken. Eventueel zouden we deze ook kunnen implementeren.

Status

Vorige semester bestudeerden we 3D particle based fluid simulations voor interactieve toepassingen. Hiertoe gebruikten we de grafische kaart als performante rekeneenheid. Dit semester gingen we opzoek naar een koppeling tussen 3D en 2D fluid simulation technieken.

2D technieken
Deze zijn beperkt maar kunnen wel op een efficiënte manier wateroppervlaktes simuleren. Om te beginnen zochten we een goede 2D representatie. We onderzochten volgende alternatieven:

• Wave Particles: hier worden ook particles gebruikt wat de interactie tussen onze particles wel eenvoudig kon maken. Maar de golven die hiermee gemodelleerd worden zijn niet het type naarwaar wij opzoek waren en boden weinig flexibiliteit.
• Heightmap: deze structuur is heel eenvoudig te implementeren. Het water wordt voorgesteld uit kolommen, waar makkelijk water aan toegevoegd kan worden, of uit weggenomen kan worden.
• Deep water: hierbij modelleren we golven. De interactie met particles is niet vanzelfsprekend.
• Shallow Water: een techniek gelijkaardig aan de heightmap, alleen is deze meer fysisch correct en wordt er ook rekening gehouden met de zeebodem. We opteerden voor deze techniek

Shallow Water
De Navier-Stokes vergelijkingen worden voor Shallow Water vereenvoudigd met de volgende veronderstellingen:

• Het water kan voorgesteld worden door een heightfield (dus geen splash of overhellende golven)
• De verticale snelheidscomponent van de watermoleculen kan genegeerd worden.
• De horizontale snelheidscomponent is in een verticale kolom overal gelijk

Om de bekomen vergelijkingen op de GPU op te lossen kunnen we Jacobi Iteration gebruiken.

2D-3D koppeling
Om toch kleinschaligere effecten te bereiken zullen we het wateroppervlakte van onze heighfield representatie samplen met SPH particles. Deze zorgen voor overhellende golven en splash. De particles liggen als het ware op het heightfield. Wanneer er veel water aanwezig is worden de particles geabsorbeerd door het heightfield. Zo zal de simulatie geen performantie verliezen bij grote volumes water.



Enkele mogelijkheden die deze koppeling zou bieden:

• Het vullen van een groot volume water
• Het simuleren van een rondvarende boot, met realistische golven en toch kleinschalige effecten zoals splash
• Breaking waves op een strand

Op het ogenblik zijn we de koppeling tussen beide technieken aan het uitwerken. Hiervoor moeten we nog twee belangrijke dingen realiseren:

• Particles opnemen in het heightfield wanneer er teveel zijn of wanneer ze praktisch stil liggen op het oppervlak
• Een vaste band van particles voorzien rond het heightfield

Heightmaps

Hier een eerste filmpje van een heightmap implementatie. De heightmap is 512x512 en draait heel vlot op mijn laptop zelfs.

De implementatie is geinspireerd door de siggraph course notes van Robert Bridson (2007). Het is een heel eenvoudige manier om een willekeurige heightmap om te vormen tot een realtime water simulatie.

In onze implementatie wordt elke pixel in een aparte thread berekend. Er wordt per pixel een snelheid berekend afhankelijk van de buren van die pixel. Aangezien er geen interactie is tussen de threads onderling, zijn er ook geen synchronisatieproblemen, en draait het geheel redelijk vlot op de gpu.

Het voordeel van dit systeem is dat het heel makkelijk is om massa toe te voegen of weg te nemen. In de video worden er willekeurige sinussen toegevoegd (en soms weggenomen) op bepaalde punten.

Het nadeel is missschien dat het geheel nogal visceus lijkt te zijn. Ook is men heel beperkt in het realisme dat men kan voorstellen, maar dit is inherent aan heightmap methodes.

Over het algemeen zou de simulatie iets realistischer kunnen gemaakt worden door betere initiele golffuncties. Ik denk dan meteen aan de paper Wave Particles van Cem Yuksel, waarbij men een wave particle om doet gaan tot een golffront met zowel transversaal als longitudinale component.

De interactie tussen heightmap en particles moet nog uitgedacht worden. Ook de visualisatie van heightmaps moet een beetje beter. Momenteel wordt er een mesh gemaakt van de oppervlak, en worden hoogtes aangegeven door de kleur. Blauw is laag, en wit is hoog. Ik zal proberen in de nabije toekomst een kleine upgrade te doen door phong shading te gebruiken met eventueel enviroment maps voor de reflectie van het water. Realtime caustics zou natuurlijk ook mooi zijn, maar daar heb ik voorlopig nog geen onderzoek naar verricht.

Uniform grid met sorting

Het is hier even stil geweest maar ondertussen is er toch wat werk gebeurd. In de eerste plaats heb ik eens onderzocht hoe de particle demo van CUDA zijn uniform grid opstelt en dit ook geïmplementeerd.

Het uniform grid wordt bijgehouden met buckets. Elke bucket kan een onbeperkte grootte hebben. Het principe wordt hier duidelijk geïllustreerd:

Door de particles te sorteren per gridcell kunnen we ook terugvinden waar er net een nieuwe gridcell begint. Het sorteren verloopt volgens het parallelle algoritme "radix sort" dat door Scott Le Grand (NVIDIA Corporation) gebruikt wordt in "Broad-Phase Collision Detection with CUDA" voor GPU GEMS 3. Dit algoritme werd reeds geïmplementeerd in de CUDA SDK dus dat konden we gewoon hergebruiken.

De implementatie is efficiënter dan wat we al hadden met de techniek van Harada. We moeten minder geheugen gebruiken, en als we geheugen uitlezen kan dit ook 1D gecached worden vermits de particles gesorteerd voorkomen (met andere woorden, buren liggen ook naast elkaar in het geheugen). Een simulatie van 65k particles loopt bijvoorbeeld aan 60 FPS op een GTX 260 waar we vroeger 50 FPS haalden op de GTX 280. Een vergelijking van beide methodes op de benchmark machine volgt nog ten gepaste tijden.

Ik heb ook nog naar de pic&flip implementatie gekeken. Door deze radix sort techniek met buckets zou de eerste stap van het algoritme geen synchronisatie problemen meer hebben en dus veel sneller zijn.
Maar de pic&flip methode zelf blijkt toch onstabiel. Wanneer de viscositeit te laag is bewegen de deeltjes gewoon de hele tijd door elkaar, wanneer de viscositeit hoger is dan aligneren de particles zich tussen gridcellen en vertonen ze hoge affiniteit voor cellen aan de rand van het oppervlak. Deze cellen hebben een lage snelheid en de viscositeit zorgt ervoor dat naburige cellen ook meer aan deze valse snelheid gaan bewegen. We kunnen cellen met lage eigenschappen filteren, maar een onderscheid maken tussen cellen met artificiële eigenschappen en cellen met echte eigenschappen is niet evident.

Van mijn CUDA uniform grid implementatie heb ik ook nog een filmpje gemaakt:

In de helft van het filmpje laat ik de zwaartekracht meebewegen met mijn view rotatie waardoor het water naar de zijkant stroomt. Dit zorgt voor een mooie golf.

Presentatie

De presentatie staat online, je kan hem hier terugvinden. Tussentijdse presentatie.

De belangrijkste waarden voor uniform grid, zijn de volgende:
De interpretatie van deze tabel is het volgende: Als men 30720 particles simuleert op een Nvidia GTX 280, dan heeft men een simulatie die 13.75 trager loopt dan de werkelijkheid.

Op naar wave particles!

PIC & FLIP met SPH

Hoe hebben we PIC & FLIP geïntegreerd met SPH? Even een algemeen overzicht van het algoritme dat uit 3 stappen bestaat:

Stap 1: mapping
Alle particles hebben een massa en een velocity. Elke particle mapped deze waarden op een gridstructuur:

• Density: met behulp van SPH kernel. Te kleine waarden (bijvoorbeeld alles kleiner dan 1) en te grote waarden (groter dan de reference density of dus de densiteit van water) worden geclamped. Extreme waarden zorgen namelijk voor instabiliteit bij het berekenen van de krachten in stap 2.
• Velocity [x y z]: met behulp van een gewogen gemiddelde van particles die binnen de smoothing range van het grid vallen. Eventueel zou dit ook met kernels kunnen maar dan moeten extreme waarden ook gefilterd worden.

Mapping van een particle density op het grid.

Stap 2: krachten tussen gridcellen uitrekenen
Deze waarden kunnen we nu gebruiken om de SPH krachten (druk, viscositeit en externe krachten) op het grid te berekenen. We gebruiken dezelfde smoothing length die we hadden voor de particles. De afstand tussen twee gridcellen nemen we namelijk (ongeveer) gelijk aan de afstand tussen twee particles. Een hogere gridresolutie is ook mogelijk, maar dan moeten we wel opletten dat de waarden die we in deze stap berekenen geschaald worden. We willen eigenlijk dat het aantal cellen die een kracht op een andere cel uitoefenen gelijk blijft ongeacht de gridresolutie. Hiertoe kunnen we een kleinere smoothing length nemen, maar dan smoothen we alles teveel uit. Het is beter om meer cellen te beschouwen en daarna te herschalen.
Om de SPH krachten uit te kunnen rekenen hebben we de massa nodig van een gridcell. Deze kunnen we berekenen aan de hand van het volume van een cell en de densiteit op het grid.

Stap 3: interpolatie naar particles
We interpoleren de gridkrachten trilineair terug naar de particles, waarna deze krachten geïntegreerd worden om de nieuwe snelheid te bepalen. In de FLIP methode moeten we enkel het verschil in snelheid terug interpoleren naar de particles om smoothing tussen het grid en de particles te beperken. Met SPH is dit inherent: via de krachten berekenen we eigenlijk de versnelling, of dus het verschil in snelheid.

Filmpje

Een kleine update:

Hier een filmpje van 245.760 particles in verschillende snelheden. Het filmpje is jammer genoeg niet real-time: Het heeft ongeveer 35 minuten geduurd om te berekenen (met uniform grid op de gpu). Maar toch, het effect is redelijk mooi.

Grid resoluties

Status update: onze implementaties werken nog niet helemaal, alles lijkt trager te gaan dan we verwachten. Het programmeren op de GPU is niet altijd even eenvoudig - hoewel debuggen met CUDA al redelijk wat makkelijker gaat, blijft het een uitgebreid werk.

De PIC & FLIP implementatie begint stilaan vormen te krijgen. De meeste code werd reeds vorig week geïmplementeerd en getest. De interpolatietechniek van vorige blogpost werd ook vervangen door een trilineaire interpolatie. De eerste stappen in het algoritme werken al (het mappen van particle eigenschappen op een grid en terug), alsook het framework (simulatieparameters, visualisatie met VBO's, ..). We wilden dan de nauwkeurigheid testen van het grid. Hiervoor gebruiken we 1536 particles (dat kwam goed uit met het aantal cores op mijn GPU) en twee intialisatie manieren: in damvorm (zoals tijdens de simulaties) en random. De resultaten voor verschillende gridwaarden kun je in onderstaande grafiek terugvinden:

Kleinere resoluties dan 32x32x32 hadden een nog grotere fout. De echte waarde zou 1000.0 moeten zijn, dus voor 32 als gridsize ondervinden we een fout van 9% bij random intialisatie en 3.5% bij de dam. Omdat de berekentijden tochwat oplopen (en er behoorlijk wat meer interpolaties nodig zijn: namelijk nog 3 voor elke component van de snelheid) zullen we toch met een zekere afrondingsfout moeten leven. Of dit echt slechte resultaten gaat opleveren is nog af te wachten.

Bij absurd lage gridresoluties (zoals 4) krijg ik nogwel vreemde resultaten (NaN), maar zelfs na enkele uren debuggen kon ik de oorzaak niet vinden. Daar moet zeker nog eens naar gekeken worden, het zou kunnen dat er nog fouten in de implementatie zitten. Daarnaast moet ik ook eens kijken naar de mogelijkheid van 3D textures om reeds zelf te interpoleren bij het uitlezen. Dat gaat wellicht sneller gaan dan mijn handmatige methode, maar voor nog grotere afrondingsfouten zorgen (die interpolaties zijn niet exact). Beide zaken laat ik even liggen om zo snel mogelijk een werkende simulatie af te hebben.

Kleine Update Uniform Grid

Hier een kleine update langs mijn kant:

De Uniform Grid implementatie op de GPU is nog niet volledig. Er zijn momenteel onstabiliteitsproblemen. De manier waarop ik momenteel mijn uniform grid opstel is door gebruik te maken van RGBA coordinaten van een texture. Elke pixel in die texture stelt een gridcel voor, en er kunnen bijgevolg maximaal 4 particles in een gridcel (cfr Smoothed Particle Hydrodynamics on gpus van Harada). De grootte van de gridcel is momenteel gelijk met de kernel radius, zodat er telkens 27 gridcellen per particle moeten worden overlopen.

Het eerste probleem is dat er van het begin al ongeveer 10% van de particles worden gedropt omdat ze gewoonweg niet in de grid passen. Bij volgende iteraties (na botsing met een wand) gaat het water zich nog harder samendrukken, waardoor er nog meer particles moeten worden gedropt. Sommige gridcellen hebben van het begin al 6 plaatsen nodig.

Bij de thesismeeting van vandaag hebben we een aantal oplossingen voorgesteld:

1. Vergroten van gridresolutie. Dus een gridcel moet kleiner worden als de kernel radius. Dit zorgt er wel voor dat er meer gridcellen per particle moeten worden overlopen, maar het kan zorgen voor stabiliteit.
2. Een tweede of derde textuur gebruiken om extra particles in op te slaan. Het aantal texturen zou experimenteel moeten worden bepaald.
3. De textuur vergroten om zo meerdere pixels per gridcel te nemen. Men kan zelfs eventueel halve pixels nemen om zo 2 particles op te slaan.
4. Tait equation gebruiken om zo te kunnen bepalen hoe hard de vloeistof samendrukbaar mag zijn.
5. Particles die niet in de gridcel zitten, zijn onzichtbaar voor alle particles (inclusief hunzelf). Een kleine verbetering zou kunnen zijn dat bij de densiteitberekening elke particle zichzelf automatisch als apart geval beschouwt, waardoor dat de onzichtbare particles op zijn minst een densiteit van zichzelf krijgen als ze geen buren hebben.
   

Dus volgende week ga ik er voor proberen te zorgen om alles stabieler te krijgen. Ook ga ik Harada een mail sturen om te vragen hoe het precies zit met het droppen van particles die niet in gridcellen passen.

Meer info van zodra ik resultaten heb!

Pic & Flip

Na onze presentatie zijn we begonnen met onderzoek naar de pic en flip methode. Na een kleine poging tot implementatie ontdekten we toch wat moeilijkheden die eerst aangepakt moesten worden. Concreet ontbraken we een manier om parallel cellwaarden up te daten zonder synchronisatie conflicten en hadden we nood aan een efficiënte interpolatiemethode.

In de pic en flip methode moeten we eigenschappen van onze deeltjes overzetten naar een gridstructuur. Cellen uit dit grid krijgen een bijdrage per deeltje, waardoor ze simultaan geupdate moeten worden (bij niet CUDA implementaties wordt hier texture splatting voor gebruikt, maar dit is in CUDA niet mogelijk). Gelukkig voorziet CUDA atomische operaties, maar jammer genoeg alleen met integers (wij gebruiken floats). In een eerste poging gebruiken we een conversion factor om floats om te zetten in integers en weer terug. De resultaten waren bruikbaar, we kunnen een factor bepalen die voor een aantal optellingen een zekere nauwkeurigheid garandeert.

Dit leek ons eerder een omslachtige methode (die niet altijd even nauwkeurig gaat werken) dus probeerden we met een soort locking systeem de synchronisatie te garanderen. Al snel bleek de atomische "compare and store" (CAS) functie het nuttigst en we vonden een werkende implementatie van atomicFloatAdd op de NVIDIA forums. In onze testsituatie bleek deze manier van werken wel 100 keer trager, maar er is natuurlijk geen verlies aan nauwkeurigheid. Er zal moeten getest worden met beide mogelijkheden om uit te maken welke de beste resultaten zal opleveren.

Als tweede puntje werd de interpolatiemethode bekeken. We kunnen op 3 manieren interpoleren:

• lineair
• met behulp van de kernel functions
• 3D Texture interpolatie van CUDA

We implementeerde reeds de eerste twee methodes en kwamen tot de conclusie dat er geen verschil is tussen beide, althans niet in nauwkeurigheid. De nauwkeurigheid hangt wel sterk af van het aantal gridcellen en de smoothing length van de deeltjes. Bij een dubbel zo kleine smoothing length heb je ongeveer dubbel zoveel gridcellen nodig, zoals je kan zien in deze grafiek:

De 3DTexture interpolatie hebben we nog niet uitgeprobeerd. Deze gaat minder nauwkeurige resultaten geven (er wordt niet exact geïnterpoleerd), maar wel sneller werken.

Met deze resultaten kunnen we verder met de implementatie van pic & flip. Deze taak zal ik op mij nemen, terwijl Daniel een uniform grid implementeert. Dit grid kent dan wel de beperking dat er af en toe deeltjes niet in gaan passen, waardoor de simulatie minder goed gaat zijn. We nemen een harde deadline: volgende week vrijdag moeten de implementaties af zijn zodat we kunnen vergelijken welke methode de beste resultaten geeft.

Als extraatje heb ik nog een linkje naar gratis Microsoft software voor studenten, bijvoorbeeld Visual Studio 2008 professional. (wie weet kan iemand dit goed gebruiken, wij blijven verder werken met kate en linux)

Eerste thesispresentatie

Vorige week zijn we vooral bezig geweest met onze tussentijdsepresentatie. Deze is hier terug te vinden. Buiten enkele nieuwe ideeën die we kunnen gebruiken in de thesis zelf kregen we nog wat tips om de presentatie te verbeteren:

• een betere aansluiting tussen de verschillende delen voorzien
• tijdens het presenteren niet naar de slides kijken
• onze grafiekjes beter verzorgen en hiermee proberen aan te tonen wat de schaalbaarheid nu net is (dit is belangrijk voor het interactieve aspect)

Vermits de grafiek met resultaten niet zo duidelijk was, geven we hier een tabel met de berekende waarden. De timestep bepaald hoe lang we moeten simuleren in een seconde en hebben we niet constant genomen, wat voor een bespreking van de schaalbaarheid eigenlijk wel zou moeten. (de reden dat we de timestep lieten variëren is natuurlijk de performantie: hoe groter we de timestep krijgen, hoe meer frames per seconde we zullen halen)


We hebben ook nog de bug uit onze CPU implementatie gehaald en het programma op de grafische kaart werkend gekregen (zonder optimisaties). Tenslotte hebben we ook een mogelijke visualisatie techniek uitgeprobeerd, namelijk blobs. Hierbij laten we de deeltjes samensmelten, zodat ze samen een volume vormen. Deze techniek is niet real time maar geeft wel meer de illusie van water weer:

PIC and FLIP

We hebben deze week een drukke week gehad, met onder andere de presentaties van CG-thesissen. Interessant was om feedback te krijgen na de presentatie zelf: De gegeven suggesties geven ons nieuw inzicht over de toekomst van onze thesis.

Zo was er een voorstel om de PIC and FLIP methode te overwegen bij het berekenen van de buren van alle particles. We hebben vandaag afgesproken met Toon, die ons een paper heeft voorgesteld waarbij dat PIC & FLIP worden gebruikt om de neighbour search algoritme te omzeilen. (Animating Sand as Fluid, Zhu & Bridson (2005))

De Particle in Cell methode (PIC) maakt gebruik van een grid om benaderde waarden bij te houden. Ik heb hier een tekening van gemaakt om te illustreren:

Figuur 1: Particle in Cell methode

Wat men normaal gezien verwacht in een SPH implementatie is dat men op zoek gaat naar buren van alle particles, om zo de krachten tussen de particles onderling te bepalen. Met de PIC methode gaat men anders redeneren. Men maakt gebruik van het feit dat het in SPH mogelijk is om op elk puntje van uw volume te gaan evalueren aan de hand van de kernel functies, zelfs al valt dat puntje dus niet op een particle.

Om hier gebruik van te maken stelt men een uniforme grid op, en gaat men in die grid op welbepaalde punten de druk, dichtheid etc... analyseren. Elke particle heeft invloed op alle gridcellen die binnen zijn kernel straal liggen. Zo kan men dus bijvoorbeeld de dichtheid bepalen van elke gridcel door alle particles exact één keer te overlopen. In de volgende stap gaat men dan de gevonden waarden terug mappen op de particles, door opnieuw gebruik te maken van de kernels.

In figuur 1 is geillustreerd hoe dat een particle invloed heeft op naburige gridcellen. De paarse cirkel stelt een kernel straal voor, en de kruisjes zijn gridcellen. De groene gridcellen zijn degene die binnen de straal vallen, en die dus zullen worden geincrementeerd door de desbetreffende particle. Het geheel kan geparalleliseerd worden in cuda, aangezien men voor alle particles tegelijk de gridcellen kan gaan uitrekenen.

Het stukje FLIP moet er dan weer voor zorgen dat kleine details zo goed mogelijk bewaard blijven tijdens de PIC methode. Het probleem met de PIC methode is dat de eigenschappen van de particles heel hard uitgemiddeld worden. De details zijn voorlopig nog niet duidelijk, maar met FLIP gaat men er voor zorgen dat men de berekeningen incrementeel toevoegt aan de particles aan de hand van de vorige berekeningen, in plaats van het volledig absoluut te bepalen. Deze methode schijnt veel beter te werken om kleinschalige effecten te bereiken.

Alleszins, we gaan het nu onderzoeken in hoeverre dat het de moeite is om deze methode te gebruiken in onze SPH implementatie.

Update

Ook een kleine updateje van mijn kant.
Zelf ben ik ook bezig aan een SPH implementatie in CUDA, maar vanuit een bestaande applicatie. In de SDK van CUDA bevond zich namelijk een voorbeeld van een particle systeem met onderlinge interactie tussen de particles. Geen SPH, maar wel een mass-spring based systeem. Het framework ziet er redelijk goed uit, met slides om waarden in te stellen en nogwat extra functionaliteiten. Er wordt een optimisatiestructuur (uniform grid) opgebouwd en in de handleiding wordt er zelfs gezegd dat deze applicatie uitgebreid kan worden om SPH te simuleren. Door het toevoegen van SPH leer ik CUDA gebruiken en zie ik ook een uitwerking van een ander type particle system.

Eerste ervaring met cuda

De voorbije week zijn we bezig geweest met het leren werken met cuda. Ik persoonlijk heb gepoogd om een naïeve implementatie van sph te maken, gebruikmakend van de parallelisme van cuda.

Aangezien men met cuda rechtstreeks naar opengl vertex buffers kan schrijven, worden posities en snelheden van particles bijgehouden als vbo op de gpu. Vervolgens worden de particles op het scherm getekend met behulp van pointsprites. Dit is een beetje een anders ten opzichte van de cpu implementatie, waarbij solid spheres werden gerenderd. De reden is gewoonweg omdat dit voorlopig de makkelijkste en efficiëntste manier was om particles weer te geven, zonder eventuele CPU<-> GPU overhead.

Het geheel ziet er nu ongeveer zo uit:

Point Sprites in OpenGL

Momenteel is enkel zwaartekracht geïmplementeerd, en hebben de particles onderling nog geen interacties. Om de zwaartekracht toe te passen op de particles, wordt er gebruik gemaakt van cuda kernels. (Niet te verwarren met SPH kernels). De kernels zijn feitelijk een soort van C methoden die uitgevoerd worden op de GPU. Het leuke aan kernels is dat men op moet geven hoeveel threads men tegelijk wilt gebruiken om zo een kernel uit te voeren.

Kernels worden uitgevoerd in threads, die zich bevinden in blocks:

CUDA Thread mangement

Een block bestaat uit meerdere threads, die gebruik kunnen maken van maximum 3 coördinaten om de threadID te bepalen. De blocks zelf worden geïdentificeerd door maximum 2 coordinaten om de blockID te bepalen in de grid. Een blok kan maximum 512 threads hebben, maar meestal minder (afhankelijk van het aantal shared en private memory dat gebruikt wordt door de threads). Op het aantal blocks zit geen limiet.

Men moet er wel rekening mee houden hoe dat het geheel wordt gescheduled door de GPU. Op hoog niveau kan elk microprocessor (afhankelijk van de gebruikte hardware) een X aantal blocks tegelijk berekenen. De blocks hebben geen vaste sequentiele volgorde, en men mag er niet van uit gaan dat sommige blocks eerder als andere blocks worden berekend. Blocks die niet actief zijn worden in een queue gestoken.

In mijn naïeve implementatie hou ik voorlopig nog geen rekening met optimale block en grid size. Ook hou ik geen rekening met de optimale particle count (die afhankelijk is van de gpu), deze wordt voorlopig op 1000 gezet vóór het compileren.

Om de neighbours bij te houden van de partices gebruik ik een lineaire array die ik uitschrijf naar een texture op de GPU. De reden hiervoor is omdat cuda enkel naar arrays kan schrijven (en niet rechtstreeks naar textures). Het uitschrijven van een lineaire array naar een texture geeft volgens geeft relatief weinig overhead, en heeft veel voordelen bij het uitlezen: Textures worden namelijk gecached, waardoor dat het heel voordelig is om te gebruiken voor constante geheugen met hoge lokaliteit.

De lineaire array wordt opnieuw parallel beschreven door alle particles onderling. Elke thread heeft een manueel afgebakend gebied waar het in kan schrijven. Uiteindelijk zullen de waarden in de texture worden gebruikt om de krachten (dus feitelijk de velocity vbo) te wijzigen voor elke tijdstap.

Toegegeven, het is een heel naïeve methode (voor 1000 particles is een gigantische texture van 1000x1000 vereist met veel redundante data), maar het is bedoeld om vertrouwd te geraken met cuda. Optimalisaties zullen nog volgen.