rozdzial2.tex

\chapter{Przeniesienie obliczeń fizycznych na procesor karty graficznej}
\label{cha:oblGPU}

Niniejszy rozdział prezentuje techniki zastosowane w celu przeniesienia
głównych obliczeń fizycznych na kartę graficzną w symulatorze \ow{Energy2D}
będącym przedmiotem tej pracy. Na początku opisane są tradycyjne podejścia do
tego problemu dla aplikacji działających w natywnym środowisku systemu
operacyjnego oraz ich odniesienie do środowiska oferowanego przez przeglądarki
internetowe. Następnie zaprezentowane są najważniejsze zagadnienia dotyczące
implementacji równoległych silników fizycznych \ow{Energy2D} działających na
procesorze karty graficznej. Informacje te mogą być szczególnie użyteczne przy
próbach podobnych optymalizacji innych aplikacji.

Dzięki przeniesieniu obliczeń fizycznych na GPU uzyskano istotny wzrost
wydajności. Dokładna analiza zysków ze zrównoleglenia symulacji jest przedstawiona
w rozdziale \ref{cha:ocena}.

\section{Typowe metody przenoszenia obliczeń na GPU, a przeglądarka internetowa}

Współczesne karty graficzne posiadają ogromną moc obliczeniową -- przy
założeniu, że obliczenia da się wykonywać w sposób równoległy, jest ona
wielokrotnie większa od mocy obliczeniowej centralnego procesora. Pomysł aby
przenieść część obliczeń ogólnego zastosowania na kartę graficzną pojawił się
wraz z dynamicznym rozwojem procesorów graficznych. Szczególnie istotnym
momentem było wprowadzenie programowalnych jednostek cieniujących
(specyfikacja \ow{DirectX 8}). Dalszy rozwój obliczeń ogólnego zastosowania na
kartach graficznych (ang. \emph{General-Purpose Computing on Graphics
Processing Units}, w skrócie GPGPU) miał miejsce wraz z wprowadzeniem
technologii, które ukryły złożoność dostępu do zasobów karty graficznej i
udostępniły interfejs wysokiego poziomu. Wiodącymi technologiami tego typu są
\ow{OpenCL} (rozwiązanie otwarte) oraz \ow{CUDA} (zamknięte rozwiązanie firmy
NVIDIA, działające wyłącznie na sprzęcie tego producenta).

Wspomniane powyżej technologie dotyczą oczywiście aplikacji pisanych w
natywnym środowisku systemu operacyjnego. Aby programować jednostki cieniujące
wystarczy podstawowy dostęp do standardowego interfejsu \ow{OpenGL} bądź
\ow{Direct3D}. Implementacje tych interfejsów można znaleźć w prawie każdym
współczesnym języku programowania. Interfejsy wyższego poziomu (\ow{OpenCL},
\ow{CUDA}) również posiadają implementacje w różnych językach programowania,
choć najczęstszym środowiskiem ich działania są aplikacje napisane w C bądź
C++.

Poniżej krótko przedstawione są sposoby przeprowadzania obliczeń ogólnego
zastosowania na kartach graficznych oraz ich związek ze specyficznym
środowiskiem przeglądarki internetowej.

\subsection{Niskopoziomowe programowanie jednostek cieniujących}
\label{subsec:niskProgJedn}

Jest to najstarsze podejście do przeprowadzania obliczeń ogólnego zastosowania
na procesorze karty graficznej. Programista w swoisty sposób ,,oszukuje'' kartę
graficzną, przeprowadzając renderowanie prostej geometrii wyłącznie w celu
uruchomienia własnych programów jednostek cieniujących, które wykonują
obliczenia często nie mające nic wspólnego z generowaniem obrazu.

W przypadku technologii \ow{WebGL}, która posłużyła do zaimplementowania
silników fizycznych \en, diagram potoku renderowania prezentuje rysunek
\ref{fig:WebGLPipeline}. Zielonym kolorem zostały oznaczone procesy, które są
programowalne. Są to dwa rodzaje programów jednostek cieniujących --
wierzchołków oraz fragmentów. Nie ma możliwości bezpośredniego sterowania
pozostałymi procesami. Odbywają się one automatycznie, ewentualnie możliwa jest
konfiguracja pewnych parametrów. Dlatego też cały algorytm przeznaczony do
zrównoleglenia musi zostać zapisany wyłącznie przy użyciu programów wierzchołków
oraz fragmentów. Ponadto dane powinny się znajdować w pamięci karty graficznej,
najczęściej w postaci dwuwymiarowych tekstur.

\begin{figure}[!h]
\centering
\includegraphics[width=\textwidth]{img/WebGLPipeline}
\caption{Diagram potoku renderowania \ow{WebGL} / \ow{OpenGL ES 2.0}}
\label{fig:WebGLPipeline}
\end{figure}

