Skip to content
Extraits de code Groupes Projets
main.py 11 ko
Newer Older
  • Learn to ignore specific revisions
  • Simon Collignon's avatar
    Simon Collignon a validé
    import time
    from numpy import *
    import matplotlib.pyplot as plt
    
    Simon Collignon's avatar
    Simon Collignon a validé
    from scipy.fftpack import dctn, idctn
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
    # 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."""
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        nn = len(c[0, :])
        dx = L / (nn - 1)
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        cp = zeros((nn + 2, nn + 2))
        f = zeros((nn + 2, nn + 2))
        cp[1:-1, 1:-1] = c
    
        for _ in range(it):
            
    
    Simon Collignon's avatar
    Simon Collignon a validé
            periodize(cp)
    
    Simon Collignon's avatar
    Simon Collignon a validé
            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])
    
    Simon Collignon's avatar
    Simon Collignon a validé
            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]
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
    
    def cahn_hilliard_spectral(it, c):
        """Cahn-Hilliard integration for it iterations using fft spectral method."""
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        nn = len(c[0, :])
        dx = L / (nn - 1)
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        # wavenumbers
    
        kr = power(2 * pi * fft.rfftfreq(nn, d=dx), 2)
        kc = power(2 * pi * fft.fftfreq(nn, d=dx), 2)
    
    Simon Collignon's avatar
    Simon Collignon a validé
        kk = array([kr,]*nn) + array([kc,]*(nn//2+1)).transpose()
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
        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
    
    Simon Collignon's avatar
    Simon Collignon a validé
    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')
        
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
    if __name__ == "__main__":
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        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()
    
    Simon Collignon's avatar
    Simon Collignon a validé
        tic = time.time()
    
    Simon Collignon's avatar
    Simon Collignon a validé
    
        # xy, dx = linspace(0, 1, nn, retstep=True)
        
        # c1 = cahn_hilliard_spectral_lss(it, c_1)
        # c2 = cahn_hilliard_spectral_lss(it, c_0)
    
    Simon Collignon's avatar
    Simon Collignon a validé
        # tt, ddt, fef, c2 = cahn_hilliard_spectral_lss_ad(1, c)
        """ 
    
    Simon Collignon's avatar
    Simon Collignon a validé
        toc = time.time() - tic
    
    Simon Collignon's avatar
    Simon Collignon a validé
        print(allclose(c1, c2))
    
    Simon Collignon's avatar
    Simon Collignon a validé
        print(f"Elapsed time is {toc} seconds.")
    
        fig, (ax1, ax2) = plt.subplots(1, 2)
    
    Simon Collignon's avatar
    Simon Collignon a validé
        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")
    
    Simon Collignon's avatar
    Simon Collignon a validé
        ax2.contourf(xy, xy, c2, cmap="turbo", vmin=-1, vmax=1)
        ax2.set_aspect('equal', 'box')
    
        ax2.set_title("Spectral method")
    
    Simon Collignon's avatar
    Simon Collignon a validé
        """ 
        """ 
        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')
        """ 
    
    
    Simon Collignon's avatar
    Simon Collignon a validé
        plt.show()