Skip to content
Extraits de code Groupes Projets
Valider 3250edd0 rédigé par Simon Collignon's avatar Simon Collignon
Parcourir les fichiers

using the hermitian symmetry of the spectrum with rfft

parent 5b815526
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
Fichier ajouté
......@@ -11,10 +11,8 @@ xy, dx = linspace(0, 1, nn, retstep=True)
c = 2 * random.rand(nn, nn) - 1
def cahn_hilliard_fdiff(it):
"""Cahn-Hilliard integration for it iterations using finite differences mehtod."""
global c
def cahn_hilliard_fdiff(it, c):
"""Cahn-Hilliard integration for it iterations using finite differences method."""
cp = zeros((nn + 2, nn + 2))
f = zeros((nn + 2, nn + 2))
......@@ -31,32 +29,40 @@ def cahn_hilliard_fdiff(it):
cp[1:-1, 1:-1] = cp[1:-1, 1:-1] + (dt / dx**2) * (f[0:-2, 1:-1] + f[2:, 1:-1] + f[1:-1, 0:-2] + f[1:-1, 2:] - 4 * f[1:-1, 1:-1])
c = cp[1:-1, 1:-1]
def cahn_hilliard_spectral(it):
"""Cahn-Hilliard integration for it iterations using fft spectral mehtod."""
return c
global c
def cahn_hilliard_spectral(it, c):
"""Cahn-Hilliard integration for it iterations using fft spectral method."""
# wavenumbers
k = power(2 * pi * fft.fftfreq(nn, d=dx), 2)
kk = array([k,]*nn) + array([k,]*nn).transpose()
kr = power(2 * pi * fft.rfftfreq(nn, d=dx), 2)
kc = power(2 * pi * fft.fftfreq(nn, d=dx), 2)
kk = array([kr,]*nn) + array([kc,]*int((nn/2)+1)) .transpose()
for _ in range(it):
c = c + dt * fft.ifft2(-kk * fft.fft2( c**3 - c - a**2 * fft.ifft2(-kk * fft.fft2(c)).real)).real
c = c + dt * fft.irfft2(-kk * fft.rfft2( c**3 - c - a**2 * fft.irfft2(-kk * fft.rfft2(c)).real)).real
return c
if __name__ == "__main__":
tic = time.time()
# cahn_hilliard_fdiff(it)
cahn_hilliard_spectral(it)
c1 = cahn_hilliard_fdiff(it, c)
c2 = cahn_hilliard_spectral(it, c)
toc = time.time() - tic
print(f"Elapsed time is {toc} seconds.")
plt.contourf(xy, xy, c, cmap="jet")
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.suptitle(f"Cahn-Hilliard integration for {it} iterations ")
ax1.contourf(xy, xy, c1, cmap="jet")
ax1.set_title("Finite differences method")
ax2.contourf(xy, xy, c2 , cmap="jet")
ax2.set_title("Spectral method")
plt.show()
0% Chargement en cours ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter