Iteracyjny algorytm rozwiązywania równań liniowych w Abaqus/Standard

Iteracyjny algorytm rozwiązywania równań liniowych w Abaqus/Standard

Iteracyjne solwery układów równań liniowych odgrywają kluczową rolę w obliczeniach inżynierskich, zwłaszcza w kontekście dużych i złożonych modeli. Spotykamu je najczęściej w solwerach CFD (Computational Fluid Dynamics), są jednak stosowane coraz częściej również w obliczeniach strukturalnych. W obu przypadkach dyskretyzacja ciągłych równań różniczkowych prowadzi do układów liniowych, które często zawierają miliony niewiadomych. Metody iteracyjne, w przeciwieństwie do bezpośrednich, pozwalają na przybliżone rozwiązanie tych układów przy znacznie mniejszym zużyciu pamięci i czasu obliczeń, co jest szczególnie istotne w przypadku problemów wielkoskalowych. W tym artykule skupimy się na charakterystyce solwera iteracyjnego dostępnego w Abaqus/Standard oraz jego zastosowaniu w praktycznych obliczeniach inżynierskich. W wielu przypadkach iteracyjny algorytm rozwiązywania równań liniowych w Abaqus/Standard jest wydajną alternatywą rozwiązywania układów równań liniowych, szczególnie przydatną w przypadku dużych modeli, o gęstej macierzy sztywności.

Jak działa iteracyjny algorytm rozwiązywania równań liniowych

W analizie nieliniowej Abaqus/Standard dzieli krok symulacji na szereg przyrostów czasu i znajduje przybliżoną konfigurację równowagi na końcu każdego przyrostu czasu. W ramach przyrostu czasu Abaqus/Standard wykorzystuje metodę Newtona, wykonując iteracje w celu znalezienia akceptowalnego rozwiązania. W każdej iteracji Newtona Abaqus/Standard używa sztywności stycznej konstrukcji, K, i wartości resztowej siły, f, do obliczenia korekcji przemieszczenia c. Odbywa się to poprzez rozwiązanie układu równań liniowych (Ku = f) za pomocą metody bezpośredniej LDU lub za pomocą solwera iteracyjnego. Metody bezpośrednie, takie jak LDU, rozwiązują układy równań w określonej liczbie kroków, zapewniając dokładne rozwiązanie z dokładnością do błędów zaokrągleń. Natomiast metody iteracyjne, stosowane do rozwiązywania układów równań liniowych, iteracyjnie dążą do uzyskania rozwiązania przybliżonego. Polegają one na iteracyjnym doskonaleniu rozwiązania początkowego, aż do momentu osiągnięcia założonej dokładności. W Abaqus/Standard stosuje się kilka iteracyjnych algorytów rozwiązywania równań liniowych:

• SQMR (Symmetric Quasi-Minimal Residual)

• FGMRES (Flexible Generalized Minimal RESidual),

które są metodami przestrzeni Kryłowa.

Włączenie solwera iteracyjnego w Abaqus

Aby włączyć solwer iteracyjnego w obliczeniach Abaqus, należy dodać w pliku wejściowym Abaqus (.inp), w definicji kroku analizy, w której chcesz użyć solwera iteracyjnego, paramter SOLVER=ITERATIVE słowa kluczowego *STEP

W Abaqus/CAE należy przejść do modułu Step, wybrać krok analizy, w którym chcesz użyć solwera iteracyjnego, i w w sekcji Other okna dialogowego Create/Edit Step wybrać Method -> Iterative:

Kiedy można i warto stosować solwer iteracyjny

