# Réduction de modèle a priori par séparation de variables espace-temps pour le problème d'élastodynamique transitoire

Lucas Boucinha, Anthony Gravouil, Amine Ammar

► To cite this version:

Lucas Boucinha, Anthony Gravouil, Amine Ammar. Réduction de modèle a priori par séparation de variables espace-temps pour le problème d'élastodynamique transitoire. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722078

HAL Id: hal-01722078

https://hal.science/hal-01722078v1

Submitted on 2 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Réduction de modèle *a priori* par séparation de variables espace-temps pour le problème d'élastodynamique transitoire

Lucas BOUCINHA<fnref n="1" />, Anthony GRAVOUIL<fnref n="1" />,<fnref n="3" />, Amine AMMAR<fnref n="3" />

**Résumé** — Les méthodes de représentation à variables séparées permettent de réduire de façon drastique les coûts numériques associés à la résolution d'une grande variété de problèmes. Cependant, l'efficacité des décompositions espace-temps pour les problèmes d'évolution de type hyperbolique reste à démontrer. Dans cette exposé, nous évaluons la capacité de la décomposition généralisée propre (PGD) à construire des représentations par variables séparées espace-temps du problème d'élastodynamique transitoire à deux dimensions d'espace. Ce problème étant non-symmétrique, la décomposition est classiquement construite par minimisation du résidu espace-temps dans la métrique euclidienne. Cependant, lorsque le problème est mal conditionné, cette décomposition n'est pas optimale. Aussi, une nouvelle définition de la PGD par minimisation du résidu dans une métrique "idéale" a été récemment introduite et permet de calculer la décomposition optimale *a-priori*. Nous appliquons cette nouvelle définition au problème d'élastodynamique formulé sur le domaine espace-temps en déplacement-vitesse. Nous comparons alors deux types de décomposition, à savoir espace-2D × temps et espace-1D × espace-1D × temps.

**Mots clés** — PGD, élastodynamique, séparation espace-temps, formulation idéale en minimum de résidu

# 1 Introduction

Les techniques classiques de discrétisation sont désormais suffisamment robustes pour être utilisées dans l'industrie. Cependant, leur application à des problématiques industrielles nécessite des ressources informatiques toujours plus importantes. Bien que les moyens de calcul actuels permettent de résoudre des modèles à plusieurs centaines de millions d'inconnues, il s'avère aujourd'hui nécessaire de proposer de nouvelles méthodes afin de mieux exploiter les ressources informatiques disponibles. Les méthodes de représentation à variables séparées sont un bon moyen pour atteindre cet objectif. Introduites dans les années 1980 pour la résolution des problèmes non-linéaires de mécanique [16], elles ont été récemment généralisées pour les problèmes multi-dimensionnels [5, 1, 10, 17]. Depuis une dizaine d'années, de nombreuses contributions ont été publiées [11], démontrant que les coûts associés à la résolution d'une grande variété de problèmes peuvent être réduits de plusieurs ordres de grandeur avec ces méthodes.

Pour les problèmes d'évolution [2, 18], la séparation s'opère sur les variables de l'espace physique et du temps. Le principe est alors d'approximer la solution espace-temps par une somme de produits de fonctions de l'espace et du temps, chaque produit pouvant être interprété comme un mode espace-temps. Si le nombre de mode espace-temps permettant de représenter la solution complète avec une bonne précision est négligeable devant la dimension des espaces de discrétisation, la mémoire nécessaire au stockage de la solution complète peut être réduite de plusieurs ordres de grandeur. Il s'agit alors de savoir combien de modes espace-temps sont nécessaires pour représenter de façon précise la solution d'un problème donné ? Pour la dynamique transitoire, cette question reste ouverte. Une première réponse a été donnée dans [7] montrant qu'un nombre raisonnable de modes espace-temps est suffisant pour représenter précisément la propagation d'ondes élastiques dans une structure mono-dimensionnelle. Nous considérons ici le problème à deux dimensions d'espace et des chargements de type impact avec des fréquences de sollicitations variant sur une large gamme de fréquences (de 1 kHz à 100 kHz).

Afin d'évaluer l'efficacité des méthodes de séparation de variables, il est tout d'abord nécessaire de donner une définition de la décomposition optimale. Celle-ci est définie *a-posteriori* (c'est-à-dire après

<footnote n="1">¹ Université de Lyon, INSA-Lyon, LaMCoS, CNRS UMR5259, Villeurbanne, France, {lucas.boucinha,anthony.gravouil}@insa-lyon.fr</footnote>

