first commit

main.py

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