TP2 Intro Modelisation

tp2.py

@@ -0,0 +1,253 @@
+
+#######################################################################
+##                           TP2 Modélisation MIC3 - GMM             ##
+#######################################################################
+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from math import *
+from copy import deepcopy
+import tp2_lib as tp
+import tp2_methodes_lib as meth
+
+
+######################################################################
+## Finite Difference Method                                         ##
+######################################################################
+
+##Initialize Data
+L = 1. 
+List_K = [10, 20, 30]
+List_tmax = [0.5, 1., 1.5, 2.]
+num1 = 0
+num2 = 0
+T0DF = np.zeros([4,11])
+T1DF = np.zeros([4,11])
+T2DF = np.zeros([4,11])
+yDF  = np.zeros(11)
+
+for tmax in List_tmax:
+    num1 = num1 + 1
+    num3 = -1
+    for K in List_K :
+        num2 = num2 + 1
+        num3 = num3 + 1
+        
+        ## Compute T
+        T = meth.methodeDF(K,tmax,L)
+
+        ## Convert vector T into a (K+1)x(K+1) Matrix
+        T2D = tp.matT(T,K,tp.indk)
+        ## print(T2D)
+
+        ## Maillage
+        x = np.linspace(0,L,K+1)
+        y = np.linspace(0,L,K+1)
+        X, Y = np.meshgrid(x, y)
+
+        ## Temperature profile extraction at x = L/4
+        i0m = floor(K/4.)
+        i0p = i0m+1
+        omega = K/4. - i0m
+        T0 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+        if num3 == 0: 
+            T00 = deepcopy(T0)
+            y0 = deepcopy(y)
+            ## Save T0 in T0DF to serve as a reference solution for MMC 
+            T0DF[num1-1,:] = T0
+            yDF = y0
+        elif num3 == 1:
+            T01 = deepcopy(T0)
+            y1 = deepcopy(y)
+        else :
+            T02 = deepcopy(T0)
+            y2 = deepcopy(y)
+
+        ## Temperature profile extraction at x = L/2
+        i0m = floor(K/2.)
+        i0p = i0m+1
+        omega = K/2. - i0m
+        T1 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+        if num3 == 0: 
+            T10 = deepcopy(T1)
+            ## Save T1 in T1DF to serve as a reference solution for MMC 
+            T1DF[num1-1,:] = T1
+        elif num3 == 1:
+            T11 = deepcopy(T1)
+        else :
+            T12 = deepcopy(T1)
+
+        ## Temperature profile extraction at x = 3L/4
+        i0m = floor(3*K/4.)
+        i0p = i0m+1
+        omega = 3*K/4. - i0m
+        T2 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+        if num3 == 0: 
+            T20 = deepcopy(T2)
+            ## Save T2 in T2DF to serve as a reference solution for MMC 
+            T2DF[num1-1,:] = T2
+        elif num3 == 1:
+            T21 = deepcopy(T2)
+        else :
+            T22 = deepcopy(T2)
+
+        ## Temperature field  
+        fig=plt.figure("ChampT_DF_"+str(num2))
+        NbIso = 25
+        h = L/K
+        CF = plt.contourf(X, Y, np.transpose(T2D), NbIso)
+        ##plt.clabel(CF, colors = 'k', fmt = '%2.1f', fontsize=12)
+        plt.colorbar(CF)
+        plt.title('Temperature field at t = ' + str(tmax) + ' for h = ' +  str(h))
+        plt.xlabel('x')
+        plt.ylabel('y')
+        #plt.savefig("ChampT_DF_"+str(num2)+".png")
+        plt.show()
+
+    fig=plt.figure("ProfilT_DF_025_"+str(num1))
+    plt.plot(y0,T00, '-bo', y1, T01, '-rs', y2, T02, '-gv')
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = L/4')
+    plt.legend(['h=L/10', 'h=L/20', 'h=L/30'], loc='best' )
+    #plt.savefig("ProfilT_DF_025_"+str(num1)+".png")
+    plt.show()
+
+    fig=plt.figure("ProfilT_DF_050_"+str(num1))
+    plt.plot(y0,T10, '-bo', y1, T11, '-rs', y2, T12, '-gv')
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = L/2')
+    plt.legend(['h=L/10', 'h=L/20', 'h=L/30'], loc='best' )
+    #plt.savefig("ProfilT_DF_050_"+str(num1)+".png")
+    plt.show()
+    
+    fig=plt.figure("ProfilT_DF_075_"+str(num1))
+    plt.plot(y0,T20, '-bo', y1, T21, '-rs', y2, T22, '-gv')
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = 3L/4')
+    plt.legend(['h=L/10', 'h=L/20', 'h=L/30'], loc='best' )
+    #plt.savefig("ProfilT_DF_075_"+str(num1)+".png")
+    plt.show()
+
+######################################################################
+## Monte Carlo Method                                               ##
+######################################################################
+
+## Initialize Data
+L = 1. 
+#List_K = [10000, 100000]
+#List_M = [10,20]
+#List_tmax = [0.5, 1., 1.5]
+List_K=[100000]
+List_M=[20]
+List_tmax = [4,5,6]
+num1 = 0
+num2 = 0
+
+for tmax in List_tmax:
+    num1 = num1 + 1
+    num3 = -1
+    for K in List_K :
+        for M in List_M : 
+            num2 = num2 + 1
+            num3 = num3 + 1
+            
+            ## Compute T 
+            #T2D = meth.methodeMC(K,M,tmax,L)
+            T2D = meth.methodeMC2(K,M,tmax,L)
+
+            ## Maillage
+            eps = L/M
+            x = np.linspace(eps/2,L-eps/2,M)
+            y = np.linspace(eps/2,L-eps/2,M)
+            X, Y = np.meshgrid(x, y)
+            
+            ## Temperature profile extraction at x = L/4
+            i0m = floor(M/4.)
+            i0p = i0m+1
+            omega = M/4. - i0m
+            T0 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+            if num3 == 0: 
+                T00mc = deepcopy(T0)
+                y0mc = deepcopy(y)
+            elif num3 == 1:
+                T01mc = deepcopy(T0)
+                y1mc = deepcopy(y)
+            elif num3 == 2:
+                T02mc = deepcopy(T0)
+                y2mc = deepcopy(y)
+            else :
+                T03mc = deepcopy(T0)
+                y3mc = deepcopy(y)
+
+            ## Temperature profile extraction at x = L/2
+            i0m = floor(M/2.)
+            i0p = i0m+1
+            omega = M/2. - i0m
+            T1 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+            if num3 == 0: 
+                T10mc = deepcopy(T1)
+            elif num3 == 1:
+                T11mc = deepcopy(T1)
+            elif num3 == 2:
+                T12mc = deepcopy(T1)
+            else :
+                T13mc = deepcopy(T1)
+
+            ## Temperature profile extraction at x = 3L/4
+            i0m = floor(3*M/4.)
+            i0p = i0m+1
+            omega = 3*M/4. - i0m
+            T2 = (1-omega)*T2D[i0m,:] + omega*T2D[i0p,:]
+            if num3 == 0: 
+                T20mc = deepcopy(T2)
+            elif num3 == 1:
+                T21mc = deepcopy(T2)
+            elif num3 == 2:
+                T22mc = deepcopy(T2)
+            else :
+                T23mc = deepcopy(T2)
+
+            ## Temperature field  
+            fig=plt.figure("ChampT_MMC_"+str(num2))
+            NbIso = 25
+            h = L/M
+            CF = plt.contourf(X, Y, np.transpose(T2D), NbIso)
+            plt.colorbar(CF)
+            #plt.title('Temperature field at t = ' + str(tmax) + ' for M = ' +  str(M) + ' and K = ' + str(K))
+            #plt.title('Temperature field at t = ' + str(tmax) + ' for $\Delta t = 0.1\epsilon^2/D$')
+            plt.title('Temperature field at t = ' + str(tmax) + ' for optimize Monte-Carlo')
+            plt.xlabel('x')
+            plt.ylabel('y')
+            #plt.savefig("ChampT_MMC_"+str(num2)+".png")
+'''
+    fig=plt.figure("ProfilT_MMC_025"+str(num1))
+    plt.plot(yDF,T0DF[num1-1,:], '--k', y0mc, T00mc, '-ro', y1mc, T01mc, '-rv', y2mc, T02mc, '-go',  y3mc, T03mc, '-gv' )
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = L/4')
+    plt.legend(['Finite Difference', 'K=10000 , $\epsilon$=1/10', 'K=10000 , $\epsilon$=1/20', 'K=100000 , $\epsilon$=1/10', 'K=100000 , $\epsilon$=1/20'], loc='best' )
+    #plt.savefig("ProfilT_MMC_025_"+str(num1)+".png")
+    plt.show()
+    
+    fig=plt.figure("ProfilT_MMC_050"+str(num1))
+    plt.plot(yDF,T1DF[num1-1,:], '--k', y0mc, T10mc, '-ro', y1mc, T11mc, '-rv', y2mc, T12mc, '-go',  y3mc, T13mc, '-gv' )
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = L/2')
+    plt.legend(['Finite Difference', 'K=10000 , $\epsilon$=1/10', 'K=10000 , $\epsilon$=1/20', 'K=100000 , $\epsilon$=1/10', 'K=100000 , $\epsilon$=1/20'], loc='best' )
+    #plt.savefig("ProfilT_MMC_050_"+str(num1)+".png")
+    plt.show()
+    
+    fig=plt.figure("ProfilT_MMC_075"+str(num1))
+    plt.plot(yDF,T2DF[num1-1,:], '--k', y0mc, T20mc, '-ro', y1mc, T21mc, '-rv', y2mc, T22mc, '-go',  y3mc, T23mc, '-gv' )
+    plt.xlabel('y')
+    plt.ylabel('T')
+    plt.title('Temperature Profile at t = ' + str(tmax) + ' and x = 3L/4')
+    plt.legend(['Finite Difference', 'K=10000 , $\epsilon$=1/10', 'K=10000 , $\epsilon$=1/20', 'K=100000 , $\epsilon$=1/10', 'K=100000 , $\epsilon$=1/20'], loc='best' )
+    #plt.savefig("ProfilT_MMC_075_"+str(num1)+".png")
+    plt.show()
+'''

