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.

Nieuwe waarden

Voor de CPU versie van SPH hebben we nog wat zitten optimaliseren om een goede vergelijking te kunnen bekomen met de uiteindelijke versie op de GPU. De optimalisaties die we hebben doorgevoerd waren voornamelijk het slechts 1 keer berekenen van de buren, en de krachten bereken in slechts 1 functie. Zo kunnen we onder andere de kernel berekeningen herbruiken bij het berekenen van de surface tension. Voor de rest hebben we het geheel gecompiled met de -ffast-math en -O1 optie, om zo nog de FPS wat te boosten.

Wel heerst er nog een bug in onze code, dat vooral bij hoge particle-count zichtbaar is. Momenteel krijgen we een onstabiele gedrag bij meer als 1000 particles. Zelfs voor 2197 particles moeten we de timestep verlagen naar 0.001. Voor de rest van de tabel gebruiken we als timestep 0.0014.

We hebben ook de werkelijke fps berekend aan de hand van de timestep en de opengl frames per second. We nemen 30 fps als referentiewaarde voor een interactieve applicatie. Als onze timestep bijvoorbeeld 0.001 milliseconden is hebben we 1000 berekeningen nodig om 1 seconde te kunnen simuleren. Uit deze 1000 berekeningen zullen we dan 30 frames op het scherm tonen. Als die 1000 berekeningen niet binnen 1 seconden lukken is onze applicatie niet meer interactief. In onderstaande tabel zien we dat onze applicatie interactief is tot 729 particles (= 30 fps). Wanneer we zouden spreken over een simulatie aan 15 fps, zouden we maar om de 2 seconden een volledige seconde gesimuleerd hebben. Deze FPS notatie is vooral handig om mee te vergelijken: als we dus in het vervolg over FPS spreken zal voor de duidelijkheid dit systeem gevolgd worden.

Particles	OpenGL Frames Per Second	100 frames met visualisatie (s)	100 frames zonder visualisatie (s)	Werkelijke Frames Per Second
125	4500-5500	0.14	0.05	231
343	1100-1300	0.36	0.11	50.4
729	450	0.75	0.23	30
1000	330	1.07	0.29	13.8
2197	140	2.11	0.63	4.2

Afspraak Toon Lenaerts

Onze excuses voor de verwarring omwille van de vorige post. Blijkbaar kunnen er twee verschillende dingen worden verstaan onder het aantal Frames Per Second. De hoeveelheden die wij hadden gepost sloegen op het aantal "OpenGL Frames Per Second" en niet op de "werkelijke" Frames Per Second die men bekomt indien men de tijdstap per frame meerekent in de berekeningen. We zullen volgende week nog onze code lichtjes optimaliseren en dan ook meteen deze werkelijke FPS berekenen.

Afspraak Toon Lenaerts 10/10:

• Deze week beginnen met CUDA goed door te nemen en al experimenteren met de code.
• Volgende week beginnen aan een basis implementatie van SPH op de GPU met CUDA.
• Zoeken naar papers waarbij men een geschikte datastructuur gebruikt om optimaal en volledig parallel SPH te berekenen.
• Voor de komende weken: Vergelijking tussen CPU en GPU statistieken van naïeve implementatie SPH.
  

SPH met grid

Vandaag dan eens een grid structuur geimplementeerd om SPH wat te versnellen. De structuur wordt opgebouwd met kubussen waarvan de zijden gelijk zijn aan de support range van de kernels. Zo moet men voor elke particle de inhoud van slechts 27 kubussen in rekening brengen ipv alles binnen de boundary.

SPH zonder grid



SPH met grid

Eerste Video

Hier een eerste video om te demonstreren hoe ver we momenteel staan. We hebben een simpele en naive SPH geimplementeerd op de CPU. Het filmpje is een versnelde capture van hoe het geheel er momenteel uit ziet.

Er zijn 4 filmpjes, waarbij dat er wordt gevarieerd tussen het wel of niet reflecteren van de snelheid van de particles bij het botsen tegen een wand, en de viscositeit.

Het is duidelijk te zien dat een hogere viscositeit leidt tot een slijmerige vloeistof dat moeilijk kan rondbewegen, terwijl die met lage viscositeit veel splash veroorzaakt.

