Osobliwy błąd przy rozwiązywaniu równania Poissona w metodzie nieskończonej objętości siatki (tylko 1D)

9

W ciągu ostatnich kilku dni próbowałem debugować ten błąd. Zastanawiałem się, czy ktoś ma porady, jak postępować.

Rozwiązuję równanie Poissona dla rozkładu ładunku krokowego (powszechny problem w elektrostatyce / fizyce półprzewodników) na niejednorodnej siatce o skończonej objętości, w której nieznane są zdefiniowane w centrach komórek i strumienie na powierzchniach komórek.

0=(ϕx)x+ρ(x)

profil opłat (termin źródłowy) podaje:

ρ(x)={1,if 1x01,if 0x10,otherwise

a warunki brzegowe są

ϕ(xL)=0ϕx|xR=0

a domena to .[10,10]

Korzystam z kodu opracowanego do rozwiązania równania dyfuzja-doradztwo-reakcja (sam napisałem, aby zobaczyć moje notatki tutaj, http://danieljfarrell.github.io/FVM ). Równanie rada-dyfuzja-reakcja jest bardziej ogólnym przypadkiem równania Poissona. Rzeczywiście równanie Poissona można odzyskać, ustawiając prędkość doradczą na zero i usuwając składnik przejściowy.

Kod został przetestowany pod kątem różnych sytuacji dla jednolitych, niejednorodnych i losowych siatek i zawsze daje rozsądne rozwiązania ( http://danieljfarrell.github.io/FVM/examples.html ) dla równania dyfuzja-doradztwo-reakcja.

Aby pokazać, gdzie kod się psuje, podałem następujący przykład. Ustawiam jednolitą siatkę złożoną z 20 komórek, a następnie robię to niejednorodną , usuwając pojedynczą komórkę. Na lewym rysunku usunąłem komórkę a po prawej zostało usunięte. Dziewiąta komórka obejmuje region, w którym termin źródłowy (tj. Ładunek) zmienia znak. Błąd pojawia się, gdy siatka jest nierównomierna w regionie, w którym znak reakcji zmienia znak . Jak widać poniżej.Ω8Ω9

Wszelkie pomysły, co może być przyczyną tego problemu? Daj mi znać, jeśli pomocne byłyby dodatkowe informacje dotyczące dyskretyzacji (nie chciałem umieszczać zbyt wielu szczegółów w tym pytaniu).

Osobliwy błąd przy rozwiązywaniu równania Poissona

boyfarrell
źródło
Czy możesz określić sposób nałożenia warunku Dirichleta na i co rozumiesz przez jako warunek początkowy (czy to równanie, które określiłeś jako stan ustalony)? x=0ρ=1
Jesse Chan,
Jak wygląda termin reakcji?
Jan
Jakiego schematu używasz do przybliżania całek terminu źródłowego? To zachowanie może być również spowodowane niewystarczającym próbkowaniem źródła. (Co, prawdopodobnie, jest taki sam mechanizm wspomniano w odpowiedzi @JLC „s.)
Jan
Zaktualizowałem pytanie, aby użyć standardowej terminologii. Mam termin źródłowy ( ), a nie termin reakcji, ponieważ, jak zauważyłeś, potrzebujemy tylko wartości stanu ustalonego. Podano teraz poprawną zależność przestrzenną (początkowa wartość była niepoprawna). ρρ
boyfarrell,
@JLC BC Dirichleta są nakładane przy użyciu metody komórki-widma (moje notatki online są nieaktualne w odniesieniu do tego szczegółu implementacji), zobacz tutaj jak to zrobić, scicomp.stackexchange.com/questions/8538/…
boyfarrell

Odpowiedzi:

9

Nawiasem mówiąc, twoja dokumentacja github jest fantastyczna.

To tylko przypuszczenie z metod DG, które mogą mieć podobne problemy, jeśli strumienie numeryczne nie zostaną starannie wybrane (sądzę, że metody FV są podzbiorem metod DG). Jeśli używasz interpolacji od centrów komórek do zdefiniowania swoich strumieni, to powinno to być równoważne z użyciem średniej jako strumienia liczbowego w DG i z zastosowaniem stałej częściowej. W przypadku standardowych metod DG Poissona prowadzi to do numerycznie nieunikalnych rozwiązań - można uzyskać nietrywialne miejsce zerowe dla operatora dyskretnego, co, jak sądzę, jest przyczyną problemów w drugim przykładzie. Zobacz ten dokument DG dla ich teorii na ten temat od strony DG.

Spróbuję wykpić przykład FV, który pokazuje, jak to się dzieje.

Edycja: oto mały przykład tego, co się dzieje. Rozważ komórki 1-9 i 11-20, w którychρ(x)=0. Z prawej strony (11-20) mamyf(x20)=0 z powodu stanu Neumanna, który mówi nam o ochronie tej komórki f(x19)==f(x11)=0. Ponieważ strumień jest średnią wartości komórek, mówi nam toϕ(x) jest stały dla wszystkich tych komórek.

Z lewej strony (1-9) mamy f(xi+1)f(xi)=0. Gdybyf(10)=0 i wtedy używamy komórek duchów f(10)=ϕ9.5ϕghost=ϕ9.5. Zapewnia to konserwacja w ciągu kilku następnych komórekf(xi)=f(10)=ϕ9.5(tj. stałe nachylenie). Pamiętaj jednak, że może to być dowolne nachylenie, tylko stałe.

Problem pojawia się w środkowej komórce. Tak jak wspomniał Jan, podkreślacie wymuszenie w drugiej siatce. To odrzuca równania równowagi w tym punkcie, daje błądf(10), który następnie propaguje się do tyłu i miesza zarówno nachylenie w lewej połowie domeny, jak i wartość ϕ(9.5).

Ta wrażliwość na błędy w wymuszaniu jest problematyczna - w przeciwieństwie do metod FEM lub FD, które wyraźnie wymuszają stan Dirchleta w x=10, FV wymusza go słabo za pomocą węzłów-widm. Intuicyjnie nałożenie słabych węzłów duchów jest jak ustawienie warunku Neumanna na twojej lewej granicy. Jeśli masz dwa warunki Neumanna dla problemu dyfuzji, twój problem jest źle postawiony i ma nietypowe rozwiązanie (możesz dodać dowolną stałą do tego problemu i nadal mieć rozwiązanie). Na poziomie dyskretnym nie ma się tego dość dobrze, ale zachowuje się bardzo wrażliwe i zależne od siatki zachowanie, jak widać w eksperymentach.

Jesse Chan
źródło
Z eksperymentów mogę wykazać, że metoda FVM jest stabilna tylko wtedy, gdy komórki po obu stronach nieciągłości (zmiany znaku) funkcji źródłowej mają równe objętości. Czy Twoja analiza się z tym zgodzi? Oznacza to, że muszę zwracać większą uwagę na generowanie rozsądnej siatki moich problemów, które zrobiłem wcześniej. Może powinienem rozważyć naukę metody MES?
boyfarrell
Odpowiedni artykuł, chociaż nie do końca podążam za wszystkimi szczegółami, jstor.org/discover/10.2307/2157873
boyfarrell 24.09.2013
Metoda FVM jest stabilna tylko w tym przypadku, gdy siatka jest w jakiś sposób wyrównana z funkcją źródła. Jeśli zmieni się funkcja źródła, musisz ponownie dostroić siatkę. Nie sądzę, aby generowanie rozsądnej siatki było właściwym podejściem do tego problemu - masz niestabilną metodę.
Jesse Chan
To dobre znalezisko. Suli jest solidnym analitykiem. Powiedziałbym, że nauka FEM może być świetną zabawą, ale FD powinna również działać na wszelkie problemy eliptyczne 1D. Możesz także zobaczyć, co robią ludzie FV (może zwiększyć swoje przepływy o kary), aby uzyskać zbieżność dla problemów eliptycznych drugiego rzędu na ogólnych polach. Matematyczna mądrość ludowa zwykle mówi, że FV / odwrócony FD jest świetny dla problemów hiperbolicznych, podczas gdy FEM / centralna różnicowa FD jest świetna dla eliptycznych.
Jesse Chan
Sprawdzam ten problem. Ponownie przeczytaj twoją odpowiedź Muszę powiedzieć, że to fantastyczne! Widzę twój punkt widzenia, że ​​metoda powinna się zmienić, ponieważ to jest podstawa problemu (nie siatki). Czy masz jakieś sugestie lub rzeczy, które mógłbym śledzić (które są dostępne dla osób niebędących ekspertami), w jaki sposób lepiej przybliżyć przepływ w tym przypadku? Tj. W sposób, który mógłby uczynić go bardziej stabilnym. Jeśli to możliwe, chciałbym znaleźć lepszą FVM dla tego równania.
boyfarrell
0

Pierwszą rzeczą, na którą należy zwrócić uwagę, są warunki brzegowe. Ponieważ możesz zmienić nachylenie i wartość, nie masz ani warunków Dirichleta, ani Neumanna.

Zatem każda linia prosta jest rozwiązaniem, w którym prawa strona ma zero. Masz tę część.

Twoje strumienie są prawdopodobnie zależne h. Czy używasz prawidłowego?h gdzie eliminujesz komórkę?

Guido Kanschat
źródło
1
Nie, to nieprawda. Problem jest dobrze postawiony. W przypadkuρ0 tylko ϕ0jest rozwiązaniem, nie ma innej funkcji liniowej, która ma zero w jednym punkcie i ma zerowe nachylenie w drugim punkcie.
Jan