<footnote n="2">² Arts et Métiers ParisTech, ENSAM Angers, amine.ammar@ensam.eu</footnote>

<footnote n="3">³ Institut Universitaire de France</footnote>
avoir calculé la solution complète sur le domaine espace-temps) par minimisation de l'erreur entre la
solution complète et l'ensemble des décompositions possibles. Le réel enjeu est alors de construire cette
décomposition *a-priori*, c'est-à-dire sans calculer (et donc stocker) la solution complète avant de calculer
sa représentation séparée. Ceci est l'objet des méthodes de décomposition généralisée propre (PGD). Le
problème d'élastodynamique étant non-symmétrique, la décomposition *a priori* est classiquement cal-
culée par minimisation du résidu espace-temps dans la métrique euclidienne [18]. Cependant, lorsque le
problème est mal conditionné, cette décomposition ne converge pas vers la décomposition optimale. Une
solution, introduite dans [6] et également évoquée dans [9], consiste à trouver une norme "idéale" dans
laquelle minimiser le résidu de façon à ce que la décomposition obtenue minimise également l'erreur.

Bien que cette nouvelle définition de la PGD permette théoriquement d'obtenir la décomposition
optimale *a-priori*, les algorithmes utilisées pour sa construction ne convergent pas nécessairement vers
la solution optimale. Notamment pour les décompositions d'ordre élevé (typiquement supérieur ou égale
à trois), une construction progressive de la décomposition ne permet pas d'obtenir la décomposition
optimale. Aussi, un algorithme optimal a été introduit dans [8] et permet de construire la décomposition
optimale également pour les décompositions d'ordre élevé. Nous illustrons cet aspect en comparant deux
types de décomposition, à savoir espace-2D × temps et espace-1D × espace-1D × temps.

Cet exposé est construit comme suit : après avoir décrit le problème d'élastodynamique modèle dans la section 2, nous introduisons la décomposition optimale *a-posteriori* dans la section 3. Nous montrons ensuite dans la section 4 comment construire la décomposition optimale *a-priori*, à l'aide de la PGD définie par minimisation du résidu dans une métrique idéale.

## 2 Description du problème