Het reflecteren van de snelheid indien er gebotst wordt tegen de rand is een naive oplossing. In "M.Becker, M.Teshner: Weakly compressible SPH for free surface flows (2007)" wordt een betere methode voorgesteld waarbij een botsing met een wand als een externe kracht wordt beschouwd.



Misschien nog belangrijk om te vermelden dat de filmpjes niet real-time zijn. Ik zal morgen de beeldjes per seconde voor en na optimalisatie berekenen. Optimalisatie zal gebeuren door betere geheugenbeheer en een eenvoudige grid systeem.

Planning komende week

We hebben net een afspraak gehad met Toon, wat ons toch weer meer inzicht bijbracht over SPH. We kunnen deze week niet verder werken, maar hebben wel al een planning voor de komende weken:

volgende week:

• SPH implementatie verder uitwerken (tegen donderdag)
• al wat optimalisaties uitvoeren
• wat testen met de parameters en resultaten bekomen
• al wat beginnen uitzoeken hoe cuda werkt
• nog wat literatuur over SPH raadplegen

De weken erna doen we verder aan onze literatuurstudie, we onderzoeken waar we met de thesis heen willen, leren CUDA om op de GPU te programmeren en maken de presentatie. Hopelijk kunnen we al een correcte animatie tonen tegen eind oktober. We kunnen ook eens nakijken wat de alternatieven zijn van SPH (bijvoorbeeld grid based technieken) en verduidelijken waarom wij SPH gebruiken.

In november beginnen we ieder apart aan de implementatie met CUDA. Als dit vlot verloopt kunnen we eens kijken naar mogelijke uitbreidingen. Een combinatie van SPH met height fields lijkt een aantrekkelijke optie.

SPH zonder zwaartekracht

Na lang zwoegen met constanten heb ik snel door dat het niet zo evident is om een stabiel systeem te bekomen. Zoals staat beschreven in de SPH Survival Kit:

Although one would like the SPH parameters to have direct correspondance to the real world, this is in general not the case.

Wat ik vooral heb geprobeerd is om het gemiddelde pressure en density op redelijke niveau's te houden. Ook wanneer men de (momenteel 3) krachten met elkaar vergelijkt, moet men ervoor zorgen dat ze ongeveer van dezelfde grootte orde zijn. Als uw surface tension te hoog is zullen de particles te snel ineen krimpen, en voor een hoop onstabiliteit zorgen.

Ook de volume speelt een grote rol, aangezien je massa per deeltje afhankelijk is van de hoeveelheid vloeistof dat je voorstelt.

Er zijn veel afhankelijkheden van constanten in de gebruikte formules. Daarom dat het niet zo evident is om ze optimaal te bepalen.

De 3 krachten die we momenteel voorstellen zijn druk, viscositeit en oppervlaktespanning. Hopelijk kunnen we morgen nog een beetje tijd vinden om zwaartekracht en collision forces er nog bij te krijgen.

begintoestand

evenwichtstoestand

Eerste screenshot

We hebben om te beginnen een miniprogrammatje gebouwd om particles te tonen op het scherm. Veel interactie is er nog niet: de particles vallen gewoon omlaag.

Next up: kernel testen..

Snel een samenvatting

Aangezien we redelijk laat met een blog zijn begonnen, hier een samenvatting van wat we tot nu toe hebben gedaan.

Afspraak met Toon Lenaerts op 4 juli 2008:

• Algemene bespreking van verloop thesis.
• Aanbevolen boeken over c++
• Afweging tussen Cg (C for graphics) en de nieuwe Nvidia CUDA.
• Algemeen bespreken van eerder werk over interactieve vloeistofsimulatie.
1. Algemene SPH: Dit is een particle model van vloeistofsimulatie.
2. SPH op de GPU.
3. Heightmap model van vloeistofsimulatie.
• Één van de grotere streef doelen is om heightmaps en SPH op de GPU te combineren om zo een nog efficiëntere simulatie te bekomen.

Afspraak met Toon Lenaerts op 20 september 2008:

Op deze afspraak hebben we in detail het SPH model besproken. Na een drukke eerste week zijn we klaar om in actie te schieten met ons eerste doel: Een naïeve implementatie van SPH op de CPU in C++. Het is een experiment beter inzicht te krijgen in de werking van SPH, en het effect van verschillende constanten die vaak "met de natte vinger" worden bepaald.