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

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.

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

a warunki brzegowe są

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.$\Omega_8$$\Omega_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).

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$\rho(x) = 0$. Z prawej strony (11-20) mamy$f(x_{20}) = 0$ z powodu stanu Neumanna, który mówi nam o ochronie tej komórki $f(x_{19}) = \ldots = f(x_{11}) = 0$. Ponieważ strumień jest średnią wartości komórek, mówi nam to$\phi(x)$ jest stały dla wszystkich tych komórek.

Z lewej strony (1-9) mamy $f(x_{i+1})-f(x_i) = 0$. Gdyby$f(-10) = 0$ i wtedy używamy komórek duchów $f(-10) = \phi_{9.5} - \phi_{\rm ghost} = \phi_{9.5}$. Zapewnia to konserwacja w ciągu kilku następnych komórek$f(x_i) = f(-10) = \phi_{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łąd$f(10)$, który następnie propaguje się do tyłu i miesza zarówno nachylenie w lewej połowie domeny, jak i wartość $\phi(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.

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ę?