ZASTOSOWANIE TECHNOLOGII GPGPU DO PRZYSPIESZENIA OBLICZEŃ
W ZAGADNIENIACH BRZEGOWYCH
ROZWIĄZYWANYCH ZA POMOCĄ PURC
Andrzej Kużelewski
1a, Eugeniusz Zieniuk
1b1Zakład Metod Numerycznych, Wydział Matematyki i Informatyki, Uniwersytet w Białymstoku
aakuzel@ii.uwb.edu.pl, bezieniuk@ii.uwb.edu.pl
Streszczenie
Celem niniejszej pracy jest zbadanie możliwości zastosowania nowoczesnych rozwiązań równoległego przetwa- rzania danych w celu przyspieszenia procesu obliczeniowego pojawiającego się w trakcie numerycznego rozwiązy- wania PURC. Autorzy wzięli pod uwagę możliwość wykorzystania procesorów kart graficznych do obliczeń nume- rycznych w zastosowaniach ogólnych (general-purpose GPU – GPGPU), a w szczególności technologię CUDA (ComputeUnified Device Architecture) firmy NVidia®.
Słowa kluczowe: PURC, CUDA, zagadnienia brzegowe
APPLICATION OF GPGPU FOR ACCELERATION
OF CALCULATIONS IN BOUNDARY VALUE PROBLEMS SOLVING USING PIES
Summary
The main purpose of this paper is examination of an application of modern parallel computing solutions to speed up the calculation in the numerical solution of parametric integral equations systems (PIES). The authors considered use of graphics cards programming in general-purpose applications (GPGPU), particularly NVidia®
CUDA(Compute Unified Device Architecture).
Keywords: PIES, CUDA, boundary value problems
1. WSTĘP
Autorzy pracy od wielu lat stosują parametryczne układy równań całkowych (PURC) do rozwiązywania zagadnień brzegowych. Dotychczas PURC wykorzysty- wane były do rozwiązywania zagadnień brzegowych 2D i 3D modelowanych równaniami różniczkowymi cząst- kowymi takimi jak: Laplace’a, Poissona, Helmholtza czy Naviera-Lamego [6-12]. PURC, stosowane do numerycz- nego rozwiązywania wspomnianych równań, w swoim formalizmie matematycznym uwzględniają kształt brzegu rozpatrywanego zagadnienia brzegowego. W celu uwzględnienia brzegu w PURC wykorzystano krzywe (np. Béziera, B-spline i inne) i płaty powierzchni (np.
płaty Coonsa, Béziera i inne) stosowane w grafice kom- puterowej. Do praktycznego zdefiniowania kształtu brzegu zadawano niewielką liczbę punktów kontrolnych.
Dotychczas autorzy w swoich pracach skupiali się na badaniu dokładności wyników uzyskiwanych za pomocą PURC w porównaniu z klasycznymi już metodami takimi jak MES czy MEB lub metodami analitycznymi.
Jednakże należy zwrócić także uwagę na inny problem dotyczący każdego algorytmu numerycznego - jego czasochłonność obliczeniową. Zwiększanie liczby danych wejściowych skutkuje wydłużeniem czasu pracy algo- rytmu. Także w przypadku PURC zostało to zaobser-
wowane, szczególnie przy bardziej złożonych równaniach (np. Naviera-Lamego) i zwiększeniu wymiarowości problemu z 2D na 3D. Zwiększenie szybkości działania algorytmu może być uzyskane na różne sposoby: wyko- rzystanie najwydajniejszych dostępnych maszyn, wyko- rzystanie maszyn wieloprocesorowych i klastrów lub procesorów graficznych.
W ostatnich latach szczególnie zauważalny jest wzrost zainteresowania naukowców możliwością wyko- rzystania procesorów kart graficznych do obliczeń nume- rycznych w zastosowaniach ogólnych (GPGPU), ze względu na możliwość uzyskania bardzo dużych przyro- stów wydajności aplikacji [3-5]. Jest to związane z architekturą kart graficznych (wieloprocesorowa i wielowątkowa), bardzo szybkimi zmiennoprzecinkowy- mi jednostkami arytmetycznymi oraz wykorzystaniem szybkich pamięci. Dlatego autorzy postanowili przyjrzeć się możliwości zastosowania nowoczesnych rozwiązań równoległego przetwarzania danych w celu przyspiesze- nia procesu obliczeniowego pojawiającego się w trakcie numerycznego rozwiązywania PURC. Wzięto pod uwagę technologię CUDA (ComputeUnified Device Architectu- re) firmy NVidia® [1].
Celem niniejszej pracy jest przedstawienie możliwości przyspieszenia obliczeń numerycznych, niezbędnych do rozwiązywania zagadnień brzegowych 3D, modelowa- nych równaniami Naviera-Lamego i rozwiązywanych za pomocą PURC z wykorzystaniem technologii CUDA.
2. PURC
DLA TRÓJWYMIAROWYCH RÓWNAŃ NAVIERA-LAMEGO
PURC dla trójwymiarowych równań Naviera- Lamego jest wynikiem analitycznej modyfikacji klasycz- nych BRC. Szczegółowy opis metodologii takiej modyfi- kacji dla zagadnień 2D modelowanych różnymi równa- niami różniczkowymi przedstawiono w pracach [6-9].
Uogólniając wspomnianą metodologię na zagadnienia 3D, uzyskano następującą postać PURC [7]:
12 , , , , ,
, , , , , ! ! (1)
przy czym: " # # , " # # , " # # , " # # , %&, '( 1,2,… , *,gdzie n – jest liczbą płatów powierzchniowych tworzących zamknięty trój- wymiarowy obszar.
Funkcje podcałkowe , , , , , , , w (1) są przedstawiane w następującej postaci macie- rzowej [7]:
, , , +, "- ./01 12 1 3
12 122 123
13 132 133
4 (2)
, , , 5, "- /607 72 73
72 722 723
73 732 733
4, (3)
gdzie: 8 2 :-9 .
Jawną postać elementów macierzy (2) i (3) można odnaleźć w [7].
Funkcja , jest jakobianem, natomiast jądra (2) i (3) uwzględniają w sposób analityczny w swoim formalizmie matematycznym zamknięty kształt brzegu.
Jest on kreowany za pomocą odpowiednich zależności pomiędzy płatami &, ' 1,2,. . . , * zdefiniowanymi we współrzędnych kartezjańskich za pomocą następujących zależności [7]:
< 7 , 7 , ,
<2 72 , 72 , ,
<3 73 , 73 , ,
< =< > <2> <3,
(4)
gdzie 7 , , 72 , , 73 , są skalarnymi składowymi wektorowego płata powierzchni ,
?7 , 72 , 73 , @A, zależnego od parametrów v, w. Powyższy zapis jest też słuszny dla płata oznaczonego indeksem l z parametrami v1, w1, tzn.
dla j = l oraz dla parametrów v = v1 i w = w1.
W zagadnieniach przestrzennych modelowanych
z wykorzystaniem PURC, wektorowe
cje , przyjmują postać znanych z grafiki kompute- rowej parametrycznych płatów powierzchni. Możliwość uwzględnienia opisu brzegu bezpośrednio w formule matematycznej jest podstawową zaletą metody oblicze- niowej PURC w porównaniu z tradycyjnymi brzegowymi równaniami całkowymi (BRC). W klasycznych BRC opis brzegu nie jest uwzględniony w formalizmie mate- matycznym tego równania, lecz zdefiniowany bardzo ogólnie za pomocą całki brzegowej. Powstaje konieczność podziału brzegu na elementy, jak to ma miejsce w klasycznej MEB. Zalety modelowania brzegu bezpo- średnio w równaniach matematycznych (PURC) były przedstawiane na zagadnieniach dwuwymiarowych [6-9]
oraz trójwymiarowych [7,11,12].
Zastosowanie PURC do rozwiązywania zagadnień płaskich oraz przestrzennych pozwoliło na wyeliminowa- nie konieczności przeprowadzenia dyskretyzacji na poziomie zarówno wspomnianej już deklaracji brzegu, a także funkcji brzegowych. W ramach dotychczasowych prac badawczych funkcje brzegowe zarówno zadawane jako warunki brzegowe, jak również poszukiwane w wyniku rozwiązania PURC przyjmowały postać szeregów aproksymujących z funkcjami bazowymi Cze- byszewa [6,7]. W odróżnieniu od zagadnień dwuwymia- rowych (modelowanych równaniem Laplace’a) problemy opisywane równaniem Naviera-Lamego wymagają uwzględnienia funkcji brzegowych w postaci wektorowej.
W związku z tym uogólnione zostały wykorzystywane
dotychczas szeregi aproksymujące tak, aby reprezento- wały one skalarne składowe wektora czeń , oraz naprężeń , [7]. Dla każdego płata powierzchni modelujących segment brzegu j zostały one przedstawione w następującej postaci:
, ∑H ∑FD G CDEC ED
C G , (5)
, ∑HC G∑FD GICDEC ED , (6) gdzie CD, ICD są poszukiwanymi współczynnikami, zaśEC ,ED to funkcje bazowe, które mają postać wielomianów interpolacyjnych Lagrange’a:
EJ K L"LM L"L … L"LN L"LNO … L"LP QL "LMRQL "L R…QL "LN RQL "LNOR…QL "LPR, (7) gdzie k={p, r}, l={M, N}, x={v, w}.
Po podstawieniu (5) i (6) do PURC (1) i po zapisa- niu go w tzw. punktach kolokacji [2] otrzymuje się układ równań algebraicznych względem niewiadomych współ- czynników, po rozwiązaniu którego otrzymamy współ- czynniki CD(5) lubICD(6). Należy mieć na uwadze też to, że zawsze jedne z tych współczynników na poszcze- gólnych segmentach będą otrzymywane w wyniku aproksymacji zadanych warunków brzegowych za pomo- cą szeregów aproksymujących (5) lub (6).
Po rozwiązaniu PURC otrzymano rozwiązanie za- gadnienia brzegowego tylko na jego brzegu, reprezento- wane w postaci szeregów (5) lub (6). W celu znalezienia rozwiązań w obszarze klasyczna tożsamość całkowa znanaz BRC została analitycznie zmodyfikowana i przekształcona analogicznie jak w przypadku zagadnień dwuwymiarowych [6,7].Otrzymana została tożsamość całkowa, która wykorzystuje rozwiązania na brzegu (5) i (6), uzyskane uprzednio w wyniku rozwiązania PURC.
Zmodyfikowana tożsamość całkowa przyjmuje następu- jąca postać [7]:
S TU S, , ,
U S, , , V , ! ! (8) Funkcje podcałkowe w tożsamości (8) przedstawiane są w następującej postaci:
U S, , +, "- .DWX1U 1U2 1U3
1U2 1U22 1U23 1U3 1U32 1U33
Y, (9) U S, , 5, "- DW" 6X7Z 7Z2 7Z3
7Z2 7Z22 7Z23 7Z3 7Z32 7Z33
Y, (10)
gdzie [W 7 , K , [W2 72 , K2, [W3
73 , K3, [W =[W > [W2> [W3.
Jawna postać elementów macierzy (9) i (10) jest przedstawiona w [7].
Funkcje podcałkowe w tożsamości całkowej (8) są wizualnie bardzo podobne do funkcji (2) i (3). Różnica polega na tym, że w funkcjach (9) i (10) poza płatami powierzchni definiującymi kształt brzegu występują współrzędne punktów w obszarze w którym interesuje nas rozwiązanie.
3. PROGRAMOWANIE KART GRAFICZNYCH W OGÓLNYCH ZASTOSOWANIACH –
TECHNOLOGIA CUDA
Architektura procesorów graficznych (GPU) zdecy- dowanie odbiega od klasycznych procesorów (CPU).
GPU składa się z wielu jednostek zmiennoprzecinkowych (FPU) i układów arytmetyczno-logicznych (ALU).
Związane jest to z charakterem wykonywanych operacji – te same operacje wykonywane są równolegle na dużej ilości danych (w grafice typowo: pikselach, tekselach lub wierzchołkach), stąd GPU są klasyfikowane jako SIMD (single instruction, multiple data).
Technologia CUDA firmy Nvidia® [1] jest dużo prostsza i bardziej intuicyjna w programowaniu w porównaniu do tradycyjnych metod programowania GPGPU. Jest rozszerzeniem klasycznego języka C posiadającym łatwe do przyswojenia odniesienia do paradygmatów programowania równoległego. Kod programu dzielony jest na dwie, mogące się przeplatać, części: obliczenia równoległe wykonywane na jednostce GPU (nazywanym device) oraz pozostałe operacje wykonywane na CPU (nazywanym host). Schemat blokowy działania technologii CUDA przedstawia rys. 1.
Procedurę obliczeń można podzielić na cztery podsta- wowe kroki: 1. inicjacja programu (host), 2. skopiowanie danych z hosta na device, 3. wykonywanie obliczeń na GPU, 4. skopiowanie danych z device na host.
Technologia CUDA opiera się na wykorzystaniu ska- lowalnej tablicy wielowątkowych procesorów strumie- niowych (Streaming Multiprocessor – SM). Najpopular- niejsze podejście do programowania odbywa się na zasadzie podziału problemu na mniejsze zagadnienia, których instrukcje skupione są w funkcjach nazywanych jądrami (z ang. kernels). Jądra są wykonywane na pewnej części wątków (z ang. threads) podzielonych na bloki (z ang. blocks). Każde 32 wątki w bloku stanowią tzw. warp. Wszystkie wątki w obrębie bloku mogą komunikować się ze sobą i być synchronizowane za pomocą funkcji synchronizujących. W celu wykorzysta- nia większej liczby wątków w ramach jednego jądra, bloki wątków grupowane są w kraty (z ang. grids).
Wątki mogą pobierać dane z różnych rodzajów pa- mięci GPU. Każdy wątek ma dostęp do pamięci global- nej GPU (globalmemory), która ma dużą pojemność, ale niską przepustowość. Wątki w ramach jednego bloku mogą korzystać z pamięci współdzielonej (sharedmemo-
ry) umieszczonej na mikroprocesorze (z ang. on-chip memory). Jest ona wspólna dla wszystkich wątków w ramach jednego bloku i ma rozmiar 16 kB. Jest nawet do 100 razy szybsza od pamięci globalnej. Kolejnym typem jest pamięć lokalna – prywatna dla każdego wątku. Ma rozmiar 8192 bajty na multiprocesor. Istnieją jeszcze pamięć stała (constantmemory) oraz pamięć tekstur.
Rys. 1. Schemat przepływu danych w CUDA.
Pomimo że CUDA jest bardzo wygodnym narzę- dziem do zrównoleglenia obliczeń z wykorzystaniem mocy obliczeniowej kart graficznych, jednak optymaliza- cja wydajności programów jest dość trudna. Przede wszystkim należało skupić się na trzech zasadniczych problemach: optymalizacja dostępu do pamięci, konfigu- racja parametrów uruchamiania jąder i optymalizacja kodu (instrukcji).
Jeszcze jednym, ważnym z punktu widzenia obliczeń naukowych, problemem jest precyzja arytmetyki zmien- noprzecinkowej. Najnowsze karty graficzne pozwalają na operacje zmiennoprzecinkowe podwójnej precyzji, starsze – tylko na pojedynczą precyzję. Decydując się na zasto- sowanie technologii CUDA, należy mieć ten fakt na uwadze, szczególnie w przypadku obliczeń wymagają- cych bardzo wysokiej dokładności.
4. PRZYKŁAD NUMERYCZNY
W celu zbadania przyspieszenia zrównoleglonego PURC w stosunku do jego niezrównoleglonej wersji wzięto pod uwagę kształt obszaru przedstawiony na rys. 2a, ale matematycznie modelowany równaniami Naviera-Lamego z warunkami brzegowymi zadanymi za pomocą funkcji analitycznych będących dokładnym rozwiązaniem równania:
(11) W obliczeniach przyjęto moduł Younga E=1MPa oraz stałą Poissona ν=0.25. Wykonano testy porównaw- cze dla różnej liczby punktów kolokacji na płatach.
Kształt brzegu zdefiniowany został za pomocą 24 prostokątnych oraz 6 trójkątnych płatów Béziera i 6 prostokątnych płatów Coonsa (rys. 2b). Pełna deklara- cja brzegu złożonego w sumie z 36 płatów wymagała 236 punktów (kontrolnych dla płatów Béziera i narożnych dla Coonsa).
W badaniach użyto karty NVidiaQuadro FX 580 z 512 MB pamięci, 4 multiprocesorami (SM) po 8 rdzeni CUDA każdy. Dopuszczalny typ operacji zmiennoprze- cinkowych – pojedynczej precyzji.
Rys. 2. a) Kształt rozważanej geometrii brzegu, b) widok z podziałem na płaty
Autorzy w niniejszej pracy przedstawiają wyniki wstęp- nych badań nad problemem przyspieszenia obliczeń numerycznych w PURC za pomocą GPGPU i technolo- gii CUDA. W celu oceny przydatności tej technologii zdecydowano się na przyspieszenie najbardziej czaso- chłonnej operacji jaką jest numeryczne obliczania całek w PURC. Nie dokonano optymalizacji obliczeń na GPU – dane przechowywane są w pamięci globalnej.
W związku z tym problem badawczy podzielono na dwie części: przyspieszenie generowania układu równań w PURC (wyniki zestawiono w tabeli 1) oraz przyspie- szenie generowania rozwiązań (wyniki zestawiono w tabeli 2).
Liczba punktów kolokacji ma wpływ na dokładność uzyskiwanych rozwiązań, a także na czas obliczeniowy.
Więcej informacji na ten temat można odnaleźć w pracy [7]. Uzyskane przyspieszenie działania procedury gene- rowania układu równań w PURC z wykorzystaniem CUDA w stosunku do wersji niezrównoleglonej jest niezależne od stosowanej liczby punktów kolokacji i bliskie 4.
Tabela 1. Przyspieszenie generowania układu równań Liczba
punktów kolokacji*
Czas realizacji obliczeń [s]
Przyspie- szenie Wersja
niezrównole- glona (C++)
Wersja zrównoleglo-
na (CUDA)
4 9.206 2.348 3.92
9 40.510 8.502 4.76
16 128.766 29.999 4.29
25 321.206 89.491 3.56
*w dziedzinie płata
Badania przyspieszenia generowania rozwiązań wy- konano na podstawie obliczeń w 24 punktach we wnę- trzu rozmieszczonych jak na rys 3.
Tabela 2. Przyspieszenie generowania rozwiązań Liczba
punktów kolokacji*
Czas realizacji obliczeń [s]
Przyspie- szenie Wersja
niezrównole- glona (C++)
Wersja zrównoleglo-
na (CUDA)
4 3.179 0.552 5.76
9 7.259 1.232 5.89
16 13.036 2.355 5.53
25 20.456 3.433 5.96
*w dziedzinie płata
Podobnie jak w przypadku generowania układu rów- nań, uzyskane przyspieszenie działania procedury gene- rowania rozwiązań z wykorzystaniem CUDA w stosunku do wersji niezrównoleglonej jest niezależne od stosowanej liczby punktów kolokacji i bliskie 6.
5. WNIOSKI
W pracy przedstawiono możliwości przyspieszenia obliczeń numerycznych pojawiających się przy rozwią- zywaniu zagadnień brzegowych 3D modelowanych równaniami Naviera-Lamego i rozwiązywanych za pomocą metody obliczeniowej PURC z jednoczesnym wykorzystaniem technologii CUDA.
Zrównoleglenie aplikacji służącej do rozwiązywania zagadnień brzegowych za pomocą metody PURC i wykorzystaniu procesorów kart graficznych pozwala na uzyskanie blisko czterokrotnego przyspieszenia szybkości generowania układu równań oraz blisko sześciokrotnego przyspieszenia generowania rozwiązań.
Rys. 3. Rozmieszczenie punktów wewnętrznych wykorzystanych w badaniach przyspieszenia generowania rozwiązań
Należy przy tym mieć na uwadze fakt, że nie wykonano praktycznie żadnej optymalizacji wykorzystania pamięci GPU.
Niniejsza praca stanowi jedną z pierwszych prób zmierzenia się autorów z zastosowaniem technologii CUDA w celu skrócenia czasu obliczeniowego aplikacji
metody obliczeniowej PURC.W kolejnych badaniach autorzy zamierzają wykorzystać szybką pamięć współ- dzieloną (shared memory) oraz zoptymalizować dostęp do pamięci globalnej. Kolejnym kierunkiem rozwoju równoległej aplikacji PURC jest zrównoleglenie procedu- ry rozwiązywania układu równań algebraicznych.
Praca finansowana ze środków na naukę w latach 2010-2013 jako projekt badawczy.
Literatura
1. CUDA C Programming Guide. http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html [dostęp 10.10.2012].
2. Gottlieb D., Orszag S.A.: Numerical analysis of spectral methods: theory and applications. Philadelphia: SIAM, 1977.
3. Kiss I., Gyimóthy S., Badics Z., Pávó J.: Parallel realization of the element-by-element FEM technique by CUDA. “IEEE Transactions on Magnetics” 2012, Vol. 48,p. 507 - 510.
4. Owens J. D., Luebke D., Govindaraju N., Harris M., Kruger J., Lefohn A. E., Purcell T. J.: A survey of general- purpose computation on graphics hardware. “Eurographics 2005 State of the Art Reports”, Dublin, 2005, p. 21 - 51.
5. Takahashi T., Hamada T.: GPU-accelerated boundary element method for Helmholtz’ equation in three dimen- sions. “International Journal of Numerical Methods in Engineering” 2009, Vol. 80, p. 1295 - 1321.
6. Zieniuk E.: Bézier curves in the modification of boundary integral equations (BIE) for potential boundary-values problems. “International Journal of Solids and Structures” 2003, Vol. 40, p. 2301 - 2320.
7. Zieniuk E.: Metoda obliczeniowa PURC w rozwiązywaniu zagadnień brzegowych. Warszawa: PWN, 2013.
8. Zieniuk E., Bołtuć A.: Bézier curves in the modeling of boundary geometries for 2D boundary problems defined by Helmholtz equation. “Journal of Computational Acoustics” 2006, Vol. 14, p. 1 - 15.
9. Zieniuk E., Bołtuć A.: Non-element method of solving 2D boundary problems defined on polygonal domains modeled by Navier equation. “International Journal of Solids and Structures” 2006, Vol. 43, p. 7939 - 7958.
10. Zieniuk E., Szerszeń K., Bołtuć A.: Globalne obliczanie całek po obszarze w PURCdla dwuwymiarowych zagad- nień brzegowych modelowanych równaniem Naviera-Lamego i Poissona. „Modelowanie Inżynierskie” 2007,nr 33,s. 181 - 186.
11. Zieniuk E., Szerszeń K., Bołtuć A.: PURC w rozwiązywaniu trójwymiarowych zagadnień brzegowych modelowa- nych równaniami Naviera-Lamego w obszarach wielokątnych. „ModelowanieInżynierskie” 2011,nr 42,s. 487-495.
12. Zieniuk E., Szerszeń K.: Triangular Bézier patches in modelling smooth boundary surface in exterior Helmholtz problems solved by PIES. “Archives of Acoustics”2009, Vol. 34,p. 1 - 11.
