commit 57c9b1e011f4aac7057eace9f0b31bb2271517e1 Author: mbui Date: Tue May 30 00:11:10 2023 +0200 TP2 Intro Modelisation diff --git a/Fig/ChampT_DF_1.png b/Fig/ChampT_DF_1.png new file mode 100644 index 0000000..598bb0d Binary files /dev/null and b/Fig/ChampT_DF_1.png differ diff --git a/Fig/ChampT_DF_10.png b/Fig/ChampT_DF_10.png new file mode 100644 index 0000000..4b549d3 Binary files /dev/null and b/Fig/ChampT_DF_10.png differ diff --git a/Fig/ChampT_DF_11.png b/Fig/ChampT_DF_11.png new file mode 100644 index 0000000..34868aa Binary files /dev/null and b/Fig/ChampT_DF_11.png differ diff --git a/Fig/ChampT_DF_12.png b/Fig/ChampT_DF_12.png new file mode 100644 index 0000000..8ce3e35 Binary files /dev/null and b/Fig/ChampT_DF_12.png differ diff --git a/Fig/ChampT_DF_2.png b/Fig/ChampT_DF_2.png new file mode 100644 index 0000000..ccd2616 Binary files /dev/null and b/Fig/ChampT_DF_2.png differ diff --git a/Fig/ChampT_DF_3.png b/Fig/ChampT_DF_3.png new file mode 100644 index 0000000..5d11269 Binary files /dev/null and b/Fig/ChampT_DF_3.png differ diff --git a/Fig/ChampT_DF_4.png b/Fig/ChampT_DF_4.png new file mode 100644 index 0000000..dda38c8 Binary files /dev/null and b/Fig/ChampT_DF_4.png differ diff --git a/Fig/ChampT_DF_5.png b/Fig/ChampT_DF_5.png new file mode 100644 index 0000000..e99bbad Binary files /dev/null and b/Fig/ChampT_DF_5.png differ diff --git a/Fig/ChampT_DF_6.png b/Fig/ChampT_DF_6.png new file mode 100644 index 0000000..97230c6 Binary files /dev/null and b/Fig/ChampT_DF_6.png differ diff --git a/Fig/ChampT_DF_7.png b/Fig/ChampT_DF_7.png new file mode 100644 index 0000000..4056bd8 Binary files /dev/null and b/Fig/ChampT_DF_7.png differ diff --git a/Fig/ChampT_DF_8.png b/Fig/ChampT_DF_8.png new file mode 100644 index 0000000..d192d4a Binary files /dev/null and b/Fig/ChampT_DF_8.png differ diff --git a/Fig/ChampT_DF_9.png b/Fig/ChampT_DF_9.png new file mode 100644 index 0000000..f9ffde2 Binary files /dev/null and b/Fig/ChampT_DF_9.png differ diff --git a/Fig/ChampT_MMC_1.png b/Fig/ChampT_MMC_1.png new file mode 100644 index 0000000..1b45cae Binary files /dev/null and b/Fig/ChampT_MMC_1.png differ diff --git a/Fig/ChampT_MMC_10.png b/Fig/ChampT_MMC_10.png new file mode 100644 index 0000000..27cb845 Binary files /dev/null and b/Fig/ChampT_MMC_10.png differ diff --git a/Fig/ChampT_MMC_11.png b/Fig/ChampT_MMC_11.png new file mode 100644 index 0000000..aab7b6a Binary files /dev/null and b/Fig/ChampT_MMC_11.png differ diff --git a/Fig/ChampT_MMC_12.png b/Fig/ChampT_MMC_12.png new file mode 100644 index 0000000..45a4f05 Binary files /dev/null and b/Fig/ChampT_MMC_12.png differ diff --git a/Fig/ChampT_MMC_2.png b/Fig/ChampT_MMC_2.png new file mode 100644 index 0000000..6391d09 Binary files /dev/null and b/Fig/ChampT_MMC_2.png differ diff --git a/Fig/ChampT_MMC_3.png b/Fig/ChampT_MMC_3.png new file mode 100644 index 0000000..d6e07e2 Binary files /dev/null and b/Fig/ChampT_MMC_3.png differ diff --git a/Fig/ChampT_MMC_4.png b/Fig/ChampT_MMC_4.png new file mode 100644 index 0000000..81f320d Binary files /dev/null and b/Fig/ChampT_MMC_4.png differ diff --git a/Fig/ChampT_MMC_5.png b/Fig/ChampT_MMC_5.png new file mode 100644 index 0000000..5e0c113 Binary files /dev/null and b/Fig/ChampT_MMC_5.png differ diff --git a/Fig/ChampT_MMC_6.png b/Fig/ChampT_MMC_6.png new file mode 100644 index 0000000..92ae557 Binary files /dev/null and b/Fig/ChampT_MMC_6.png differ diff --git a/Fig/ChampT_MMC_7.png b/Fig/ChampT_MMC_7.png new file mode 100644 index 0000000..e6e96d6 Binary files /dev/null and b/Fig/ChampT_MMC_7.png differ diff --git a/Fig/ChampT_MMC_8.png b/Fig/ChampT_MMC_8.png new file mode 100644 index 0000000..b203126 Binary files /dev/null and b/Fig/ChampT_MMC_8.png differ diff --git a/Fig/ChampT_MMC_9.png b/Fig/ChampT_MMC_9.png new file mode 100644 index 0000000..e50d6a2 Binary files /dev/null and b/Fig/ChampT_MMC_9.png differ diff --git a/Fig/ProfilT_DF_025_1.png b/Fig/ProfilT_DF_025_1.png new file mode 100644 index 0000000..3e2d62c Binary files /dev/null and b/Fig/ProfilT_DF_025_1.png differ diff --git a/Fig/ProfilT_DF_025_2.png b/Fig/ProfilT_DF_025_2.png new file mode 100644 index 0000000..dac2e6e Binary files /dev/null and b/Fig/ProfilT_DF_025_2.png differ diff --git a/Fig/ProfilT_DF_025_3.png b/Fig/ProfilT_DF_025_3.png new file mode 100644 index 0000000..656dafa Binary files /dev/null and b/Fig/ProfilT_DF_025_3.png differ diff --git a/Fig/ProfilT_DF_025_4.png b/Fig/ProfilT_DF_025_4.png new file mode 100644 index 0000000..36a4e19 Binary files /dev/null and b/Fig/ProfilT_DF_025_4.png differ diff --git a/Fig/ProfilT_DF_050_1.png b/Fig/ProfilT_DF_050_1.png new file mode 100644 index 0000000..80d8105 Binary files /dev/null and b/Fig/ProfilT_DF_050_1.png differ diff --git a/Fig/ProfilT_DF_050_2.png b/Fig/ProfilT_DF_050_2.png new file mode 100644 index 0000000..82dde76 Binary files /dev/null and b/Fig/ProfilT_DF_050_2.png differ diff --git a/Fig/ProfilT_DF_050_3.png b/Fig/ProfilT_DF_050_3.png new file mode 100644 index 0000000..cdf1ca2 Binary files /dev/null and b/Fig/ProfilT_DF_050_3.png differ diff --git a/Fig/ProfilT_DF_050_4.png b/Fig/ProfilT_DF_050_4.png new file mode 100644 index 0000000..dbdd213 Binary files /dev/null and b/Fig/ProfilT_DF_050_4.png differ diff --git a/Fig/ProfilT_DF_075_1.png b/Fig/ProfilT_DF_075_1.png new file mode 100644 index 0000000..902321f Binary files /dev/null and b/Fig/ProfilT_DF_075_1.png differ diff --git a/Fig/ProfilT_DF_075_2.png b/Fig/ProfilT_DF_075_2.png new file mode 100644 index 0000000..22525ac Binary files /dev/null and b/Fig/ProfilT_DF_075_2.png differ diff --git a/Fig/ProfilT_DF_075_3.png b/Fig/ProfilT_DF_075_3.png new file mode 100644 index 0000000..94862f0 Binary files /dev/null and b/Fig/ProfilT_DF_075_3.png differ diff --git a/Fig/ProfilT_DF_075_4.png b/Fig/ProfilT_DF_075_4.png new file mode 100644 index 0000000..0b321bd Binary files /dev/null and b/Fig/ProfilT_DF_075_4.png differ diff --git a/Fig/ProfilT_MMC_025_1.png b/Fig/ProfilT_MMC_025_1.png new file mode 100644 index 0000000..33fa7f7 Binary files /dev/null and b/Fig/ProfilT_MMC_025_1.png differ diff --git a/Fig/ProfilT_MMC_025_2.png b/Fig/ProfilT_MMC_025_2.png new file mode 100644 index 0000000..c93b554 Binary files /dev/null and b/Fig/ProfilT_MMC_025_2.png differ diff --git a/Fig/ProfilT_MMC_025_3.png b/Fig/ProfilT_MMC_025_3.png new file mode 100644 index 0000000..6abacab Binary files /dev/null and b/Fig/ProfilT_MMC_025_3.png differ diff --git a/Fig/ProfilT_MMC_050_1.png b/Fig/ProfilT_MMC_050_1.png new file mode 100644 index 0000000..286cc05 Binary files /dev/null and b/Fig/ProfilT_MMC_050_1.png differ diff --git a/Fig/ProfilT_MMC_050_2.png b/Fig/ProfilT_MMC_050_2.png new file mode 100644 index 0000000..4f5b369 Binary files /dev/null and b/Fig/ProfilT_MMC_050_2.png differ diff --git a/Fig/ProfilT_MMC_050_3.png b/Fig/ProfilT_MMC_050_3.png new file mode 100644 index 0000000..061eb87 Binary files /dev/null and b/Fig/ProfilT_MMC_050_3.png differ diff --git a/Fig/ProfilT_MMC_075_1.png b/Fig/ProfilT_MMC_075_1.png new file mode 100644 index 0000000..de073d8 Binary files /dev/null and b/Fig/ProfilT_MMC_075_1.png differ diff --git a/Fig/ProfilT_MMC_075_2.png b/Fig/ProfilT_MMC_075_2.png new file mode 100644 index 0000000..2e839d5 Binary files /dev/null and b/Fig/ProfilT_MMC_075_2.png differ diff --git a/Fig/ProfilT_MMC_075_3.png b/Fig/ProfilT_MMC_075_3.png new file mode 100644 index 0000000..cd70776 Binary files /dev/null and b/Fig/ProfilT_MMC_075_3.png differ diff --git a/__pycache__/tp2_lib.cpython-39.pyc b/__pycache__/tp2_lib.cpython-39.pyc new file mode 100644 index 0000000..755dc80 Binary files /dev/null and b/__pycache__/tp2_lib.cpython-39.pyc differ diff --git a/__pycache__/tp2_methodes_lib.cpython-39.pyc b/__pycache__/tp2_methodes_lib.cpython-39.pyc new file mode 100644 index 0000000..c8a8f91 Binary files /dev/null and b/__pycache__/tp2_methodes_lib.cpython-39.pyc differ diff --git a/tp2.py b/tp2.py new file mode 100644 index 0000000..8525ca1 --- /dev/null +++ b/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() +''' diff --git a/tp2_lib.py b/tp2_lib.py new file mode 100644 index 0000000..a422022 --- /dev/null +++ b/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 + + + + + + + diff --git a/tp2_methodes_lib.py b/tp2_methodes_lib.py new file mode 100644 index 0000000..a64db96 --- /dev/null +++ b/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 \ No newline at end of file