tp2_lib.py

@@ -0,0 +1,156 @@
+import numpy as np
+import scipy as sp
+import scipy.sparse as spsp    ## Sparse matrix 
+import matplotlib as mpl        
+import matplotlib.pyplot as plt
+from math import * 
+
+def indk(i,j, K): #indice k corres à (i,j)
+    ## to be completed 
+    k=i+(K+1)*j
+    return k
+
+def vit(x,y,V0,L): #composanté de la vitesse
+    V1 = -V0*np.sin(np.pi*x/L)*np.cos(np.pi*y/L)
+    V2 = V0*np.cos(np.pi*x/L)*np.sin(np.pi*y/L)
+    V = np.array([V1,V2]) # shape (2,n)
+    ## to be completed  
+    return V
+
+def source(x,y,L): #fonction f
+    S = 256*(x/L)**2*(y/L)**2*(1-x/L)**2*(1-y/L)**2
+    ## to be completed  
+    return S
+
+def MatA(L, V0, D, h, dt, K, indk, vit):
+    K2 = (K+1)**2
+    A = spsp.lil_matrix((K2, K2))   ## Declaration of A as a sparse matrix             
+## Loop on the internal nodes
+    for i in range(1,K):
+        for j in  range(1,K):
+            k = indk(i,j,K)
+            ke = indk(i+1,j,K)
+            ko = indk(i-1,j,K)
+            kn = indk(i,j+1,K)
+            ks = indk(i,j-1,K)
+            ##to be completed 
+            V1,V2 = vit(i*h,j*h,V0,L)[:]
+            A[k,k] = 1 - 4*dt*D/h**2
+            A[k,ke] = -dt*V1/(2*h) + dt*D/h**2 
+            A[k,ko] = dt*V1/(2*h) + dt*D/h**2
+            A[k,kn] = -dt*V2/(2*h) + dt*D/h**2 
+            A[k,ks] = dt*V2/(2*h) + dt*D/h**2
+ ## Loop on the nodes located in x = 0 (Neumann Boundary Condition) 
+    for j in range(1,K):
+        ## to be completed
+        k = indk(0,j,K)
+        ke = indk(1,j,K)
+        #i=0 => non ko
+        kn = indk(0,j+1,K)
+        ks = indk(0,j-1,K)
+        V1,V2 = vit(0,j*h,V0,L)[:]
+        A[k,k] = 1 - 3*dt*D/h**2
+        A[k,ke] = dt*D/h**2 
+        A[k,kn] = -dt*V2/(2*h) + dt*D/h**2 
+        A[k,ks] = dt*V2/(2*h) + dt*D/h**2
+ ## Loop on the nodes located in x = L (Neumann Boundary Condition) 
+    for j in range(1,K):
+        ## to be completed 
+        k = indk(K,j,K)
+        #i=K => non ke
+        ko = indk(K-1,j,K)
+        kn = indk(K,j+1,K)
+        ks = indk(K,j-1,K)
+        V1,V2 = vit(K*h,j*h,V0,L)[:]
+        A[k,k] = 1 - 3*dt*D/h**2
+        A[k,ko] = dt*D/h**2 
+        A[k,kn] = -dt*V2/(2*h) + dt*D/h**2 
+        A[k,ks] = dt*V2/(2*h) + dt*D/h**2
+    return A
+def secmb(K, dt, h, L, Tbas, Thaut, indk, source): #second member S
+    ## to be completed
+    K2 = (K+1)**2
+    S = np.zeros(K2)
+    for i in range(K+1):
+        for j in range(1,K): #S[indk(i,0,K) = 0]
+            k = indk(i,j,K)
+            S[k] = dt*source(i*h,j*h,L)
+    ## On the nodes located in y = L <=> j=K  
+        S[indk(i,K,K)] = 1
+    return S 
+           
+    
+def matT(T, K, indk):
+    T2D = np.zeros((K+1,K+1))
+    ## to be completed
+    '''for i in range (K+1):
+        for j in range (K+1):
+            T2D[j,i] = T[indk(i,j,K)]
+    '''
+    T2D = (T.reshape((K+1,K+1))).T
+    #print(T2D)
+    return T2D
+
+def vecT(T2D,K):
+    '''T = np.zeros((K+1)**2)
+    for i in range (K+1):
+        for j in range (K+1):
+            T[indk(i,j,K)] = T2D[i,j]
+    '''
+    T = T2D.reshape(-1)
+    #print(T)
+    return T2D
+
+def posinit(K, L):
+    ## to be completed
+    X0 = np.random.uniform(0,L,K)
+    Y0 = np.random.uniform(0,L,K)
+    return X0,Y0
+
+def evolution(X, Y, Theta, K, dt, L, D, V0, Tbas, Thaut, source, vit):
+    Vit = np.zeros([K,2])
+    S = np.zeros(K)
+    alpha = np.random.randn(K,2)
+    
+    """for k in range(0, K):
+        V = vit(X[k],Y[k],V0,L)
+        Vit[k,0]= V[0]
+        Vit[k,1]= V[1]
+        S[k] = source(X[k],Y[k],L)"""
+    ## to be completed
+    Vit = vit(X,Y,V0,L).T
+    S = source(X,Y,L)
+    
+    Xetoile = X + dt*Vit[:,0] + np.sqrt(2*D*dt)*alpha[:,0]
+    Yetoile = Y + dt*Vit[:,1] + np.sqrt(2*D*dt)*alpha[:,1]
+    ThetaEtoile = Theta + dt*S
+    
+    X = np.minimum(L,np.maximum(Xetoile,0))
+    Y = np.minimum(L,np.maximum(Yetoile,0))
+    
+    Theta = (0 < Yetoile)*(Yetoile < L)*(ThetaEtoile) + (Yetoile>=L)
+    
+    return X,Y,Theta
+
+
+def tmoy(X, Y, Theta, K, M, L):
+
+    Tm = np.zeros([M,M])    ## mean temperature in a cell 
+    nbp = np.zeros([M,M])   ## number of particles in a cell 
+    e = L/M
+    
+    for k in range(0, K):
+        ## to be completed 
+        i = min(floor(X[k]/e),M-1)
+        j = min(floor(Y[k]/e),M-1)
+        Tm[i,j] += Theta[k]
+        nbp[i,j] += 1
+    Tm = (nbp != 0)*Tm/nbp
+    return Tm
+
+    
+    
+
+
+    
+    

