Nowoczesne Technologie Obliczeniowe - 2019

Projekt CUDA - propagacja światła za przeźroczem - model 3D

Autorzy: Szymon Baczyński & Łukasz Szeląg

Projekt opierający się o komunikację MPI i bibliotekę CUDA, by wykorzystać pełny potencjał obliczeniowy klastra DWARF wydziału Fizyki PW.

Program w skrócie wykonuje się w następujących krokach:

1. Stworzenie dwuwymiarowej tablicy o wymiarach zadanych przez użytkownika (tablica ta jest przezroczem, przez które będzie propagować się światło).
2. Wpisanie tablicy utworzonej w punkcie 1 do dwa razy większej tablicy (mniejsza tablica umieszczona w środku większej).
3. Nowo utworzona tablica (od teraz nazywana tablicą wejściową) kopiowana jest na kartę graficzną (GPU).
4. Na GPU liczona jest szybka transformata Fouriera (FFT) tablicy wejściowej.
5. Wynik wysyłany jest do procesów MPI, który następnie kopiowany jest na GPU.
6. Każdy proces oblicza dwuwymiarową tablicę odpowiedzi impulsowej (h(z) - zależna od odległości propagacji z) i wysyła ją na GPU.
7. GPU obliczają FFT odpowiedzi impulsowej, a następnie mnożą transformatę tablicy wejściowej z transformatą odpowiedzi impulsowej.
8. Wynik mnożenia zostaje poddany odwrotnej transformacie Fouriera. Po wykonaniu odwrotnej transformaty, jej wynik kopiowany jest z GPU na hosta.
9. Otrzymana tablica danych zostaje poddana operacji ROLL - zamiana miejscami ćwiartek, dla prawidłowego wyniku.
10. Ostateczny wynik zostaje zapisany do pliku w postaci amplitudy.

Pełny kod dostępny jest w repozytorium pod linkiem: prop.cu

Propagacja światła koherentnego przez dowolne przeźrocze - dyfrakcja - jest ciągle ciekawym tematem dla naukowców - pomaga zrozumieć naturę światła i jest szeroko wykorzystywana podczas symulacji układów optycznych. Istnieje kilka metod symulacji propagacji światła m.in. liczenie całki tablicy wejściowej oraz metoda splotowa.
Jako metodę symulacji propagacji światła wykorzystano metodę splotową. Bardzo dobrą literaturą ukazującą algorytm, jak i zalety i wady tej metody jest artykuł profesora Macieja Sypka [1]. Ze względu na czasochłonne obliczenia związane ze splotem dwóch funkcji, w tym przypadku tablicy wejściowej (u₁) oraz odpowiedzi impulsowej (h) - PSF (Point Spread Function), wykorzystano własności transformacji Fouriera. Szybszą i wydajniejszą metodą będzie zrobienie FFT{u₁} oraz FFT{h} i wymnożenie obu tablic ze sobą, U₂ = FFT{u₁} x FFT{h}. Wynikiem odwrotnej transormaty Fouriera tablicy U₂ , będzie tablica zespolonych wartości, u₂ , które zawierają informację o amplitudzie i fazie w danym miejscu w przestrzeni.

Program przedstawia rozkład amplitudy fali świetlnej za przezroczem w zadanych przez użytkownika odległościach (z) od tego przezrocza. Rozkłady w różnych odległościach liczone są równolegle przez procesy MPI. Jedną z wad tego algorytmu są krawędzie tablic - wprowadzające zniekształcenia i wysokoczęstotliwościowy szum. Aby pominąć szkodliwe efekty FFT na krawędziach wprowadza się modyfikację algorytmu w postaci powiększenia tablic obliczeniowych do wymiarów 2Nx2N. Tablica wejściowa, NxN , w tym przypadku przepisywana jest na środek większej tablicy. Tablicą wyjściową jest środkowa część NxN.

Stworzenie tablicy wejściowej, wpisanie jej do większej tablicy oraz wykonanie transformaty tej tablicy robione jest wyłącznie przez proces "0". Reszta procesów oblicza odpowiedź impulsową, zależną od odległości z od przezrocza, różnej dla każdego procesu. Następnie każdy proces mnoży tranformatę tablicy wejściowej z transformatą odpowiedzi impulsowej. otrzymany wynik następnie poddawany jest odwrotnej transformacie Fouriera i przesyłany z GPU na hosta.

