POŚREDNI NARASTANIA FAZY STAŁEJ KRZEPNIĘCIA WYKORZYSTUJĄCEJ MODEL OBLICZENIA RÓWNOLEGŁE W SYMULACJI

Pełen tekst

(1)MODELOWANIE INŻYNIERSKIE nr 55, ISSN 1896-771X. OBLICZENIA RÓWNOLEGŁE W SYMULACJI KRZEPNIĘCIA WYKORZYSTUJĄCEJ MODEL POŚREDNI NARASTANIA FAZY STAŁEJ Robert Dyja1a, Andrzej Grosser1b 1 a. Instytut Informatyki Teoretycznej i Stosowanej, Politechnika Częstochowska robert.dyja@icis.pcz.pl, bandrzej.grosser@icis.pcz.pl. Streszczenie W pracy zaprezentowano wyniki symulacji numerycznej krzepnięcia. Prezentowany model obliczeniowy zbudowany jest na podstawie równania przewodzenia ciepła z członem źródłowym w jawnym sformułowaniu pojemnościowym. Model narastania fazy stałej jest zgodny z pośrednim modelem, w którym zakładana jest skończona dyfuzja domieszki w fazie ciekłej oraz brak dyfuzji domieszki w fazie stałej. Wynikające z modelu równania różniczkowe rozwiązane są za pomocą metody elementów skończonych. Obliczenia były realizowane na komputerach o pamięci rozproszonej. Wykonano testy dla siatek elementów skończonych składających się z maksymalnie 25 milionów elementów. Testy skalowalności były przeprowadzone dla maksymalnie 2048 procesorów. Słowa kluczowe: metoda elementów skończonych, krzepnięcie, komputery dużej mocy. PARALLEL COMPUTING OF SOLIDIFICATION SIMULATIONS WITH INDIRECT SOLID PHASE GROWTH MODEL Summary The results of numerical simulations are presented in this paper. Computational model presented in this work uses heat transfer equation with heat source term in explicit enthalpy formulation. An indirect model for a solid phase growth was used. The indirect model assumes finite diffusion of inclusion in liquid phase and no diffusion in solid state. Resulting differential equations are solved with use of Bubnov-Galerkin Finite Element Method (space discretization). The calculations were performed on computer cluster with distributed memory positioned 7th on TOP500 list from November 2014. We carried out tests using meshes of up to 25 million of tetrahedral elements. The scalability tests were conducted for up to 2048 processor cores. Keywords: finite element method, solidification, high performance computing. 1. WSTĘP Modelowanie procesu krzepnięcia jest trudne do wykonania, ponieważ krzepnięcie obejmuje takie zjawiska jak: przepływ ciepła z uwzględnieniem przemiany fazowej, transport domieszki w mikro- i makroskali oraz występowanie co najmniej dwóch obszarów o różnych własnościach fizycznych. Z tych powodów oraz z faktu, że krzepnięcie jest procesem zachodzącym w czasie, wynika duża czasochłonność symulacji. Problem czasochłonności symulacji może być jednak rozwiązany przez obliczenia wykonywane równolegle.. Prezentowany model obliczeniowy zbudowany jest na podstawie równania przewodzenia ciepła z członem źródłowym w jawnym sformułowaniu pojemnościowym. Model narastania fazy stałej jest zgodny z pośrednim modelem, w którym zakładana jest skończona dyfuzja domieszki w fazie ciekłej oraz brak dyfuzji domieszki w fazie stałej. Wynikające z modelu równania różniczkowe rozwiązane są przy za pomocą metody elementów skończonych w sformułowaniu Bubnova-Galerkina (dyskretyzacja przestrzenna) oraz wstecznego schematu Eulera (dyskretyzacja czasu). Proces oddawania ciepła. 21.