Takie podejście jednak jest wymagające i obarczone pewnymi problemami. Łatwo
popełnić błędy, wymagana jest też przynajmniej podstawowa wiedza o programowaniu
grafiki trójwymiarowej. Najważniejsze koncepcje niskopoziomowego programowania
jednostek cieniujących zostały doskonale przedstawione przez Marka Harrisa w
jednym z jego artykułów do popularnej serii \emph{GPU Gems} \cite{GPUConcepts}.

\subsection{Technologie wyższego poziomu}
\label{subsec:techWyzPoz}

Technologie, które przyczyniły się do gwałtownego wzrostu popularności obliczeń
na kartach graficznych to szczególnie \ow{OpenCL} oraz \ow{CUDA}.
Udostępniają  one znacznie wyższy poziom abstrakcji -- programista jest
zwolniony z obowiązku przyswojenia sobie niskopoziomowych mechanizmów rządzących
działaniem kart graficznych (choć ta wiedza pozwala tworzyć aplikacje
efektywniejsze). Udostępniony jest specjalny interfejs oraz składnia, co
znacznie ułatwia pracę, szczególnie programistom bez doświadczenia w pracy z
programowaniem grafiki trójwymiarowej. Najważniejsze dokumenty opisujące te
technologie to OpenCL Specification (\cite{OpenCLSpec}) oraz 
CUDA Reference Manual (\cite{CUDARef}).

\subsection{Możliwości w środowisku przeglądarki internetowej}
\label{subsec:srodPrzegInt}

Należy uściślić, że ,,środowisko przeglądarki internetowej'' to zbiór
możliwości nowoczesnych, wiodących przeglądarek internetowych,
rozpowszechniony na tyle, aby był dostępny dla większości użytkowników
internetu. Równocześnie możliwości te nie mogą wymagać instalowania żadnych
dodatków i rozszerzeń, gdyż znacznie ogranicza to ich osiągalność. Jest to
szczególnie istotne dla aplikacji edukacyjnych, które wymagają łatwego i
powszechnego dostępu.

Przeglądarka internetowa z założenia oferuje środowisko bardzo ograniczone,
przede wszystkim ze względów bezpieczeństwa. Jednak niedawno wraz z
nadejściem standardu HTML5, został dodany podstawowy dostęp do zasobów karty
graficznej. Realizuje go technologia WebGL, która jest implementacją standardu
\ow{OpenGL ES 2.0}. Tym samym pojawiła się możliwość programowania jednostek
cieniujących kart graficznych, a więc i przeprowadzania na nich obliczeń
ogólnego zastosowania (por. \ref{subsec:niskProgJedn}).

