Współczynnik konwergencji solwera Poissona FFT

Jaki jest teoretyczny współczynnik zbieżności dla narzędzia FFT Poison?

Rozwiązuję równanie Poissona: przy

\nabla^{2)} {V.}_{H.} (x, y, z) = - 4 π n (x, y, z)

$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$

w domenie

z okresowymi warunkami brzegowymi. Ta gęstość ładunku jest neutralna netto. Rozwiązanie daje:

n (x, y, z) = \frac{3)}{π} ((x - 1)^{2)} + (y - 1)^{2)} + (z - 1)^{2)} - 1)

$n(x, y, z) = {3\over\pi} ((x-1)^2 + (y-1)^2 + (z-1)^2 - 1)$

[0, 2] \times [0, 2] \times [0, 2]

$[0, 2] \times [0, 2] \times [0, 2]$

gdzie

. W przestrzeni wzajemnej

{V.}_{H.} (x) = \int \frac{n (y)}{| x - y |} {re}^{3)} y

$V_H({\bf x}) = \int {n({\bf y})\over |{\bf x}-{\bf y}|} d^3 y$

x = (x, y, z)

${\bf x}=(x, y, z)$

gdzie

są wzajemnymi wektorami kosmicznymi. Interesuje mnie energia Hartree:

{V.}_{H.} (sol) = 4 π \frac{n (sol)}{{sol}^{2)}}

$V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}$

G

${\bf G}$

W przestrzeni wzajemnej staje się (po dyskretyzacji):

{mi}_{H.} = \frac{1}{2)} \int \frac{n (x) n (y)}{| x - y |} {re}^{3)} x {re}^{3)} y = \frac{1}{2)} \int {V.}_{H.} (x) n (x) {re}^{3)} x

$E_H = {1\over 2} \int {n({\bf x}) n({\bf y})\over |{\bf x}-{\bf y}|} d^3 x d^3 y = {1\over 2} \int V_H({\bf x}) n({\bf x}) d^3 x$

Termin pominięty, który sprawia, że gęstość ładunku netto neutralnego (a ponieważ jest już neutralne, a wszystko jest zgodna).

{mi}_{H.} = 2) π \sum_{sol \neq 0} \frac{| n (sol) |^{2)}}{{sol}^{2)}}

$E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2}$

G = 0

${\bf G}=0$

{mi}_{H.} = \frac{128}{35 π} = 1.16410 ...

$E_H = {128\over 35\pi} = 1.16410...$

Oto program korzystający z NumPy, który wykonuje obliczenia.

from numpy import empty, pi, meshgrid, linspace, sum
from numpy.fft import fftn, fftfreq
E_exact = 128/(35*pi)
print "Hartree Energy (exact):      %.15f" % E_exact
f = open("conv.txt", "w")
for N in range(3, 384, 10):
    print "N =", N
    L = 2.
    x1d = linspace(0, L, N)
    x, y, z = meshgrid(x1d, x1d, x1d)

    nr = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
    ng = fftn(nr) / N**3

    G1d = N * fftfreq(N) * 2*pi/L
    kx, ky, kz = meshgrid(G1d, G1d, G1d)
    G2 = kx**2+ky**2+kz**2
    G2[0, 0, 0] = 1  # omit the G=0 term

    tmp = 2*pi*abs(ng)**2 / G2
    tmp[0, 0, 0] = 0  # omit the G=0 term
    E = sum(tmp) * L**3
    print "Hartree Energy (calculated): %.15f" % E
    f.write("%d %.15f\n" % (N, E))
f.close()

A oto wykres konwergencji (wystarczy wykreślić conv.txtpowyższy skrypt, oto notatnik, który robi to, jeśli chcesz się tym bawić):

Wykres konwergencji FFT

Jak widać, zbieżność jest liniowa, co było dla mnie zaskoczeniem, pomyślałem, że FFT zbiega się znacznie szybciej.

Aktualizacja :