(2) OBLICZENIA RÓWNOLEGŁE W SYMULACJI KRZEPNIĘCIA WYKORZYSTUJĄCEJ MODEL… (1) wraz z równaniem (5) otrzymuje się podstawowe sformułowanie entalpowe krzepnięcia, które ma postać:. do otoczenia był modelowany zgodnie z warunkiem brzegowym Newtona. Obliczenia zostały wykonane w autorskim programie, który został zaimplementowany z wykorzystaniem języka C++ i biblioteki zgodnej ze standardem MPI. Opracowany program umożliwia wczytanie obszaru o dowolnym kształcie złożonego z czworościennych elementów skończonych.. ∇ ∙ ∇ =. Równanie przewodzenia ciepła wraz ze źródłem ciepła ma postać [4, 5]: డ் డ௧. − . డு(்) డ௧. డ௙ೞ. డு డ். డு(்) డ௧. డ் డ௧. − ௦ . డ௙ೞ డ௧. ் ೝ೐೑. + ௦ (1 −

(3) ௦ ()). డ௧. = . డ் డ௧. − ௦ . డ௙ೞ డ௧. డு డ். (7). డ் డ௧. = − ௦ . డ௙ೞ. (8). డ். = − ௦ . . డ௙ೞ డ் డ௧. (9). డ௧. . డ௙ೞ డ் డ௧. డ௧. (10). Wyrażenie w nawiasach okrągłych oznaczane jest w sformułowaniu pojemnościowym przez symbol c*, który oznacza zastępczą pojemność cieplną, stąd też wzór (10) przybiera postać: ∇ ∙ ∇ = ∗ (). డ் డ௧. (11). gdzie c* obliczane jest z odpowiednich wzorów aproksymacyjnych [10]. Powyższe równanie jest sformułowaniem pojemnościowym krzepnięcia. Sformułowanie pojemnościowe ma tę zaletę w stosunku do podstawowego entalpowego, że zmienną w tym równaniu jest temperatura, stąd nie zachodzi konieczność przechodzenia z rozwiązania równania w postaci entalpii do interesującej wartości w postaci temperatury. Do rozwiązania równania (11) wykorzystana została metoda elementów skończonych [13].. (3). 3. MODEL NARASTANIA FAZY STAŁEJ Przebieg krzepnięcia stopów metali analizowany może być za pomocą wykresów równowagi fazowej. Przykład takiego wykresu dla stopu dwuskładnikowego. (4). przedstawiony jest na rys. 1 [8]. Na wspomnianym wykresie na osi odciętych oznacza się stężenie domieszki C, natomiast na osi rzędnych temperaturę T. Na rys. 1 przedstawione są także przebiegi krzepnięcia stopu o składzie oznaczonym jako C0 dla różnych modeli krzepnięcia. Jak można zauważyć, różne modele zakładają różny stosunek fazy stałej. w którym to wzorze Tref jest temperaturą odniesienia. Wyrażenie (4) po zróżniczkowaniu względem t przybiera postać: డு. =. ∇ ∙ ∇ = − ௦ . Równanie w tej postaci jest równaniem nieliniowym, znacznie trudniejszym do rozwiązania od problemów liniowych. Uzyskanie rozwiązania metodami numerycznymi równania (3) wymagałoby stosowania metod iteracyjnych, które to oprócz wydłużenia czasu obliczeń powodują szereg problemów natury numerycznej. W celu uproszczenia równania (3) wykorzystuje się sformułowania entalpowe krzepnięcia. Przy wyprowadzeniu sformułowania entalpowego korzysta się z definicji entalpii w postaci: HT = ். డ௧. która w wyniku podstawienia do (6) daje:. gdzie L jest utajonym ciepłem krzepnięcia, a f udziałem odpowiedniej fazy w odlewie, indeks dolny s (lub l) oznacza, że dany parametr odnosi się odpowiednio do fazy stałej (lub ciekłej), w tym przypadku ρs oraz fs oznaczają odpowiednio gęstość fazy stałej oraz udział fazy stałej. Udział fazy stałej przyjmuje wartości z zakresu od 0 do 1. Sposób obliczania udziału fazy stałej jest ściśle uzależniony od wyboru modelu krzepnięcia i zostanie szerzej omówiony dalej. Podstawiając zależność (2) do równania (1), otrzymuje się: ∇ ∙ ∇ = . డு(்ሺ௫,௧ሻ). a sama formuła (7) przybiera postać:. (2). డ௧. =. Obliczenie ∂H/∂T jest możliwe na podstawie zależności (4). W tym przypadku, po zróżniczkowaniu definicji entalpii względem T otrzymuje się:. (1). W równaniu tym T oznacza temperaturę, c - ciepło właściwe, ρ - gęstość, t - czas, natomiast q jest strumieniem ciepła ze źródła wewnętrznego, które w przypadku krzepnięcia związane jest z przemianami fazowymi zachodzącymi w krzepnącym odlewie. W przypadku krzepnięcia przyjmuje się, że źródło ciepła ma postać: q = ௦ . (6). డ௧. Sformułowanie pojemnościowe krzepnięcia, które było wykorzystywane w tej pracy, wyprowadza się na podstawie podstawowego sformułowania entalpowego [9]. Przede wszystkim należy zauważyć, że pochodna entalpii względem czasu jest w istocie pochodną funkcji złożonej:. 2. MODEL KRZEPNIĘCIA. ∇ ∙ ∇ = . డு. (5). Prawa strona równania (5) jest również prawą stroną równania (1). Stąd też, w wyniku połączenia równania. 22.