Technologie programowania kart graficznych wyższego poziomu (por.
\ref{subsec:techWyzPoz}) nie są (jeszcze) dostępne w przeglądarce
internetowej. Aktualnie trwają intensywne prace nad implementacją standardu
OpenCL o roboczej nazwie WebCL. Jednak ta technologia jest jeszcze na bardzo
wczesnym etapie rozwoju. Więcej informacji można znaleźć na stronie
internetowej: \url{http://www.khronos.org/webcl/}.

Dlatego też jedyną metodą na przeprowadzenie obliczeń ogólnego zastosowania na
procesorze karty graficznej w środowisku przeglądarki internetowej jest sposób
zaprezentowany w sekcji \ref{subsec:niskProgJedn}. Właśnie taka,
niskopoziomowa metoda została zastosowana dla silników fizycznych symulatora
\ow{Energy2D}.

\section{Implementacja silników fizycznych Energy2D przy użyciu WebGL}
\label{sec:implSilFizWebGL}

Symulator \ow{Energy2D} składa się z dwóch kluczowych silników fizycznych -
przewodnictwa cieplnego oraz dynamiki płynów. Są one dokładniej przybliżone w
sekcji \ref{sec:silnikiFizyczne}. Oba te silniki są niezwykle wymagające
obliczeniowo. Dlatego też ich optymalizacja była zadaniem kluczowym, aby
stworzyć aplikację wartościową edukacyjnie. Zbyt wolny przebieg symulacji może
skutecznie zniechęcić większość potencjalnych użytkowników. Jednym z
podstawowych wymagań było, aby wirtualne laboratorium było aplikacją w pełni
interaktywną czyli również działającą jak najpłynniej.

Poniżej przedstawione są najważniejsze zagadnienia związanie z przeniesieniem
obliczeń fizycznych \ow{Energy2D} na kartę graficzną.

\subsection{Analiza algorytmów silników fizycznych pod kątem przetwarzania
równoległego}  

Wszystkie algorytmy rozwiązujące równania fizyczne w symulatorze \ow{Energy2D}
zostały poddane analizie pod kątem możliwości równoległego przetwarzania.
Zostały wyróżnione następujące kryteria, które algorytm musi spełniać:

\begin{itemize}

\item Możliwość reprezentowana danych w pamięci karty graficznej.

\item Niezależność wyniku algorytmu od sekwencji wykonywania obliczeń.

\end{itemize}

Jeżeli kryteria te są spełnione, powinna istnieć możliwość zaimplementowania
algorytmu w sposób równoległy. Oczywiście zmienia się strona techniczna, użyte
rozwiązania, język programowania oraz techniki, jednak sama koncepcja
algorytmu i główne kroki powinny pozostać bez większych zmian. Problem pojawia
się wtedy, gdy któryś z algorytmów nie spełnia jednego z tych kryteriów. W
takim przypadku konieczne są gruntowne modyfikacje lub rezygnacja z
implementacji takiego algorytmu.

Pierwszy warunek spełniają wszystkie algorytmy, jako że obliczenia są wykonywane
na prostokątnych siatkach symulacyjnych. Można je reprezentować w pamięci karty
graficznej przy użyciu dwuwymiarowych tekstur. Temat organizacji danych w
pamięci GPU oraz związane z tym problemy porusza sekcja \ref{sec:orgDanychWGPU}.

Drugi warunek nie został spełniony przez wszystkie zastosowane algorytmy.
Problematyczne okazały się metody rozwiązywania układów równań liniowych metodą
relaksacji \ow{Gaussa-Seidela} (por. \cite{GaussSeidel}). Podczas jednego
przebiegu przez wszystkie komórki macierzy, nowa wartość komórki $(i, j)$ jest
zależna od wartości komórek sąsiednich. Następnie komórka macierzy jest
niezwłocznie aktualizowana. Powoduje to, iż przy sekwencyjnym przetwarzaniu
wszystkich komórek, nowa wartość dla każdej komórki jest zależna od wartości
obliczonych zarówno w poprzednim kroku relaksacji jak i kroku aktualnym. W
sposób uproszczony schemat takiej relaksacji prezentuje algorytm
\ref{alg:GaussCPU}.

\begin{algorithm}[!h]
  \caption{Relaksacja metodą Gaussa-Seidela na CPU}
  \label{alg:GaussCPU}
\begin{algorithmic}
\For {$0 \to relaxation\_steps$}
  \ForAll {$ \textrm{grid cells} $}
    \State $new\_val\gets f(x_{i,j}, x_{i+1,j}, x_{i-1,j}, x_{i,j+1}, x_{i,j-1})$
    \State $x_{i,j}\gets new\_val$
  \EndFor
\EndFor
\end{algorithmic}
\end{algorithm}

Niestety, przy obliczeniach równoległych nie można polegać na takiej zależności.
Każda komórka zostanie zaktualizowana na podstawie wartości komórek sąsiednich
wyłącznie z poprzedniego kroku relaksacji. Co więcej, ograniczeniem są też
kwestie czysto technologiczne - tekstury w których przechowywane są dane mogą
być podczas obliczeń wyłącznie przeznaczone do odczytu lub do zapisu. Tak więc
niemożliwy jest jednoczesny odczyt z danej tekstury oraz natychmiastowy zapis.

Problem ten został rozwiązany przez zmianę algorytmu rozwiązywania układów
równań liniowych na metodę relaksacji \ow{Jacobiego}. Główną różnicą jest moment
aktualizowania komórek siatki symulacji. W przeciwieństwie do metody \ow{Gaussa-
Seidela} nie następuje to natychmiast po obliczeniu nowej wartości dla danej
komórki, ale dopiero obliczeniu nowych wartości dla wszystkich komórek.
Uproszczony schemat tej relaksacji prezentuje algorytm \ref{alg:JacobiCPU}.

\begin{algorithm}[!h]
  \caption{Relaksacja metodą Jacobiego na CPU}
  \label{alg:JacobiCPU}
\begin{algorithmic}
\For {$0 \to relaxation\_steps$}
  \ForAll {$ \textrm{grid cells} $}
    \State $temp_{i,j}\gets f(x_{i,j}, x_{i+1,j}, x_{i-1,j}, x_{i,j+1}, x_{i,j-1})$
  \EndFor
  \ForAll {$ \textrm{grid cells} $}
    \State $x_{i,j}\gets temp_{i,j}$
  \EndFor
\EndFor
\end{algorithmic}
\end{algorithm}

Taki algorytm można już przetworzyć na wersję równoległą. Macierze $x$ oraz
$temp$  zostają zamienione na odpowiednie tekstury. Również w celach
wydajnościowych zawartość tekstur nie jest przepisywana tylko zamieniane są
ich referencje. Uproszczony schemat relaksacji metodą \ow{Jacobiego} na GPU
prezentuje algorytm \ref{alg:JacobiGPU}.

\begin{algorithm}[!h]
  \caption{Relaksacja metodą Jacobiego na GPU}
  \label{alg:JacobiGPU}
\begin{algorithmic}
\For {$0 \to relaxation\_steps$}
  \ForAll {$ \textrm{grid cells} $} [IN PARALLEL]
    \State $temp_{i,j}\gets f(x_{i,j}, x_{i+1,j}, x_{i-1,j}, x_{i,j+1}, x_{i,j-1})$
  \EndFor
  \State swap $x$ with $temp$
\EndFor
\end{algorithmic}
\end{algorithm}

\begin{figure}[!p]
\centering

\includegraphics[width=.45\textwidth]{img/linSolvCPU5}
\caption{Symulacja przy użyciu metody Gaussa-Seidela na CPU, 
5 kroków relaksacji}
\label{fig:linSolvCPU5}

\includegraphics[width=.45\textwidth]{img/linSolvGPU5}
\caption{Symulacja przy użyciu metody Jacobiego na GPU, 
5 kroków relaksacji}
\label{fig:linSolvGPU5}

\includegraphics[width=.45\textwidth]{img/linSolvGPU10}
\caption{Symulacja przy użyciu metody Jacobiego na GPU, 
10 kroków relaksacji}
\label{fig:linSolvGPU10}

\end{figure}

Oczywiście zmiana algorytmu rozwiązywania równań liniowych niesie ze sobą
pewne konsekwencje. Inna jest konwergencja tych dwóch algorytmów. Metoda \ow
{Gaussa-Seidela} jest szybciej zbieżna niż metoda \ow{Jacobiego}. Efekty
symulacji dla różnej konfiguracji algorytmów rozwiązywania układów równań
liniowych prezentują rysunki \ref{fig:linSolvCPU5}, \ref{fig:linSolvGPU5} oraz
\ref{fig:linSolvGPU10}.

Na tych rysunkach można zaobserwować wyraźne różnice w rezultatach symulacji. W
przypadku przedstawionej symulacji oczekiwanym wynikiem było powstanie ścieżki
wirowej von Kármána dla większej przeszkody. Przy implementacji równoległej metody
\ow{Jacobiego} widać, iż pożądany efekt symulacji zostaje utracony. Dlatego
też należy wykonać więcej kroków relaksacji. Okazało się, że wartością
wystarczającą jest dziesięć. Wartość ta została ustalona empirycznie, tak aby
wyniki symulacji odpowiadały oczekiwaniom oraz były jak najbardziej zbliżone
do rezultatów sekwencyjnych silników fizycznych. Jest to niezbędne, ponieważ
wersja równoległa aplikacji \en powinna być w pełni zgodnym i kompatybilnym
rozszerzeniem wersji podstawowej (sekwencyjnej, napisanej bez użycia
technologii \ow{WebGL}). Oczywiście wiąże się to ze spowolnieniem symulacji,
jednak w żadnym wypadku nie neguje opłacalności przeniesienia obliczeń na
kartę graficzną -- relaksacja na GPU jest wciąż dużo szybszych niż na CPU mimo
konieczności wykonania dwukrotnie większej liczby kroków.

Istnieją również implementacje algorytmu \ow{Gaussa-Seidela} na maszyny
równoległe. Nie są to dokładne kopie wersji sekwencyjnej, lecz imitują
obliczanie nowych wartości komórek na podstawie wartości z kroku relaksacji
poprzedniego i aktualnego. Doskonałe opracowanie i porównanie algorytmów
rozwiązujących układy równań liniowych na CPU oraz GPU, które okazało się
niezwykle przydatne podczas implementacji równoległych algorytmów dla \en,
zostało przygotowane przez G. Amadora oraz A. Gomesa \cite{LinSolvers}.
Podobną tematykę podejmuje praca magisterska Filipa Vesely
\cite{masterSolvers}. Niestety niskopoziomowa implementacja algorytmu \ow
{Gaussa-Seidela} na GPU wymaga wykonania dwóch procesów renderowania do
tekstury dla jednego kroku relaksacji -- powoduje to narzut czasowy, który w
efekcie niweluje zysk z szybszej konwergencji algorytmu.

\subsection{Podstawowy zarys implementacji}

Zgodnie z wnioskami sekcji \ref{subsec:srodPrzegInt}, ze względu na ograniczenia
środowiska przeglądarki internetowej wymuszona została implementacja
niskopoziomowa przy użyciu technologii \ow{WebGL}. Opiera się ona na
programowaniu jednostek cieniujących karty graficznej, głównie wykorzystując
programy fragmentów (ang. fragment programs). Taka metoda przeprowadzania
obliczeń ogólnego przeznaczenia na karcie graficznej wymaga zrozumienia
podstawowych zagadnień związanych z programowaniem grafiki trójwymiarowej (por.
sekcja \ref{subsec:niskProgJedn}).

W przypadku implementacji \ow{Energy2D} podstawowy schemat przeprowadzania
obliczeń fizycznych na GPU z wykorzystaniem technologii \ow{WebGL} wygląda
następująco:

\begin{itemize}

\item Przygotowanie oraz kompilacja programów wierzchołków oraz fragmentów,
które zawierają implementacje algorytmów silników fizycznych.

\item Przygotowanie tekstur, które przechowują dane symulacyjne. 

\item Przygotowanie danych geometrii płaszczyzny pokrywającej cały zakres
współrzędnych, które mieszczą się w obszarze renderowania.

\item Wykonywanie kolejnych kroków algorytmów fizycznych, co sprowadza się do:
	
	\begin{itemize}
	\item Renderowania wcześniej przygotowanej geometrii płaszczyzny przy użyciu
	wcześniej przygotowanych programów wierzchołków i fragmentów.

  \item Kopiowania danych z bufora ramki do wybranej tekstury przechowującej
	dane symulacji \footnote{Odbywa się to automatycznie z użyciem obiektu bufora 
  ramki (ang. FrameBuffer Object) z dołączoną do niego teksturą}.
	\end{itemize}

\end{itemize}

Poszczególne podpunkty w szerszym zakresie przybliżają sekcje od \ref{sec:progWierzFrag} do
\ref{sec:wykonNaGPU}.

\subsection{Programy wierzchołków oraz fragmentów}
\label{sec:progWierzFrag}

Właściwa implementacja algorytmów w programach wierzchołków i fragmentów
decyduje o poprawności oraz wydajności całej aplikacji. Bardzo często w
przypadku obliczeń ogólnego zastosowania na GPU, a także w przypadku
symulatora \ow{Energy2D}, większość pracy przypada na programy fragmentów.
Program wierzchołków zwykle sprowadza się do skopiowania wejściowych
współrzędnych wierzchołka i tekstury. Nie ma potrzeby jakiejkolwiek
modyfikacji geometrii renderowanej płaszczyzny. Dlatego też większość
programów wierzchołków w symulatorze Energy2D odpowiada poniższej
implementacji:

\begin{lstlisting}[language=GLSL, caption=Typowa implementacja programu
wierzchołków w symulatorze \ow{Energy2D} (język GLSL)]
attribute vec4 vertexPos;
attribute vec4 texCoord;

varying vec2 coord;
void main() {
  coord = texCoord;
  gl_Position = vertexPos;
}
\end{lstlisting}

Programy fragmentów nie są już tak trywialne, ponieważ wykonują właściwe
obliczenia. Bardzo ważne jest poprawne określenie współrzędnych tekseli
tekstury. Posiadając siatkę symulacji o wymiarach NxN, należy ją zrzutować na
zakres domyślnych współrzędnych tekstury, które zawierają się w przedziale [0,
1]. Jeżeli zrobi się to nieprawidłowo trudno będzie wykryć taki błąd, gdyż
domyślnie tekstury interpolują wartości leżące pomiędzy rzeczywistymi danymi.
Może to prowadzić do nieoczekiwanych rezultatów, stąd niezwykle istotne jest
precyzyjne określanie współrzędnych. W tym celu, praktycznie każdy program
fragmentów zawierał wektor o nazwie \ow{grid} równy (1.0 / N, 1.0 / N), gdzie
NxN to wymiary siatki symulacyjnej. Dodając lub odejmując odpowiedni jego
komponent można uzyskać dokładną wartość komórki sąsiedniej. Warto też
pamiętać o fakcie, iż pierwsza kolumna tekseli nie ma współrzędnej X równej
0.0, lecz 0.5 / N. To samo dotyczy pierwszego rzędu tekseli oraz współrzędnej
Y równej 0.5 / N. Błędne założenie, iż ta kolumna lub rząd tekseli mają
współrzędne 0.0 prowadzi do problemów przy wymuszaniu warunków brzegowych
podczas symulacji. Ostatecznie, typowy szkielet programów fragmentów wygląda
następująco:

\begin{lstlisting}[language=GLSL, caption=Szkielet implementacji programu
fragmentów w symulatorze \ow{Energy2D} (język GLSL)]
uniform sampler2D simulationData;
uniform vec2 grid;
varying vec2 coord;

vec4 F(vec4 data) {
  // Funkcja wykonująca właściwe obliczenia dla danego kroku symulacji.
  // ...
}

void main() {
  vec4 data = texture2D(simulationData, coord);
  // Instrukcja warunkowa sprawdzająca czy nie są przetwarzane brzegi siatki.
  if (coord.x > grid.x && coord.x < 1.0 - grid.x &&
      coord.y > grid.y && coord.y < 1.0 - grid.y) {
    data = F(data);
  }
  gl_FragColor = data;
}
\end{lstlisting}

Oczywiście program, który wymusza warunki brzegowe posiada odwrotny warunek w
liniach 13 oraz 14. Implementacja funkcji $F$ nie jest przytoczona, gdyż nie
da się wyróżnić jakiegoś ogólnego jej schematu czy wzoru. Można powiedzieć, że
przenosząc dany krok algorytmu do języka GLSL, funkcja $F$ stanowi
implementację ,,wewnętrznych'' instrukcji zagnieżdżonych pętli iterujących po
wszystkich komórkach siatki symulacyjnej.


\subsection{Organizacja danych w pamięci karty graficznej}
\label{sec:orgDanychWGPU}
Dane symulacji (takie jak np. macierz temperatury czy macierz prędkości płynu)
przechowywane są w dwuwymiarowych teksturach zmiennoprzecinkowych. Tego typu
tekstury nie wchodzą w skład podstawowej specyfikacji \ow{WebGL 1.0}
(\cite{WebGLSpec}). W związku z tym wymagane jest użycie rozszerzenia
\ow{OES\_texture\_float}, które jest dostępne na zdecydowanej większości
współczesnych urządzeń.

Niezwykle ważną kwestią jest organizacja danych w teksturze. Jest sporo
możliwości ponieważ tekstura z zasady nie jest wierną kopią tablicy
\js, a obiektem przystosowanym do przechowywania obrazów. Przede
wszystkim posiada kanały kolorów. Dlatego też można rozważyć kilka
potencjalnych sposobów na rozmieszczenie danych z tablic \js w teksturach
\ow{WebGL}:

\begin{itemize}

\item Schemat 1 -- wymaga użycia tekstur jednokanałowych (format ALPHA lub
LUMINANCE). Dane znajdują się w jedynym dostępnym kanale koloru (Alpha lub
Luminance), stąd jedna tekstura odpowiada jednej tablicy \js. Dalej nazywany
schematem \textbf{A1} na potrzeby niniejszego opracowania.

\item Schemat 2 -- wymaga użycia tekstur czterokanałowych (format RGBA). Dane
znajdują się tylko w jednym kanale (dowolnym -- R, G, B, lub Alpha), stąd
jedna tekstura odpowiada jednej tablicy \js. Dalej nazywany schematem
\textbf{RGBA1} na potrzeby niniejszego opracowania.

\item Schemat 3 -- wymaga użycia tekstur czterokanałowych (format RGBA). Dane
znajdują się w każdym z kanałów, stąd jedna tekstura odpowiada czterem
tablicom \js. Dalej nazywany schematem \textbf{RGBA4} na potrzeby niniejszego
opracowania.

\item Schemat 4 -- wymaga użycia tekstur czterokanałowych (format RGBA). Dane
znajdują się w każdym z kanałów, natomiast rozmiar tekstury jest zredukowany
czterokrotnie, gdyż każdy kanał odpowiada jednej ćwiartce tablicy. Dlatego też
jedna tekstura odpowiada jednej tablicy \js. Dalej nazywany schematem
\textbf{RGBA1/4} na potrzeby niniejszego opracowania.

\end{itemize}

Każdy z powyższych schematów został przetestowany podczas implementacji
symulatora \ow{Energy2D}. Wyniki testów wydajnościowych przedstawia wykres
\ref{fig:texPerf}\footnote{Szczegółowa metodologia testów oraz konfiguracja
sprzętowa nie zostaje tutaj przedstawiona, gdyż test ten ma tylko na celu
pokazanie względnych różnic między schematami. Szczegółowe testy
wydajnościowe całej aplikacji, wraz z dokładnym opisem metodologii prezentuje
rozdział \ref{sec:testyWydajnosciowe}}.

\begin{figure}[!h]
\centering
\includegraphics[width=.9\textwidth]{img/texPerf}
\caption{Wpływ organizacji danych w teksturach na wydajność symulacji 
\ow{Energy2D}}
\label{fig:texPerf}
\end{figure}

Rozwiązaniem optymalnym ze względu na wydajność oraz czytelność kodu wydaje się
schemat A1. Umożliwia on przeniesienie danych z tablic w relacji 1:1 do tekstur.
Do tego, każdy odczyt czy zapis do tekstury dotyczy tylko jednego kanału.
Niestety na większości urządzeń nie ma możliwości dołączenia tekstury typu
ALPHA lub LUMINANCE do własnego obiektu bufora ramki (ang. FrameBuffer Object)
-- czyli renderowania do takiej tekstury\footnote{Testy zostały przeprowadzone
na systemie Linux (Ubuntu 12.04) i przeglądarce Google Chrome w wersji 22, która
korzysta z natywnego sterownika \ow{OpenGL}. Natomiast na tej samej konfiguracji
sprzętowej, ale działającej pod kontrolą systemu Windows nie udałoby się
uruchomić aplikacji. Podobnie w przypadku równie popularnego systemu
operacyjnego Mac OS X.}. Wyklucza to możliwość użycia tego rozwiązania mimo
oczywistych zalet. Być może w przyszłości rozwój specyfikacji \ow{WebGL} na to
pozwoli.

Warto tutaj nadmienić, że specyfikacja \ow{WebGL} (\cite{WebGLSpec}) nie
gwarantuje, iż jakikolwiek z formatów tekstur zmiennoprzecinkowych będzie
zaakceptowany jako cel renderowania. Programista powinien wykonać test, aby
sprawdzić czy maszyna użytkownika wspiera daną konfigurację. Można to zrobić w
następujący sposób:

\begin{lstlisting}[language=JavaScript, caption=Weryfikacja poprawności formatu
i typu tekstury używanej jako cel renderowania]
var gl 		= getWebGLContext(),
	texture = gl.createTexture(),
	fbo 	= gl.createFramebuffer();

if (!gl.getExtension('OES_texture_float')) {
	throw new Error("Rozszerzenie OES_texture_float niedostępne.");
}
gl.bindTexture(gl.TEXTURE_2D, texture);
gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA, 128, 128, 0, gl.RGBA, gl.FLOAT, null);
gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, texture, 0);
if (gl.checkFramebufferStatus(gl.FRAMEBUFFER) !== gl.FRAMEBUFFER_COMPLETE) {
	throw new Error("Dana tekstura nie jest wspierana jako cel renderowania.");
}
\end{lstlisting}

W praktyce tekstury zmiennoprzecinkowe posiadające cztery kanały kolorów są
najczęściej akceptowanym formatem do którego można zapisywać dane podczas
renderowania. Dlatego też schematy organizacji danych RGBA1, RGBA4 oraz RGBA1/4
używają takiego formatu tekstury.

Schemat RGBA1 posiada te same zalety co A1 jeśli chodzi o organizacje i
czytelność kodu źródłowego, jednak w tym przypadku dochodzi do dużego narzutu
wydajności związanego z odczytem i zapisem tekstur. Przy każdej z tych operacji
karta graficzna musi odczytać cztery kanały, jednak praktycznie wykorzystywany
jest tylko jeden z nich. Operacje dostępu do pamięci są czasochłonne, dlatego
też taka organizacja danych nie jest korzystna ze względów wydajnościowych.

Rozwiązaniem tego problemu są schematy RGBA4 oraz RGBA1/4. Organizacja danych
wg. trzeciego schematu pozwala zredukować narzut związany z odczytem oraz
zapisem pod warunkiem dobrej organizacji danych w teksturach. Pomysł ten
obrazuje diagram \ref{fig:rgba4Tex}. 

\begin{figure}[!p]
\centering
\includegraphics[width=.65\textwidth]{img/rgba4Tex}
\caption{Schemat organizacji danych z czterech macierzy w jednej teksturze}
\label{fig:rgba4Tex}
\end{figure}