Wybór między iteracyjnym a bezpośrednim solverem równań liniowych w Abaqus/Standard zależy od wielu czynników, takich jak: rozmiar i struktura modelu, typ analizy, oczekiwana dokładność i dostępne zasoby obliczeniowe. Iteracyjny solver jest szczególnie zalecany w przypadku modeli o wysokiej spójności siatki, takich jak jednolite struktury blokowe, czyli przypominają bardziej sześcian niż płytę, lub złożenia z niewielką liczbą części. W przypadku modeli o wielkości przekraczającej 10 milionów stopni swobody (DOF) i wymagających ponad 1.0E+12 operacji zmiennoprzecinkowych (FLOP) na pojedynczy przebieg solwera, solwer iteracyjny powinien zapewnić lepszą wydajność niż solwer bezpośredni. Solwer iteracyjny ma też znacznie mniejsze zużycie pamięci niż w przypadku solwera bezpośredniego a ponadto koszt rozwiązania rośnie w jego przypadku proporcjonalnie do rozmiaru modelu.

Solwera iteracyjnego należy unikać w przypadku modeli, które są źle uwarunkowane, ponieważ mogą powodować problemy ze zbieżnością. Są to np. modele z dużą liczbą elementów o niekorzystnym współczynniku kształtu, np. z bardzo dużymi lub małymi kątami między krawędziami. Źle uwarunkowane modele często generują ostrzeżenia o ujemnych wartościach własnych, zerowych pivotach lub osobliwościach numerycznych podczas analizy z użyciem solwera bezpośredniego. Podobnie, w przypadku modeli, które mają duże silnie nieliniowe a tym bardziej nieciągłe właściwości materiałowych, mogą zbiegać się wolno lub wcale. Również w przypadku specyficznych elementów lub obciążeń, jak np. elementy kohezyjnych zastosowanie solwera iteracyjnego prawdopodobnie doprowadzi będzie do braku zbieżności. Tak samo analizy, w których mamy do czynienia z wieloma kontaktami zdefiniowanymi za pomocą metody bezpośredniej Lagrange'a możemy oczekiwać problemów ze zbieżnością. W przypadku solwera iteracyjnego zalecane jest stosowanie metod wymuszania więzów kontaktowych takich jak metoda funkcji kary (*SURFACE BEHAVIOR, PENALTY) lub metody rozszerzonego Lagrange'a (*SURFACE BEHAVIOR, AUGMENTED LAGRANGE). W przypadku analizy kontaktowej zaleca się przetestowanie modeli z kontaktem w analizie perturbacyjnej przed uruchomieniem pełnej, docelowej analizy nieliniowej, aby sprawdzić wydajność solwera iteracyjnego a tym samym ocenić korzyści jego zastosowania. Trzeba bowiem pamiętać, że sam algorym Newtona-Raphsone będzie zbiegał się co do zasady wolnie w przypadku stosowania solwera iteracyjnego. W przypadku problemów z uzyskaniem zbieżności, zawsze zaleca się przetestowanie modelu za pomocą solwera bezpośredniego, aby ustalić, czy problem ten jest spowodowany użyciem solwera iteracyjnego.

Solwer iteracyjny nie może być stosowany w analizach geostatycznych i analizach dyfuzji płynu w porach dla modeli z więzami Lagrange'a. Nie jest również obsługiwany w analizach termomechanicznych i termoelektrycznych, a także w analizach transportu w stanie ustalonym.

Co zrobić gdy iteracyjny algorytm rozwiązywania równań liniowych nie zbiega się

Brak zbieżności solwera iteracyjnego lub algorytmu Newtona-Raphsona w Abaqus może być czasami trudnym problemem do rozwiązania. Istnieje jednak kilka potencjalnych przyczyn braku zbieżności i sposobów rozwiązania tego problemu. Jedną z nich mogą być duże różnice we wartościach właściwościach materiałowych lub materiały o charakterystyce nieściśliwej (współczynniku Poissona bliskim 0,5) czy silnie zniekształcone elementy. Problemy ze zbieżnością mogą również wynikać z interakcji kontaktowych.

