# CAPACITE 2 : EQUATION DE LA DIFFUSION THERMIQUE
##################################################################################################################
import numpy as np
import matplotlib.pyplot as plt
#################################################################################################################
# Exemple 1 : Résolution de l'équation de la diffusion thermique avec températures imposées aux extrémités 
#(conditions de  Dirichlet)

# Ecriture de la fonction de résolution de l'équation de la diffusion thermique

def Equationdiffusion(M,N,Tini,T0,TL):    

# Grandeurs d'entrée :
    # M : nombre d'intervalles de la plage de temps [0,tmax] (donc M+1 valeurs du temps t)
    # N : nombre d'intervalles de la plage d'abscisse [0,L]  (donc N+1 valeurs de l'abscisses x)
    # Tini : température à t = 0 pour tout x  (condition initiale)
    # T0 : température en x = 0  (condition aux limites en x = 0) 
    # TL : température en x = L  (condition aux limites en x = L)        
  
    T = np.zeros((M+1,N+1))  # Initialisation du tableau des températures T(ti,xk)
    T[0,:] = Tini            # Première ligne (indice i = 0) : température initiale des points xk    
    
    for i in range (0,M):    # Itération en temps ti = idt  (M valeurs de 0 à M-1)
    
        T[i+1,0] = T0          # Condition aux limites en x = 0 
        T[i+1,N] = TL          # Condition aux limites en x = L
        for k in range(1,N):   # Pour chaque abscisse de x1 à xN-1 (N-1 valeurs de k de 1 à N-1)
            # Discrétisation équation de la diffusion
            T[i+1,k] = T[i,k] + C*(T[i,k+1] + T[i,k-1] - 2*T[i,k])     
   
    return T
   
# Paramètres physiques

lambd = 1.65                    # Conductivité thermique du mur        (W.m-1.K-1)
µ = 2150                        # Masse volumique du mur               (kg.m-3)  
c = 1000                        # capacité thermique massique du mur   (J.K-1.kg-1) 
D = lambd/(µ*c)                 # Diffusivité thermique du mur         (m2.s-1)
L = 0.40                        # Epaisseur du mur                     (m)
tmax = 86400                    # Durée de l'étude                     (s) 

# Paramètres numériques

N = 100                         # Nombre d'intervalles en abscisses 
dx = L/N                        # Pas en abscisse 
x = np.linspace(0,L,N+1)        # Plage d'abscisse

M = 8640                        # Nombre d'intervalles de temps 
dt = tmax/M                     # Pas de temps 

# Conditions initiales et aux limites

T0 = 20                         # Conditions aux limites en x = 0      (°C) 
TL = 10                         # Condition aux limites en x = L       (°C)
Tini= 10                        # Température initiale                 (°C) 

# Constante intervenant dans l'équation de diffusion discrétisée

C = D*dt/dx**2                 

if C < 0.5:        # Test sur la condition de convergence de l'algorithme

    T = Equationdiffusion(M,N,Tini,T0,TL)       # Résolution de l'équation de la diffusion 
                                                          
    plt.xlabel('x')
    plt.ylabel('T', rotation=0)
    plt.title('Equation de la diffusion thermique 1D - Exemple 1')

    # Tableau des itérations choisies pour le tracé de la température à différents instants
    it = np.array([0,int(M/100),int(M/10),int(2*M/10),int(4*M/10),int(6*M/10),int(8*M/10),M]) 
    t = it*dt                                 # Instants associés aux itérations choisies

    for j in range(0,len(it)):
        plt.plot(x,T[it[j],:],label='t={} s'.format(t[j]))

    plt.grid()
    plt.legend()
    plt.show()

else:
    print("Les pas sont mal choisis car la constante C est supérieure à 0,5")
##################################################################################################################