Wykonano przykładowy kod w Pythonie wykorzystujący biblioteki numpy, scipy, PIL i pathlib. Używając tablicy 1024x1024 piksele sprawdzono amplitudę i fazę po propagacji na odległość 1 metra. Wynik wydaje się być prawidłowy porównując do obrazów z programu LightSword stworzonego specjalnie do takich obliczeń. Optymalizacji może jeszcze ulec funkcja PSF by uzyskać od razu FFT{h}. Działający algorytm jest podstawą do wyboru kierunku obliczeń przyszłego programu.
Program napisany w języku Python dostępny jest tutaj: Propagacja_Python
Dodatkowy kod w języku Python wspomagający przygotowanie danych i późniejszą obróbkę wyliczonych płaszczyzn: Dodatek

Wykonane testy skalowania:

Program został odpowiednio zoptymalizowany pod kątem pamięci, by nie było niepotrzebnych tablic zajmujących miejsce. W jednym momencie na GPU znajdują się dwie tablice: FFT tablicy wejściowej oraz funkcja H(z), czyli FFT{h(z)}. Te dwie tablice następnie są mnożone w celu symulacji rozkładu na zadanej odległości. W trakcie działania programu zauważono następujące maksymalne wartości użycia pamięci RAM na GPU:

u_in	Mnożnik M	Rozmiar Tablic tymczasowych	Max. użycie RAM-u na GPU [MiB]
1024x1024	1	1024x1024	130
1024x1024	2	2048x2048	356
1024x1024	4	4096x4096	1380
1024x1024	6	6144x6144	2980
1024x1024	8	8192x8192	5220
1024x1024	10	10240x10240	8100

Karty Tesla K80 posiadają 11439 MiB możliwej do użycia pamięci RAM. Jak widać tymczasowe użycie pamięci GPU przy mnożniku M=10, sięga ponad 8 GiB. Przy wyższych mnożnikach, np. M=12, program wykonuje część zadań a podczas próby alokowania zbyt dużej ilości danych na GPU - funkcje zwracają problem z tworzeniem/alokowaniem pamięci.

Do testowania propagacji posłużyła nam tablica przechowująca informacje o przezroczu o wymiarach 1024x1024 (wcześniej przetworzonej z BMP do txt).
Przezrocze (jako jasne punkty w tablicy) wyglądało następująco:


Podstawowe parametry do obliczenia rozkładu amplitudy w zadanej odległości za przezroczem wyglądały następująco:

Sampling = 10 [mikrometrów] (odległość od kolejnych punktów w przestrzeni x,y)
Wielkość tablicy = 1024x1024 = 10,24 x 10,24 [mm^2]
Długość fali = 633 [nm]
Wektor falowy = 2*pi/lambda

Kompilacja i uruchamianie programu:

Cały kod został napisany w pliku *.cu - co wymuszało do kompilacji zastosowania dosyć złożonej komendy linkowania i ścieżek dla kompilatora:

compile: nvcc -o prop.x prop.cu -O3 -gencode=arch=compute_35,code=sm_35 -gencode=arch=compute_37,code=sm_37 -gencode=arch=compute_60,code=sm_60 -I/usr/local/cuda/inc -L/usr/local/cuda/lib -lcufft -I/opt/openmpi-gcc721-Cuda90/3.1.1/include -Xcompiler "-pthread -fPIC" -L/opt/openmpi-gcc721-Cuda90/3.1.1/lib -lmpi

Aby korzystać z programu należy uruchomić go za pomocą danej komendy:

[dteam001@nodeXXX Projekt_Propagacja]$ mpirun -n NP prop.x TABLICA_WEJSCIOWA NX NY M

NP - liczba procesów MPI, której chcemy użyć
TABLICA_WEJSCIOWA - tablica zawierająca rozkład przezrocza
NX, NY - rozmiary tablicy wejściowej
M - mnożnik tablic użytych do obliczeń (ile razy większe mają być NX i NY)

Przykładowe tablice (1024x1024) po propagacji na zadane odległości:

Z = 500 mm	Z = 600 mm	Z = 700 mm	Z = 800 mm
Z = 900 mm	Z = 1000 mm	Z = 1500 mm	Z = 1700 mm