Błąd rozwiązania przybliżonego jest mierzony za pomocą względnej wartości resztowej układu równań liniowych, zdefiniowanej jako || Ku – f || / || f ||. Zbieżność solwera iteracyjnego zostaje osiągnięta, gdy ta względna wartość resztowa spadnie poniżej określonej tolerancji (domyślnie tolerancja wynosi 1.0E-06) w określonej liczbie iteracji (domyślnie maksymalna liczba iteracji to 500). Możliwa jest zmiana tolerancji i maksymalnej liczby iteracji. W przypadku problemów ze zbieżnością algorytmu Newtona-Raphsona skuteczne może być zmniejszenie tolerancji solwera iteracyjnego (np. do 1.0E-08) w celu poprawy dokładności. Celowe będzie również zwiększenie maksymalna liczba iteracji do 1000. Tolerancję zbieżności i maksymalną liczbę iteracji solwera iteracyjnego można zmienić za pomocą słów kluczowych w pliku wejściowym dla danego zadania:

Można również ustawić je globalnie, za pomocą zmiennych środowiskowych - dla tolerancji zbieżności solwera iteracyjnego jest to zmienna ABA_EQSITR_RTOL:

Dla maksymalnej liczby iteracji solwera iteracyjnego jest to zmienna ABA_EQSITR_MAX_ITER:

Dla maksymalnej liczby iteracji solwera iteracyjnego jest to zmienna ABA_EQSITR_MAX_ITER:

Należy pamiętać, że modyfikacja zmiennych środowiskowych Abaqus, zwłaszcza tych wpływających na działanie solwera, wymaga dokładnego zapoznanie się z dokumentacją Abaqus. Niewłaściwe ustawienia mogą prowadzić do niestabilności obliczeń, błędów lub nieprawidłowych wyników.

Co to jest prekondycjoner

Algorytm iteracyjnego rozwiązywania równań liniowych w Abaqus/Standard wykorzystuje prekondycjoner do przyspieszenia zbieżności rozwiązania. Prekondycjoner odgrywa kluczową rolę w metodach przestrzeni Kryłowa, przyspieszając zbieżność rozwiązania układu równań. Szybkość zbieżności metody Kryłowa jest silnie uzależniona od tzw. uwarunkowania układu równań. Układy dobrze uwarunkowane zbiegają szybciej niż układy źle uwarunkowane. Prekondycjoner przekształca oryginalny układ równań w równoważny układ o lepszym uwarunkowaniu. W uproszczeniu, prekondycjoner P jest macierzą, która przybliża macierz odwrotną do macierzy K. Zamiast rozwiązywać oryginalny układ równań Ku = f, metoda Kryłowa z prekondycjonerem rozwiązuje układ przekształcony:

P⁻¹Ku = P⁻¹f.

Jeśli prekondycjoner P jest dobrym przybliżeniem macierzy K , to macierz P⁻¹K będzie miała lepsze uwarunkowanie niż macierz K. Prekondycjoner jest obliczany tylko raz, na początku procesu rozwiązywania układu równań liniowych.

W Abaqus/Standard jako prekondycjoner stosowana jest algebraicznych metod wielosiatkowych (AMG), której kluczowym elementem jest funkcja wygładzająca, zwana w języku angielskim "smoother". Metoda AMG działa na zasadzie rozwiązywania układu równań na różnych poziomach "siatek", gdzie siatka o rzadszej strukturze reprezentuje aproksymację problemu na poziomie globalnym, a siatki o drobniejszej strukturze reprezentują szczegóły lokalne. Rola funkcji wygładzającej polega na redukcji tzw. oscylacji o wysokiej częstotliwości w błędzie rozwiązania na danej siatce. W uproszczeniu, funkcja wygładzająca działa jak filtr, który wygładza szybkie zmiany w rozwiązaniu, pozostawiając wolnozmienne składowe. Jest to ważne, ponieważ metody iteracyjne, takie jak metody przestrzeni Kryłowa, są bardziej efektywne w redukcji błędu o niskiej częstotliwości.

