-------------------------------------------------------------
COURS M2 Bijan Mohammadi Universite de Montpellier
-------------------------------------------------------------
https://umontpellier-fr.zoom.us/my/bijan.mohammadi.zoom
--------------------------------------------------------------

contenus sur https://imag.umontpellier.fr/~mohamadi/ADAPT/

-----------------------------
FAIRE VOTRE PAGE WEB (sur GIT)
Creer un compte github

Creer votre page web + d'autre repo avec vos projets passes

Y DEPOSER:
-VOTRE CV
-Vos realisations dans ce cours
-TOUS LES AUTRES PROJETS
m'envoyer le lien par email: la qualite de votre page sera prise en compte dans la note à l'examen.

-----------------------------

SEANCE 1 
-Prise en main de spyder ou jupyter notebook

Faire un code Python pour resoudre l'equation differentielle suivante:

u,t = - lambda u, u(0)=1,  lambda = 1 

en utilisant un schema d'Euler explicit et comparer avec la solution exacte: uex(t)= u0 exp(-lambda t)

Pour Dt = 1s, tracer les solution exacte et numerique et l'erreur en temps sur 1 minute.

Tracer la courbe de l'erreur L2 de la fonction et de sa derivee "en fonction du pas de temps".
Mettre figures cote a cote. 

------------------------------------------------------------------------

LIRE JUSQU'A 4/ poly_adapopt_bm.pdf

----------------------------------------------------------------------

Voir Annexe B du POLY.pdf pour rappel des definitions

V=(v1, v2)
u,t + <V , grad u > -  nu div( grad u) = -lambda u + f
u,t + v1 u,x + v2 u,y - nu ( u,xx + u,yy ) = -lambda u + f

Conditions aux limites de Dirichlet uniquement lorsque V.n(s)<0 (bords entrants)

Ouvrir un compte sur une IA generative.

Utiliser une IA generative pour generer un code python pour resoudre cette equation 
dans un domaine rectangulaire avec 
f(t,s) = Tc exp(-k d(s,sc)**2)   avec d(s,sc)**2 = (s1 - sc1)**2 + (s2 - sc2)**2 

----------------------------------------------------------------------
passage 1D:
u,t + v1 u,x - nu  u,xx = -lambda u + f
+ condition initiale u(t=0,x)=u0(x)
+ conditions aux limites:   u(t,x=0)=ul,  u,x(t,x=L)=g
Dirichlet non homogene a gauche et Neuman non homogene a droite.

Demander a l'IA generative u0 compatible avec les CL a gauche et a droite.

------------------------------------------------------------------------

SEANCE 2

Soit:
-comprendre le code Python adrs.py pour resoudre:

-pour u(s)=uex(s) = exp(-10*(s-L/2)**2), trouver f, la porter dans le code
u,t + v u,s - nu  u,ss  + lambda u = f(s)
v=1, nu=0.01, lambda=1   L=1

-Implementer f(s) dans le code et verifier si on retrouve uex(s) ?

maillage uniforme en N points = discretisation de [0,L] en (N-1) intervalles s_i = (i-1)*h avec h=L/(N-1) 

-s'assurer de la convergence vers la solution stationnaire pour un maillage de 100 points
-tracer la convergence ||un+1 - un||_L2  normalisee : on a la solution stationnaire
-calculer apres convergence les normes L2 et H1 de (u-uexact) pour differents maillages, partant de 3 points.
-tracer l'erreur L2 et H1 en fonction de h=dx pour ces cinq maillages
-Placer les images cote a cote.

Pour le schema numerique:
Identifier ou sont les elements suivants dans le code vu en seance
-decentrage = centre + viscosite numerique  (chap 12)
-condition de stabilite CFL (chap 10)
-quelle condition en sortie est implementee actuellement ? 
u1(n)=u1(n-1)-ux*(xmax-x(n-1))+uxx(n-1)*(xmax-x(n-1))**2/2 
u1(n)=u1(n-1)+ux*h+uxx(n-1)*h**2/2  = u1(n)=u1(n-1)+uxx(n-1)*h**2/2 
-implementer u,s(s=L)=0

-------------------------------------------------------------

-4  -> interpolation Ph (lineaire)
-regarder thm 18.1.3 (relation entre l'erreur d'interpolation et la derivee seconde)

-Evaluer en discret les normes des differences du poly sur les differents maillages 
(L2 difference solutions et derivees approchee et exacte, 
L2 difference solutions exacte et son interpole P1, derivee seconde).

