#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 18 10:48:03 2022

@author: adomps
"""

######################################################################
#
##  Résolution de l'éq de la chaleur 1d par méthode explicite
##  conformément aux capacités numériques du programme officiel de 2022
#
#######################################################################
#
#
#  Le terme source, les conditions de bords et tous les paramètres
# du problème sont traités comme des variables globales.
#  Les conditions de bord sont exprimés par diverses fonctions.
#
#  Une fonction d'animation est importée depuis un fichier séparé.
# 
#
##############################################################################


#cd /home/adomps/cours/chaleur/capacite_numerique_prog_22


import numpy as np
import matplotlib.pyplot as plt
import time
from matplotlib import animation
from math import erf
from fonction_animation import studio_poinca
v_erf = np.vectorize(erf)

#########################################################################
#### Paramètres numériques du problème 
#### traités comme des variables globales
#########################################################################


L = 1.
N = 51
Lambda = 400 # conductivité du cuivre
D = 1.2e-4  # diffusivité thermique du cuivre
x = np.linspace(0, L, N)
delta_x = x[1] - x[0]
delta_t = 0.1
rrr = 0.4 # en dimension 1, la condition de stabilité est rrr < 1/2
#rrr = D * delta_t / delta_x**2
print('rrr', rrr)
delta_t = rrr * delta_x**2 / D
q = D * delta_t / Lambda



def pdt_naif(T_anc) :
    ''' Fait évoluer l'intérieur du domaine. Les conditions de bord
    sont traitées à part. '''
    T_nouv = np.zeros(N)
    for i in range(1, N-1) :
        T_nouv[i] = T_anc[i] + rrr * (T_anc[i+1] + T_anc[i-1] - 2 * T_anc[i]) \
                    + q * p_source[i]
    return T_nouv


def pdt_vecto(T_anc) :
    T_nouv = np.zeros(N)
    T_nouv[1:-1] = T_anc[1:-1] + rrr * (T_anc[0:-2] + T_anc[2:] - 2 * T_anc[1:-1]) \
                    + q * p_source[1:-1]
    return T_nouv
                    
                    
def construit_matrice(N, rrr) :
    M = np.eye(N)  # avec des 1 dans les coins, les bords sont fixés
    for i in range(1, N-1) :
        M[i, i-1] = rrr
        M[i, i+1] = rrr
        M[i, i]  -= 2 * rrr
    return M

# Pour l'utilisation éventuelle de la méthode matricielle
Mat = construit_matrice(N, rrr)

def pdt_matrice(T_anc) :
    T_nouv = np.dot(Mat, T_anc)
    T_nouv[1:-1] += q * p_source[1:-1]
    return T_nouv


def bords_fixes(T_anc, T_nouv) :
    T_nouv[0] = T_anc[0]
    T_nouv[-1] = T_anc[-1]
    return
    

def impose_j_gauche(T_anc, T_nouv) :
    T_nouv[0] = T_nouv[1] + j0 * delta_x / Lambda    
    T_nouv[-1] = T_anc[-1]
    return
    

def evolution(T_init, N_pas, fonction_pdt, fonction_bords) :
    champ_T = np.zeros((N_pas+1, N))
    champ_T[0] = T_init
    for i in range(N_pas) :
        champ_T[i+1] = fonction_pdt(champ_T[i])
        fonction_bords(champ_T[i], champ_T[i+1] )
    return champ_T

def representation(x, champ_T, modulo, num_fig) :
    ''' Représente le champ de température en fonction de x, 
    à diverses dates multiplies de modulo. '''
    plt.figure(num_fig)
    N_pas = len(champ_T)
    plt.xlabel('$x$ (m)')
    plt.ylabel('T (°C)')    
    for j in range(0, N_pas, modulo) :
        plt.plot(x, champ_T[j])
    return


# Champ de température initial : plusieurs variantes au choix

def T_unif(x, T_interieur, T_gauche, T_droite) :
    ''' Température uniforme à l'intérieur, valeurs à choisir aux
    deux bords. '''
    N = len(x)
    T = T_interieur * np.ones(N)
    T[0] = T_gauche
    T[-1] = T_droite
    return T

def T_mode(x, t, n, A ):
    ''' Solution en mode propre. Pour l'utiliser comme condition initiale,
    prendre t=0. On peut bien sûr faire la somme de plusieurs modes.'''
    return  A * np.sin(n * np.pi * x / L) * np.exp(- D * t * n** 2 * np.pi**2 / L**2)


def source_uniforme(x, p_0) :
    ''' Définit une source uniforme de puissance volumique p_0. '''
    N = len(x)
    return p_0 * np.ones(N)

def source_portion(x, p_0, x_g, x_d):
    ''' Défini une source valant p_0 entre les abscisse x_g et x_d,
    et 0 ailleurs. '''
    N = len(x)
    p = np.zeros(N)
    masque = (x_g <= x)  * (x <= x_d) # multiplication de tableaux booléens, équivaut à AND
    p[masque] = p_0
    return p


def solution_1(x, t, T_paroi, T_avant) :
    ''' Solution exacte dans un milieu semi-infin initialement
    à T_avant dont la paroi x=0 est porté à T_paroi. '''
    return T_paroi + (T_avant - T_paroi) * v_erf(x / (2 * np.sqrt(D * t)))
    

    
    



N_pas = 100000






#p_source = source_portion(x, 100, 0.3, 0.8)






#####################################################################
###  Applications
################################################################


## Température imposée à une extrémité.
p_source = source_uniforme(x, 0)
T_init = T_unif(x, 0, 10, 0)
t1 = time.time()
champ_T = evolution(T_init, N_pas, pdt_naif, bords_fixes)
print('durée', time.time() - t1)
representation(x, champ_T, 100, 1)

studio_poinca(x, champ_T, delta_t)



# Température imposée à une extrémité et comparaison
# avec la solution en erf.
T_init = T_unif(x, 0, 10, 0)
p_source = source_uniforme(x, 0)
t1 = time.time()
champ_T = evolution(T_init, N_pas, pdt_vecto, bords_fixes)
print('durée', time.time() - t1)
representation(x, champ_T, 100, 2)
i_choisi = 200 # choisir un des indices de temps affichés
t_choisi =  i_choisi * delta_t # en secondes, prendre une date 
T_instant = solution_1(x, t_choisi, 10, 0) # solution en erf
plt.plot(x, T_instant, linestyle = 'dashed', linewidth = 2, color = 'red')
studio_poinca(x, champ_T, delta_t, intervalle = 100)



# Température identique fixée aux deux bords, et condition
# initiale en combinaison de modes propres.
t1 = time.time()
T_init = 10 + T_mode(x, 0, 1, 3) + T_mode(x, 0, 3, 3)
champ_T = evolution(T_init, N_pas, pdt_matrice, bords_fixes)
print('durée', time.time() - t1)
representation(x, champ_T, 100, 3)
i_choisi = 100
t_choisi = delta_t * i_choisi # en secondes
T_instant = 10 + T_mode(x, t_choisi, 1, 3) + T_mode(x, t_choisi, 3, 3)
plt.plot(x, T_instant, linestyle = 'dashed', linewidth = 2, color = 'red')

studio_poinca(x, champ_T, delta_t)

###################################################################
##  Exemple avec flux constant imposé à l'extrémité gauche
S = 4.6e-4 # cf TP sur la barre chauffante
Phi_0 = 20
j0 = Phi_0 / S  
T_init = T_unif(x, 0, 0, 0)
champ_T = evolution(T_init, N_pas, pdt_vecto, impose_j_gauche)
representation(x, champ_T, 100, 4)
studio_poinca(x, champ_T, delta_t, intervalle = 1)


#####################################################################
#### Exemple avec terme source
## avec source uniforme et bords fixes, on attend Delta T_max =  p0 * L**3 / (8 * Lambda)
T_init = T_unif(x, 0, 0, 0)
p0 = 1000
print('Delta T max : ', p0 * L**2 / (8 * Lambda))
p_source = source_uniforme(x, p0)
champ_T = evolution(T_init, N_pas, pdt_vecto, bords_fixes)
representation(x, champ_T, 100, 5)
studio_poinca(x, champ_T, delta_t)



##########################################################################
#### Réacteur chimique avec possible explosion, exo cf Landau
# Selon la valeur du paramètre gamma, un régime stationnaire existe ou
# non (cas d'explosion). Valeur critique gamma = 0.88
# Quand le régime stationnaire existe, on peut prédire T_max au centre.
## Se reporter à l'exo et au fichier << resolution_explosion_thermique >>.
##########################################################################
def source_reaction_chi(T, p0, alpha, T0) :
    return  p0 * np.exp(alpha * (T - T0))


def evolution_reacteur(T_init) :
    global p_source
    champ_temp = np.zeros((N_pas+1, N))
    champ_temp[0] = T_init
    for i in range(N_pas) :
        p_source = source_reaction_chi(champ_temp[i], p0, alpha, T0)
        champ_temp[i+1] = pdt_vecto(champ_temp[i])
        bords_fixes(champ_temp[i], champ_temp[i+1])
    return champ_temp



p0 = 40000
alpha = 0.01
L = 1.
ell = L/2 #  ATTENTION : demi-longueur du réacteur
Lambda = 400
T0 = 10
gamma = p0 * alpha * ell**2 / Lambda  # valeur critique gamma = 0.88
T_init = T0 + T_mode(x, 0, 1, 0.6) + T_mode(x, 0, 3, -0.2)
N_pas = 30000

champ_T = evolution_reacteur(T_init)

from resolution_explosion_thermique_Landau import resolution_theta1
print('gamma', gamma)
if gamma > 0.88 : 
    print("Il n'existe pas de régime stationnaire.")
else :
    theta1_attendu = resolution_theta1(gamma) # calculé avec le code << resolution_explosiont_thermique >>
    T_max_attendu = T0 + theta1_attendu / alpha
    print('Régime staionnaire avec T_max attendu ', T_max_attendu)


representation(x, champ_T, 100, 5)
studio_poinca(x, champ_T, delta_t)
#studio_poinca(x, champ_T[0:500], delta_t) # Si gamma > 0.88, on prend seulement le début, après on a des NaN

##########################################################################
###  Chauffage alternatif comme dans le TP sur la barre
##########################################################################
def evolution_onde_flux(T_init) :
    global j0
    periode = 300 # en secondes

    
    champ_T = np.zeros((N_pas+1, N))
    champ_T[0] = T_init
    for i in range(N_pas) :
        date = i * delta_t
        j0 = jmax * (1 + np.cos(2 * np.pi / periode * date)) / 2
        champ_T[i+1] = pdt_vecto(champ_T[i])
        impose_j_gauche(champ_T[i], champ_T[i+1])
    return champ_T  



S = 4.6e-4 # section de la barre chauffante du TP
Phi_0 = 20 # puissance en watts.
jmax = Phi_0 / S

L = .3
N = 51
Lambda = 400 # conductivité du cuivre
D = 1.2e-4  # diffusivité thermique du cuivre
x = np.linspace(0, L, N)
delta_x = x[1] - x[0]
delta_t = 0.1
rrr = 0.4 # en dimension 1, la condition de stabilité est rrr < 1/2
#rrr = D * delta_t / delta_x**2
delta_t = rrr * delta_x**2 / D
q = D * delta_t / Lambda
print('rrr', rrr)
N_pas = 20000


T_init = T_unif(x, 20, 20, 20)
champ_T_onde = evolution_onde_flux(T_init)
representation(x, champ_T_onde, 123, 6)
t = [i * delta_t for i in range(N_pas + 1)]

studio_poinca(x, champ_T_onde, delta_t, intervalle = 1)

plt.figure(7)
plt.clf()
plt.plot(t, champ_T_onde[:, 1], color = 'black')
plt.plot(t, champ_T_onde[:, 5], color = 'red')
plt.plot(t, champ_T_onde[:, 9], color = 'blue')
plt.xlim(0, 1500)
plt.ylabel('T (°C)')
plt.xlabel('t (s)')