Abaqus/Standard oferuje trzy typy funkcji wygładzających dla obu rodzajów problemów, z symetryczną i niesymetryczną macierzą sztywności: "weak_sparse" (domyślna), "jacobi_chebyshev" i "strong_sparse", które mogą przynieść lepszą wydajność lub stabilność dla określonych typów modeli. Zmiana typu funkcji wygładzającej możliwa jest za pomocą zmiennej ABA_PAMG_SMOOTHER_TYPE:

Typ funkcji wygładzającej (smoother) w prekondycjonerze AMG może znacząco wpłynąć tak na wydajność jak i na stabilność solwera iteracyjnego. Domyślny typ smoothera (tzw. "weak_sparse") jest zazwyczaj dobrym wyborem, ale w niektórych przypadkach warto eksperymentować z pozostałymi typami, aby znaleźć optymalne ustawienie dla danego modelu.

Szybko, szybciej, na wielu rdzeniach najszybciej

Wszystkie etapy iteracyjnego procesu rozwiązywania, w tym agregacja macierzy, konfiguracja prekondycjonera i samo rozwiązywanie metodą Kryłowa, są obsługiwane w trybie obliczeń równoległych. Abaqus/Standard oferuje trzy główne tryby zrównoleglania obliczeń: SMP (Shared Memory Parallelization), DMP (Distributed Memory Parallelization) oraz tryb hybrydowy. SMP wykorzystuje wielowątkowość na pojedynczym komputerze wielordzeniowym, gdzie wszystkie wątki współdzielą dostęp do tej samej pamięci RAM. Ten tryb jest domyślny, gdy praca jest przesyłana do jednego hosta. DMP dzieli model na mniejsze domeny i uruchamia każdy proces MPI (Message Passing Interface) na osobnym rdzeniu, z dedykowaną pulą pamięci. W tym trybie komunikacja między rdzeniami odbywa się poprzez sieć, co może generować znaczne jej obciążenie a jej wydajność znacząco wpływać na czas obliczeń. Tryb hybrydowy łączy zalety obu podejść, uruchamiając wiele procesów MPI, z których każdy wykorzystuje wielowątkowość na przydzielonych mu rdzeniach. Pozwala to na zrównoważenie wydajności i zużycia pamięci.  Na przykład, polecenie:

uruchomi 4 procesy MPI, z których każdy będzie wykorzystywał 4 wątki. Tryb DMP oferuje najlepszą wydajność dla solwera iteracyjnego, ale może prowadzić do znacznego wzrostu zużycia pamięci. Tryb hybrydowy z 4 wątkami na proces MPI jest rekomendowany w przypadku ograniczeń pamięci, ponieważ redukuje ilość procesów MPI i wykorzystuje wielowątkowość w obrębie każdego procesu, co zmniejsza zużycie pamięci i pozwala na optymalne wykorzystanie zasobów obliczeniowych.

Nie zawsze domyślne ustawienia gwarantują maksymalną wydajności obliczeń. Dzięki kilku  zmiennym środowiskowych, które można dodać do pliku konfiguracyjnego Abaqus (abaqus_v6.env), można w określonych przypadkach zwiększyć wydajność obliczeń jak również monitorować zużycie zasobów. Domyślnie Abaqus wykorzystuje 80% dostępnej pamięci fizycznej. Można użyć zmiennej:

która, zwiększa ten limit do 90%, co może być korzystne w przypadku dużych modeli. Podobnie włączenie równoległego przetwarzania wstępnego może przyspieszyć analizę, szczególnie w przypadku ponownego uruchamiania analizy z pliku restartu (*RESTART, WRITE). Pomaga to uniknąć ograniczeń rozmiaru globalnej struktury bazy danych dla dużych modeli:

W celu monitorowanie zużycia zasobów można użyć zmiennej ABA_RESOURCE_MONITOR:

