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, c): 
    """Cahn-Hilliard integration for it iterations using finite differences method."""

    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]

    return c
    

def cahn_hilliard_spectral(it, c):
    """Cahn-Hilliard integration for it iterations using fft spectral method."""

    # wavenumbers
    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.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()
    
    c1 = cahn_hilliard_fdiff(it, c)
    c2 = cahn_hilliard_spectral(it, c)

    toc = time.time() - tic
    
    print(f"Elapsed time is {toc} seconds.")
    
    fig, (ax1, ax2) = plt.subplots(1, 2)
    fig.suptitle(f"Cahn-Hilliard integration for {it} iterations ")
    ax1.imshow(c1, cmap="jet", vmin=-1, vmax=1)
    ax1.set_title("Finite differences method")
    ax2.imshow(c2, cmap="jet", vmin=-1, vmax=1)
    ax2.set_title("Spectral method")
    plt.show()