Computer Simulations in Physics. Numerical modelling of natural phenomena. All written by me in C/C++/HTML5. Most of them are with the code available from me,
1) Fill the array with values 0/1 with probability 0.5.
2)
Loop over entire array, calculate the sum from the 8 neighbors (plus
the given node) and change the value according to the rule:
sum = 0 1 2 3 4 5 6 7 8 9 new value = 0 0 0 0 1 0 1 1 1
3)
Apply the ping-pong rule, that is, for two buffers in which we keep the
data A and B and creating new values we write into the second one so as
not to overwrite those from which we are to draw data in a given
iteration. Then we replace A with B and so on over and over again.
1) Wypełnij tablicę wartościami 0/1 z prawdopodobieństwem 0.5
2) Przejdź całą tablicę, oblicz sumę z 8 sąsiadów (plus dany węzeł) i zmień wartość wg reguły:
suma = 0 1 2 3 4 5 6 7 8 9
nowa wartość = 0 0 0 0 1 0 1 1 1 1
3) Zastosuj regułę ping-pong, to znaczy dla dwóch buforów, w których trzymamy dane A i B i tworząc nowe wartości z jednego z nich odczytujemy dane, a wpisujemy do drugiego po to, żeby nie nadpisywać tych, z których mamy czerpać dane w danej iteracji. Potem zamieniamy A z B i tak w kółko.
Przedstawiam kilka symulacji, które mogą pomóc zrozumieć podstawowe
zjawiska z fizyki fal. To, co je wyróżnia (tak mi się przynajmniej
wydaje ;) to, że są to symulacje trójwymiarowe z wizualizacją 4D (x,y,z i
kolor dla amplitudy. Zjawiska o których mówię, to:
Fala w pudełku 3D z pojedynczego źródła
Dyfrakcja fali na pojedynczej szczelinie
Dyfrakcja na podwójnej szczelinie + interferencja = eksperyment Younga (powstają prążki i tu je dobrze widać)
Ten
projekt to wynik mojej pracy z shaderami obliczeniowymi. Jest to
trójwymiarowa wersja różnic skonczonych, gdzie rozwiązywane jest
równanie falowe. W sumie to nic trudnego, sam model jest dobrze opisany w
mojej książce o symulacjach (Helion 2021). Najwięcej problemów sprawiła
tu wizualizacja, bo nie jest łatwo wizualizować dane wolumetryczne.
Porzuciłem jednak pomysł raymarchingu i wykonałem po prostu wizualizację
na cząsteczkach z alpha-blendingiem. Nie jest idealnie, ale coś
wreszcie widać.
Jeśli ktoś woli - można obejrzeć to samo w formie kompilacji z muzyką.
Zainteresowanych techniczną stroną programowania shaderów zapraszam do wysłuchania wykładu na ich temat - świeży, z tego tygodnia.
This sound small but takes me some time to find out.
My colleague asked me how to access (read or write) GPU buffer on the CPU side.
This is simple, you have to use mapping from OpenGL, so we basically do:
float *data = A.map<float>(GL_READ_ONLY);
// now you can read from data[] table on CPU side
This was trivial, what was not trivial (and I found it out from some examples as my program crashed all the time) - you have to.. unmap this as well after you use it, so you basically do:
This tutorial is the first part of the mini compute shader serie. See part 2 on continuation with more specific fluid flow 3D simulation in compute shaders.
See lecture in which I describe the tutorial in details.
I was asked by my friend to prepare short lecture for students about doing LBM (fluid flow simulation) in compute shaders. I've done it in the past using OpenGL low level code and I must say I didn't like it. It worked but adding anything was a headache (a lot of binding, buffers, and other machine state related unpleasant stuff). Thus, although it was quite cool to see large simulations run in realtime with proper visualization, I didn't continue this.
Now, as we have 2022 and OpenGL has progressed a bit as well as myself learned to use third-party libraries I found another way. This is Open Frameworks (OF), really cool and nice C++ library that mainain much of the work that is unrelated to the core simulation algorithm but must be done anyway (graphics, window, shader loading etc.). Although it is now easier to work with (you will see, compute shaders are working like a charm) I found it rather difficult to start with OF from scratch, even though I learned compute shaders before and I knew 80% of the stuff here. This is mostly because OF documentation is weak and very often point to OF forum. This is why I decided to write this tutorial - I write it hoping that it will be useful to beginners who want to make
first steps into compute shaders.
The aim of this tutorial is to provide a quick introduction to compute shaders in Open Frameworks. I will not introduce OF itself (see forum and documentation for that - it covers it nicely) nor I will help you with initial configuration (again forum + documentation) but you will get actually everything to set up your first compute shader. In the background you will learn how to simulate Gray Scott reaction-diffusion system as I use it to ilustrate the thing.
Gray Scott Reaction-Diffusion using Compute Shader with Open Frameworks
Requirements
To start with you need to have at least OpenGL 4.3 installed on your machine. I tested all in Windows and Linux and all examples work here. Please download the latest OpenFrameworks from the openframeworks.cc website for your system and make sure you can compile the gl/computeShaderTextureExample project from the library. If you can't run it - that means your system is not compute shader ready (too low opengl, drivers problem, etc.). Here I used of_v0.9.8_vs_release for this tutorial.
What is compute shader?
In simple words - this is a program that runs entirely on GPU. GPU is a vector processor means it can run thousands of same processes at the same time. We will use that and run some simple procedures in the NX x NY image. This may be i.e. simple simulation.
Fig: So, this is GPU - another procesor, separated from main CPU.
Starting point
We start using three basic and empty OF files: main.cpp, ofApp.cpp and ofApp.h. Note I even removed some ususal OF standard methods to make it clear. Our ofApp class will consist of three crucial functions - setup, update and draw.
---------------------------------- ofApp.h #pragma once #include "ofMain.h"
class ofApp : public ofBaseApp { public: void setup(); void update(); void draw(); }; ---------------------------------- ofApp.h (end)
int main( ) { ofGLWindowSettings settings; //settings.setSize(W,H);// this is needed if you work in linux settings.width = W; settings.height = H; ofCreateWindow(settings); ofRunApp(new ofApp()); } ---------------------------------- main.cpp (end)
What do we want to write?
Let us write some simple computer simulation to ilustrate how shaders work. We will implement the Gray Scott model of diffusion-reaction system that may be implemented in discrete world easily. See some web resources to study the subject in more details: 1. https://en.wikipedia.org/wiki/Reaction%E2%80%93diffusion_system or 2. https://www.karlsims.com/rd.html 3. https://www.lanevol.org/resources/gray-scott
In brief: we will work on four rectangular fields A1,A2,B1 and B2 which are all of the (NX,NY) size. We need four, because we will use the ping-pong strategy where fields "1" will be source for changes in fields "2" in even and vice-versa in odd time steps. The algorithm seeks for the solution to a set of equations (by discretizing them). The equations are rather simple and you may find them elsewhere. Here I will just type them as written in the code:
where laplA and laplB are discrete versions of laplace operators, 1 and 2 interchange depends on the odd/even time step and f, k are some model parameters that you may change. It's up to you what you choose but be aware that changing too much may blow your model and simulation. Some safe values chosen by me you will find below in the code.
That is - the algorithm is simple. After you initialize your A1 and A2 by value 1.0 plus B2 by value 0.0 and B1 randomly to 0.0 and somewhere by 1.0 you simply iterate the above two equations (remember to swab 1 and 2 from step to step) and visualize A and B fields.
Prepare buffers for shader implementation
As it was said our shader will work entirely on GPU. That means one need to prepare memory to keep A1,A2,B1 and B2 there. What I do below is that I create buffers on CPU first, fill them with some initial values and then allocate (with CPU->GPU copy) the same memory on GPU using GL_SHADER_STORAGE_BUFFER structure. First, we need to declare CPU arrays and GPU buffers, which in Open Frameworks are called ofBufferObject.
---------------------------------- ofApp.h #pragma once #include "ofMain.h"
#define W 1280 #define H 720
class ofApp : public ofBaseApp { public: void setup(); void update(); void draw();
Note, that we used one dimensional array mapping of 2D array (means we recalculate idx each time we access the array).
Buffer bindings - let the shader know which buffer is which
As we are done with allocating buffers on GPU, we need to let the shader know which buffer is which.
Fig: Binding may be understood as wire to slot connections.
This is called "binding" and in simple words it is just assigning id numbers to allocated buffers. Why do we need this? In CPU world buffers have names but those are not visible on GPU side, thus, it is much easier to operate on numbers. First buffer in CPU world will be first in GPU and so on. OK, stop talking, let's bind the stuff now.
// allocate buffers on GPU side A1.allocate(W*H*sizeof(float),A1cpu,GL_STATIC_DRAW); A2.allocate(W*H*sizeof(float),A2cpu,GL_STATIC_DRAW); B1.allocate(W*H*sizeof(float),B1cpu,GL_STATIC_DRAW); B2.allocate(W*H*sizeof(float),B2cpu,GL_STATIC_DRAW);
That was easy right? I promise this is not even a bit harder, really, just stay with me.
We want to see things - visualize results as a texture
Our simulation will perform steps on A1...B2 arrays in GPU world but what we want at the end is to visualize results and put some fancy col4 results on the screen. For that we will use simply a texture which has to be allocated, binded (again!) and drawn in our program. All steps are trivial and I will just illustrate them by coloring changes in code below.
---------------------------------- ofApp.h #pragma once #include "ofMain.h"
#define W 1280 #define H 720
class ofApp : public ofBaseApp { public: void setup(); void update(); void draw();
Please note that binding differ a bit from binding of the buffers - we will now bind the texture as an image. Also allocation needs to be done such that our image will consist of colors (GL_RGBA8 for the pixel format is used here). ---------------------------------- ofApp.cpp #include "ofApp.h" #include "ofConstants.h"
// allocate buffers on GPU side A1.allocate(W*H*sizeof(float),A1cpu,GL_STATIC_DRAW); A2.allocate(W*H*sizeof(float),A2cpu,GL_STATIC_DRAW); B1.allocate(W*H*sizeof(float),B1cpu,GL_STATIC_DRAW); B2.allocate(W*H*sizeof(float),B2cpu,GL_STATIC_DRAW);
Culmination - now we actually implement the shader
To get the shader we need to, again, declare it (shader in OF is the ofShader object), load it's source and comple it on GPU side. Now, with OF this is easy. I just add the proper object in ofApp.h class definition:
---------------------------------- ofApp.h #pragma once #include "ofMain.h"
#define W 1280 #define H 720
class ofApp : public ofBaseApp { public: void setup(); void update(); void draw();
Next, I implement it's initialization. We will load the source code of the shader (setupShaderFromFile method of ofShader object) and then compile it and link to shader program (linking in OpenGL words).
// allocate buffers on GPU side A1.allocate(W*H*sizeof(float),A1cpu,GL_STATIC_DRAW); A2.allocate(W*H*sizeof(float),A2cpu,GL_STATIC_DRAW); B1.allocate(W*H*sizeof(float),B1cpu,GL_STATIC_DRAW); B2.allocate(W*H*sizeof(float),B2cpu,GL_STATIC_DRAW);
Now, we will implement the update function. The aim of this function is to run shader step by step and swap the buffer binding 1->2 and 2->1 depends on the time step. The technique is simple, we will have some semaphor kind of variable "c" which will be 0 or 1 and do proper binding depend on its value. I use c = 1-c each time step to change 0 to 1 and opposite.
Fig: We use standard ping-pong scheme in which buffers A1/B1 swap with A2/B2 to become source or destination from step to step.
Then, after binding is done, we simply run the shader. This is called "dispatchCompute" (why not run?.. don't ask me). And three parameters we provide are just to divide our taks into parallel, independent threads that will run simultaneously. This is, where all the magic happens - I mean this is where parallel computing is done.
Fig: Here we used DispatchCompute using W/20 x H/20 size for single working group. Just remember, all threads in working group run in parallel and single working group cannot be larger than allowed by your drivers. More details i.e. on Khronos pages here.
Side note: I used W/20 and H/20 for horizontal and vertical number of threads in the run grid. Details on how GPU works and organize threads may be found elsewhere (Khronos, nVidia, AMD or CUDA documentation) but what you need to know is that you can't push too many threads if you'r card doesn't allow for that, so if your simulation doesn't work or work too slow try to decrease or increase this "20" hardcoded there to some other values.
// allocate buffers on GPU side A1.allocate(W*H*sizeof(float),A1cpu,GL_STATIC_DRAW); A2.allocate(W*H*sizeof(float),A2cpu,GL_STATIC_DRAW); B1.allocate(W*H*sizeof(float),B1cpu,GL_STATIC_DRAW); B2.allocate(W*H*sizeof(float),B2cpu,GL_STATIC_DRAW);
Up to now whar we were doing was the most tedious part. Now we can eat the candy and have fun with compute shader implementation. We will use GLSL languague for that. Don't worry it is very simple. Our starting point is empty shader that says which version of OpenGL we require and empty "main" function. Btw. the ".cs" name is not the only option, I've seen people use ".glsl" as well.
Now we progress by adding binding layout indications. We need to "receive bindings from the CPU side" so that we can use them later in our code and main function. In addition layout of the execution grid is specified, I use the same integer number of grid cells as used in dispatch command while running our shader (means 20x20 block size).
I use a little "per" help function to deal with periodicity. This is a function that will be called every time our index want to point into out of domain cell. So, our solution will be actually 2D like wraped on the torus (periodic boundary conditions).
Another helpful function will be colour mapping one that I base on the article from the IQ's website (https://iquilezles.org/www/articles/palettes/palettes.htm). It is simple and returns vec4 color mapped to some nice palette based on the float 0-1 input. Adjust coltab[] if you want to get another palette. Nice examples can be found on IQ's website.
To implement the main function we need to know which cell do we compute. For this, a standard variable gl_GlobalInvocationID (which is vector with x/y components) is utilized. So we read the x and y and for this cell compute everything that is needed. Note that we use buffers binded from CPU side so this way we can control our shader (this is what we do here to do ping-pong double buffered computation using binding from the cpu side).
Note I hardcode the "W" and "H" which coud be passed by uniform variables from CPU to shaders. But I didn't want to complicate the code further.
I compute the grid index idx and then set all the model constats.
As you may notice I put 5 cents to make the Gray Scott model look better and change the "k" and "f" parameters smoothly with the x index in the grid. This is because I found it quite boring if single values were used. But it is up to you to play with them. One thing I didn't try and may turn to be nice would be to change them in time (i.e. using some sinuses or even impulses).
Next I caluclate indices for all neighbours (those idx0-idx8 variables) and then compute laplacians and final Gray Scott model that consists of... two lines of code. This is funny, so much work to prepare for those two lines.
Finally visualization is done by storing rgb color in the texture that is calculated directly from current A and B data.
For those of you that are too lazy to copy from above or for those that jumped from the title to the end section - the final project consists of 4 files listed below.
int main( ) { ofGLWindowSettings settings; //settings.setSize(W,H);// this is needed if you work in linux settings.width = W; settings.height = H; ofCreateWindow(settings); ofRunApp(new ofApp()); } ---------------------------------- main.cpp (end)
---------------------------------- ofApp.h #pragma once #include "ofMain.h"
#define W 1280 #define H 720
class ofApp : public ofBaseApp { public: void setup(); void update(); void draw();
// allocate buffers on GPU side A1.allocate(W*H*sizeof(float),A1cpu,GL_STATIC_DRAW); A2.allocate(W*H*sizeof(float),A2cpu,GL_STATIC_DRAW); B1.allocate(W*H*sizeof(float),B1cpu,GL_STATIC_DRAW); B2.allocate(W*H*sizeof(float),B2cpu,GL_STATIC_DRAW);
As an result of the above code we get the time evolution of Gray Scott model on a 1280x720 grid that runs at interactive rates on moderate GPU. Please see screenshot below or jump directly tothe youtube video that was captured while I was writing the article for you.
Please leave a comment or drop an email if you found the article useful. You may also like to subscribe to my youtube channel if you want more frequent updates from me.
This tutorial is the first part of the mini compute shader serie. See part 2 on continuation with more specific fluid flow 3D simulation in compute shaders.
Jest kwiecień 2021 i właśnie we Wrocławiu spadł po raz kolejny tej zimy, prawdziwy śnieg. Przedstawiam wykład o płatkach śniegu, o ich strukturze i modelowaniu.
Płatki śniegu powstają przez wzrost kryształków lodu, które przyczepiają się do siebie i tworzą interesujące, różnorodne bryły. Bryły te w większości mają określoną symetrię, ale bogactwo form jest ogromne (zdarzają się takie z różnymi defektami lub wzrastające w różny sposób). Przykłady różnorodności płatków śniegu można znaleźć na przykład na fotografiach Wilsona Bentleya, amerykańskiego farmera żyjącego w latach 1835-1931, który sfotografował ponad 5000 płatków śniegu opisując je jako piękne małe obiekty.
Wilson Bentley (1835-1931) i jego fotografie płatków śniegu.
Dziś można już oglądać płatki śniegu pod mikroskopami różnego typu, popatrzcie na zdjęcia wykonane np. mikroskopią elektronową.
Szczególnie ciekawe jest to, że zdjęcia płatków śniegu są wybierane na obraz roku w Wikipedii (ja znalazłem takie dwa zdjęcia). Polecam szczególnie stereogram, zdjęcie, które oglądane w odpowiedni sposób (należy zbliżyć oczy do ekranu i skupić wzrok na punkcie pomiędzy ekranem, a oczami) daje złudzenie oglądania bryły 3D. Mi się to prawie udało na kilka sekund, ale nie jestem w tym dobry. Ciekawe jak Wam to pójdzie (uwaga - polecam wyświetlić na pełnym ekranie):
Żr: Wikipedia
W 1611 r. Johannes Kepler, niemiecki astronom, autor praw dotyczących ruchu planet (słynne prawa Keplera), po raz pierwszy analizuje sześciokątną symetrię płatków śniegu twierdząc, że jest związana ze sposobem w jaki układają się układy gęsto upakowane sfery (ang. sphere packings). Publikuje on też krótki esej na ten temat [2]. Faktycznie, najgęstsze upakowanie dwy-wymiarowych kół w dwóch wymiarach realizowane jest na sieci o strukturze trójkątnej (heksagonalnej). Jeśli chodzi o sam wzrost kryształów na sieci heksagonalnej to jest on związany z budową cząsteczki wody.
Oryginalny rysunek z pracy J. Keplera [2].
Jeśli przyjrzymy się prostej symulacji zjawiska upakowania kółek dwuwymiarowych o tej samej średnicy i będziemy ich dodawać więcej i więcej, to na koniec zauważymy, że układają się na sieci o strukturze sieci trójkątnej.
Symulacja: Dwuwymiarowe dyski upakowane na tyle gęsto, że tworzą sieć trójkątną.
Jednym z najprostszych sposobów modelowania wzrostu płatków śniegu jest model automatu komórkowego, który opisany jest na sieci trójkątnej w dwóch wymiarach. Sieć trójkątna powstaje przez przesunięcie nieparzystych rzędów regularnej sieci kwadratowej o połowę długości sieci.
Rys. Sieć trójkątka i jej implementacja z użyciem tablicy kwadratowej.
Każda z komórek w tym modelu jest komórką lodu lub komórką pustą - te zapisywać będziemy odpowiednio jako zera i jedynki. W kolejnych krokach algorytmu komórki, które były lodem nie zmieniają swojego stanu. Komórki, które lodem nie były zostają lodem wtedy i tylko wtedy, kiedy mają dokładnie jednego sąsiada będącego lodem. Implementacja wydaje się zatem prosta, z poprawką na implementację sieci trójkątnej. Trzeba bowiem zaprogramować sieć trójkątną. To nie jest bardzo trudne, trzeba jedynie pamiętać o tym czym w tej sieci jest sąsiedztwo punktu, bo nie ma on już czterech sąsiadów jak w sieci regularnej. Tu sąsiadów jest sześciu.
Rysunek: Model wzrostu płatka na sieci trójkątnej.
W omawianym modelu wzrost rozpoczyna się od zarodka (mogącego przybrać dowolną formę zlepku kilku komórek). Oryginalny artykuł opisujący ten model pochodzi z roku 1986 ([5] Packard NH. Lattice models for solidification and aggregation. In: Wolfram S, editor. Theory and applications of cellular automata. Singapore: World Scientific Publishing; 1986. p. 305–10.) i został spopularyzowany przez S. Wolframa w jego książce o automatach komórkowych.
Wyniki omawianego modelu (kod C++ z olcPixelGameEngine) dla różnych parametrów.
Model ten potrafi symulować wzrost abstrakcyjnych form płatków, ale to
nie potrafi wygenerować płatków z globalnym wzrostem dendrytów i gwiazd,
dlatego dziedzina symulacji i modelowania tego typu obiektów jest dość
rozwinięta. Zachęcam do własnych poszukiwań!
Przy próbie implementacji zdarzyła mi się ciekawostka. Otóż próbując zaimplementować oryginalny model z użyciem dyfuzji temperatury (próba nieudana, w artykule [5] jest chyba za mało szczegółów dotyczących parametrów) uzyskałem taki oto artefakt w modelu. Dywan Sierpińskiego na płatku śniegu!
Niespodziewany wynik w modelu płatka śniegu z temperaturą (prawdopodobnie błąd w kodzie) - dywan Sierpińskiego.
W literaturze można spotkać wiele prac z modelowania płatków śniegu. Nie będę się tu silił na wielki przegląd, ale zachęcam do spojrzenia na model C.A. Reitera z wartościami 0-1 na sieci heksagonalnej [1]. Wartość 0-1 (autorzy wskazują, że jest to duże przybliżenie i nieformalne założenie) oznacza frakcję wody w komórce. Wartości powyżej 1 oznaczają lód. Woda (wartość < 1) może przepływać pomiędzy sąsiednimi komórkami. Algorytm ma dwa kroki i jest to dwustopniowy automat komórkowy (oba kroki to dwa oddzielne automaty pracujące na wartościach rzeczywistych w komórkach). W pierwszym kroku określane są komórki receptywne (odbiorcze). Są to komórki, które są lodem, lub mają sąsiadów, którzy są lodem. Te doznają skokowego wzrostu masy (patrz rysunek z algorytmem w omawianym artykule). Reszta komórek rozpływa się (woda) z przybliżeniem dyfuzyjnym. Następnie oba efekty są sumowane. Zachęcam, aby spojrzeć do oryginalnego artykułu - efekty i zakres uzyskanych płatków śniegu jest o wiele większy niż w oryginalnym, prostym, automacie komórkowym.
I na koniec wersja video tego wykładu:
Zapraszam do zapoznania się z resztą materiałów z kursu modelarni na kanale "Koder na strychu". Materiały są ułożone w playliście:
Literatura: [1] C.A. Reiter / Chaos, Solitons and Fractals 23 (2005) 1111–1119 [2] Kepler J. The six sided snowflake (1611) translation. Oxford University Press; 1966. [3] Bentley WA, Humphreys WJ. Snow crystals. McGraw-Hill Book Company; 1931, also, New York: Dover Publications; 1962. [4] Nakaya U. Snow crystals: natural and artificial. Harvard University Press; 1954. [5] Packard NH. Lattice models for solidification and aggregation. In: Wolfram S, editor. Theory and applications of cellular automata. Singapore: World Scientific Publishing; 1986. p. 305–10. [6] Wolfram S. A new kind of science. Champaign: Wolfram Media; 2002. [7] https://en.wikipedia.org/wiki/Timeline_of_snowflake_research