Exporter dans un fichier la norme L2, la norme H1, la semi-norme H2 pour chaque dx
Utiliser un algorithme d'optimisation pour identifier (C,k) et trouver M.  
Poly 2.3.1  (c'est un probleme de regression)
quel est l'ordre en espace de ce code ?

Apres avoir trouve x_opt=(C,k)  pour chaque probleme de moindres carres:
-Visualiser les courbes  : Q1 vs. C h^(k+1)  et Q2 vs. C h^k

|u − uh|_{0,2} / |u|_{2,2} < ? h**? 

-------------------------------------------------------------

SEANCES 3-4

Comprendre l'adaptation a posteriori en comparant l'integration de Riemann et Lebesgues
pour calculer l'integrale de :
f(x)=exp(-a*(x-L/3)**2)-2*exp(-a*(x-2*L/3)**2), a=100, L=1
1. integrer en Riemann cette fonction (pas uniforme en x)
2. adapter le pas d'integration (Lebesgues) pour avoir un pas uniforme en y=f(x)

--------------------------------------------------------------

code adrs_multiple_mesh_adap.py

-lire 18.4.2
-Faire les etapes de 5/
-identifier la definition de la metrique dans le code (quelle loi est en place ?)
-implementer la loi 3
-tracer NX(err) pour err=0.04, 0.02, 0.01, 0.005, 0.0025
-comment evolue NX en fonction de err ?
-Quels sont les criteres d'arret pour l'iteration d'adaptation. 
Introduire un critere mixte portant sur le nombre de points de maillage et l'erreur L2: ne pas s'arreter 
tant que les 2 ne sont pas realises.
-Identifier la contraction dans ce code

-Restructurer le code en introduisant 3 fonctions (adrs, metric, mesh)

-------------------------------------------------------------

SEANCES 4-5 

Code adrs_insta.py : solution instationnaire avec solution exacte cible instationnaire a creer.
-Visualiser l'erreur a T/2 and T(fin de calcul) pour differents maillages uniformes.
-Visualiser l'evolution de l'erreur au point milieu du domaine pour differents Runge-Kutta (1-4).

----

Code adrs_insta_multiple_mesh_adap.py a partir de adrs_multiple_mesh_adap.py

AVEC UN TERME SOURCE DEPENDANT DU TEMPS calcule pour uex(t,s)=u(t) v(s)
--> f(s,t)

-modifier le code avec ce terme source (ou bien verifier ce qui est deja implemente).
-le residu ne converge pas vers 0 car instationnaire:
u(t)v(s)=sin(4*pi*t) v(s) => u'*v=4*pi*cos(4*pi*t)*v(s)

-visualiser la solution a differents instants pour Time=1s (deja en place ?)

-Quels sont les criteres d'arret pour l'iteration d'adaptation. 
Introduire un critere mixte portant sur le nombre de points de maillage et l'erreur L2: ne pas s'arreter 
tant que les 2 ne sont pas realises (fait plus haut).

-que donne le maillage adaptatif avec controle de metrique stationnaire (base sur 
la solution finale en t=Time).

- introduire les etapes de l'adaptation en instationnaire (intersection
des metriques en temps) dans le code (moyenne en temps des metriques).
Lire : 18.7.1 (instationnaire) --> 427

-------------------------------------------------------------
SEANCE 6 :  Partie 9

code optim_adrs.py 

--> Optimisation en utilisant la linearite de l'equation d'etat par rapport aux controles.

ADRS avec N controles sources:

Tracer la surface J(x1,x2,....) en echantillonnant les deux premiers controles (les autres fixes).

Si l'adaptation de maillage est utilisee lors des calcul des solutions et seconds membres elementaires, 
ces quantites ne seront pas disponibles  sur le meme maillage.
Comment calculer alors numeriquement les expressions Aij et Bi ?

Quelle doit etre la precision de l'integration numerique approchee (Aij, Bi)
pour ne pas introduire d'erreur supplementaire a ce niveau ? voir poly_cours

