Metody rozwiązywania nieliniowych układów doradczo-dyfuzyjnych poza Newton-Raphson?

9

Pracuję nad projektem, w którym mam dwie domeny sprzężone za pomocą adv-diff za pomocą odpowiednich terminów źródłowych (jedna domena dodaje masę, druga odejmuje masę). Dla zwięzłości modeluję je w stanie ustalonym. Równania to standardowe równanie transportowo-dyfuzyjne z terminem źródłowym wyglądającym tak:

c1t=0=F1+Q1(c1,c2)c2t=0=F2+Q2(c1,c2)

Gdzie Fi jest dyfuzyjnym i korzystnym strumieniem dla gatunków i, i Qi jest terminem źródłowym dotyczącym gatunku i.

Udało mi się napisać solver dla mojego problemu za pomocą metody Newtona-Raphsona i całkowicie sparowałem dwie domeny za pomocą macierzy mas blokowych, tj .:

Fcoupled=[A100A2][c1,ic2,i]xi[b1(c1,i,c2,i)b2(c1,i,c2,i)]

Termin Fcoupled służy do określania macierzy jakobiańskiej i aktualizacji obu c1 i c2:

J(xi)[xi+1xi]=Fcoupled

lub

xi+1=xi(J(xi))1Fcoupled

Aby przyspieszyć, nie obliczam jakobianu przy każdej iteracji - w tej chwili bawię się co pięć iteracji, które wydają się działać wystarczająco dobrze i utrzymywać rozwiązanie na stałym poziomie.

Problem polega na tym: przeniosę się do większego systemu, w którym obie domeny są w 2D / 2.5D, a obliczenie macierzy jakobińskiej szybko wyczerpie moje dostępne zasoby komputerowe. Buduję ten model do późniejszego wykorzystania w ustawieniach optymalizacyjnych, więc nie mogę też być za kierownicą przy każdej iteracji, dostosowując współczynnik tłumienia itp.

Czy mam rację, szukając gdzie indziej bardziej niezawodnego i algorytmu dla mojego problemu, czy jest to tak dobre, jak to możliwe? Zajrzałem trochę do quasi-linearyzacji, ale nie jestem pewien, jak ma zastosowanie w moim systemie.

Czy są jakieś inne zręczne algorytmy, które mogłem przeoczyć, które mogą rozwiązać układ równań nieliniowych bez uciekania się do ponownego obliczania jakobianów jako przestępstwa?

cbcoutinho
źródło
2
Czy zastanawiałeś się nad rozwiązaniami iteracyjnymi, takimi jak AMG - algebraiczne metody wielosiatkowe. Może być konieczne wymyślenie dobrych warunków wstępnych opartych na fizyce.
NameRakes 27.04.16
1
Czy można uzyskać dostęp do klastra obliczeniowego, w którym można rozpowszechniać formację i rozwiązanie jakobianów za pomocą równoległego pakietu algebry liniowej?
Bill Barth
Nie, nie zastanawiałem się nad AMG, myślałem, że są one przeznaczone tylko dla systemów symetrycznych i nie mogą być stosowane w problemach zdominowanych przez konwekcję. Jeszcze raz spojrzę w literaturze na AMG.
cbcoutinho 27.04.16
Równoległe obliczenia są trudne, ponieważ projekt ten jest opracowywany jako samodzielna aplikacja dla współpracowników, którzy nie mają dostępu do tego rodzaju zasobów. Rozważyłem rozwinięcie MPI w projekt dla samego siebie, ale to zwiększyłoby barierę wejścia dla innych, co było
sednem tego wszystkiego
3
Dlaczego obliczanie jakobianów jest tak problematyczne? Jeśli używasz różnic skończonych / objętości / elementów, powinna ona mieć rzadką część, która jest zawsze taka sama i część ukośną, która się zmienia, ale jest łatwa do obliczenia.
David Ketcheson

Odpowiedzi:

4

Zakładam, że ograniczeniem w 2D i 3D jest przechowywanie jakobianów.

Jedną z opcji jest zachowanie pochodnych czasu i użycie wyraźnego „pseudo” krokowego czasu, aby przejść do stanu ustalonego. Zwykle liczba CFL potrzebna do systemów dyfuzyjnych i reaktywnych może być zbyt mała. Możesz wypróbować nieliniową multigrid (jest to również nazywane multigrid z pamięcią o pełnej aproksymacji) i lokalnym krokiem w celu przyspieszenia konwergencji.

Inną opcją jest użycie w pełni niejawnego schematu, tak jak teraz, ale nie przechowywanie globalnego jakobianu. Możesz użyć niejawnego schematu wolnego od macierzy.

DF(un)δun=F(un)
(gdzie DF jest jakobianem) można rozwiązać za pomocą solvera do podprzestrzeni Kryłowa, takiego jak GMRES i BiCGStab, wykorzystując fakt, że
DF(un)δuF(un+ϵδuδu)F(un)ϵ.
Wynika to z faktu, że GMRES i BiCGStab nie wymagają matrycy LHS A, muszą tylko móc obliczyć swój produkt Ax dany wektor x.

Teraz z odpowiednią wartością ϵ (zwykle około 107dla pływaków podwójnej precyzji) możesz wykonać pętlę Newtona bez obliczania lub przechowywania jakobianów. Wiem na pewno, że ta technika jest używana do rozwiązywania niektórych nietrywialnych przypadków w obliczeniowej dynamice płynów. Należy jednak pamiętać, że liczba ocen funkcjiF będzie więcej niż w technice przechowywania macierzy, zamiast wymagać produktu macierz-wektor.

