updates

.gitignore

+data
+figs
+__pycache__

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

free_energy.py

+from numpy import *
+from matplotlib import pyplot as plt
+import matplotlib
+
+
+matplotlib.rcParams.update({
+        "pgf.texsystem": "pdflatex",
+        'font.family': 'serif',
+        'toolbar': 'None',
+        'text.usetex': True,
+        'pgf.rcfonts': False,
+        'legend.fancybox': False,
+        'legend.shadow': False,
+        'figure.figsize': [3.4, 3.6],
+    })
+cc = linspace(-1.5, 1.5, 1000)
+
+fig, ax = plt.subplots(1)
+ax.axvline(-1, color='k', linestyle='dashed', linewidth=.75)
+ax.axvline(1, color='k', linestyle='dashed', linewidth=.75)
+plt.title(r"Mixture Free Energy $\Delta_{mix}G_m \approx F(c)$", fontsize=10)
+ax.set_ylabel(r"$F(c) = 0.25(c^2 - 1)^2$")
+ax.set_xlabel(r"Concentration $c = c_{B} - c_{A}$")
+ax.set_yticks([0])
+ax.set_yticklabels(["0"])
+ax.set_xticks([-1, 0, 1])
+ax.set_xticklabels(["$-c_{A} = -1$", "0", "$c_{B} = 1$"])
+ax.plot(cc, .25 * (cc**2 - 1)**2, 'k', linewidth=1)
+plt.savefig('homogeneous_free_energy.pgf')
+# plt.show()

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     
+from scipy.fftpack import dctn, idctn
 
 
+# c = loadtxt('data/initial_c_128.txt').reshape((nn, nn))
+    
 def cahn_hilliard_fdiff(it, c): 
     """Cahn-Hilliard integration for it iterations using finite differences method."""
 
+    nn = len(c[0, :])
+    dx = L / (nn - 1)
+
     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]
-        
+        periodize(cp)
         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])
-
-        f[0, :]  = f[-2, :]
-        f[-1, :] = f[1, :] 
-        f[:, 0]  = f[:, -2]
-        f[:, -1] = f[:, 1]
-
+        periodize(f)
         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]) 
 
     return cp[1:-1, 1:-1]
@@ -40,10 +29,13 @@ def cahn_hilliard_fdiff(it, c):
 def cahn_hilliard_spectral(it, c):
     """Cahn-Hilliard integration for it iterations using fft spectral method."""
 
+    nn = len(c[0, :])
+    dx = L / (nn - 1)
+
     # 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()