\begin{figure}[!p]
\centering
\includegraphics[width=.65\textwidth]{img/rgba14Tex}
\caption{Schemat organizacji danych z macierzy $NxN$ w teksturze $(N/2)x(N/2)$}
\label{fig:rgba14Tex}
\end{figure}

Przy korzystnym ułożeniu danych jest możliwe praktycznie całkowite
zredukowanie narzutu związanego z odczytem wszystkich czterech kanałów
tekstury, jednak w praktyce jest to często niewykonalne. Przez korzystne
rozmieszczenie danych przyjmuje się takie ich ułożenie, żeby program jednostek
cieniujących odczytując teksturę faktycznie korzystał z danych zawartych w
każdym z kanałów. Podobnie przy zapisie, program renderujący powinien
modyfikować wszystkie cztery kanały. W przypadku symulatora \ow{Energy2D}
udało się uzyskać taką organizację danych, żeby odczyt był w znacznym stopniu
zoptymalizowany. Jednak podczas zapisu modyfikowany jest tylko jeden lub dwa
kanały (kanał zawierający dane o temperaturze lub kanały zawierające
komponenty wektorów prędkości). Mimo tego schemat RGBA4 okazał się około 54\%
wydajniejszy od schematu RGBA1.

Ciekawą organizacją danych może wydawać się również pomysł przedstawiony w
schemacie RGBA1/4. Pozwala przechowywać tablicę JavaScript o wymiarach
\ow{NxN} w teksturze o wymiarach \ow{(N/2)x(N/2)}. Pomysł ten obrazuje diagram
\ref{fig:rgba14Tex}. Taki układ danych teoretycznie posiada sporo zalet
cechujących użycie tekstur jednokanałowych. Jednak znacznie zmniejsza on
czytelność kodu i komplikuje implementacje. Utrudnione zostaje przede
wszystkim kontrolowanie warunków brzegowych, programy jednostek cieniujących
przetwarzają cztery pola jednocześnie i stają się bardzo skomplikowane. W
przypadku \ow{Energy2D} komplikacje programów fragmentów były tak znaczne, że
w efekcie odnotowano bardzo słabe rezultaty pod względem wydajności. Symulacja
okazała się być wolniejsza około 34\% od symulacji korzystającej ze schematu
RGBA4.