Włącza ona generowanie plików użycia (.use), które zawierają informacje o zużyciu czasu procesora i pamięci przez Abaqus.

Benchmark - porównanie solwera iteracyjnego z solwerem bezpośrednim

W benchmarku porównano wydajność solwera iteracyjnego i bezpośredniego (direct) w Abaqus. Szczególną uwagę poświęcono wymaganiom dotyczącym pamięci oraz czasu przetwarzania. Analizowane zadanie to benchmark wydajnościowy ABAQUS/Standard S4, symulujący przykręcanie głowicy cylindrów do bloku silnika. Model ten składał się z 23 608 804 elementów, 5 750 962 węzłów oraz 16 776 241 zmiennych. Do modelowania głowicy i bloku zastosowano elementy czworościenne typu C3D4, do śrub elementy C3D8I, a uszczelkę zamodelowano za pomocą dedykowanych elementów GK3D8. Materiały dla bloku, głowicy i śrub zdefiniowano jako liniowo sprężyste, podczas gdy uszczelka miała nieliniową charakterystykę z uwzględnieniem plastyczności. Nieliniowość problemu wynikała z warunków kontaktu oraz uplastyczniania się materiału uszczelki podczas dokręcania śrub. Wszystkie zadania uruchomiono na 16 rdzeniach (8 rdzeni na jeden węzeł NUMA).

Porównanie różnych konfiguracji uruchomienia solwerów wykazało istotne różnice w wymaganiach pamięciowych i czasie obliczeń. Solwer bezpośredni (direct) uruchomiony domyślnie na 16 rdzeniem z wykorzystaniem 16 wątków miał największe wymagania co do ilości pamięci wymaganej do minimalizacji operacje wejścia/wyjścia (Memory Min I/O: 218 825 MB), co pozwala na uzyskanie maksymalnej wydajności. W tej konfiguracji czas obliczeń wyniósł 31 647 s, co stanowiło punkt odniesienia (znormalizowany czas = 1).

W przypadku solwera iteracyjnego, Abaqus domyślnie (odpowiada za to specjalny moduł JLE - ) uruchamiany jest w trybie hybrydowym (MPI+wątki) z użyciem PlatformMPI, w tym przypadku w konfiguracja z 4 rankami MPI wymagała znacznie mniej pamięci wymaganej do minimalizacji operacje wejścia/wyjścia (Memory Min I/O: 108 470 MB). Solwer iteracyjny pozwolił skrócić czas obliczeń do 11 310 s (ok. 35,7% czasu bazowego). Jeszcze lepsze wyniki uzyskano przy zastosowaniu solwera iteracyjnego uruchominionego w trybie hybrydowym jednak z 2 rankami MPI z użyciem OpenMPI oraz nie domyślną koligacją (affinity) do węzłów NUMA, gdzie czas obliczeń wyniósł jedynie 9 210 s (29,1% czasu bazowego), przy minimalnie niższych wymaganiach pamięci na poziomie 108 171 MB.

Podsumowanie

Iteracyjny algorytm rozwiązywania równań liniowych w Abaqus/Standard to alternatywa dla domyślnego solwera bezpośredniego, szczególnie przydatna w przypadku dużych modeli o wysokiej spójności siatki, takich jak "blokowe" struktury, gdzie solwer bezpośredni wymagałby nadmiernej liczby operacji zmiennoprzecinkowych. Algorytm ten, może znacznie skrócić czas obliczeń i obniżyć zużycie pamięci w porównaniu do solwera bezpośredniego. Należy jednak pamiętać, że jego wydajność i zbieżność zależą od specyfiki modelu, takich jak typy elementów, definicje kontaktu, nieliniowości materiałowe i geometryczne, a w przypadku modeli źle uwarunkowanych, solwer iteracyjny może zbiegać bardzo wolno lub w ogóle nie osiągnąć zbieżności.