Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
import time
from numpy import *
import matplotlib.pyplot as plt
nn = 128
it = 12000
a = 0.01
dt = 1E-7
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
cp = zeros((nn + 2, nn + 2))
f = zeros((nn + 2, nn + 2))
cp[1:-1, 1:-1] = c
for _ in range(it):
cp[0, :] = cp[-2, :]
cp[-1, :] = cp[1, :]
cp[:, 0] = cp[:, -2]
cp[:, -1] = cp[:, 1]
f[1:-1, 1:-1] = cp[1:-1, 1:-1]**3 - cp[1:-1, 1:-1] - (a / dx)**2 * (cp[0:-2, 1:-1] + cp[2:, 1:-1] + cp[1:-1, 0:-2] + cp[1:-1, 2:] - 4 * cp[1:-1, 1:-1])
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."""
global c
# wavenumbers
k = power(2 * pi * fft.fftfreq(nn, d=dx), 2)
kk = array([k,]*nn) + array([k,]*nn).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
if __name__ == "__main__":
tic = time.time()
# cahn_hilliard_fdiff(it)
cahn_hilliard_spectral(it)
toc = time.time() - tic
print(f"Elapsed time is {toc} seconds.")
plt.contourf(xy, xy, c, cmap="jet")
plt.show()