Ostatecznie w \ow{Energy2D} zdecydowano się zastosować schemat RGBA4. Zapewnia
on najlepszy kompromis pomiędzy dobrą wydajnością oraz wsparciem przez
większość dostępnych obecnie urządzeń.

\subsection{Geometria}

Podczas wykonywania obliczeń ogólnego przeznaczenia przy pomocy bezpośredniego
programowania jednostek cieniujących karty graficznej niezbędne jest
stworzenie obiektu (a właściwie jego geometrii), który będzie renderowany. W
teorii może być to dowolna bryła. Jednak najczęściej pożądane jest, aby
program fragmentów wykonał się dla każdej komórki siatki symulacyjnej, którą
stanowią piksele tekstury. Można to osiągnąć renderując płaszczyznę, która
pokrywa całą dostępną przestrzeń renderowania.

Również w przypadku \en wykorzystywane zostało takie podejście. Do
przechowywania własności płaszczyzny użyte zostały bufory wierzchołków
(ang. vertex buffers) oraz indeksów (ang. index buffers) znajdujące się w
pamięci karty graficznej. Dzięki temu podczas wielokrotnego renderowania tej
samej płaszczyzny, nie jest konieczne nieustanne przesyłanie atrybutów
wierzchołków do pamięci GPU.

W finalnej implementacji \en zostały przygotowane klasy pomocnicze
zarządzające geometrią oraz buforami. Jednak najprostszy sposób na stworzenie
płaszczyzny pokrywającej cały obszar renderowania oraz umieszczenie jej
atrybutów w pamięci karty graficznej zaprezentowany jest poniżej:

\begin{lstlisting}[language=JavaScript, caption=Definicja geometrii
płaszczyzny pokrywającej cały obszar renderowania]
var gl           = getWebGLContext(),
    vertexBuffer = gl.createBuffer(),
    indexBuffer  = gl.createBuffer(),
    vertexData,
    indexData;

// Współrzędne wierzchołków.
vertexData = new Float32Array([
    -1, -1, 
     1, -1,
    -1,  1,
     1,  1
]);
// Transfer do bufora wierzchołków.
gl.bindBuffer(gl.ARRAY_BUFFER, vertexBuffer);
gl.bufferData(gl.ARRAY_BUFFER, vertexData, gl.STATIC_DRAW);
// Indeksy trojkątów płaszczyzny.
indexData = new Uint16Array([
    0, 1, 2,
    2, 1, 3
]);
// Transfer do bufora indeksów.
gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, indexBuffer);
gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, indexData, gl.STATIC_DRAW);
\end{lstlisting}

\subsection{Wykonywanie kroków algorytmów na GPU}
\label{sec:wykonNaGPU}

Dysponując skompilowanymi programami wierzchołków oraz fragmentów, danymi
symulacji zapisanymi w teksturach dwuwymiarowych oraz niezbędną geometrią,
można przejść do faktycznej realizacji kroków algorytmów
zapisanych w programach jednostek cieniujących.

W przypadku technologii WebGL narzuca się schemat związany z techniką
renderowania do tekstury przy użyciu obiektu bufora ramki (ang. \emph{Frame
Buffer Object}). Związanie takiego obiektu z wybraną teksturą, a następnie
uaktywnienie przed właściwym renderowaniem, powoduje iż karta graficzna
automatycznie kopiuje zawartość bufora ramki do tekstury po zakończonym
renderowaniu.

Technika ta ma jednak kilka ograniczeń. Tekstura związana z obiektem bufora
ramki (czyli przeznaczona do zapisu) nie może być używana jednocześnie do
odczytu wartości. W związku z tym zawsze trzeba używać minimalnie jednej
tekstury tymczasowej. Następnie należy kopiować jej zawartość do tekstury
docelowej lub podmienić referencje. Oczywiście modyfikacja wyłącznie
referencji  jest znacznie efektywniejsza przez co i częściej stosowana.
Całościowo taki schemat nazywany jest ,,ping-pong rendering''.

Implementacja tego schematu w języku \ow{JavaScript} przy użyciu technologi
\ow{WebGL} nie odbiega znacząco od implementacji w innych językach przy użyciu
,,tradycyjnego'' API OpenGL. Mark Harris zaprezentował najważniejsze aspekty
tej techniki w swoich artykułach opublikowanych w serii \emph{GPU Gems}
(\cite{GPUConcepts} oraz \cite{GPUFluid}).