#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov  3 20:01:44 2021

@author: maximebouchereau
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

# Résolution de l'équation de la chaleur en une dimension d'espace: dyke & encaissant

# Paramètres physiques [à ajuster]

L = 20 # Longueur du domaine

Ldyke = 2 # Longueur du dyke
Tdyke = 1200 # Température du dyke
kdyke = 3 # Conductivité thermique du dyke
rhodyke = 2900 # Masse volumique du dyke
cdyke = 1000 # Capacité thermique massique du dyke

T = 300 # Température de l'encaissant
k = 1.5 # Conductivité thermique de l'encaissant
rho = 2400 # Masse volumique de l'encaissant
c = 1000 # Capacité thermique massique de l'encaissant

tau = 1.577E7 # Durée de la simulation (en s - égal à 6 mois)



# Paramètres de la simulation numérique [à ajuster]

J = 200 # Nombre de subdivisions en espace
N = 200 # Nombre de subdivisions en temps



# Paramètres calculés

domaine = np.linspace(0,L,J+1) # Discrétisation du domaine
Ddyke = kdyke/(rhodyke*cdyke) # Coefficient de diffusion thermique du dyke
D = k/(rho*c) # Coefficient de diffusion thermique de l'encaissant

Diff = np.zeros(J+1) # Vecteur contenant les coefficients de diffusion en chaque point du domaine

for j in range(J+1):
    if abs(domaine[j]-L/2) <= Ldyke/2 : # Test de vérification si on se trouve sur le dyke
        Diff[j] = Ddyke
    else: # C'est que l'on se trouve sur l'encaissant
        Diff[j] = D
        
MatDiff = np.diag(Diff) # Matrice contenant les coefficients de diffusion en fonction de l'endroit

hx = L/J # Pas de subdivision en espace
ht = tau/N # Pas de subdivision en temps



# Construction des matrices impliquées dans le schéma numérique

# A est la matrice qui va imiter le laplacien en 1D (avec conditions de Neumann)
A = np.diag(np.ones(J),-1) + np.diag(np.ones(J),1) - np.diag(2*np.ones(J+1),0)
A[0,0] = -1
A[J,J] = -1

# M est la matrice qui sera utilisée à chaque itération dans le schéma numérique

I = np.eye(J+1,J+1) # Matrice identité
M = np.linalg.inv(I-(ht/hx**2)*MatDiff@A) # Inversion à cause du schéma implicite



# Calcul de la solution approchée de l'équation de la chaleur

U = np.zeros((J+1,N+1)) # Collection de vecteurs contenant les solutions aux temps t_n (attention aux doubles parenthèses !!!)

# Construction de la condition initiale

for j in range(J+1):
    if abs(domaine[j]-10) <= Ldyke/2: # Test de vérification si on se trouve sur le dyke
        U[j,0] = Tdyke
    else:
        U[j,0] = T
    
    

# Calcul de la solution à chaque temps

for n in range(N):
    U[:,n+1] = M@(U[:,n])
    
    

# Calcul de la longueur ayant dépassé les 600°C durant le processus

nb600 = 0 # Compte le nombre de points où la température a dépassé les 600°C dans le domaine discrétisé

for j in range(J+1):
    if max(U[j,:]) >= 600:
        nb600 = nb600+1

print("Longueur de l'auréole métamorphique:",format(nb600/(J+1)*L - Ldyke,'.2f'),"m") # On retire Ldyke, car ce domaine n'est pas contitué d'encaissant



# Tracé de la solution au bout de six mois

plt.figure()
plt.title("Profil de température au bout de six mois [°C]")
plt.xlabel("Domaine")
plt.ylabel("Température")
plt.plot(domaine,U[:,N],color="red") # Tracé de la solution au bout de six mois
plt.grid()
plt.show()


# Tracé de la solution au cours du temps (bonus)

plt.figure()
plt.title("Evolution de la température au cours du temps [°C]")
plt.xlabel("temps (en mois)")
plt.ylabel("Domaine")
plt.imshow(U,cmap="jet",aspect="auto",extent=(0,6,0,L))
plt.colorbar()
plt.show()

# Animation du tracé de la solution au cours du temps (bonus++)

deltat = 5000/N # Durée entre deux frames (en ms)

fig = plt.figure()
line, = plt.plot([],[],color="red")
plt.xlabel("Domaine")
plt.ylabel("température")
plt.grid()
plt.xlim(0,L)
plt.ylim(T,Tdyke)
plt.title("Evolution de la température")

def animate(n):
    line.set_data(domaine,U[:,n])
    return line,

ani = animation.FuncAnimation(fig, animate, frames=N, blit=True, interval=deltat, repeat=True)
#%matplotlib qt
plt.show()