Rozwiązanie ma granicę na granicy (wcześniej nie zdawałem sobie z tego sprawy). Aby FFT szybko się zjednoczył, rozwiązanie musi mieć gładkie wszystkie pochodne. Więc spróbowałem również następującej prawej strony:

nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

$V_H=\sin(\pi x)\sin(\pi y)\sin(\pi z)$ $E_H = \frac{3\pi}{8}$

Czy ktoś zna jakiś punkt odniesienia w 3D, dzięki czemu mogę zobaczyć szybszą zbieżność niż liniową?

— Ondřej Čertík
źródło

Ondrej, czy transformacja Fouriera twojej gładkiej gęstości nie jest funkcją delty? Przyznaję, że jestem zbyt leniwy, aby go uruchomić, ale powinien uzyskać dokładną odpowiedź przy pierwszej próbie.

— Matt Knepley,

Myślę, że to jest. Ale nie zbiega się w jednej iteracji, co można zobaczyć na podstawie wykresów zeszytów. Nie mam pojęcia o co chodzi.

— Ondřej Čertík

Ondrej, czy na pewno Twoja implementacja jest poprawna? Pamiętam, jak próbowałem wdrożyć solwery spektralne do zadań domowych w szkole gradowej i całkowicie flubować stałe. Zauważam, że mierzysz błąd, patrząc na bezwzględną odległość między obliczoną a dokładną energią. Jak wygląda twoja zbieżność z rzeczywistym rozwiązaniem problemu? Powinno to być łatwe do obliczenia, a nawet nakreślenia 2-wymiarowego fragmentu problemu.

— Aron Ahmadia

Aron --- sprawdziłem moją implementację pod kątem innego kodu, ale sprawdziłem ją pod kątem mojego błędnego początkowego próbkowania, więc miałem ten sam błąd w obu kodach. Matt miał rację, teraz zbiega się za pierwszym razem. Zobacz moją odpowiedź poniżej.

— Ondřej Čertík

Pozwól mi najpierw odpowiedzieć na wszystkie pytania:

Jaki jest teoretyczny współczynnik zbieżności dla narzędzia FFT Poison?

Teoretyczna zbieżność jest wykładnicza, o ile rozwiązanie jest wystarczająco gładkie.

Jak szybko ta energia powinna się zbiegać?

$E_H$

Czy ktoś zna jakiś punkt odniesienia w 3D, dzięki czemu mogę zobaczyć szybszą zbieżność niż liniową?

Każda prawa strona, która wytwarza rozwiązanie, które jest okresowe i nieskończenie zróżnicowane (w tym przez granicę okresową), powinna zbiegać się wykładniczo.

W powyższym kodzie występuje błąd, który powoduje, że konwergencja jest wolniejsza niż wykładnicza. Za pomocą kodu gładkiej gęstości ( https://gist.github.com/certik/5549650/ ) następująca łatka naprawia błąd:

@@ -6,7 +6,7 @@ f = open("conv.txt", "w")
 for N in range(3, 180, 10):
     print "N =", N
     L = 2.
-    x1d = linspace(0, L, N)
+    x1d = linspace(0, L, N+1)[:-1]
     x, y, z = meshgrid(x1d, x1d, x1d)

     nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

Problem polegał na tym, że próbkowanie przestrzeni rzeczywistej nie może powtórzyć pierwszego i ostatniego punktu (który ma tę samą wartość ze względu na okresowe warunki brzegowe). Innymi słowy, problemem było ustawienie początkowego próbkowania.

Po tej poprawce gęstość zbiega się w jednej iteracji, jak powiedział Matt powyżej. Więc nawet nie wykreśliłem wykresów konwergencji.

Można jednak spróbować trudniejszej gęstości, na przykład:

     nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4

następnie konwergencja jest wykładnicza, zgodnie z oczekiwaniami. Oto wykresy konwergencji dla tej gęstości: wprowadź opis zdjęcia tutaj

— Ondřej Čertík
źródło

Niesamowite. Przepraszam, że nie byłam więcej pomocą!

— Aron Ahmadia