Taux de convergence du solveur de Poisson FFT

16

Quel est le taux de convergence théorique d'un solveur FFT Poison?

Je une équation de Poisson: avec sur le domaine avec périodique condition limite. Cette densité de charge est nette neutre. La solution est donnée par: où . Dans l'espace réciproque où sont les vecteurs spatiaux réciproques. Je m'intéresse à l'énergie Hartree: n ( x , y , z ) = 3

2VH(X,y,z)=-4πn(X,y,z)
[0,2]×[0,2]×[0,2]VH(x)=n( y )
n(X,y,z)=3π((X-1)2+(y-1)2+(z-1)2-1)
[0,2]×[0,2]×[0,2]x=(x,y,z)VH(G)=4πn(G)
VH(X)=n(y)|X-y|3y
X=(X,y,z) GEH=1
VH(g)=4πn(g)g2
g
EH=12n(X)n(y)|X-y|3X3y=12VH(X)n(X)3X
Dans l'espace réciproque cela devient (après discrétisation): Le terme est omis, ce qui rend effectivement la densité de charge nette neutre (et puisqu'il est déjà neutre, alors tout est cohérent).
EH=2πg0|n(g)|2g2
g=0

Pour le problème de test ci-dessus, cela peut être évalué analytiquement et on obtient: quelle vitesse cette énergie devrait-elle converger?

EH=12835π=1.16410 ...

Voici un programme utilisant NumPy qui fait le calcul.

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()

Et voici un graphique de convergence (en traçant simplement le conv.txtscript ci-dessus, voici un cahier qui le fait si vous voulez jouer avec cela vous-même):

Graphique de convergence FFT

Comme vous pouvez le voir, la convergence est linéaire, ce qui m'a surpris, j'ai pensé que la FFT converge beaucoup plus vite que cela.

Mise à jour :

La solution a un point culminant à la frontière (je ne m'en étais pas rendu compte auparavant). Pour que la FFT converge rapidement, la solution doit avoir tous les dérivés lisses. J'ai donc également essayé le côté droit suivant:

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

Que vous pouvez simplement mettre dans le script ci-dessus ( script mis à jour ). La solution exacte est , qui devrait être infiniment différentiable. L'intégrale exacte dans ce cas est . Pourtant, le solveur FFT ne converge toujours que linéairement vers cette solution exacte, comme cela peut être vérifié en exécutant le script ci-dessus et en traçant la convergence ( bloc-notes mis à jour avec des tracés).VH=péché(πX)péché(πy)péché(πz)EH=3π8

Quelqu'un connaît-il une référence en 3D pour que je puisse voir une convergence plus rapide que linéaire?

Ondřej Čertík
la source
Ondrej, la transformée de Fourier de votre densité lisse n'est-elle pas une fonction delta? J'avoue être trop paresseux pour l'exécuter, mais il devrait obtenir la réponse exacte du premier coup.
Matt Knepley
Je pense que c'est. Mais il ne converge pas en une seule itération, comme on peut le voir sur les graphiques du carnet. Je n'ai aucune idée de ce qu'il se passe.
Ondřej Čertík
Ondrej, êtes-vous sûr que votre implémentation est correcte? Je me souviens avoir essayé d'implémenter des solveurs spectraux pour une affectation aux devoirs dans une école doctorale et avoir complètement effacé les constantes. Je remarque que vous mesurez l'erreur en regardant la distance absolue entre l'énergie calculée et l'énergie exacte. À quoi ressemble votre convergence avec la solution réelle du problème? Cela devrait être facile à calculer et même à tracer sur une tranche 2D du problème.
Aron Ahmadia
Aron --- J'ai vérifié mon implémentation par rapport à un autre code, mais je le vérifiais pour mon mauvais échantillonnage initial, donc j'avais le même bug dans les deux codes. Matt avait raison, maintenant il converge du premier coup. Voir ma réponse ci-dessous.
Ondřej Čertík

Réponses:

10

Permettez-moi d'abord de répondre à toutes les questions:

Quel est le taux de convergence théorique d'un solveur FFT Poison?

La convergence théorique est exponentielle tant que la solution est suffisamment lisse.

À quelle vitesse cette énergie devrait-elle converger?

EH

Quelqu'un connaît-il une référence en 3D pour que je puisse voir une convergence plus rapide que linéaire?

Tout côté droit qui produit une solution périodique et différentiable à l'infini (y compris à travers la frontière périodique) doit converger de façon exponentielle.


Dans le code ci-dessus, il se trouve qu'il y a un bogue, qui ralentit la convergence par rapport à l'exponentielle. À l'aide du code de densité lisse ( https://gist.github.com/certik/5549650/ ), le correctif suivant corrige le bogue:

@@ -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

Le problème était que l'échantillonnage dans l'espace réel ne peut pas répéter le premier et le dernier point (qui ont la même valeur en raison de la condition aux limites périodique). En d'autres termes, le problème résidait dans la mise en place de l'échantillonnage initial.

Après ce correctif, la densité converge en une seule itération, comme Matt l'a dit ci-dessus. Je n'ai donc même pas tracé les graphiques de convergence.

Cependant, on peut essayer une densité plus difficile, par exemple:

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

alors la convergence est exponentielle, comme prévu. Voici les graphiques de convergence pour cette densité: entrez la description de l'image ici entrez la description de l'image ici

Ondřej Čertík
la source
Impressionnant. Désolé, je n'ai pas été plus aidé!
Aron Ahmadia