Inną rzeczą do odnotowania jest to, że jeśli twój system jest taki, że potrzebny jest potężny warunek wstępny (tj. Jacobi lub blok-Jacobi nie wystarczą), możesz spróbować użyć wyżej wspomnianej metody jako płynniejszej w schemacie wielosieciowym. Jeśli chcesz wypróbować punktowy lub blokowy warunek wstępny Jacobiego, możesz obliczyć i zapisać tylko ukośne elementy lub ukośne bloki jakobianu, co nie jest wiele. Chciałbym również wspomnieć, że wdrożenie warunku Gaussa-Seidla lub SSOR może być możliwe bez wyraźnego przechowywania jakobianu. W tym artykule opisano implementację wolnego od macierzy GMRES wstępnie wyposażonego w pozbawiony macierzy symetryczny Gaussa-Seidela w kontekście obliczeniowej dynamiki płynów.

Aditya Kashi
źródło
1

Z mojego doświadczenia z równaniami Naviera-Stokesa można sobie radzić bardzo dobrze bez całkowicie ukrytych schematów.

Jeśli potrzebujesz tylko szybkiego schematu numerycznego dla rozwiązania ewolucji czasu, spójrz na schematy IMEX (niejawne); patrz np. ten artykuł Aschera, Ruuth, Spiteri Implicit-Explicit Metody Runge-Kutta dla zależnych od czasu równań częściowych różniczkowych .

Możesz nawet spróbować po prostu użyć jawnego schematu integracji czasowej wysokiego poziomu z kontrolą wielkości kroku (jak Matlab ODE45). Możesz jednak wpaść w kłopoty z powodu sztywności twojego systemu, która pochodzi z części dyfuzyjnej. Na szczęście część dyfuzyjna jest liniowa, dzięki czemu można ją leczyć niejawnie (taka jest idea schematów IMEX).

Jan
źródło
0

Na początku zastanawiałem się nad dodaniem tylko uwagi, ale miejsca było za mało, dlatego dodam krótki opis moich doświadczeń z tym tematem.

Po pierwsze, spójrz na swoją notację Fcoupled Przypuszczam, że nie widzę sprzężonej formy b1 i b2 będą zależne od c1,i i c2,i. Ponadto jeśliA1 i A2 są reprezentacjami macierzowymi przybliżeń F1 i F2 to nie powinny zależeć tylko od ci, ale także wartości sąsiada, ale może to być tylko błędne zrozumienie notacji.

Jako komentarz ogólny chciałbym dodać, że użycie analitycznego jakobianu wydaje się być jedynym sposobem na uzyskanie kwadratowej zbieżności nieliniowego solwera iteracyjnego (tzn. W twoim przypadku solvera Newtona-Raphsona). Czy zaobserwowałeś to w swoim przypadku? Jest to dość ważne, ponieważ w przeciwnym razie może dojść do nieporozumień w przybliżeniach (linearyzacja).

We wszystkich aplikacjach, w które byłem zaangażowany (niektóre z nich obejmowały obliczenia na dużą skalę) nigdy nie mieliśmy problemów z czasochłonnością składania jakobianów, najbardziej czasochłonnym problemem zawsze było zastosowanie liniowego solvera. Analityczny jakobian (jeśli jest dostępny) zawsze był w aplikacjach, nad którymi pracowałem nad preferowanym wyborem ze względu na kwadratową zbieżność. W nielicznych przypadkach taki nieliniowy solver wytwarza macierz, która powoduje problemy z konwergencją iteracyjnego solwera liniowego, więc próbowaliśmy zastosować prostszą linearyzację niż analityczna metoda jakobowska, aby pomóc liniowemu solverowi. Taki kompromis między zachowaniem nieliniowego i liniowego rozwiązania algebraicznego w zależności od linearyzacji nieliniowego układu algebraicznego zawsze był trudny i nie mogłem dać ogólnej rekomendacji.

Ale masz rację, że wadą (lub właściwością) analitycznego jakobianu dla systemu PDE jest to, że wytwarza on sprzężony system algebraiczny, więc jeśli oddzielisz taki system w jakikolwiek sposób, np. Rozwiązując osobno przybliżenie każdego PDE, powiedzmy, iteracyjny podział następnie tracisz kwadratową zbieżność globalnego solvera. Ale przynajmniej, jeśli rozwiążesz osobno każdy dyskretny (oddzielony) PDE, możesz ponownie przyspieszyć rozwiązanie tego konkretnego problemu, stosując metodę Newtona-Raphsona.

Peter Frolkovič
źródło
Cześć @Peter, masz rację co do sprzężenia, edytowałem główne równanie do pokazania i b1 i b2 są obie funkcje c1 i c2. MatryceA1 i A2w tym przypadku są to macierze sztywności obu układów, które opracowuje się metodą elementów skończonych. Są to tylko funkcje współrzędnych węzłów, a nie zmiennych stanu.F1 i F2są wektorami, więc są funkcjami wektora zmiennych stanu, a nie tylko jednej zmiennej. Obliczam jakobian numerycznie, używając różnic skończonych. Do tej pory nie badałem analitycznego jakobianu.
cbcoutinho