+    kk = array([kr,]*nn) + array([kc,]*(nn//2+1)).transpose()
 
     for _ in range(it): 
         
@@ -52,21 +44,338 @@ def cahn_hilliard_spectral(it, c):
     return c
 
 
+def cahn_hilliard_spectral_lss(it, c):
+    """Cahn-Hilliard integration for it iterations using fft spectral method."""
+
+    nn = len(c[0, :])
+    dx = L / (nn - 1)
+
+    # 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): 
+       
+        s = fft.rfft2(c) - dt * kk * fft.rfft2( c**3 - 3 * c) 
+        v = s / (1 + dt * (2 * kk + a**2 * kk**2))
+        c = fft.irfft2(v).real 
+
+    return c
+
+def cahn_hilliard_spectral_lss_ad(t, c, adapt=True):
+    """Cahn-Hilliard integration for it iterations using fft spectral method with adaptative time stepping."""
+
+    nn = len(c[0, :])
+    dx = L / (nn - 1)
+
+    # 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()
+
+    # adaptative time stepping
+    c_n  = zeros((nn, nn))
+    tt = zeros((1,))
+    ddt  = array([dt]) 
+
+    # free energy functional
+    fef0 = free_energy_functional(c) 
+    fef = array([1])
+
+    while(tt[-1] < t): 
+       
+        s = fft.rfft2(c) - ddt[-1] * kk * fft.rfft2( c **3 - 3 * c) 
+        v = s / (1 + ddt[-1] * (2 * kk + a**2 * kk**2))
+        c_n  = fft.irfft2(v).real 
+
+        tt = append(tt, tt[-1] + ddt[-1])
+        ddt = append(ddt, adaptative_time_step((c_n - c)/ddt[-1]))
+
+        fef = append(fef, free_energy_functional(c_n))
+
+        c = c_n 
+
+    # normalize fef
+    fef[1:] = fef[1:] / fef0
+
+    return tt, ddt, fef, c
+
+
+def cahn_hilliard_spectral_dct(it, c):
+    """Cahn-Hilliard integration for it iterations using fft/dct spectral method."""
+
+    nn = len(c[0, :])
+    dx = L / (nn - 1)
+
+    # wavenumbers
+    kr1 = power(2 * pi * fft.rfftfreq(nn, d=dx), 2)
+    kc1 = power(2 * pi * fft.fftfreq(nn, d=dx), 2)
+    kk1 = array([kr1,]*nn) + array([kc1,]*(nn//2+1)).transpose()
+    kc2 = (pi * arange(0, nn))**2
+    kk2 = array([kc2,]*(nn)) + array([kc2,]*(nn)).transpose()
+
+    for _ in range(it): 
+
+        c = c + dt * idctn(-kk2 * dctn(c**3 - c - a**2 * fft.irfft2(-kk1 * fft.rfft2(c)).real, 1, norm="ortho"), 1, norm="ortho")
+
+    return c
+
+
+def periodize(cp): 
+
+    cp[0, :]  = cp[-2, :]
+    cp[-1, :] = cp[1, :] 
+    cp[:, 0]  = cp[:, -2]
+    cp[:, -1] = cp[:, 1]
+
+    return cp
+
+
+def free_energy_functional(c):
+
+    nn = len(c[0, :])
+
+    cp = zeros((nn + 2, nn + 2))
+    cp[1:-1, 1:-1] = c
+    
+    cp = periodize(cp)
+
+    I = 0
+
+    for ii in range(1, nn + 1):
+        for jj in range(1, nn + 1):
+
+            I += dx**2 * (0.25 * (cp[ii, jj]**2 - 1)**2) + (a**2/2) * ((cp[ii + 1, jj] - cp[ii, jj])**2 + (cp[ii, jj + 1] - cp[ii, jj])**2)
+
+    return I
+
+
+def adaptative_time_step(delta):
+    
+    tau =  5e-2
+    # verifier que ça marche avec ord=inf
+    return maximum(tau / linalg.norm(delta, ord=inf), dt)
+
+
+def inf_norm_error(c1, c2):
+    """
+        accurcy: if c1 is n x n, then c2 is 2n x 2n 
+    """
+
+    err = zeros((len(c1), len(c1)))
+    err[:, :] = c1[:, :] - c2[::2, ::2]
+    
+    return linalg.norm(err, ord=inf) / linalg.norm(c2[::2, ::2], ord=inf) 
+
+def inf_norm_error(c1, c2):
+    """
+        accurcy: if c1 is n x n, then c2 is 2n x 2n 
+    """
+
+    err = zeros((len(c1), len(c1)))
+    err[:, :] = c1[:, :] - c2[::2, ::2]
+    
+    return linalg.norm(err, ord=inf) / linalg.norm(c2[::2, ::2], ord=inf) 
+
+def accuracy_test(smooth=True): 
+    """
+        accuracy_test: function which runs tests to compute the relative errors between different grid sizes.
+    """ 
+    a = 0.001
+
+    xy_2048 = linspace(0, 1, 2048).reshape((2048, 1)); 
+    xy_1024 = linspace(0, 1, 1024).reshape((1024, 1)); 
+    xy_512 = linspace(0, 1, 512).reshape((512, 1)); 
+    xy_256 = linspace(0, 1, 256).reshape((256, 1)); 
+    xy_128 = linspace(0, 1, 128).reshape((128, 1)); 
+    xy_64 = linspace(0, 1, 64).reshape((64, 1)); 
+    xy_32 = linspace(0, 1, 32).reshape((32, 1)); 
+
+    if(smooth):
+
+        c_2048 = zeros((2048, 2048)) 
+        c_1024  = zeros((1024, 1024)) 
+        c_512  = zeros((512, 512)) 
+        c_256  = zeros((256, 256))
+        c_128  = zeros((128, 128))
+        c_64   = zeros((64, 64))
+        c_32   = zeros((32, 32))
+
+        c_2048[:, :] = 2 * exp(- ((xy_2048 - .5)**2 + (xy_2048.T - .5)**2) * 20 ) - 1
+        c_1024[:, :] = 2 * exp(- ((xy_1024 - .5)**2 + (xy_1024.T - .5)**2) * 20 ) - 1
+        c_512[:, :]  = 2 * exp(- ((xy_512  - .5)**2 + (xy_512.T  - .5)**2) * 20 ) - 1
+        c_256[:, :]  = 2 * exp(- ((xy_256  - .5)**2 + (xy_256.T  - .5)**2) * 20 ) - 1
+        c_128[:, :]  = 2 * exp(- ((xy_128  - .5)**2 + (xy_128.T  - .5)**2) * 20 ) - 1
+        c_64[:, :]   = 2 * exp(- ((xy_64   - .5)**2 + (xy_64.T   - .5)**2) * 20 ) - 1
+        c_32[:, :]   = 2 * exp(- ((xy_32   - .5)**2 + (xy_32.T   - .5)**2) * 20 ) - 1
+
+    else: 
+
+        c_2048 = zeros((2048, 2048)) - 1 
+        c_1024  = zeros((1024, 1024)) - 1 
+        c_512  = zeros((512, 512)) - 1 
+        c_256  = zeros((256, 256)) - 1
+        c_128  = zeros((128, 128)) - 1
+        c_64   = zeros((64, 64)) - 1
+        c_32   = zeros((32, 32)) - 1
+
+        c_2048[:, :] += 2 * (sqrt((xy_2048 - L/2)**2 + (xy_2048.T - L/2)**2) <= .25) 
+        c_1024[:, :] += 2 * (sqrt((xy_1024 - L/2)**2 + (xy_1024.T - L/2)**2) <= .25) 
+        c_512[:, :] += 2 * (sqrt((xy_512 - L/2)**2 + (xy_512.T - L/2)**2) <= .25) 
+        c_256[:, :] += 2 * (sqrt((xy_256 - L/2)**2 + (xy_256.T - L/2)**2) <= .25) 
+        c_128[:, :] += 2 * (sqrt((xy_128 - L/2)**2 + (xy_128.T - L/2)**2) <= .25) 
+        c_64[:, :] += 2 * (sqrt((xy_64 - L/2)**2 + (xy_64.T - L/2)**2) <= .25) 
+        c_32[:, :] += 2 * (sqrt((xy_32 - L/2)**2 + (xy_32.T - L/2)**2) <= .25) 
+
+
+
+    c_2048= cahn_hilliard_spectral_lss(120, c_2048)
+    c_1024 = cahn_hilliard_spectral_lss(120, c_1024)
+    c_512 = cahn_hilliard_spectral_lss(120, c_512)
+    c_256 = cahn_hilliard_spectral_lss(120, c_256)
+    c_128 = cahn_hilliard_spectral_lss(120, c_128)
+    c_64 = cahn_hilliard_spectral_lss(120, c_64)
+    c_32 = cahn_hilliard_spectral_lss(120, c_32)
+
+    plt.figure()
+    plt.plot(xy_2048[:, 0], c_2048[2048//2, :])
+    plt.plot(xy_1024[:, 0], c_1024[1024//2, :])
+    plt.plot(xy_512[:, 0], c_512[512//2, :])
+    plt.plot(xy_256[:, 0], c_256[256//2, :])
+    plt.plot(xy_128[:, 0], c_128[128//2, :])
+    plt.plot(xy_64[:, 0], c_64[64//2, :])
+    plt.plot(xy_32[:, 0], c_32[32//2, :])
+
+
+    err = zeros((6,))
+    err[5] = inf_norm_error(c_1024, c_2048) 
+    err[4] = inf_norm_error(c_512, c_1024) 
+    err[3] = inf_norm_error(c_256, c_512) 
+    err[2] = inf_norm_error(c_128, c_256) 
+    err[1] = inf_norm_error(c_64, c_128) 
+    err[0] = inf_norm_error(c_32, c_64) 
+
+    print("=== accuracy test ===")
+    print(f"inf-norm_2048_512 = {inf_norm_error(c_1024, c_2048)}")
+    print(f"inf-norm_1024_512 = {inf_norm_error(c_512, c_1024)}")
+    print(f"inf-norm_512_256 = {inf_norm_error(c_256, c_512)}")
+    print(f"inf-norm_256_128 = {inf_norm_error(c_128, c_256)}")
+    print(f"inf-norm_128_64  = {inf_norm_error(c_64, c_128)}")
+    print(f"inf-norm_64_32   = {inf_norm_error(c_32, c_64)}")
+
+    fig, (ax1, ax2) = plt.subplots(1, 2)
+    fig.suptitle(f"Accuracy and Convergence")
+    ax1.semilogy([2**x for x in range(5, 11)], err, 'ko', base=10)
+    ax1.grid('on')
+    ax1.set_title("Accuracy")
+    ax2.semilogy([2**x for x in range(5, 10)], [err[i]/err[i+1] for i in range(5)], 'ko', base=2)
+    ax2.set_title("Convergence")
+    ax2.grid('on')
+
+    plt.figure()
+    plt.contourf(xy_2048[:, 0], xy_2048[:, 0], c_2048, cmap="jet");
+
+
+def efficiency_test(method, max_size):
+
+    i_max = int(log2(max_size))
+    c_max = 0.5 * random.rand(max_size, max_size) - 0.25
+    toc = zeros((i_max - 4,))
+
+    print("=== efficiency test ===")
+
+    for i in range(5, i_max + 1):
+        
+        tic = time.time() 
+        method(100, c_max[::2**(i_max - i), ::2**(i_max - i)])
+        toc[i - 5] = time.time() - tic
+
+        print(f"Elpased time for N = {2**i} is {toc[i - 5]} seconds")
+    
+    fig, (ax1) = plt.subplots()
+
+    nn = linspace(32**2, 1024**2, 1000)
+    ax1.plot([2**(2*x) for x in range(5, i_max + 1)], toc * 300000,  'k+',)
+    ax1.plot(nn, nn * log10(nn), '--k')
+    ax1.grid('on')
+        
+    
+def cahn_hilliard_1D_test(): 
+
+    n = 128;
+    dx = L / (n - 1)
+    it = 12000;
+    a = 0.01
+
+    x = linspace(0, 1, n);
+
+    c = zeros((n,)) - 1
+    c += 2 * (abs(x - L/2) <= .25) 
+
+    # wavenumbers
+    kr = power(2 * pi * fft.rfftfreq(n, d=dx), 2)
+
+    for _ in range(it): 
+       
+        s = fft.rfft(c) - dt * kr * fft.rfft(c**3 - 3 * c) 
+        v = s / (1 + dt * (2 * kr + a**2 * kr**2))
+        c = fft.irfft(v).real 
+
+    
+    plt.figure(figsize=(4, 4))
+    plt.title("test")
+    plt.plot(x, c, 'k')
+    
+
 if __name__ == "__main__":
 
+    nn = 512 
+    it = 1200 
+    a = 0.01
+    L = 1
+    dt = 1E-6
+
+    c_0 = 0.5 * random.rand(nn, nn) - 0.25
+    c_1 = c_0[::2, ::2]
+    
+    # accuracy_test(smooth=True)    
+    efficiency_test(cahn_hilliard_spectral_lss, 1024)
+    # cahn_hilliard_1D_test()
     tic = time.time()
+
+    # xy, dx = linspace(0, 1, nn, retstep=True)
+    
+    # c1 = cahn_hilliard_spectral_lss(it, c_1)
+    # c2 = cahn_hilliard_spectral_lss(it, c_0)
     
-    c1 = cahn_hilliard_fdiff(it, c)
-    c2 = cahn_hilliard_spectral(it, c)
 
+    # tt, ddt, fef, c2 = cahn_hilliard_spectral_lss_ad(1, c)
+    """ 
     toc = time.time() - tic
-
+    print(allclose(c1, c2))
     print(f"Elapsed time is {toc} seconds.")
-    
     fig, (ax1, ax2) = plt.subplots(1, 2)
-    fig.suptitle(f"Cahn-Hilliard integration for {it} iterations ")
-    ax1.contourf(xy, xy, c1, cmap="jet", vmin=-1, vmax=1)
+    fig.suptitle(f"Cahn-Hilliard integration for {it} iterations")
+    
+    ax1.contourf(xy, xy, c2, cmap="turbo", vmin=-1, vmax=1)
+    ax1.set_aspect('equal', 'box')
     ax1.set_title("Finite differences method")
-    ax2.contourf(xy, xy, c2, cmap="jet", vmin=-1, vmax=1)
+    ax2.contourf(xy, xy, c2, cmap="turbo", vmin=-1, vmax=1)
+    ax2.set_aspect('equal', 'box')
     ax2.set_title("Spectral method")
+    """ 
+    """ 
+    plt.figure(1)
+    plt.contourf(xy, xy, c2, cmap="turbo", vmin=-1, vmax=1)
+
+    plt.figure(2)
+    plt.plot(tt, ddt, 'k')
+    plt.grid('on')
+
+    plt.figure(3)
+    plt.plot(tt, fef, 'k')
+    plt.grid('on')
+    """ 
+
     plt.show()

matrix_generator.py

+from numpy import *
+
+nn = 512 
+L = 1
+
+xy = linspace(0, 1, 512); 
+x = xy.reshape((512, 1))
+
+c = zeros((nn, nn)) - 1
+# c = 0.5 * random.rand(nn, nn) - 0.45     
+c[:, :] += 2 * (sqrt((x - L/2)**2 + (x.T - L/2)**2) <= 0.25) 
+
+c = c.flatten()
+
+
+print(shape(c))
+
+savetxt("initial_c_512_disc.txt", c, '% .18e', delimiter="", newline="\n")
+