Algorytm 9: Algorytm jednorodnej p adaptacji wykonywany apropri
Algorytm 9: Algorytm jednorodnej p adaptacji wykonywany apropri
Algorytm p-adaptacji w ogólnej postaci zapisać można w sposób następujący:
1. jednorodna p adaptacja wykonywana apropri (siatka początkowa, dokładność, maksymalna liczba kroków) 2. aktualna siatka = siatka początkowa
3. rozwiąż swój problem obliczeniowy na aktualnej siatce (na przykład projekcję bitmapy lub transport ciepła), obliczając aktualne rozwiązanie
4. pętla dla i = 1 do maksymalnej liczby kroków 5.
6. dla każdego elementu aktualnej siatki wykonaj następujące czynności
7. oblicz bład aproksymacji na elemencie obliczając seminormę gradientu rozwiązania na elemencie , oraz normy z pochodnych kierunkowych oraz
8. jeśli wtedy
9. jeśli wtedy
10. jeśli wtedy
11. koniec pętli po elementach
12. jeśli wymagana dokładność, to koniec 13. pętla po elementach siatki
14. jeśli wówczas podnieś stopień aproksymacji wielomianowej w dwóch kierunkach, kontynuuj
pętlę po elementach
15. jeśli wówczas podnieś stopień aproksymacji wielomianowej w kierunku osi , kontynuuj pętlę po elementach
16. jeśli wówczas podnieś stopień aproksymacji wielomianowej w kierunku osi , kontynuuj pętlę po elementach
17. koniec pętli po elementach 18. konieć pętli adaptacji
Jeśli dany element wymaga modyfikcji stopnia wielomianów w celu polepszenia jakości aproksymacji, możliwy jest jeden z trzech opisanych scenariuszy, podniesenie stopnia funkcji bazowych w kierunku osi , podniesienie stopnia funkcji bazowych w kierunku osi , lub podniesnie stopnia w obydwu kierunkach.
W algorytmie adaptacyjnym zazwyczaj nie wykonuje się adaptacji na wszystkich elementach siatki, ale jedynie na tych elementach na których błąd jest większy niż 33 procent błedu maksymalnego (liczonego po wszystkich elementach). Jest to oczywiście arbitralnie przyjęta wartość. W zaproponowanej wersji algorytmu próbujemy najpierw wykonać adaptację w obydwu kierunkach (podnieść stopień wielomianu o 1 w obydwu kierunkach), a następnie (jeśli nie jest to wskazane) próbujemy podnieść stopień wielomianów w jednym z dwóch kierunków. Jest to oczywiście pzykładowa wersja algorytmu, możliwe są różne modyfikacje. Najbardziej zaawansowaną wersją algorytmu adaptacyjnego jest algorytm hp adaptacji opisany w innym module. Scalając ze sobą wiele patchów (grup elementów) na których rozpięte są różne wielomiany B-spline, uzyskamy siatkę obliczeniową na której możliwe będą lokalne modyfikacje stopnia funkcji B-spline, na przykład stosując opisany algorytm p-adaptacji.
Gdy patch elementowy składa się z jednego elementu, funkcje bazowe rozpięte na tym elemencie odpowiadają klasycznemu algorytmowi p-adaptacji pracującemu na wielomianach Lagrange'a. Dzieje się tak dlatego, iż wektory węzłów w których każdy element oddzielony jest separatem, tak jak ma to miejsce w klasycznych wielomianach Lagrange'a. Obie bazy generują wówczas identyczną przestrzeń wielomianów.
Matematyczne właściwości algorytmu p-adaptacji opisał po raz pierwszy prof. Ivo Babuśka z Uniwersytetu Teksańskiego w Austin [35]
Algorytm hp-adaptacji
Algorytm hp-adaptacji
Algorytm automatycznej hp adaptacji został po raz pierwszy szczegółowo opisany w książkach prof. Leszka Demkowicza, matematyka polskiego pochodzenia, pracującego na Uniwersytecie Teksańskim w Austin. Analizę właściwości matematycznych algorytmu hp-adaptacji dokonał prof. Ivo Babuška, amerykanin czeskiego pochodzenia, również pracujący na Uniwersytecie Teksańskim w Austin [6] [7] [36] [37].
ALGORYTM
Algorytm 10: Algorytm automatycznej hp adaptacji
Algorytm 10: Algorytm automatycznej hp adaptacji
1. Algorytm automatycznej hp adaptacji (siatka początkowa, wymagana dokładność, maksymalna liczba iteracji) 2. siatka rzadka = siatka początkowa
3. pętla i = 1 do maksymalnej liczby iteracji
4. Rozwiąż problem obliczeniowy na aktualnej siatce rzadkiej (na przykład problem projekcji bitmapy lub problem transportu ciepła) uzyskując rozwiązanie
5. siatka gęsta = siatka rzadka
6. Połam każdy element siatki rzadkiej na 4 elementy (w przypadku elementów prostokątnych w dwóch wymiarach) lub 8 elementów (w przypadku elementów sześciennych w trzech wymiarach) (możliwe również łamania elementów trójkątnych w dwóch wymiarach i elementów czworościennych, pryzmatycznych i piramid w trzech wymiarach) oraz zwiększ stopień wielomianów o jeden w kążdym kierunku na każdym elemencie (co jest naturalne wprzypadku elementów prostokątnych i sześciennych, natomiast wymaga doprecyzowania w przypadku elementów innego rodzaju)
7. Rozwiąż problem obliczeniowy na aktualnej siatce gęstej uzyskując rozwiązanie 8. Bład maksymalny
9. Pętla po elementach siatki rzadkiej
10. Dla każdego elementu siatki rzadkiej oszacuj błąd względny (normę różnicy pomiędzy rozwiązaniem na siatce rzadkiej a gęstej)
11. Jeśli wówczas
12. Koniec pętli po elementach
13. Jeśli wymaganej dokładności, wówczas koniec. 14. Pętla po elementach siatki rzadkiej
15. Jeśli wówczas wybierz optymalny sposób adaptacji elementu z siatki rzadkiej i zaaplikuj go do elementu siatki rzadkiej
16. Koniec pętli po elementach 17. Koniec iteracji
Rysunek 28: Siatka początkowa, wszystkie elementy oddzielone są C0 separatorami, na każdym elemencie siatki rozpietęto funkcje bazowe powstałe przez iloczyn tensorowy wielomianów kwadratowych w kierunku poziomym oraz pionowym.
Rysunek 29: Rozwiązanie (pole skalarne temperatury) na siatce rzadkiej.
Rysunek 30: Siatka gęsta powstała przez połamanie każdego elementu siatki rzadkiej na cztery elementy oraz podniesnie stopnia wielomianów w kierunku poziomym i pionowym o jeden, na każdym elemencie.
Rysunek 31: Rozwiązanie (pole skalarne temperatury) na siatce gęstej.
Rysunek 32: Wybór optymalnej strategii adaptacji dla każdego elementu siatki rzadkiej, poprzez porównanie rozwiązań na siatkach żadkiej i gęstej. .
Algorytm hp-adaptacji zilustrowany jest na rysunkach Rys. 28 - Rys. 32.
Konwencja przyjęta podczas rysowania siatek obliczeniowych jest następująca. Kolory na rysunkach Rys. 28, Rys. 30 i Rys. 32 oznaczają stopnie wielomianów rozpiętych na krawędziach i wnętrzach elementów. Na każdej krawędzi ustalony jest jeden stopień wielomianu, natomiast na każdym wnętrzu elementu ustalone są stopień w kierunku poziomym oraz stopień w kierunku pionowym.
Na przykład, na rysunku Rys. 28 wszystkie wielomiany na wszystkich krawędziach mają stopień 2, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 2 w kierunku pionowym i poziomym. Na rysunku Rys. 30 wszystkie wielomiany na wszystkich krawędziach mają stopień 3, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 3 w kierunku pionowym i poziomym. Z koleji na rysunku Rys. 32 przeglądając elementy wierszami od lewej do prawej strony, w pierwszym wierszu na pierwszym elemencie wszystkie wielomiany na krawędziach mają stopień 2, z wyjątkiem wielomianu na dolnej krawędzi elementu, który ma stopień 1. Na pierwszym elemencie wielomiany we wnętrzu mają stopień 2 w kierunku pionowym i poziomym. Na drugim elemencie, na krawędziach lewej, górnej i dolnej, wielomiany mają stopień 2, natomiast we wnętrzu stopień 3 w kierunku pionowym i stopień 2 w kierunku poziomym. Na trzecim i czwartym elemencie wszystkie wielomiany na krawędziach mają stopień 3. Podobnie wielomiany we wnętrzach w kierunkach pionowym i poziomym. W drugim wierszu elementów, pierwszy element (czyli piąty globalnie) posiada wielomiany stopnia 3 na lewej i prawej krawędzi, wielomian stopnia 1 na górnej i dolnej krawędzi. We wnętrzu element ten posiada wielomiany stopnia 1 w kierunku poziomym i wielomiany stopnia 3 w kierunku pionowym. Drugi element w drugim wierszu posiada stopień wielomianu 2 w kierunku poziomym oraz stopień wielomianu 3 w kierunku pionowym. Na trzecim i czwartym elemncie wszystkie wielomiany na wszystkich krawędziach mają stopień 3, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 3 w kierunku pionowym i poziomym. Pozostałe elementy w trzecim i czwartym wierszu mają wielomiany rozłożone symetrycznie w stosunku do wielomianów w pierwszym i drugim wierszu na pierwszym i drugim elemencie.
W jaki sposób dokonywany jest wybór optymalnego sposobu adaptacji pojedynczego elementu? Zauważmy że w przypadku algorytmu hp adaptacji mamy wiele możliwości modyfikacji pojedynczego elementu:
1. Pozostawienie elementu bez zmian
2. Możliwe jest połamanie elementu na dwa nowe elementy w kierunku pionowym 3. Możliwe jest połamanie elementu na dwa nowe elementy w kierunku poziomym
4. Możliwe jest połamanie elementu na cztery nowe elementy (dwa w kierunku poziomym i dwa w kierunku pionowym) Ponadto dla każdego połamanego elementu możliwa jest modyfikacja stopnia wielomianów we wnętrzu elementu
1. Pozostawienie stopni elementu bez zmian
2. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku poziomym 3. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku pionowym
4. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku poziomym i pionowym
Stopnie na krawędziach zmodyfikowanych elementów ustalane są na podstawie reguły minimum. Innymi słowy stopień wielomianu na krawędzi dzielonej przez dwa sąsiadujące elementy ustalany jest jako minimum stopni we wnętrzu elementów w odpowiednim kierunku.
Dla krawędzi pionowej otoczonej przez dwa elementy oraz ustalamy stopień .
Dla krawędzi poziomej otoczonej przez dwa elementy oraz ustalamy stopień .
Mamy więc bardzo dużo możliwości adaptacji pojedynczego elementu (szacując na podstawie wymienionych możliwych rodzajów adaptacji mamy 4 (modyfikacje stopnia elementu bez łamania elementu) + 4*4 (połamanie elementu na dwa elementy w kierunku poziomym i po cztery możliwości modyfikacji stopnia każdego elementu) + 4*4 (połamanie elementu na dwa elementy w kierunku pionowym i po cztery możliwości modyfikacji stopnia każdego elementu) + 4*4*4*4 (połamanie elementu na cztery elementy i po cztery możliwości modyfikacji stopnia każdego elementu). W sumie 4+4*4+4*4+4*4*4*4=292 możliwości. W jaki sposób podejmowana jest decyzja o wyborze rodzaju adaptacji pojedynczego elementu?
ALGORYTM
Algorytm 11: Algorytm wyboru optymalnej strategii adaptacji elementu
Algorytm 11: Algorytm wyboru optymalnej strategii adaptacji elementu
1. Wybór optymalnej strategii adaptacji elementu (rozwiązanie na siatce żadkiej zawężone do elementu , rozwiązanie na siatce gęstej zawężone do elementu )
2. Pętla po możliwych rodzajach adaptacji elementu od 1 do 292 3. Największa prędkość spadku błędu = 0
4. Optymalna adaptacja = 0
5. Zrób kopię , elementu , i wykonaj na nim rozważaną adaptację
6. Oblicz projekcję rozwiązania na siatce gęstej na adaptowany element
7. Oblicz o ile spadnie błąd względny jeśli wykonamy adaptacje reprezentowaną przez kod, spadek błędu(kod)= 8. Oblicz ile niewiadomych (ile funkcji bazowych) musimy dodać żeby zrealizować adaptację reprezentowaną przez kod,
koszt adaptacji (kod)
9. Oblicz i zapamiętaj prędkość spadku błędu (kod) = spadek błędu (kod) / koszt adaptacji(kod)
10. Jeśli prędkość spadku błędu (kod) jest większa niż największa prędkość spadku błędu, to zapamiętaj największa wartość spadku błędu = pedkość spadku błędu (kod), optymalna adaptacja=kod
11. Koniec pętli po możliwych rodzajach adaptacji
12. Wykonaj na elemencie znalezioną optymalną adaptację
Innymi słowy, wybieramy taki rodzaj adaptacji dla elementu, który pozwala uzyskać największy spadek błędu najmniejszym kosztem. Wielkość ta reprezentowana jest przez prędkość spadku błędu, która rośnie wraz ze spadkiem błędu, ale maleje wraz z poniesionym kosztym (z liczbą funkcji dodaną do elemetu, ponieważ koszt obliczenia na danym elemencie zależy od liczby niewiadomych, które są współczynnikami funkcji bazowych). Algorytm ten zilustrowany jest na rysunku Rys. 33
Rysunek 33: Wybór optymalnej strategii adaptacji elementu spośród wielu możliwości, na podstawie obliczonej prędkości spadku błędu.
Dla przykładowego problemu transportu ciepła na obszarze w kształcie litery L możliwe jest uzyskanie siatki hp adaptacyjnej która pozwala na rozwiązanie problemu z błędem względnym 0.0001. Siatka taka przedstawiona jest na rysunkach Rys. 34 - Rys. 55
Rysunek 34: Siatka hp adaptacyjna umożliwiająca rozwiązanie problemu transportu ciepła na obszarze w kształcie litery L z dokładnością 0.001 w sensie błędu względnego (różnicy w normie H1 pomiędzy rozwiązaniem na siatce żadkiej a siatce gęstej)
Rysunek 35: Zoom 10 na siatkę optymalną
Rysunek 36: Zoom 100 razy na siatkę optymalną.
Rysunek 37: Zoom 1000 razy na siatkę optymalną.
Rysunek 38: Zoom 10000 razy na siatkę optymalną.
Rysunek 39: Zoom 100000 razy na siatkę optymalną.
Rysunek 40: Zoom 1000000 razy na siatkę optymalną.
Algorytm automatycznej hp adaptacji jest więc bardzo skomplikowany i trudny do zaimplementowania. Dlaczego więc wart jest on naszego zainteresowania? Rozważmy możliwe różne algorytmy adaptacyjne dla przykładowego problemu transportu ciepła na obszarze w kształcie litery L, czyli
1. Algorytm jednorodnej h adaptacji, w którym równomiernie na całym obszarze zwiększana jest liczba elementów 2. Algorytm jednorodnej p adaptacji, używający siatki zbudowanej z 3 elementów i zwiększający równomiernie stopień
wielomianów we wnętrzach wszystkich elementów w kierunku poziomym oraz pionowym, oraz stopień wielomianów na wszystkich krawędziach
3. Algorytm jednorodnej p adaptacji, używający siatki zbudowanej z 12 elementów i zwiększający równomiernie stopień wielomianów we wnętrzach wszystkich elementów w kierunku poziomym oraz pionowym, oraz stopień wielomianów na wszystkich krawędziach
4. Algorytm automatycznej p adaptacji, zwiększający stopień wielomianów we wnętrzach wybranych elementów w wybranych kierunkach, i modyfikujący stopnie na krawędziach zgodnie z regułą minimum
5. Algorytm automatycznej h adaptacji, łamiący wybrane elementy w wybranych kierunkach 6. Algorytm automatycznej hp adaptacji, opisany w tym rozdziale
Jeśli narysujemy wykres zbieżności tych algorytmów w taki sposób, że na osi poziomej umieścimy rozmiar siatki rozumiany jako liczbę funkcji bazowych na całej siatce, co odpowiada liczbie niewiadomych współczynników tych funkcji bazowych czyli rozmiarowi macierzy układu równań do rozwiązania, a na osi pionowej błąd wzgłędny liczony w normie H1 pomiędzy rozwiązaniem na siatce rzadkiej i siatce gęstej
wówczas okaże się że jedynie algorytm automatycznej hp adaptacji posiada cechę zbieżności eksponencjalnej. Jest on zdecydowanie najszybszy, i potrafi rozwiązać zadany problem obliczeniowy z dokładnością niemożliwą do uzyskania przez inne algorytmy adaptacyjne.
Rysunek 41: Zbieżność różnych algorytmów adaptacyjnych.
W praktyce większość problemów inżynieryjnych wymaga dokładności z zakresu 5 procent i wielomianów kwadratowych. Są jednak takie zadania obliczeniowe, gdzie wysoka dokładność rozwiązania jest niezbędna. Przykładem takich zadań obliczeniowych są na przykład symulacje przepływów wokół skrzydeł samolotów, które wymagają bardzo bardzo dużej dokładności w warstwie przyściennej na skrzydłach samolotu, lub symulacje propagacji fal elektromagnetzcznych w warstwach górotworu, w celu identyfikacji złóż ropy naftowej, które wymagają bardzo dużej dokładności w okolicy anteny odbiorczej.