tp2_methodes_lib.py

@@ -0,0 +1,112 @@
+import numpy as np
+import scipy as sp
+import scipy.sparse as spsp  
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from math import * 
+import tp2_lib as tp
+
+def methodeDF(K,tmax,L):
+#######################################################################
+##                 Finite Difference Method                          ##
+#######################################################################
+
+    ## Continuous Problem Data
+    V0 = 1.0
+    D = 0.2
+    Tbas = 0.   ## (BC en y = 0)
+    Thaut = 1.  ## (BC en y = L)
+
+
+    ## Numerical parameters
+    K2 = (K+1)**2 
+    h = L/K
+    dt = 0.8 * h*h / (4*D)
+    
+    # Initialisation of T
+    T = np.zeros(K2) #(K2,)
+
+    # Assembling of the Rigth Hand Side
+    S = tp.secmb(K, dt, h, L, Tbas, Thaut, tp.indk, tp.source) #(K2,)
+    #S = S.reshape((K2,1))
+    #print(np.shape(S))
+
+    # Assembling of the Scheme Matrix 
+    A = tp.MatA(L, V0, D, h, dt, K, tp.indk, tp.vit) #(K2,K2)
+
+    # Time Loop 
+    time = 0.0
+    while time < tmax :
+        ## to be completed
+        time = time+dt
+        #print(np.shape(A),np.shape(T),np.shape(S))
+        T = A@T + S
+    return T
+
+
+
+def methodeMC(K,M,tmax,L):
+#######################################################################
+##                       Monte Carlo Method                          ##
+#######################################################################
+
+    ## Continuous Problem Data
+    V0 = 1.0
+    D = 0.2
+    Tbas = 0.   ## (BC en y = 0)
+    Thaut = 1.  ## (BC en y = L)
+    
+    # Initialisation of Theta and X,Y
+    X,Y = tp.posinit(K,L)
+    Theta = np.zeros(K)
+    ## Numerical parameters
+    eps = L/M
+
+    dt = 0.25 * eps * eps / D 
+    #dt = 10 * eps * eps / D
+    #dt = 1 * eps * eps / D 
+    #dt = 0.1 * eps * eps / D 
+    # Time Loop 
+    time = 0.0
+    while time < tmax :
+        ## to be completed
+        time = time+dt
+        (X,Y,Theta) = tp.evolution(X, Y, Theta, K, dt, L, D, V0, Tbas, Thaut, tp.source, tp.vit)
+    # Compute Temperaure field
+    T = tp.tmoy(X, Y, Theta, K, M, L)
+    return T
+
+
+def methodeMC2(K,M,tmax,L):
+#######################################################################
+##                      Optimize Monte Carlo Method                          ##
+#######################################################################
+
+    ## Continuous Problem Data
+    V0 = 1.0
+    D = 0.2
+    Tbas = 0.   ## (BC en y = 0)
+    Thaut = 1.  ## (BC en y = L)
+    
+    # Initialisation of Theta and X,Y
+    X,Y = tp.posinit(K,L)
+    Theta = np.zeros(K)
+    ## Numerical parameters
+    eps = L/M
+    dt = 0.25 * eps * eps / D 
+    # Time Loop 
+    time = 0.0
+    i = 0 #times for T >= 3
+    TSum = np.zeros([M,M])
+    while time < tmax :
+        time = time+dt
+        (X,Y,Theta) = tp.evolution(X, Y, Theta, K, dt, L, D, V0, Tbas, Thaut, tp.source, tp.vit)
+        if time >= 3:
+            i +=1
+            TSum += tp.tmoy(X, Y, Theta, K, M, L)
+    # Compute Temperaure field
+    if i==0:
+        T = tp.tmoy(X, Y, Theta, K, M, L)
+    else:
+        T = TSum/i
+    return T