archive.py

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

    k1 = zeros((nn, nn))
    k2 = zeros((nn, nn))
    k3 = zeros((nn, nn))
    k4 = zeros((nn, nn))

    # 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,]*(nn//2+1)).transpose()

    for _ in range(it): 
        
        k1 = fft.irfft2(-kk * fft.rfft2( c**3 - c - a**2 * fft.irfft2(-kk * fft.rfft2(c)).real)).real
        k2 = fft.irfft2(-kk * fft.rfft2( (c+dt*k1/2)**3 - (c+dt*k1/2) - a**2 * fft.irfft2(-kk * fft.rfft2((c+dt*k1/2))).real)).real
        k3 = fft.irfft2(-kk * fft.rfft2( (c+dt*k2/2)**3 - (c+dt*k2/2) - a**2 * fft.irfft2(-kk * fft.rfft2((c+dt*k2/2))).real)).real
        k4 = fft.irfft2(-kk * fft.rfft2( (c+dt*k3)**3 - (c+dt*k3) - a**2 * fft.irfft2(-kk * fft.rfft2((c+dt*k3))).real)).real

        c = c + 1/6 * dt * (k1 + 2*k2 + 2*k3 + k4) 

    return c