(4) Robert Dyja, Andrzej Grosser i ciekłej w temperaturowym przedziale krzepnięcia znajdującym się pomiędzy liniami likwidus i solidus.. Rys. 1 Wykres równowagi układu z przemianą eutektyczną; przebieg krzepnięcia: 1 - model równowagowy, 2 - model Scheila, 3 - model pośredni. Rys. 2 Wykres stężenia domieszki w krzepnącym ziarnie dla modelu pośredniego. Model pośredni zakłada dyfuzję domieszki w fazie stałej i pełne rozprowadzenie domieszki w fazie ciekłej, ciekłej czego ilustracją jest rys. 2: t ξ, t ξ, (12) ξ, t t. Na wykresie oznaczone są temperatury charakterycharakter styczne dla rozważanego stopu, takie jak: • TM - temperatura krzepnięcia czystego czys składnika, • TL - temperatura likwidus, czyli temperatura, w której pojawiają się pierwsze kryształy krzepnącego stopu (niezależna od przyjętego modelu), zako • TS - temperatura solidus (temperatura zakończenia krzepnięcia) modelu równowagowego, • TE - temperatura eutektyczna, • TSE - temperatura zakończenia krzepnięcia dla modelu pośredniego. Zgodnie z tym wykresem, ten sam stop w wyniku stosowania różnych modeli krzepnięcia kończy krzepniękrzepni cie w różnych temperaturach. W przypadku modelu równowagowego (linia 1) temperatura końca krzepnięcia jest zależna tylko od składu chemicznego stopu. W przypadku modelu nierównowagowego (linia 2) temperatura końca krzepnięcia jest zawsze temperaturą eutektyczną. Natomiast model pośredni (linia 3) uzależuzale nia przebieg krzepnięcia od prędkości krzepnięcia odlewu (tj. od wielkości powstających ziaren). CharakterystyczCharakterystyc ną cechą modelu pośredniego jest to, iż model równowarównow gowy oraz model nierównowagowy są skrajnymi przyprz padkami tego modelu (odpowiednio: dla bardzo wolnego stygnięcia - model odel równowagowy, a dla bardzo szybkiego - model nierównowagowy). Z tego wynika również kolejna zaleta modelu pośredniego: uwzględnienie szybszy kości krzepnięcia w modelu i jej wpływ na powstającą mikrostrukturę odlewu. Istnieje możliwość uzyskania równań analitycznych ana opisujących udział fazy stałej dla poszczególnych modeli krzepnięcia. W tym celu należy zapisać równanie bilansu masy dla krzepnącego ziarna [2, 5, 6].. Równanie bilansu su masy dla pośredniego modelu nan rastania fazy stałej dla ziarna płaskiego ma postać [8]: ηt. . , . 1 . ೗ . ೥.

(5) . . . ೞ , . 0. (13). W równaniu tym CS i CL są odpowiednio stężeniem domieszki w fazie stałej oraz ciekłej, natomiast ξ jest bieżącą współrzędną.. Tak zapisane równanie r odpowiada sytuacji przedstawionej na rys. 2, na którym kt wyróżnić można trzy obszary: Ω1 jest obszarem, który już zaz krzepł, w obszarze Ω2 znajduje się tylko ciekły metal, natomiast Ω3 to obszar, w którym odbywa się krzepniękrzepni cie. Wielkość η(t) (t) wyznacza rozmiar zakrzepniętego zakrzepni materiału. u. Maksymalny rozmiar ziarna jest znormalizoznormaliz wany i wynosi 1. W tym przypadku pojawiają się w równaniu dodatkowe kowe parametry: rz jest grubością ziarna, natomiast Ds jest współczynnikiem dyfuzji domieszki w fazie stałej. Przekształcenia równania (13), przy wykorzystaniu zależności: k. ೞ. (14). ೗. oraz założenia, że narastanie fazy stałej odbywa się zgodnie z prawem em pierwiastkowym, pierwiastkowym prowadzą do uzyskania rozwiązania, którym jest:. బ. . ೖషభ. (15).

