# Exercices école d'été CRM
-----------------------------------------------------------

Nous présenterons le problème physique que nous allons utiliser pour expérimenter les différents préconditionneurs!

---------------------
## __Le problème physique__
--------------------

Soit le problème de Boussinesq est constitué d'un système couplant les équations de Navier-Stokes et de la diffusion.

$-\nabla \cdot (\boldsymbol{K}(T) \nabla T) = 0$

$-\nabla \cdot (2\eta\dot{\boldsymbol{\gamma}}(\boldsymbol{u})) + \nabla p = \rho {\bf r}- \rho \left(\frac{\partial {\bf u}}{\partial t} + ({\bf u}\cdot \nabla){\bf u}\right)$

où $\boldsymbol{K}(\boldsymbol{x})$ est le tenseur (matrice) de
conductivité thermique du matériau (voir Fortin et Garon, chapitre 6,
bien que les conditions limites diffèrent).

Nous allons nous intéresser à un problème couplé d'un système de refroidissement d'un transformateur d'électricité. Le transformateur simplifié est représenté par un noyau de cuivre qui génère de la chaleur par induction. Ce noyau baigne dans une huile qui circule à basse vitesse pour le refroidir. Voici un schéma des principales composantes:

<center>
</center>

Pour la partie fluide, on impose un profile de vitesse au fluide à l'entrée du domaine (représentée par
l’entité “inlet”) et une vitesse nulle sur les parois, sauf la sortie ()"outlet") pour laquelle la vitesse est laisée libre. Pour la partie thermique, on impose une température à l'entrée de 298.15 degrés Kelvin et une source thermique de 500W/m^2.

<center>
$\displaystyle \int_\Omega (\boldsymbol{K}(\boldsymbol{x})
\nabla T) \cdot \nabla w\ T dA = \int_{\Gamma_{sortie}}-100\,
w\ T ds$
</center>

Notons qu’il faudrait relever les conditions de Dirichlet non nulles
pour être plus rigoureux. Ce relèvement est toutefois fait
automatiquement dans MEF.

Voici les fichiers d'entrées associés au problème sont dans les répertoires d'exercice:
> [feracheval.mail](./feracheval.mail)
> [feracheval.enti](./feracheval.enti)
> [feracheval_stationnaire.champs](./feracheval.champs)

Nous allons maintenant lancer plusieurs fois le calcul en changeant les options pour le préconditionnement

-------------
## __Étapes à suivre__
-------------

Il faut éditer le fichier .champs :

#### Étape 1 : Lancer le calcul

Dans la VM... lancer...

vous aurez:

coucou

#### Étape : Visualiser les résultats

Paraview

#### Étape : Afficher la matrice

Éditer le fichier .champs

Éditer le fichier .petscrc:

#Pour faire une pause à l'affichage de la matrice:
#-mat_draw_pause -1

Relancer le calcul

On voit la matrice avant assemblage, à l'assemblage avec 3 couleurs: significations
Important: clique droit pour continuer