Modifier le code fourni pour repondre aux questions precedentes (en introduisant l'adaptation 
et l'interpolation entre maillages).

comparer le controle optimal obtenu avec la solution sur maillage fixe (assez fin pour servir de
reference).

///////////////////////////////////////////////////
0. ADRS (xcible = (1,2,3,4) --> udes=u(xcible) --> J(xcible)= 1/2 ||u(xcible)-udes||^2 
1. ADRS (x) --> u(x) --> J(x)= 1/2 ||u(x)-udes||^2
2. Utiliser un minimiseur python pour trouver xopt 
3. Utiliser la linéarité de l'équation ADRS pour simplifier les calculs (voir plus bas)

////////////////////////////////////////////////////

Algorithm

Probleme inverse
define Xopt=(1,2,3,4)  in R^4
define udes=u(Xopt) par ex

J(x) = min_x 1/2 ||u(x)-udes||^2  -> X* = Xopt ?

X* such that Grad_x J(X*) = 0

---
Refinement loop Th2 

find U0
find U1, U2, U3, U4 <- U(ei)

Calculate :  Grad_x J(u(x)) = Ax-b = 0  (to show)
Aij= <ui,uj>L2
bi=<ui,(u(x)-udes)>L2

Invert A X* = b with x=(alpha_1,...,alpha_4)
 
check if u(X*)=udes? calculating J(Xopt) = ||u(X*)-udes||^2 and if Xopt = X*
//////////////////////////////////////////////////////

Introduce a refinement loop and see the impact on the convergence to Xopt (plot each component of X*(h)). 
Plot J(X*(h)) and error=||X*(h) - Xopt||

Now look for udes=1 (we do not know Xopt) and see how close you can get with this control ?
Plot J(X*(h)) and the evolution of each component of X*(h)

Gradient calculation
Grad_x J(u(x)) = Ax-b = 0 
J=1/2 h sum_i=1^NX (u0 + (sum_j=1^nbc x_j u_j) - udes)**2   
Grad_xj J(u(x)) = h sum_i=1^NX (u0(i) + (sum_j=1^nbc x_j u_j(i)) - udes) u_k(i),  for j,k=1,...,nbc

h sum_i=1^NX ((sum_j=1^nbc x_j u_j) u_k(i)) = - h sum_i=1^NX (u0(i) - udes(i)) u_k(i)  for k=1,...,nbc
sum_j=1^nbc x_j (h sum_i=1^NX ( u_j(i) u_k(i)) = h sum_i=1^NX (udes(i) - u0(i) ) u_k(i)

A x = b

Aij= <ui,uj>L2
bi=<ui,(udes - u0(x))>L2

-------------------------------------------------------------
-------------------------------------------------------------
SEANCE 7 : EXAMEN

Examen ecrit + releve de compteur sur les sites web

-------------------------------------------------------------

OPTIONNEL:
-------------------------------------------------------------
-------------------------------------------------------------
-------------------------------------------------------------

On considere u,t=sin(2 pi t), u(0)=a et J(a)=(1/2) u^2(t=1,a)

-Calculer dJ/da = u u,a
-Ecrire un petit programme pour resoudre l'equation par Euler explicit
-Calculer dJ/da par differences finies et variables complexes
-Comparer au gradient exact

u(t)=-1/(2 pi) (cos(2 pi t)-1)+a 
u,a=1
dJ/da = u u,a

-------------------------------------------------------------

Ex 10.2

-Quelle est la complexite en terme de stockage de l'information pour
pouvoir calculer l'adjoint par integration retrograde ?

-Refaire les calculs avec l' equation nonlineaire suivante:
ut + uu,s - nu u,ss = lambda u + f(x):
Quelles nouvelles consequences sur le stockage (besoin d'utilisation des etats intermediaires) ?

Le programme advdiffadapadj propose une implementation de l'ensemble de cette discussion.

Quelle est la fonctionnelle consideree ?

Quelle est l'approximation aux differences finies utilisee ?

comparer le gradient par differences finies et l'adjoint: illustrer les deux gradients sur une meme figure.

Que se passe-t-il au niveau du maillage quand l'adjoint est pour une fonctionnelle donnee par la norme H^1.

------------

Utiliser ce script pour generer une BD d'apprentissage contenant pour chaque scenario:
(largeur entree, vitesse entree, viscosite, Temperature entree, coupe a 1/3 domaine, coupe a 2/3 domaine)
les coupes seront en 10 points et releveront les valeurs des champs de vitesse et temperature : 
30 valeurs (vecteur vitesse + champ temperature)x10 points

Pour plus de facilite, commencer avec une seule entree.

Apprendre cette BD, et donner l'erreur en crossfold avec 
un modele lineaire donne dans le scipt python skl_... pour deux problemes:

R4 -- > R30 : connaissant les parametres de fonctionnement (largeur entree, vittesse entree, viscosite, Temperature entree)
trouver les champs sur la coupe en 2/3
 
R30 -- > R30: connaissant le champs a 1/3 trouver le champ a 2/3 du domaine.

Le premier probleme permet d'utiliser la ML pour trouver le champ pour de nouveaux points de fonctionnement.

Le second permet de trouver apres une mesure non intrusive, le champs dans une partie non visible du domaine.


Le script skl... permet de lire des fichiers csv, de visualiser les donnees, et
d'apprendre la base de donnees pour identifier le modele lineaire adapte, 
en utilisant differentes regularisations (voir exo 39 poly cours_optim_bm.pdf)

-------------------------------------------------------------

C++  

1/advect.cpp est un petit code simple en C++ pour la solution de l'equation de transport. 
Comprendre le code. Compiler (g++ -O advect.cpp -o advect.exe) / executer et visualiser le resultat avec 
gnuplot. 

2/adap.cpp est une version C++ de notre travail pour la solution de l'equation ADRS et l'adaptation de maillage par controle de metrique.
Compilation: g++ -O adap.cpp -o adap.exe
La solution est dans results.dat et peut etre visualise avec gnuplot.

Ce sont des codes statiques et auto-suffisants. Il y a un exemple d'insertion de librairie externe en statique
dans advect.cpp

- decouper adap.cpp programme en routines/librairies avec un programme principal.

- creer un makefile 

- completer cet ensemble avec la partie optimisation qu'on a vu en Fortran et le calcul d'adjoint.

Ce travail vous preparera au cours C++ a suivre.

-------------------------------------------------------------

--------------------------
exemple sortie gnuplot
set title 'TITRE'
set xlabel 'x(M)'
set ylabel 'u(kg/m3)'
plot [0:1][] 1/x w l lw 4 title 'Model',1/x**2 lw 3 title 'Model 2'
--------------------------
plot 1:(log10($2))  w l
plot [0:1][] 1/x**2 w linespoints lw 3