(6) మೖഀషభ. gdzie α to współczynnik czynnik Brody'ego-Flemmingsa: Brody'ego α. 23. ೞ ೑ ೥మ. (16).

(7) OBLICZENIA RÓWNOLEGŁE W SYMULACJI KRZEPNIĘCIA WYKORZYSTUJĄCEJ MODEL… Wielkość oznaczona jako tf jest lokalnym czasem krzepnięcia. Końcową grubość ziarna można obliczyć z zależności: ௭ = r௕ 1 − exp ሶ ିଵ. Stampede jest nazwą klastra komputerowego, który powstał z myślą o wykorzystaniu go przez naukowców w symulacjach intensywnych obliczeniowo. W listopadzie 2014 roku komputer ten znajdował się na 7. pozycji listy TOP500, będącej listą rankingową najpotężniejszych komputerów na świecie [16]. Stampede składa się z 6400 węzłów, z których każdy jest w istocie stacją obliczeniową, zawierającą dwa ośmiordzeniowe procesory Intel Xeon E5 oraz 32GB pamięci. Teoretycznie, pełna moc obliczeniowa Stampede to 102400 procesorów oraz ponad 200TB pamięci RAM.. (17). ். w której to rb jest maksymalną dopuszczalną grubością ziarna dla danego materiału, a szybkością stygnięcia. Na podstawie równania (15) oraz zależności wynikających z wykresu równowagi fazowej można uzyskać zależność udziału fazy stałej od temperatury dla modelu pośredniego:

(8) ௦ =. ଵ ଵିଶ௞ఈ. 1 − . ்ಾ ି் ்ಾ ି்ಽ. . భషమೖഀ ೖషభ. . (18). 5. PRZYKŁAD Implementację przetestowano na przykładzie odlewu, który posiadał skomplikowany kształt, jaki można spotkać często przy symulacjach obiektów technicznych. Plik geometrii pobrany został z repozytorium INRIA [14] i przedstawiony jest na rys. 3.. 4. OPIS IMPLEMENTACJI 4.1 IMPLEMENTACJA Przedstawione wcześniej wyprowadzenia zostały zaimplementowane we własnym programie, który korzystał z przygotowanyych wcześniej i dobrze przetestowanych bibliotek numerycznych. Istotną cechą implementacji była konieczność uwzględnienia w projekcie wykonywania symulacji na komputerach z pamięcią rozproszoną. Do implementacji komunikacji pomiędzy poszczególnymi zadaniami w pojedynczej symulacji wykorzystywany był standard MPI. Na komputerze Stampede, na którym wykonano symulacje, dostępna była implementacja MPI firmy Intel. Implementacja struktur do przechowywanie macierzy głównej, jak również metod rozwiązywania układów równań liniowych, została zaczerpnięta z biblioteki PETSc [1]. W programie wykorzystywana była metoda iteracyjna rozwiązywania układów równań liniowych GMRES [7] również zaimplementowana w bibliotece PETSc. Ułatwieniem w procesie tworzenia oprogramowania była również możliwość wykorzystania biblioteki TalyFEM, która zawiera implementacje zależności wykorzystywanych w metodzie elementów skończonych. Biblioteka ta była dotychczas używana głównie w symulacjach dotyczących technologii organicznych paneli słonecznych [3, 11].. 5.1 WŁASNOŚCI FIZYCZNE Dla odlewu prezentowanego na rys. 3 założono jego wykonanie ze stopu dwuskładnikowego Al z domieszką 2% Cu. Własności fizyczne stopu, niezbędne do przeprowadzenia symulacji, są zgromadzone w tabeli 1. Tabela 1 Własności fizyczne stopu Al2%Cu [8]. Wielkość Wsp. przewodzenia fazy stałej [W/mK]. Wartość 262. Wsp. przewodzenia fazy ciekłej [W/mK]. 104. Gęstość fazy stałej [kg/m3]. 2824. Gęstość fazy ciekłej [kg/m3]. 2498. Ciepło właściwe fazy stałej [J/kgK]. 1077. Ciepło właściwe fazy ciekłej [J/kgK]. 1275. Temperatura solidus [K]. 853. Temperatura likwidus [K]. 926. Temperatura eutektyczna [K]. 821. Temperatura krzepnięcia czystego składnika [K] Maks. dopuszczalna grubość ziarna [m]. 4.2 STAMPEDE Ze względu na swoją czasochłonność, zaawansowane obliczenia inżynierskie wymagają stosowania komputerów o dużych mocach obliczeniowych. Obecnie najpotężniejszymi komputerami są komputery złożone z wielu niezależnych jednostek, takie jak Stampede [15]. Komputery takie (nazywane również klastrami) wymuszają na wykorzystującym je oprogramowaniu stosowanie modelu programowania z przekazywaniem komunikatów [12].. 24. 933 10-4. Wsp. Brody-Flemminga [m2]. 4,6·10-8. Ciepło krzepnięcia [J/kgK]. 390000.