Nous considérons le problème de dynamique transitoire décrit sur la Figure 1a). Ce problème mod-
elise la propagation des ondes élastiques planes dans une structure bi-dimensionnelle Ω au cours de
l'intervalle de temps I = [0, T]. Le domaine spatial Ω = Ω₁ × Ω₂ est choisi de forme rectangulaire (de
dimensions 1 m × 0.5 m). On considère un comportement linéaire élastique et isotrope du matériau (mod-
ule d'Young E = 2.10¹<fnref n="1" /> N/m², coefficient de poisson ν = 0.3, masse volumique ρ = 8000 kg/m³). La
structure est encastrée à l'une de ses extrémités et un chargement sinusoïdal de type impact σ̄ est appliqué
au milieu de l'extrémité opposée durant une période ΔTload (voir Figure 1a)). Les conditions initiales en
déplacement et vitesse sont imposées nulles. La simulation est effectuée jusqu'à T = 10⁻³ s.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – a) Géométrie du cas test et conditions aux limites dans le domaine espace-temps $\Omega \times I$, et b) Décompositions espace-temps considérées</figcaption></figure>

Le problème d'élastodynamique est formulé de façon faible sur le domaine espace-temps avec la
méthode de Galerkin discontinue en temps en déplacement-vitesse [15]. Cette formulation à deux champs
permet de contrôler les hautes fréquences parasites introduites par les schémas classiques d'intégra-
tion en temps<fnref n="1" />. Les inconnues du problème sont les champs de vecteur déplacement $\underline{u}(x_1, x_2, t) = \sum_{i=1}^{2} u_i(x_1, x_2, t) e_i$ et de vitesse $\underline{v}(x_1, x_2, t) = \sum_{i=1}^{2} v_i(x_1, x_2, t) e_i$, définis sur le domaine espace-temps $\Omega_1 \times$

<footnote n="1">1. Une comparaison numérique de l'influence des perturbations hautes fréquences sur l'efficacité des décompositions espace-temps est effectuée dans [7] et montre notamment que l'utilisation de schémas contrôlant les oscillations parasites permet de réduire le nombre de modes espace-temps nécessaires pour représenter la solution complète à précision fixée.</footnote>
$\Omega_2 \times I$. La discrétisation du problème est effectuée avec les mêmes espaces d'approximation pour les champs de déplacement et de vitesse, à savoir des éléments finis linéaires (continus en espace et discontinus en temps). Chaque composante des champs de vecteur déplacement et vitesse est considéré comme un champ à part entière. Ainsi après discrétisation, les inconnues du problème sont les représentations discrètes des composantes des champs de déplacement et de vitesse $\mathbf{u}_1, \mathbf{u}_2$ et $\mathbf{v}_1, \mathbf{v}_2$ respectivement.

La résolution du problème espace-temps est tout d'abord effectuée de façon incrémentale, pour des fréquences de sollicitation variant sur une large gamme (nous considérons les cas tests $1/\Delta T_{load} = 1-5-20-100$ kHz). Nous fixons les paramètres de discrétisation de façon à ce que l'erreur de discrétisation sur le champ de déplacement soit du même ordre de grandeur pour tous les cas tests (voir Tableau 1). L'erreur de discrétisation est calculée sur le domaine espace-temps avec l'indicateur suivant :

$$err_{DISC}(\mathbf{u}) = \frac{\sqrt{\|\mathbf{u}_1^e - \mathbf{u}_1\|_2^2 + \|\mathbf{u}_2^e - \mathbf{u}_2\|_2^2}}{\sqrt{\|\mathbf{u}_1^e\|_2^2 + \|\mathbf{u}_2^e\|_2^2}} \text{ avec } \mathbf{u}_1^e, \mathbf{u}_2^e \text{ les composantes du champ de déplacement calculées sur un maillage très fin.}$$

<table><thead><tr><th>1/ΔT<sub>load</sub> (kHz)</th><th>Δt (s)</th><th>err<sub>DISC</sub>(u) (%)</th></tr></thead><tbody><tr><td>1</td><td>2.5 × 10<sup>-5</sup></td><td>&lt; 2</td></tr><tr><td>5</td><td>6.25 × 10<sup>-6</sup></td><td>&lt; 2</td></tr><tr><td>20</td><td>1.5625 × 10<sup>-6</sup></td><td>&lt; 2</td></tr><tr><td>100</td><td>6.25 × 10<sup>-7</sup></td><td>&lt; 4</td></tr></tbody></table>

Tableau 1 – Paramètres de discrétisation et erreur de discrétisation associée (par rapport à la solution exacte) – $\Delta x_1 = \Delta x_2 = c_L \Delta t$ avec $c_L = \sqrt{E/\rho} = 5000$ m/s la célérité des ondes longitudinales.

La solution en espace obtenue à différents instants pour le cas test $1/\Delta T_{load}=100$ kHz est représentée sur la Figure 2. Les solutions espace-temps obtenues en déplacement et vitesse pour toutes les fréquences de sollicitation considérées sont représentées sur la Figure 3. En augmentant la fréquence de sollicitation on observe un comportement de plus en plus localisé des champs de déplacement et de vitesse dans le domaine espace-temps. Cette localisation de la solution physique est communément considérée comme problématique pour la construction de représentations à variables séparées espace-temps. Nous quantifions cet aspect dans la section suivante.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Amplitude du champ de déplacement sur le domaine spatial à différents instants.</figcaption></figure>

## 3 Décompositions espace-temps A-POSTERIORI

Ayant calculé les solutions complètes sur le domaine espace-temps, nous pouvons désormais construire leurs décompositions optimales *a-posteriori*. Une décomposition différente est construite pour chaque composante des champs de déplacement et de vitesse. Pour simplifier les notations, la définition de la décomposition optimale est donnée pour une seule composante, notée $\mathbf{u}$.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Amplitude des champs de déplacement et de vitesse sur le domaine espace-temps.</figcaption></figure>

## 3.1 Définition de la décomposition optimale

Considérons un tenseur $\mathbf{u} \in \mathbb{R}^{N(1) \times \cdots \times N(D)}$ d'ordre $D$ où $N(d)$ est l'opérateur donnant la taille de l'espace d'approximation dans la dimension $d$. Une décomposition possible de ce tenseur est appelée décomposition tensorielle canonique [14] et consiste à approximer le tenseur $\mathbf{u}$ sous la forme :

$$ \mathbf{u} \approx \mathbf{u}_M = \sum_{k=1}^{M} \underline{w}_k^1 \otimes \cdots \otimes \underline{w}_k^D \in S_M \quad (1) $$

où les termes $\underline{w}_k^1 \otimes \cdots \otimes \underline{w}_k^D$ sont appelés modes et le nombre de modes est le rang-$M$ de la décomposition. La décomposition $\mathbf{u}_M$ étant seulement une approximation du tenseur $\mathbf{u}$, il est nécessaire de définir la décomposition optimale dans une certaine norme. Ainsi, la décomposition optimale est classiquement définie comme un extremum du problème de minimisation suivant :

$$ \mathbf{u}_M = \arg \min_{\mathbf{u}^* \in S_M} \|\mathbf{u} - \mathbf{u}^*\|_2^2 \quad (2) $$

Lorsque $D = 2$, le tenseur $\mathbf{u}$ est une matrice et la décomposition optimale de rang-$M$ coïncide avec la décomposition en valeurs singulières (SVD) tronquée de la matrice. Pour les tenseurs d'ordre élevé ($D > 3$), la décomposition optimale peut être vue comme une généralisation de la décomposition en valeurs singulières [13]. Elle peut être construite avec un algorithme de type "Alternating Least Square" (ALS) [14]. Nous utilisons ici l'implémentation proposée dans [4].

**Remarque 1** Pour les composantes $\mathbf{u}_1$, $\mathbf{u}_2$ et $\mathbf{v}_1$, $\mathbf{v}_2$ des champs de déplacement et de vitesse, nous considérons deux types de décompositions (d'ordre $D = 2$ et $D = 3$ respectivement), illustrés sur la Figure 1b). La première décomposition correspond à une séparation des variables d'espace $(x_1, x_2)$ de la variable temporelle $t$ et s'écrit sous forme discrète pour les composantes du champ de déplacement (et de façon similaire en vitesse) :

$$ \forall i = 1..2 \quad \mathbf{u}_i \in \mathbb{R}^{N_{12} \times N_T} \approx \mathbf{u}_{iM} = \sum_{k=1}^{M} \underline{w}_{ik}^{12} \otimes \underline{\lambda}_{ik} \quad (\text{décomposition S12} \times T) \quad (3) $$

Pour la seconde décomposition, nous considérons également une décomposition du mode spatial, les variables $x_1$, $x_2$ et $t$ étant toutes séparées. Cette décomposition s'écrit avec $N_{12} = N_1 \times N_2$ comme suit :

$$ \forall i = 1..2 \quad \mathbf{u}_i \in \mathbb{R}^{N_1 \times N_2 \times N_T} \approx \mathbf{u}_{iM} = \sum_{k=1}^{M} \underline{w}_{ik}^1 \otimes \underline{w}_{ik}^2 \otimes \underline{\lambda}_{ik} \quad (\text{décomposition S1} \times S2 \times T) \quad (4) $$
## 3.2 Gain mémoire optimal

L'erreur de décomposition en fonction du rang-$M$ est calculé par rapport à la solution discrète de la façon suivante pour l'amplitude du déplacement (et de façon similaire en vitesse) :

$$ err(\mathbf{u}, M) = \frac{\sqrt{\|\mathbf{u}_1 - \mathbf{u}_{1M}\|_2^2 + \|\mathbf{u}_2 - \mathbf{u}_{2M}\|_2^2}}{\sqrt{\|\mathbf{u}_1\|_2^2 + \|\mathbf{u}_2\|_2^2}} $$

L'erreur obtenue pour les décompositions (3) et (4) et tous les cas tests est représentée sur la Figure 4. En comparant les courbes des décompositions $S12 \times T$ et $S1 \times S2 \times T$ à erreur fixée, on observe que les rangs des deux types de décomposition sont du même ordre de grandeur pour tous les cas tests. Aussi, on constate que plus de modes espace-temps sont nécessaires pour obtenir une erreur de décomposition donnée lorsque la fréquence de sollicitation augmente.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Erreur de décomposition *a-posteriori* pour tous les cas tests.</figcaption></figure>

Afin de quantifier cet aspect nous proposons de calculer le gain mémoire obtenu avec les décompositions optimales. Les solutions espace-temps ayant été calculées à erreur de discrétisation fixée, nous identifions le rang-$M$ tel que l'erreur de décomposition soit inférieure (ou égale) à l'erreur de discrétisation. Le gain mémoire optimal est alors défini comme le rapport entre la mémoire nécessaire pour stocker la solution complète et la mémoire nécessaire pour stocker sa décomposition optimale de même précision. Les gains obtenus pour les deux décompositions considérés et l'ensemble des cas tests sont reportés dans le Tableau 2. Même si plus de modes espace-temps sont nécessaires lorsque des sollications de fréquence élevée sont appliquées, le gain de mémoire est meilleur lorsque la fréquence de sollicitation augmente, ceci s'expliquant par la dimension des espaces d'approximations nécessairement plus grande pour obtenir une erreur de discrétisation fixée lorsque la solution présente un comportement de plus en plus localisé. Aussi, la séparation des modes spatiaux s'avère avantageuse dans les cas traités. Le gain obtenu pour la fréquence de sollicitation la plus haute est en effet un ordre de grandeur plus élevé avec la décomposition $S1 \times S2 \times T$ qu'avec la décomposition $S12 \times T$. Pour le problème de dynamique transitoire

<table><caption>Tableau 2 – Gains de mémoire obtenue avec la décomposition optimale du champ de déplacement ($M$ est choisi de façon à ce que l'erreur de décomposition soit égale (ou inférieure) à l'erreur de discrétisation).</caption><thead><tr><th rowspan="2">1/ΔT<sub>load</sub> (kHz)</th><th rowspan="2">2(N<sub>1</sub>N<sub>2</sub>N<sub>T</sub>)</th><th colspan="2">Décomp. S12xT</th><th colspan="2">Décomp. S1xS2xT</th></tr><tr><th>M</th><th>gain=<sup>N<sub>1</sub>N<sub>2</sub>N<sub>T</sub></sup>&frasl;<sub>M(N<sub>1</sub>N<sub>2</sub>+N<sub>T</sub>)</sub></th><th>M</th><th>gain=<sup>N<sub>1</sub>N<sub>2</sub>N<sub>T</sub></sup>&frasl;<sub>M(N<sub>1</sub>+N<sub>2</sub>+N<sub>T</sub>)</sub></th></tr></thead><tbody><tr><td>1</td><td>6.4×10<sup>3</sup></td><td>2</td><td><b>13</b></td><td>2</td><td><b>17</b></td></tr><tr><td>5</td><td>3.5×10<sup>5</sup></td><td>6</td><td><b>34</b></td><td>8</td><td><b>59</b></td></tr><tr><td>20</td><td>2.1×10<sup>7</sup></td><td>32</td><td><b>35</b></td><td>51</td><td><b>142</b></td></tr><tr><td>100</td><td>3.3×10<sup>8</sup></td><td>81</td><td><b>37</b></td><td>174</td><td><b>257</b></td></tr></tbody></table>

traité, la séparation *a-posteriori* des variables spatiales et temporelle permet donc de réduire de plusieurs ordres de grandeur la mémoire nécessaire au stockage de la solution espace-temps. La question est maintenant : peut-on calculer ces décompositions optimales sans avoir à stocker une solution espace-temps complète à aucun moment de la construction ? Ceci est l'objet des méthodes de décomposition *a-priori*.
# 4 Décompositions espace-temps A-PRIORI

Les méthodes de décomposition a-priori, basées sur la PGD, exploite le fait que le tenseur **u** que l'on cherche à décomposer n'est pas quelconque mais solution d'un problème dont on connait une modélisation. Pour les problèmes linéaires, cette modélisation peut s'écrire sous la forme suivante :

$$ \mathbf{A}^{\mathrm{D}}\mathbf{u} = \mathbf{b} \quad \text{avec} \quad \mathbf{A} \in \mathbb{R}^{N(1) \times N(1) \cdots N(D) \times N(D)} \quad \text{et} \quad \mathbf{b} \in \mathbb{R}^{N(1) \times \cdots \times N(D)} \qquad (5) $$

où le produit intérieur “D” est défini comme suit $\mathbf{a}^{\mathrm{D}}\mathbf{b} = \sum_{i_1, \dots, i_D} a_{i_1 \dots i_D} b_{i_1 \dots i_D}$ et $(\mathbf{A}^{\mathrm{D}}\mathbf{a})_{i_1 \dots i_D} = \sum_{j_1, \dots, j_D} A_{i_1 j_1 \dots i_D j_D} a_{j_1 \dots j_D}$. Le problème est alors dit séparable s'il est possible d'identifier une représentation séparée des opérateurs du système linéaire. Cette propriété de séparabilité s'écrit :

$$ \mathbf{A} = \sum_{k=1}^{M_A} \underline{\mathbf{A}}_k^1 \otimes \cdots \otimes \underline{\mathbf{A}}_k^D \quad \text{et} \quad \mathbf{b} = \sum_{k=1}^{M_B} \underline{\mathbf{b}}_k^1 \otimes \cdots \otimes \underline{\mathbf{b}}_k^D \qquad (6) $$

**Remarque 2** La première étape est donc d'identifier la représentation séparée des opérateurs associés au système linéaire représentatif du problème espace-temps. Dans le cas du problème d'élastodynamique considéré, ce système linéaire fait intervenir les quatre tenseurs inconnus $\mathbf{u}_1, \mathbf{u}_2$ et $\mathbf{v}_1, \mathbf{v}_2$. On obtient alors un système linéaire par blocs où chaque bloc vérifie les propriétés de séparabilité (6). La construction de ce système est décrite en détails dans [8] pour les décompositions (3) et (4).

## 4.1 Définition par minimisation du résidu dans la norme euclidienne

Différentes définitions peuvent ensuite être utilisées pour construire la PGD, les plus utilisées étant les constructions basées sur un critère d'orthogonalité de Galerkin ou par minimisation du résidu [18]. Pour les problèmes non symmétriques ($\mathbf{A} \neq \mathbf{A}'$), la convergence de la décomposition a été prouvée seulement pour la décomposition en minimum de résidu [3]. Le problème d'élastodynamique étant non-symmétrique, nous considérons seulement la décomposition construite par minimisation du résidu. Clas- siquement, le résidu est calculé dans la norme euclidienne et le problème de minimisation associé s'écrit :

$$ \mathbf{u}_M = \arg \min_{\forall \mathbf{u}^* \in S_M} \| \mathbf{b} - \mathbf{A}^{\mathrm{D}} \mathbf{u}^* \|_2^2 \qquad (7) $$

Cependant, lorsque le conditionnement de l'opérateur $\mathbf{A}$ est élevé, la décomposition obtenue par minimisation du résidu dans la norme euclidienne ne permet pas d'obtenir la décomposition optimale 2 (voir Tableau 3).

**Remarque 3** Une extension aux problèmes multi-champs de la décomposition par minimisation du résidu a été introduite dans [7]. Le calcul de la décomposition est réalisé en minimisant le résidu du problème multi-champs, formulé de façon monolithique. Le caractère espace-temps du problème permet alors de tirer avantages des propriétés de convergence des approches monolithiques tout en autorisant l'utilisation d'approximations différentes en espace et en temps pour chaque champ. Aussi, l'extension de la PGD par minimisation du résidu dans la norme idéale est détaillée dans [8]

## 4.2 Définition par minimisation du résidu dans la norme idéale

Afin de construire a-priori la décomposition optimale (2), l'idée introduite dans [6], consiste à trouver une norme “idéale” dans laquelle minimiser le résidu de façon à ce que la décomposition obtenue minimise également l'erreur. C'est-à-dire que l'on cherche un opérateur de norme $\mathbf{N}$ (définissant la norme $\|\mathbf{x}\|_{\mathbf{N}} = \sqrt{\mathbf{x}^{\mathrm{D}}\mathbf{N}^{\mathrm{D}}\mathbf{x}}$) tel que $\arg \min_{\forall \mathbf{u}^* \in S_M^D} \|\mathbf{u} - \mathbf{u}^*\|_2 = \arg \min_{\forall \mathbf{u}^* \in S_M^D} \|\mathbf{b} - \mathbf{A}^{\mathrm{D}}\mathbf{u}^*\|_{\mathbf{N}}$. Un choix idéal est alors $\mathbf{N} = (\mathbf{A}^{\mathrm{D}}\mathbf{A}')^{-1}$. La définition a-priori de la décomposition optimale par minimisation du résidu s'écrit alors :

$$ \mathbf{u}_M = \arg \min_{\forall \mathbf{u}^* \in S_M} \| \mathbf{b} - \mathbf{A}^{\mathrm{D}} \mathbf{u}^* \|_{(\mathbf{A}^{\mathrm{D}} \mathbf{A}')^{-1}}^2 \qquad (8) $$

<footnote n="2">2. L'optimalité de la décomposition dépend également de l'algorithme utilisé pour sa construction. Ici nous utilisons une version multi-champs de l'algorithme optimale basé sur l'itération du sous-espace noté (R)PGD-S dans [18].</footnote>
<table><caption>Tableau 3 – Comparaison des décompositions *a-posteriori* et *a-priori* de rang M=10 en fonction de la discrétisation (la discrétisation est fixée à $(\Delta x_1, \Delta x_2, \Delta t)/factor$ avec $\Delta t = 12.5 \times 10^{-6}$).</caption><thead><tr><th rowspan="3">S12xT<br>factor</th><th colspan="2">err(u, 10)</th><th colspan="4">(pour le cas test 1/ΔT<sub>load</sub> = 5 kHz)</th></tr><tr><th rowspan="2">ALS</th><th rowspan="2">(R)PGD-S</th><th colspan="4">(IMR)PGD-S</th></tr><tr><th>R=5</th><th>R=10</th><th>R=20</th><th>R=40</th></tr></thead><tbody><tr><td>1</td><td>7.10e-03</td><td>3.84e-02</td><td>8.99e-03</td><td>8.05e-03</td><td>7.17e-03</td><td>7.13e-03</td></tr><tr><td>2</td><td>9.24e-03</td><td>4.83e-02</td><td>1.06e-02</td><td>1.02e-02</td><td>9.71e-03</td><td>9.36e-03</td></tr><tr><td>4</td><td>9.21e-03</td><td>1.14e-01</td><td>1.52e-02</td><td>1.38e-02</td><td>1.40e-02</td><td>1.33e-02</td></tr><tr><td>8</td><td>9.11e-03</td><td>8.32e-01</td><td>1.72e-02</td><td>1.68e-02</td><td>1.66e-02</td><td>1.62e-02</td></tr></tbody></table>

Aussi, bien que cette définition donne théoriquement la décomposition optimale, l'algorithmne utilisé pour sa construction peut ne pas être optimal. Notamment, pour les décompositions d'ordre élevé, une construction progressive ne permet pas de converger vers la décomposition optimale. Aussi, nous avons proposé dans [8] un algorithme optimal de construction de la PGD idéale en minimum de résidu, basé sur l'itération du sous-espace (voir Algorithm 1) et permettant de calculer la décomposition optimale également pour les décompositions d'ordre élevé, ceci au même coût numérique qu'une construction progressive. Cet algorithme est basé sur l'introduction d'un problème auxiliaire symétrique, résolu à l'aide de la PGD par critère de Galerkin (étape (4)). Cette PGD auxiliaire est l'étape la plus coûteuse de la construction<fnref n="3" />. Différents algorithmes peuvent être utilisés pour la résolution de la PGD auxiliaire. Pour les tenseurs d'ordre élevé, il est également nécessaire d'utiliser un algorithme optimal pour la résolution de la PGD auxiliaire.

<caption>**Algorithm 1** Ideal Minimal Residual PGD with Subspace iteration algorithm – (IMR)PGD-S</caption>

1: $\mathbf{u}_M = \mathbf{0}$
2: **for** $k = 1$ to $k_{max}$ **do**
3:     $\mathbf{r} = \mathbf{b} - A^{\mathrm{D}}\mathbf{u}_M$
4:     $\mathbf{y}_R = (\mathbf{G})\mathrm{PGD}(A^{\mathrm{D}}A', \mathbf{r}, \mathbf{R})$ (choose between $(\mathbf{G})\mathrm{PGD-P}$ or $(\mathbf{G})\mathrm{PGD-S}$ algorithms)
5:     $\mathbf{u}_M = \mathrm{ALS}(\mathbf{u}_M + A'^{\mathrm{D}}\mathbf{y}_R, \mathbf{M})$
6:     **Check convergence** of $\mathbf{u}_M$
7: **end for**

Finalement, les résultats obtenus avec l'algorithmne (IMR)PGD-S montre que la précision de la décomposition idéale en minimum de résidu dépend peu de la discrétisation (voir Tableau 3). Ainsi, il est possible d'obtenir une bonne approximation de la décomposition optimale *a-priori*, la précision de la décomposition obtenue dépendant du rang-R de la PGD auxiliaire (voir la Figure 5).

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Erreur de décomposition *a-priori* pour le cas test 1/ΔT<sub>load</sub> = 5kHz, obtenue avec l'algorithmne (IMR)PGD-S et différentes valeurs du rang-R pour la résolution de la PGD auxiliaire.</figcaption></figure>

<footnote n="3">3. L'étape (5) de l'algorithmne (IMR)PGD-S est une projection sur l'espace des tenseurs de rang-M d'un tenseur de rang plus élevé. Ce tenseur étant donné sous forme séparée, cette étape est négligeable en terme de coût numérique par rapport à l'étape de résolution de la PGD auxiliaire.</footnote>
# 5 Conclusions et Perspectives

Dans cet exposé, nous montrons que les représentations par variables séparées espace-temps permettent de réduire de façon drastique la mémoire nécessaire au stockage des solutions espace-temps des problèmes d'élastodynamique transitoire (jusqu'à deux ordres de grandeur dans les cas considérés). Aussi, les décompositions optimales peuvent être calculées *a-priori* par minimisation du résidu dans une métrique idéale, ceci également pour les décompositions d'ordre élevé. Cependant, le temps de calcul associé à la construction *a-priori* de la décomposition optimale est supérieur au temps nécessaire pour une résolution incrémentale. Aussi, les développements actuels concernent l'accélération des algorithmes PGDs. Des résultats prometteurs ont notamment été obtenus avec la méthode de GMRES non-linéaire proposée dans [12] pour des constructions *a-posteriori*, et que nous avons testée *a-priori* avec la PGD. Ce travail est financé par l'Agence Nationale de la Recherche dans le cadre du projet SIM-DREAM (ANR-10-22198).

## Références

[1] Ammar A, Mokdad B, Chinesta F, Keunings R. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory. *J. Non-Newtonian Fluid Mech.* **139** (2006) 153-176

[2] Ammar A, Mokdad B, Chinesta F, Keunings R. A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II : Transient simulation using space-time separated representations. *J. Non-Newtonian Fluid Mech.* **144** (2007) 98-121

[3] Ammar A, Chinesta F, Falcó A. On the convergence of a Greedy rank-one update algorithm for a class of linear systems. *Arch. Comput. Methods Eng.* **17** (2010) 473-486

[4] Bader BW, Kolda TG. MATLAB Tensor Toolbox Version 2.5, Available online, January 2012. URL : http://www.sandia.gov/tgkolda/TensorToolbox/

[5] Beylkin G, Mohlenkamp MJ. Algorithms for numerical analysis in high dimensions. *SIAM J. Sci. Comput.* **26** (2005) 2133-2159

[6] Billaud-Friess M, Nouy A, Zahm O. A tensor approximation method based on ideal minimal residual formulations for the solution of high dimensional problems. (in preparation)

[7] Boucinha L, Gravouil A, Ammar A. Space-time proper generalized decompositions for the resolution of transient elastodynamic models. *Comput. Methods Appl. Mech. Engrg.* **255** (2013) 67-88

[8] Boucinha L, Ammar A, Gravouil A, Nouy A. Optimal proper generalized decompositions for non-symmetric high-dimensional multi-field models – Application to transient elastodynamics in space-time domains. (to be submitted)

[9] Cances E, Ehrlacher V, Lelievre T. Greedy algorithms for high-dimensional non-symmetric linear problems. (2012) arXiv:1210.6688

[10] Chinesta F, Ammar A, Lemarchand F, Beauchene P, Boust F. Alleviating mesh constraints : model reduction, parallel time integration and high resolution homogenization *Comput. Methods Appl. Mech. Engrg.* **197** (2008) 400-413

[11] Chinesta F, Ladevèze P, Cueto E. A short review on model order reduction based on Proper Generalized Decomposition. *Arch Comput Methods Eng* **404** (2011) 395-404

[12] De Sterck H. A Nonlinear GMRES Optimization Algorithm for Canonical Tensor Decomposition *SIAM J. Sci. Comput.* **34** (2012) 1351-1379

[13] Hackbusch W. Tensor Spaces and Numerical Tensor Calculus. Springer, 2012

[14] Harshman RA. Foundations of the PARAFAC procedure : Models and conditions for an "explanatory" multimodal factor analysis. *UCLA Working Papers in Phonetics* **16** (1970) 1-84

[15] Hulbert GM, Hughes TJR. Space-time finite element methods for second-order hyperbolic equations. *Comput. Methods Appl. Mech. Eng.* **84** (1990) 327-348

[16] Ladevèze P. Nonlinear computational structural mechanics : new approaches and non-incremental methods of calculation. Springer, 1999

[17] Nouy A. Proper generalized decompositions and separated representations for the numerical solution of high dimensional stochastic problems. *Arch Comput Methods Eng* **17** (2010) 403-434

[18] Nouy A. A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations. *Comput. Methods Appl. Mech. Engrg.* **199** (2010) 1603-1626