Jak debugować kod numeryczny, co może być źródłem tego błędu oscylacyjnego?

16

Cichą wiedzę można uzyskać z doświadczenia, zastanawiałem się tylko, czy ktoś wcześniej widział coś podobnego do tego. Wykres pokazuje początkowy warunek (zielony) dla równania rada-dyfuzja, następnie roztwór przy iteracji 200 (niebieski), a następnie ponownie przy iteracji 400 (czerwony).

Błąd oscylacyjny

Rozwiązanie równania dyfuzyjno-doradczego pojawia się po kilku iteracjach. Liczba Pécleta i warunek CFL jest spełniony, , więc równania powinny być stabilne. Przewiduję, że mam błąd w kodzie numerycznym.μ0,07do0,0015

Tło. Dyskretyzacja stanowi zasadniczą różnicę zarówno w odniesieniu do warunków doradztwa, jak i rozpowszechniania. Uważam, że jest to pierwsza kolejność porad i druga kolejność rozpowszechniania. Zaimplementowałem to, stosując podejście o skończonej objętości (po raz pierwszy), w którym wartości współczynników (prędkości i współczynnika dyfuzji) na powierzchniach komórek określa się poprzez interpolację liniową z średnich komórek. Stosuję warunek brzegowy Robina na lewej i prawej powierzchni i ustawiam strumień na granicach na zero.

Jak debugujesz swój kod numeryczny? Czy ktoś wcześniej coś takiego zrobił, gdzie byłoby dobre miejsce, aby zacząć szukać?

Aktualizacja

Aktualizacja

Rozwiązanie nie może być prostsze! Właśnie popełniłem błąd znakowy dla terminu dyfuzji. To dziwne, jestem pewien, że nie opublikowałem tego, nie znalazłbym błędu! Jeśli ktoś chce podzielić się wskazówkami na temat debugowania swojego kodu numerycznego, nadal jestem zainteresowany. Nie mam metody, jest trochę trafiona i chybiam, próbuję różnych rzeczy, aby zdobyć wskazówki, ale ten proces może zająć tygodnie (czasem).

Dowód, że to działa (należy pamiętać, że przy metodzie o skończonej objętości wszystko, co musisz zrobić, aby obliczyć obszar, to suma szerokości wysokości dla wszystkich komórek, jeśli użyjesz metody integracji, takiej jak numpy.trapz, wyniki zawierają wartości liczbowe błąd metody trapezowej). Co tu się dzieje? Istnieje stała prędkość i współczynniki dyfuzji, ale przy zamkniętych warunkach brzegowych. Dlatego na granicy widzimy równowagę między polem prędkości popychającym w prawo a dyfuzją popychaną w lewo.×

Równanie dyfuzyjne z zamkniętymi warunkami brzegowymi metodą skończonej objętości.

boyfarrell
źródło
2
Jakiego rodzaju dyskretyzacji używasz? Jaka metoda zamówienia? Jakie są twoje warunki brzegowe?
Geoff Oxberry
Dziękuję @GeoffOxberry, zaktualizowałem o więcej szczegółów. Chociaż, aby naprawdę zrozumieć, co zrobiłem, możesz przeczytać notatki z mojej książki laboratoryjnej na powyższym linku github.
boyfarrell
2
Podziwiam twoje sumienne podejście do tego problemu. Czy te wykresy pochodzą z równań zależnych od czasu? Jeśli tak, co zaobserwujesz dla ? θ=0,0,5,1
Jan
1
Próbowałem tego, daje zasadniczo to samo. Ale teraz czuję się głupio, testowałem pod kątem błędu znaku ... ustawienie i sprawia, że ​​działa (w granicy dyfuzji tak dyfuzja ujemna sprawia, że ​​działa!). Dla przyszłego odniesienia powyższe może być zgodne z błędem znaku. Jeszcze nie w pełni przetestowany, więc wkrótce się z nim skontaktuję! PS. dzięki za komplement :)d < 0za=0re<0
boyfarrell

Odpowiedzi:

9

Zebrałem trochę mojego doświadczenia w debugowaniu kodów numerycznych tutaj: deal.II FAQ: debugowanie . Nie wiem, czy to pomogłoby ci w tym konkretnym przypadku, ale może w innych.

Wolfgang Bangerth
źródło
Cześć, ten link zepsuł się, gdy projekt się przesunął - czy to teraz jest właściwy? github.com/dealii/dealii/wiki/…
hyperpallium
Tak, to jest poprawne. Dziękujemy za aktualizację w komentarzu!
Wolfgang Bangerth,