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