(Mettre une image de ce qu'ils voient)

#### Étape : Modifier l'ordre des DDLs

Déclarer le split des variables via le problemeef:

problemeef ProbGlobal [defaut, defaut, [(u*,p*),T]]

Relancer et observer la matrice:
(mettre une image ...)

#### Étape : PCFieldSplit

Block Jacobi => additive
Gauss-Seidel => multiplicative
Gauss-Seidel symmétrique => symmetric_multiplicative

Pour en savoir plus, voir la doc PETSc: lien..

{{ecritEtape(2,"Déclarer le(s) problème(s) éléments finis à résoudre")}}

problemeef probDida

Le mot clé « problemeef » déclare un problème éléments finis nommé « probDida ». Si on souhaite résoudre plusieurs problèmes différents, il faudra tous les déclarer de la même façon avec chacun leur propre nom. Une bonne pratique est de tous les déclarer au début avec des noms clairs décrivant bien le problème à résoudre (probNavierStokes ou problemeEFDiffusion par exemple). Nous verrons dans des exemples futurs ce genre de cas. Cela permet au lecteur de voir rapidement quels types de problèmes sont résolus par le .champ (thermique, fluide, contact, optimisation de forme, etc...).

### Les variables du problème physique.

{{ecritEtape(3, "Déclarer la variable T de notre problème")}}

scallin T 0

Le mot clé « scallin » déclare un champ éléments finis scalaire P1 nommé « T » et initialisé à 0. Il existe beaucoup d’autre types de champs éléments finis. Vous pouvez trouvez une liste non exhaustive ici : Fichier_Champs_CL. Ou bien utiliser l’éditeur de fichier champs [eMEF++](https://giref.ulaval.ca/eMEF++/).

{{ecritEtape(4, "Déclarer des blocs matériaux")}}

{ Acier
to2sym K [80.79, 80.79, 80.79, 0, 0, 0]
}
{ Cuivre
scalaire K1 413.66
empilementto2 K [K1,K1,K1,Zero,Zero,Zero]
}

{{ecritEtape(5, "Associer les matériaux aux parties du domaines")}}

champssurgeom mat1=Acier mat2=Cuivre

Dans le `.enti`, on a déclaré des entités mat1 et mat2 qui correspondent à la répartition sur le domaine de nos deux matériaux. On associe ici les valeurs des champs définis dans les blocs Acier et Cuivre à leurs entités.
Si votre domaine est constitué d’un seul matériau il n’est pas nécessaire de déclarer les paramètres dans un namespace ni de les associer à une entité, il suffit de les déclarer comme pour « T » . De plus ce n’est pas la seule méthode possible pour séparer un domaine en différents matériaux mais c’est la manière la plus simple et la plus élégante quand la répartition des matériaux ne change pas dans le domaine pendant la résolution.

{{ecritEtape(6, "Déclarer un champ pour le terme source")}}

scalaire F 0

Ici cette étape est facultative car le terme source est nul dans notre problème, on pourrait se contenter de ne pas l’ajouter au problème éléments finis, mais pour des raisons pédagogiques nous allons quand même le faire.
Le mot clé « scalaire » déclare un champs qui n’est pas éléments finis. Dans notre cas à nous il contient une expression algébrique qui sera évalué aux coordonnées voulues (ddls, point d’intégration etc..).
Ce genre de champs existe aussi pour des vecteurs ou des tenseurs (voir : Fichier_Champs_CL, eMEF++).

Nous avons finis de déclarer toutes les variables physiques de notre problème. Une bonne pratique de programmation est de tout déclarer au même endroit afin de ne pas avoir à chercher à différents endroits si vous souhaitez modifier votre mise en données.

Passons maintenant aux schémas d’intégration :

### Le(s) schéma(s) d’intégrations:

{{ecritEtape(7, "Définir un namespace et construire un schéma d’intégration")}}

{ ParamSchema
scalaire ptintgarete 2
scalaire ptintgface3 3
scalaire ptintgface4 2
scalaire ptintgelement1d 2
scalaire ptintgtriangle 12
scalaire ptintgquadrangle 2
scalaire ptintgtetra 5
scalaire ptintgprismetri_legendre 2
scalaire ptintgprismetri_triangle 3
scalaire ptintghexa 2
}

schemaintg schema_intg ParamSchema

Si vous n’êtes pas familier avec et que vous ne savez pas combien de points est nécessaire pour intégrer l’ordre que vous désirez, {{reqLien("eMEF++")}} est là pour vous.

Le mot clé « schemaintg » déclare un schéma d’intégration portant le nom « schema_intg » et dépendant des options « ParamSchema ». Il est tous à fait possible de construite plusieurs schémas d’intégration, chacun à partir d’un namespace différent afin d’avoir un schéma différent pour chaque problème si nécessaire. Nous verrons plus tard que chaque terme de formulation se voit attribué un schéma d’intégration.

Passons maintenant au solveur.

### Le(s) solveur(s)

Chaque problème élément fini a besoin d’un solveur. Comme nous avons un problème linéaire stationnaire, nous allons donc construire un solveur linéaire pour la matrice.

{{ecritEtape(8, "Définir un namespace pour les options de solveur et construire le solveur linéaire")}}

{ options_prefixees_petsc options_slin
ksp_type cg
ksp_atol 1e-20
ksp_rtol 1e-12
ksp_divtol 1e30
ksp_max_it 4000
ksp_monitor
pc_type gamg
}

{ solveur_lin solveur(probDida)
prefixe_options options_slin
}

Le mot clé « options_prefixees_petsc » déclare un namespace pour des options de solveurs associé au nom « options_slin ».
Ici nous avons des options pour un solveur de type gradient conjugué « cg ». Les différentes options de solveurs seront détaillées dans une leçon dédiée.

Le mot clé « solveur_lin » déclare un namespace pour le solveur « solveur(probDida) » associé au problème EF « probDida ».
Le mot clé « prefixe_options » indique les options de solveur à associer, ici « options_slin ».

### Le problème élément fini

Passons enfin au vif du sujet. La grande différence entre MEF++ et les autres mains comme probDiffusion est que le problème physique à résoudre n’est pas défini en dur dans le code. Il faut le déclarer dans le .champs. Cela offre une grande souplesse pour tester de nouvelles formulations mais demande plus de travail dans le .champs.

{{ecritEtape(9, "Définir un namespace pour le problème EF")}}

On commence par déclarer un namespace portant le nom du problème EF que l’on veut renseigner

{ probDida

On déclare le terme de formulation associé à l'intégrale

<center>
$\displaystyle \int\ \Omega (\boldsymbol{K}(\boldsymbol{x})
\nabla T) \cdot \nabla w\ T dA$
</center>

TFDiffusion TFDiffT [champgeo, T, T, K, schema_intg]

On commence par le mot clé indiquant le type de TF à ajouter au problème, ici « TFDiffusion », on lui donne ensuite un nom « TFDiffT » puis la liste des paramètres dont le TF a besoin « `[champgeo, T, T, K, schema_intg]` »
Les paramètres varient d’un TF à l’autre. Pour avoir la liste des TFs disponibles, leur mot clé associé ainsi que la liste des paramètres il est fortement recommendé d’utiliser {{reqLien("eMEF++")}}.

Nous allons détaillé les paramètres pour notre TFDiffusion. En premier on renseigne le champ géométrique (`« champgeo »`), ensuite on renseigne l’inconnu variationnelle du TF (`« T »`), puis un champ dont on veut les fonctions de base (encore `« T »` car les fonctions de base sont du même type que notre inconnu mais ce n’est pas toujours le cas), on renseigne ensuite le tenseur d’ordre 2 (`« K »`). Enfin on ajoute le schéma d’intégration que l’on désire pour le TF (`« schema_intg »`).

On fait la même chose pour le terme source

<center>
$\displaystyle \int_\Omega \boldsymbol{F} \cdot w_T dA$
</center>

TFDiffusion TFDiffT [champgeo, T, T, K, schema_intg]

Pour les paramètres, ici les paramètres sont : le champ géométrique `« champgeo »`, le champ source « F », le champ pour les fonctions de base `« T »` et le schéma d’intégration `« schema_intg »`.

Passons maintenant aux conditions aux limites : contrairement aux autres mains, dans MEF++ les conditions aux limites sont dans le .champs et non dans un fichier séparé. On ajoute les CLs dans les namespaces des problèmes finis, comme pour les termes de formulations.

On commence par une condition de dirichlet :

scalaire cond_entree 100
dirichlet scalaire "entree" "T" "cond_entree"

Les deux premiers mots clés d’une CL indiquent toujours respectivement le type de condition que l’on souhaite imposé (`« dirichlet »`) et le type de champ sur lequel on souhaite l’imposer (`« scalaire »`). Ensuite il faut indiquer l’entité géométrique où imposer la condition (`« "entree" »`) puis l’inconnu du problème (`« "T" »`) et enfin la valeur à imposer (`« "cent" »`).

Puis une neumann :

scalaire cond_sortie -100
neumann scalaire "sortie" "T" "cond_sortie" "schema_intg

Pour une neumann, il faut un paramètre de plus, en effet contrairement à une condition de type dirichlet la neumann implique une intégrale sur le bord et donc nécessite un schéma d’intégration (`« "schema_intg" »`).

Une fois fini, on ferme le namespace :

}

### Les exportations

{{ecritEtape(10, "Définir un namespace d'exportation")}}

Il faut maintenant définir un bloc exportation. On ouvre un namespace d’exportation :

{ exportation probDida

En lui donnant le nom du problème EF « "probDida" », on attache l’exporation au problème EF en question. L’exportation sera alors automatiquement exécuté à la fin de la résolution du problème. Il est tout à fait possible de créer des blocs d’exportations sans les associés à un problème EF (en leur donnant un nom différent) et de contrôler le moment de l’exportation (voir [PrePostTraitements](./../Construire fichier champs PrePostTraitements)).

On indique le type d’interpolation désirée pour l’exportation:

asgnDegreInterpolation Exportation::eMAILLAGE_LINEAIRE

Un booléen optionnel qui permet de désactiver l’exportation du maillage:

asgnExporteToujoursMaillage faux

Par défaut, il vaut mieux le garder à faux. Il est utile de le passer à vrai que si votre maillage change pendant le calcul.
On ajoute le champ K à l’exportation :

ajouteChamp K

Par défaut un problème EF exporte toujours ses inconnues variationnelles, il est alors nécessaire d’ajouter uniquement les champs qui nous intéresse et qui ne sont pas des inconnus du problèmes. Si votre bloc d’exportation n’est associé à un problème, il vous faut ajouter tous les champs qui vous intéresse.

Enfin, on donne le chemin d’exportation (par défaut « ./ ») et le préfixe pour le nom des fichiers exportés:

prefixe_sorties $CheminAccesExport$/resultats_FerachDiff

Puis on ferme le namespace:

}

Il ne reste plus qu’à lancer la résolution et visualiser les résultats dans paraview comme vu précédement [Utiliser_MEF++.](https://redmine.giref.ulaval.ca/projects/mefpp/wiki/Utiliser_MEF++)

-------------------
## __Résultats numériques__
--------------------

<center>


</center>