(9) Robert Dyja, Andrzej Grosser. Rys. 5 Średni czas w s potrzebny na rozwiązanie jednego kroku czasowego. Czas obliczeń pojedynczego kroku czasowego, czaso który jest średnią z wszystkich wykonanych kroków, kroków przedstawiony jest na rys. 5. Czas obliczeń jednego kroku czasoczas wego zawiera czas niezbędny na asemblację układu równań i na jego rozwiązanie. Na tym wykresie widać wyraźny, proporcjonalny do liczby procesorów p spadek czasu obliczeń wraz z dokładaniem kolejnych procesoproces rów. Spadek ten jest bardzo bliski teoretycznemu, idealidea nemu przypadkowi, w którym N-krotne N zwiększenie liczby procesorów powoduje N-krotne krotne zmniejszenie czasu obliczeń. Jeszcze wyraźniej widać dać to na rys. 6., na którym przedstawione jest przyspieszenie względne. Przyspieszenie względne liczone jest jako:. Rys. 3 Kształt odlewu. Na osiach skala w cm. Do symulacji oddawania ciepła przez odlew do otoot czenia wykorzystany został warunek brzegowy Newtona, w którym wartość współczynnika wnikania ciepła była równa 10 W/m2K, natomiast temperatura otoczenia wynosiła 400 K. Warunek ten został założony na pop wierzchni odlewu. W symulacji użyty został krok całkocałk wania po czasie równy 0,05 s. Przykładowy wykres rozkładu temperatur po upływie 30 s jest przedstawiony na rys. 4. Na tym rysunku maksymalna temperatura wynosi 960 K i jest oznaczona kolorem czerwonym, natomiast minimalna, równa 800 K, oznaczona jest kolorem niebieskim.. . ೝ ೛. (19). gdzie τr to czas referencyjny (tutaj jest to czas obliczeń oblicze realizowanych na 16 procesorach), natomiast τp to czas obliczeń, dla którego liczone jest przyspieszenie. W obydwu przypadkach koniec odliczania czasu wypawyp dał na moment, w którym ostatni procesor z grupy zakończył przewidziane obliczenia.. Rys. 4 Rozkład temperatury w odlewie po 30 s. 5.2 PRZETWARZANIE RÓWNOLEGŁE Do testowania wydajności przetwarzania równoległerównoległ go wykorzystano przedstawiony wcześniej przykład, w którym wykonano obliczenia dla 100 kroków czasoczas wych. Rozmiary siatek elementów skończonych, dla których wykonywane były obliczenia,, to: 3 074 501, 10 723 578, 16 739 748 oraz 24 596 008 czworościennych elementów skończonych. Z uwagi na fakt, że pojedynczy węzeł komputera Stampede zawiera 16 procesorów, taka została przyjęta najmniejsza liczba procesorów, na której uruchamiano no obliczenia. Maksymalna liczba procesorów wykorzystanych do obliczeń w przypadku największego zadania równa była 2048.. Rys. 6 Skalowalność. Szara linia wskazuje przyspieszenie idealne (teoretyczne). 6. PODSUMOWANIE W pracy zaprezentowano możliwości wykorzystania autorskiego programu w obliczeniach realizowanych na komputerach wielkiej mocy. Ponieważ komputery takie dysponują olbrzymią ią ilością pamięci operacyjnej, można było użyć siatek o rozmiarach wystarczających, a nawet większych niż były potrzebne do wiernego odwzorowania testowej geometrii.. 25.

(10) OBLICZENIA RÓWNOLEGŁE W SYMULACJI KRZEPNIĘCIA WYKORZYSTUJĄCEJ MODEL… Zaprezentowane zostały czasy potrzebne na rozwiązanie zadania. Nawet przy użyciu największej z testowanych siatek ten czas potrzebny do wykonania obliczeń jednego kroku czasowego jest niewielki, a przy użyciu większej liczby procesorów spada do wartości poniżej jednej sekundy. Prezentowane wyniki pokazują bardzo dobrą, niemal idealną, skalowalność przy używaniu aż do 1000 procesorów. Skalowalność obliczeń przy wykorzystaniu powyżej 1000 procesorów jest wciąż dobra, obserwowalne jest wyraźne przyspieszenie obliczeń. Prawdopodobnie. gdyby przygotować jeszcze większe siatki niż użyta w pracy siatka elementów skończonych o rozmiarze 24 milionów, można by było uzyskać lepsze przyspieszenia. Autorzy nie zdecydowali się na przeprowadzenie eksperymentów z większymi siatkami z powodu ograniczeń przyznanych zasobów oraz z faktu, że Stampede ogranicza maksymalny rozmiar zadania do 4096 procesorów. Pewnym problemem może być również rozmiar plików wyjściowych, który wynosi od około 150MB na jeden krok czasowy przy najmniejszej siatce do 1GB przy największej.. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources (via XSEDE allocation TG-CTS110007) that have contributed to the research results reported within this paper.. Literatura 1. 2. 3.. 4. 5. 6. 7. 8. 9. 10. 11.. 12. 13. 14. 15. 16.. Balay S. et al.: PETSc Users Manual, Argonne National Laboratory, 2014. Bokota A.: Symulacja numeryczna na poziomie makroskopowym i mikroskopowym krzepnięcia odlewu rozwiązanie MES i MEB. "Solidification of Metals and Alloys" 1994, 19, p. 31-40. Kodali H.K., Ganapathysubramanian B.: A computational framework to investigate charge transport in heterogeneous organic photovoltaic devices. „Computer Methods in Applied Mechanics and Engineering” 2012, 247, p. 113-129. Mendakiewicz J.: Identification of solidification process parameters. “Computer Assisted Mechanics and Engineering Sciences”, 2010, 17, p. 59-73. Mochnacki B., Suchy J.: Numerical methods in computations of foundry processes. Kraków: Wyd. PFTA, 1995. Parkitny R., Sczygiol N.: Równania krzepnięcia objętościowego dwuskładnikowych stopów metali. „Krzepnięcie Metali i Stopów” 1987, 12, p.29-44. Saad Y., Schultz M.H.: Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. „SIAM Journal on Scientific Computing” 1986, 3, Vol. 7, p. 856–869. Sczygiol N.: Modelowanie numeryczne zjawisk termomechanicznych w krzepnącym odlewie i formie odlewniczej, Częstochowa: Wyd. Pol. Częstoch., 2000. Sczygiol N., Szwarc G.: Application of enthalpy formulations for numerical simulation of castings solidification. „Computer Assisted Mechanics and Engineering Sciences” 2001, 8, p. 99-120. Voller V.R., Swaminathan C.R., Thomas B.G.: Fixed grid techniques for phase change problems: a review. „International Journal of Numerical Methods in Engineering” 1990, 30, p. 875-898. Wodo O., Ganapathysubramanian B.: Computationally efficient solution to the Cahn–Hilliard equation: adaptive implicit time schemes, mesh sensitivity analysis and the 3D isoperimetric problem. „Journal of Computational Physics” 2011, 230, p. 6037-6060. Wyrzykowski R.: Klastry komputerów PC i architektury wielordzeniowe: budowa i wykorzystanie. Warszawa: Akad. Ofic. Wyd. EXIT, 2009. Zienkiewicz O.C., Taylor R.L.: The finite element method. Vol. 1 – The Basis. Oxford: Butterworth-Heinemann, 2000. INRIA GAMMA Group Database: www.rocq.inria.fr/gamma/gamma/download/download.php Stampede - Texas Advanced Computing Center: www.tacc.utexas.edu/stampede/ Top 500 Supercomputer Sites: www.top500.org/system/177931. 26.

(11)