R Ésum

Géraud

Marc

Guilhem Jean-Philippe Chu

Jean-François Jérôme

A PROGRESSIVE REDUCED BASIS/EMPIRICAL INTERPOLATION METHOD FOR NONLINEAR PARABOLIC PROBLEMS

Cette thèse propose trois nouveaux développements de la méthode des bases réduites (RB) et de la méthode d'interpolation empirique (EIM) pour des problèmes non-linéaires. La première contribution est une nouvelle méthodologie, la méthode progressive RB-EIM (PREIM) dont l'objectif est de réduire le coût de la phase de construction du modèle réduit tout en maintenant une bonne approximation RB finale. L'idée est d'enrichir progressivement l'approximation EIM et l'espace RB, contrairement à l'approche standard où leurs constructions sont disjointes. La deuxième contribution concerne la RB pour les inéquations variationnelles avec contraintes non-linéaires. Nous proposons une combinaison RB-EIM pour traiter la contrainte. En outre, nous construisons une base réduite pour les multiplicateurs de Lagrange via un algorithme hiérarchique qui conserve la positivité des vecteurs de cette base. Nous appliquons cette stratégie aux problèmes de contact élastique sans frottement pour les maillages non-coïncidents. La troisième contribution concerne la réduction de modèles avec assimilation de données. Une méthode dédiée a été introduite dans la littérature pour combiner un modèle numérique avec des mesures expérimentales. Nous élargissons son cadre d'application aux problèmes instationnaires en exploitant la méthode POD-greedy afin de construire des espaces réduits pour tout le transitoire temporel. Enfin, nous proposons un nouvel algorithme qui produit des espaces réduits plus représentatifs de la solution recherchée tout en minimisant le nombre de mesures nécessaires pour le problème réduit final.

CHAPTER 1

INTRODUCTION

Ce chapitre développe plusieurs éléments permettant de situer les travaux présentés dans ce manuscrit. Tout d'abord, nous détaillons le contexte industriel ayant donné lieu à cette thèse. Ensuite, nous abordons des aspects théoriques généraux sur la réduction de modèles, en particulier sur la méthode des bases réduites. Ces éléments serviront d'introduction aux chapitres suivants pour le lecteur peu familier avec le sujet. Enfin, nous dressons un bref état de l'art suivi d'une description des principaux résultats de cette thèse. avait pour objectif d'étudier et de comparer les données obtenues lors de l'essai aux résultats des simulations numériques multi-physiques. Les accords sont des éléments de validation et les désaccords indiquent les limites des modèles et outils de calcul.

Contexte industriel

Pour la robinetterie et bien d'autres projets, plusieurs travaux sont initiés à EDF R&D sur la réduction de modèles et l'assimilation de données. Afin d'éviter des coûts d'essais prohibitifs, les constructeurs de robinets pour les centrales nucléaires ont de plus en plus recours à la simulation numérique afin d'étudier le comportement de leur matériel en conditions de qualification. EDF R&D recourt également à la simulation numérique pour évaluer les simulations des fabricants. Dans certaines configurations, même avec les ressources de calcul actuelles, il peut s'avérer difficile, voire impossible, d'utiliser les méthodes numériques standards. Ces calculs (par éléments finis) sont généralement lourds et impliquent des maillages très fins, beaucoup de pas de temps, des lois de comportement complexes et parfois la prise en compte de contraintes physiques non-linéaires comme les conditions de non-interpénétration dans les phénomènes de contact [START_REF] Meunier | On the influence of flows in clearances for thermal shocks in a globe valve[END_REF]. De plus, ces calculs sont souvent indispensables pour étudier l'influence de paramètres sur un modèle ou recaler des paramètres. De telles analyses nécessitent beaucoup d'appels à un modèle éléments finis complexe; elles sont alors rendues difficiles à cause des temps de calcul considérables qu'elles requièrent. Les méthodes de 'méta-modèles' consistent à remplacer un modèle com-plexe par un modèle approché beaucoup moins coûteux. Les méthodes dites de 'chaos polynomial' sont une famille de méthodes de méta-modèles étudiées [START_REF] Berveiller | Eléments finis stochastiques: approches intrusive et non intrusive pour des analyses de fiabilité[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF] et implémentées [START_REF] Baudin | Openturns: An industrial software for uncertainty quantification in simulation[END_REF] depuis plusieurs années à EDF R&D.

En raison de contraintes de temps et de moyens, il serait d'un grand intérêt industriel de concevoir des modélisations simplifiées, permettant des simulations à temps de calcul réduit, et dont les résultats seraient similaires (à une marge d'erreur acceptable près) à ceux des simulations plus complexes qu'EDF R&D développe en interne, et par conséquent au coût de calcul nettement plus élevé.

Réduction de modèles

Plusieurs domaines de l'ingénierie requièrent de pouvoir résoudre numériquement des équations aux dérivées partielles (EDP) modélisant des phénomènes physiques paramétriques. Deux configurations sont particulièrement récurrentes

• Les études multi-requêtes: Ce type d'études se caractérise par le besoin de résoudre plusieurs problèmes du même genre mais qui diffèrent légèrement les uns des autres; ce qui soulève naturellement la question du coût excessif de l'utilisation répétée de la même méthode numérique. Bien que ces méthodes puissent être très performantes en soi, elles n'ont pas été développées pour une utilisation répétitive.

• Les simulations en temps réel: Ce type d'études exige l'immédiateté de la résolution. Néanmoins, les méthodes standard ne parviennent pas à fournir des résultats aussi rapidement.

Les techniques de réduction de modèle [START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF][START_REF] Prud'homme | Reduced-basis output bound methods for parametrized partial differential equations[END_REF] présentent un intérêt considérable dans un tel contexte où on considère un problème dont la solution dépend de paramètres. Ainsi, pour un ensemble donné de paramètres, ces méthodes permettent d'obtenir la solution correspondante après un temps de calcul relativement court. Pour parvenir à de tels résultats, il est crucial de calculer et stocker en amont certaines quantités caractéristiques du problème qui seront utilisées lors du calcul des futures solutions. Ainsi, les méthodes de réduction de modèle se structurent en deux étapes consécutives. Pendant la première étape, dite 'hors-ligne', les calculs coûteux utilisant les modèles fins sont réalisés, permettant ainsi de créer une bibliothèque de calculs. La caractéristique principale de la phase hors-ligne est d'être effectuée une seule et unique fois. A l'opposé, la deuxième partie, dite 'en-ligne', consiste en la résolution de systèmes réduits découlant d'une reformulation du problème initial basée sur l'apprentissage hors-ligne. La phase en-ligne est donc effectuée autant de fois qu'un nouveau paramètre est choisi pour la résolution. Son coût de résolution est très faible afin d'assurer l'efficacité de l'approche. Dans cette thèse, nous nous intéressons à deux applications importantes de la réduction de modèles. Une première application est l'étude et la simulation des problèmes thermiques instationnaires. Cela est motivé conjointement par la multitude d'applications de la thermique non-linéaire instationnaire (robinetterie, générateur de vapeur, cuve, etc) et la simplicité théorique des problèmes de thermique, au moins dans le cas linéaire, par rapport aux problèmes de mécanique. La deuxième application est celle des problèmes de mécanique de contact qui sont également très importants, notamment pour les calculs chaînés thermo-mécaniques. En effet, la robinetterie a la caractéristique d'impliquer du contact unilatéral, car les robinets sont des assemblages de plusieurs pièces (cf. Figure 1.2). En outre, les quantités d'intérêt pour EDF R&D sont celles issues d'un chaînage thermo-mécanique, d'où la complémentarité des deux problèmes de réduction thermique et mécanique. Il est également à noter que les calculs mécaniques en robinetterie sont nettement plus coûteux que les calculs thermiques. En effet, les calculs thermiques se résolvent en un temps qui est de l'ordre de l'heure alors que les calculs mécaniques se résolvent en un temps de calcul qui est de l'ordre de la dizaine d'heures. Cela justifie l'intérêt porté à la réduction de modèles en mécanique de contact où les gains escomptés sont importants. De plus, il s'agit d'un domaine de recherche relativement neuf où les premiers travaux datent de 2016 [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF], alors que les premiers travaux sur la thermique non-linéaire datent de 2007 [START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF].

Assimilation de données

Plusieurs projets au sein d'EDF R&D traitent à la fois de modèles numériques et de mesures expérimentales. Ainsi, la capacité à exploiter conjointement les résultats issus des deux approches de modélisation et d'expérimentation est d'un grand intérêt.

Parmi les applications concernées à EDF R&D figurent la fatigue au département MMC, le génie civil à travers la maquette de vérification réaliste du confinement des réacteurs (VERCORS) au département MMC, la neutronique au département PER-ICLES (Performance et Prévention des Risques Industriels du Parc par la Simulation et les Études) ainsi que la robinetterie à MMC (cf. Section 1.1.1). L'application industrielle envisagée dans le futur concerne les tests liés à ce sujet. De telles campagnes d'essais sont particulièrement onéreuses et chronophages et ne peuvent donc être systématiquement envisagées. Compte tenu de ces contraintes de temps et de moyens, les constructeurs s'appuient sur des modélisations simplifiées, permettant des simulations à un temps de calcul raisonnable. Dans ce contexte, disposer à la fois d'un modèle complexe, de modèles plus ou moins simplifiés et de mesures expérimentales ouvre plusieurs perspectives. En premier lieu, la capacité d'EDF R&D à porter un regard critique sur les modèles numériques utilisés en se servant des mesures expérimentales comme d'une référence. En outre, l'assimilation de données combinée avec la réduction de modèles pourra être considérée dans l'étude d'un robinet RAMA.

Réduction de modèles

Dans cette section, nous rappelons quelques notions relatives à la réduction de modèles qui seront utiles pour la lecture de ce manuscrit. En particulier, nous détaillons la méthode des bases réduites (RBM; de l'anglais: Reduced Basis Method). La RBM est une méthode dite a priori. Elle ne nécessite aucune connaissance de la solution, car cette connaissance est construite à mesure que le calcul progresse et doit donc être réadaptée à tout nouveau type de problème. Deux méthodes de construction lors de la phase d'apprentissage hors-ligne d'une base réduite dans le cadre de la RBM seront présentées: une compression d'un ensemble de solutions précalculées par décomposition propre orthogonale (POD) [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] et un algorithme glouton -ou greedy -pour générer ces solutions de manière progressive [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF]. La POD est une méthode dite a posteriori car elle requiert des informations préalables sur la solution du problème d'intérêt.

Problème modèle

Cette section se limite aux problèmes linéaires. Nous nous intéressons à un modèle mathématique qui décrit le comportement physique d'un système en s'appuyant sur une EDP paramétrique. L'ensemble des paramètres sert à identifier une configuration particulière du système décrit par le modèle (et donc les EDPs sous-jacentes) et sera identifié comme entrée du modèle. Notons P l'ensemble de paramètres dont les éléments seront notés µ P R P , avec P ¥ 1. Ces paramètres d'entrée peuvent être des propriétés physiques de matériaux, des variables caractérisant la géométrie, des efforts variables, etc. Soit Ω un domaine borné, régulier, dans R d (d 1, 2, 3). Le domaine Ω peut dépendre de µ. Ce sera le cas lorsque nous étudierons les problèmes de contact, mais dans cette introduction, nous supposons pour simplifier que Ω ne dépend pas de µ. Soit X un espace de fonctions sur Ω muni d'un produit scalaire p¤, ¤q X et de la norme associée } ¤ } X . Nous notons X N X un sous-espace vectoriel de X de dimension N (X N est typiquement un espace d'approximation éléments finis [START_REF] Ern | Theory and practice of finite elements[END_REF]). Soient une forme bilinéaire symétrique et paramétrique a : P ¢X ¢X Ñ R et une forme linéaire paramétrique f : P ¢ X Ñ R. Les caractères bilinéaire symétrique de a et linéaire de f sont par rapport à leurs variables d'entrée dans X. Nous nous intéressons à la résolution du problème paramétrique suivant: trouver upµq X vérifiant apµ; upµq, vq f pµ; vq, dv X.

(1.1)

Le caractère bien posé du problème (1.1) est garanti par le théorème de Lax-Milgram sous les hypothèses suivantes:

(i) apµ; ¤, ¤q est coercive et continue sur X ¢ X uniformément pour tout µ P, i.e., il existe une constante positive α ¡ 0 et une constante β a V vérifiant 5 dv X, dµ P, α}v} 2 X ¤ α LB pµq}v} 2 X ¤ apµ, v, vq, dv, w X, dµ P, apµ; v, wq ¤ β a }v} X }w} X .

(1.2) (ii) f pµ; ¤q est uniformément continue pour tout µ P, i.e., il existe une constante β f V vérifiant dv X, dµ P, f pµ; vq ¤ β f }v} X .

(

Sous ces hypothèses, le problème (1.1) est bien posé. Toutefois, pour la majorité des cas pratiques, la solution exacte du problème (1.1) est inaccessible. En revanche, la projection de Galerkin de (1.1) dans X N , notée u N pµq, est accessible au calcul. Elle est solution du problème suivant: trouver u N pµq X N vérifiant apµ; u N pµq, v N q f pµ; v N q, dv N X N .

(1.4)

Ici, l'entier N désigne la dimension de l'espace X N et permet de quantifier le coût de calcul de u N pµq pour chaque µ P. Dorénavant, la solution u N pµq de (1.4) sera dite solution 'haute-fidélité'. Cette solution est habituellement calculée par des solveurs éléments finis (code_aster [START_REF] De | Finite element code aster , analysis of structures and thermomechanics for studies and research[END_REF], FreeFem++ [START_REF] Hecht | New developments in freefem++[END_REF], etc). En général, la résolution de (1.4) pour une valeur µ P est onéreuse.

La méthode des bases réduites (RBM)

Nous introduisons la notion de variété -ou manifold -de solutions M qui est l'ensemble des solutions du problème paramétrique (1.1) pour l'ensemble des paramètres P. La variété M est donc définie par M tupµq X, dµ Pu X.

(1.5)

De la même manière, nous définissons la version discrète M N de la variété de solutions M:

M N tu N pµq X N , dµ Pu X N .

(1.6)

La méthode des bases réduites (RBM) repose sur la prémisse selon laquelle la variété M N peut être approchée avec une très bonne précision par un sous-espace vectoriel 'réduit' X N X N de dimension N bien plus faible devant N . Ainsi, la RBM se base sur une projection de type Galerkin sur l'espace X N construit de telle sorte à ce qu'il approche au mieux la variété M N (cf. Figure 1.3). Ainsi, le problème (1.4) Figure 1.3 -Illustration de l'espace réduit X N et des variétés continue M et discrète M N de solutions [START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. en formulation réduite s'écrit: trouver u N pµq X N tel que apµ; u N pµq, v N q f pµ; v N q, dv N X N .

(1.7)

Soient tφ 1 , . . . , φ N u une base de l'espace éléments finis X N et tθ 1 , . . . , θ N u une base de l'espace réduit X N dont la représentation matricielle sera notée X N R N ¢N . La construction de la base tθ 1 , . . . , θ N u sera discutée dans la section 1.2.3 ci-dessous. La solution réduite u N pµq se décompose dans cette base sous la forme u N pµq °N n1 u N,n pµqθ n . Par conséquent, le problème (1.7) est reformulé sous la forme N ņ1 u n pµqapµ; θ n , θ p q f pµ, θ p q, d1 ¤ p ¤ N,

À ce stade, nous définissons la matrice réduite  et le vecteur réduit f tels que p Âpµqq np papµ; θ n , θ p qq np , et p f pµqq p pfpµ; θ p qq p , d1 ¤ n, p ¤ N, (1.9) qui sont facilement calculables via les formules

Âpµq X T N ApµqX N , et f pµq X T N f pµq, (1.10)
où la matrice A et le vecteur f sont issus du problème haute-fidélité (1.4) résolu par éléments finis et sont évalués en utilisant la base tφ 1 , . . . , φ N u. De cette manière, la RBM revient pour tout paramètre µ P à résoudre le système linéaire réduit Âpµqu N pµq f pµq.

(1.11)

Pour une résolution rapide et peu onéreuse de (1.7) ou de sa formulation matricielle (1.11), nous faisons l'hypothèse cruciale que a et f admettent une dépendance affine en le paramètre µ P, i.e. il existe Q a , Q f N ¦ petits devant N tels que apµ; v, wq Q a q1 σ a q pµqa q pv, wq, et f pµ; vq Q f q1 σ f q pµqf q pvq, (1.12) où pour tout entier q, σ a q : P Ñ R et σ f q : P Ñ R sont des fonctions du paramètre uniquement, a q : X ¢ X Ñ R est une forme bilinéaire continue et symétrique et f q : X Ñ R est une forme linéaire continue. Cette hypothèse restrictive sur a et f augmente significativement les performances des calculs numériques en évitant de recalculer ex nihilo la forme bilinéaire apµ; v, wq et la forme linéaire f pµ; vq pour chaque nouveau paramètre µ P. En effet, une fois les formes bilinéaires pa 1 , . . . , a Q a q et les formes linéaires f 1 , . . . , f Q f ¨stockées lors de la phase d'apprentissage hors-ligne, il suffit d'évaluer les quantités σ a 1 pµq, . . . , σ a Q a pµq ¨et σ f 1 pµq, . . . , σ f Q f pµq ¨pendant la phase en-ligne pour assembler ensuite les expressions (1.12) et obtenir ainsi les valeurs de apµ; v, wq et de f pµ; vq. Le raisonnement est identique pour la version matricielle (1.11) du problème. Le calcul et le stockage des matrices pA 1 , . . . , A Q a q et des vecteurs f 1 , . . . , f Q f ¨s'effectuent hors-ligne, tandis que le calcul des quantités σ a 1 pµq, . . . , σ a Q a pµq ¨et σ f 1 pµq, . . . , σ f Q f pµq ¨se fait en-ligne. Cela permet l'assemblage de l'équivalent matriciel de (1.12) donné par les formules Apµq Q a q1 σ a q pµqA q , et f pµq

Q f q1 σ f q pµqf q , (1.13)
L'hypothèse (1.12) est rarement satisfaite d'emblée. Nous verrons dans la section 1.2.5 que dans le cas général où l'EDP ne dépend pas des paramètres de manière affine, il est possible de se ramener au cadre (1.12) moyennant une méthode d'interpolation empirique (EIM); ce qui introduit une nouvelle source d'erreur d'approximation qui peut être a priori contrôlée par le biais des entiers Q a et Q f .

Génération de la base réduite

Le meilleur espace réduit X N est -sous réserve d'existence-celui qui minimise l'erreur maximale de projection, dite épaisseur de Kolmogorov d'ordre N de la variété discrète M N , et donnée par la formule

d N pM N q inf dimpW qN sup µP inf wW }u N pµq ¡ w} X .
(1.14)

Quand l'épaisseur d N est petite, unebonne approximation de la varièté discrète M N est possible par un espace X N de très petite dimension N 3 N . Toutefois, l'espace qui réalise le minimum dans la définition de d N est souvent inaccessible. Dans la littérature, on recense plusieurs stratégies de construction d'espaces réduits. Les plus largement utilisées sont la POD [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] et l'algorithme glouton -ou greedy - [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF].

Dans les deux cas, nous partons d'un ensemble d'apprentissage P tr P. À titre d'exemple, P tr peut être constitué par un échantillonnage uniforme, log-uniforme ou aléatoire de P. Cet ensemble d'apprentissage peut aussi être déterminé grâce à une connaissance physique de valeurs de paramètres pertinentes pour le problème d'intérêt.

La décomposition propre orthogonale (POD)

Supposons disposer d'un ensemble de solutions S M N de (1.4) pour un sousensemble de paramètres d'apprentissage P tr P. Si S est suffisamment riche, alors l'espace vectoriel engendré par ses éléments constitue une bonne approximation de M N . Le but de la POD est de compresser (i.e. réduire la taille) de l'ensemble S tout en conservant des bonnes propriétés d'approximation de l'espace vectoriel engendré par ses éléments. La qualité de l'approximation est quantifiée par un paramètre pod donné par l'utilisateur. Nous utilisons la notation pθ 1 , . . . , θ N q PODpS, pod q, (1.15) où S pv 1 , . . . , v R q est composé de R ¥ 1 fonctions de l'espace X N et pod est la tolérance définie par l'utilisateur. Nous adoptons une description algébrique pour plus de clarté, et nous renvoyons le lecteur à [START_REF] Gubisch | Chapter 1: Proper orthogonal decomposition for linear-quadratic optimal control[END_REF] pour un exposé détaillé de la méthode POD. Soit p 1 , . . . , N q une base de X N . Pour une fonction w X N , on note w : pw j q 1¤j¤N son vecteur de composantes dans R N , de telle sorte que w °N j1 w j j . La version algébrique de la procédure (2.29) consiste à se don- ner R vecteurs formant une matrice rectangulaire S : pv 1 , . . . , v R q R N ¢R , et à chercher N vecteurs formant la matrice rectangulaire Θ : pθ 1 , . . . , θ N q R N ¢N . Les vecteurs pθ 1 , . . . , θ N q doivent être orthonormaux par rapport à la matrice de Gram du produit scalaire de X. Dans ce contexte, on considère la matrice de Gram

C N R N ¢N telle que C N ¡ apµ 0 ; n , p q © 1¤p,n¤N , (1.16)
où la forme bilinéaire apµ 0 ; ¤, ¤q est celle utilisée dans (1.1) pour une valeur représentative µ 0 du paramètre. Ainsi, nous cherchons à vérifier θ T n C N θ p δ n,p , où δ n,p est l'indice de Kronecker, pour tous n, p t1, . . . , N u.

Soit T : pC N q 1 2 S R N ¢R et soit l'entier D minpN , Rq (en général, D R et D 3 N , ce que nous supponsons par la suite). La décomposition en valeurs singulières (SVD) [START_REF] Noble | Applied Linear Algebra[END_REF] de la matrice T retourne les réels σ 1 ¥ σ 2 ¥ ¤ ¤ ¤ ¥ σ D ¥ 0, la famille de vecteurs colonnes orthonormaux pξ n q 1¤n¤D pR N q D (tels que ξ T n ξ p δ p,n ) et la famille de vecteurs colonnes orthonormaux p ψn q 1¤n¤D pR R q D (tels que ψT n ψp δ p,n ) de telle manière que soit satisfaite la relation pθ n q 1¤n¤N sont obtenus à travers la relation

T D ņ1 σ n ξ n ψT n . ( 1
θ n pC N q ¡ 1 2 ξ n 1 σ n pC N q ¡ 1 2 T ψn 1 σ n S ψn . (1.20)
L'intérêt de la méthode POD est qu'elle permet en général de capturer un maximum d'information en norme X en utilisant très peu de modes. Par conséquent, en cherchant la solution RBM de (1.7) sous la forme

u N pµq N ņ1 u n pµqθ n , (1.21)
on obtient une troncature raisonnable en norme X puisque les vecteurs propres associés aux plus grandes valeurs propres sont conservés. Par ailleurs, une tolérance définie par l'utilisateur guide la POD dans la sélection des modes dominants. Concrètement, sont conservés les premiers modes dominants dont la somme des valeurs propres associées par rapport à la trace de la matrice d'auto-corrélation est supérieure à 1 ¡ pod pour une tolérance pod donnée (10 ¡2 par exemple).

L'algorithme glouton

Contrairement à la POD, l'algorithme glouton [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF]

µ n 1 argmax µP tr ∆ n pµq, (1.23)
calcule la solution haute-fidélité u N pµ n 1 q, et enrichit l'espace réduit pour obtenir X n 1 vecttu N pµ 1 q, . . . , u N pµ n 1 qu. Cette démarche est répétée jusqu'à atteindre une estimation d'erreur max µP tr ∆ N pµq en-dessous d'une tolérance fixée par l'utilisateur rb ¡ 0. En pratique, le critère max µP tr ∆ N pµq ¤ rb est atteint pour une valeur de N petite (N 3 N ).

Estimation d'erreur a posteriori

Un ingrédient essentiel à l'évaluation de la qualité de l'approximation d'un modèle réduit par RBM est un estimateur d'erreur garanti, précis et peu coûteux [START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF]. u n pµqapµ; θ n , vq,

Q f q1 σ f q pµqf q pvq ¡ N ņ1 u n pµq Q a q1 σ a q pµqa q pθ n , vq, (1.27)
si bien que pêpµq, vq X rpµ; vq

(1.28) Q f q1 σ f q pµqf q pvq ¡ N ņ1 Q a q1 σ a
q pµqa q pθ n , vqu n pµq .

(1.29)

Nous obtenons donc êpµq

Q f q1 σ f q pµqC q N ņ1 Q a q1 σ a q pµqL q n u n pµq , (1.30)
où C q X N et L q n X N sont les solutions éléments finis du problème 5 pC q , vq X f q pvq dv X N , pL q n , vq X ¡a q pθ n , vq dv X N .

(1.31)

La norme du résidu est alors calculée par la formule

||êpµq|| 2 X Q f q1 Q f qI 1 σ f q pµqσ f q I pµqpC q , C q I q X 2 Q f qI 1 Q a q1 N ņ1 σ f q I pµqσ a q pµqpC q I , L q n q X u n pµq (1.32) Q a q1 N ņ1 σ a q pµq 5 Q a qI 1 N ņI 1 σ a q IpµqpL q n , L q I n I q X u I n pµq C u n pµq.
Durant la phase hors-ligne, sont calculées les quantités indépendantes du paramètre, à savoir pC q , C q I q X , pC q I , L q n q X et pL q n , L q I n I q X pour 1 ¤ n, n I ¤ N et 1 ¤ q, q I ¤ Q a ou Q f . Durant la phase en-ligne, seules les quantités σ a q pµq et σ f q pµq pour 1 ¤ q ¤ Q restent à calculer; l'assemblage du résidu se fait en utilisant les calculs hors-ligne et la formule (1.32).

L'estimation précise d'une borne inférieure α LB pour la constante de coercivité est un point délicat lors de l'estimation d'erreur a posteriori. La méthode des contraintes successives (SCM) est une méthode robuste qui a été introduite pour la première fois dans [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and infsup stability constants[END_REF] afin de pallier ce problème; d'autres améliorations de cette technique sont présentées dans [START_REF] Chen | Improved successive constraint method based a posteriori error estimate for reduced basis approximation of 2D Maxwell's problem[END_REF][START_REF] Huynh | A natural-norm successive constraint method for inf-sup lower bounds[END_REF]. La méthode repose sur une stratégie hors-ligne/enligne efficace qui conduit à un problème d'optimisation linéaire. L'inconvénient de la SCM est son coût de calcul. Une stratégie alternative moins onéreuse a été récemment proposée dans [START_REF] Manzoni | Heuristic strategies for the approximation of stability factors in quadratically nonlinear parametrized PDEs[END_REF]. Cette méthode se base sur le calcul d'une approximation interpolatrice de α LB .

La méthode d'interpolation empirique (EIM)

L'objectif de la méthode d'interpolation empirique (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF] est de construire une approximation affine du type décrit dans (1.12). Comme indiqué plus haut, ce type d'approximation est crucial pour le succès d'une décomposition hors-ligne/enligne dans le cadre de la méthode des bases réduites. Dans le contexte le plus général, l'EIM a été introduite pour approcher une fonction continue bivariée γ :

P ¢ Ω Ñ R.
Cette approximation est décrite par un opérateur d'interpolation γ M qui interpole la fonction γ en des points d'interpolation px i q 1¤i¤M dans Ω tr Ω en tant que combinaison linéaire de fonctions pq j q 1¤j¤M avec q j : Ω Ñ R. Ces dernières fonctions ne sont pas des fonctions polynomiales ou trigonométriques mais sont déduites directement à partir de la famille de fonctions tγpµ; ¤q, dµ P tr u par des combinaisons linéaires de M 'échantillons' γpµ 1 ; ¤q, . . . , γpµ M ; ¤q, où les paramètres µ 1 , . . . , µ M P tr sont sélectionnés par un algorithme glouton dans un ensemble d'apprentissage P tr P (qui peut être différent de celui introduit plus haut). Ainsi, pour tout paramètre µ P et tout x Ω, l'approximation de rang M , notée γ M , est donnée par la formule γ M pµ; xq M j1 ϕ j pµqq j pxq.

(1.33) Le but de la séparation de variables dans (1.33) est de permettre le calcul des fonctions pq j q 1¤j¤M qui sont indépendantes de µ pendant la phase hors-ligne afin de n'avoir que les fonctions pϕ j q 1¤j¤M qui dépendent de µ à évaluer lors de phase en-ligne; celles-ci vérifient la propriété d'interpolation

γ M pµ, x i q M j1 ϕ j pµqq j px i q, d1 ¤ i ¤ M. (1.34)
Pour des fonctions à valeurs réelles v définies sur Ω, on définit }v} V pΩ tr q : max xΩ tr |vpxq|. Soit un compteur d'itérations m ¥ 1 et une fonction γ m¡1 définie sur P tr ¢ Ω, avec la convention γ 0 0. Une itération EIM se définit comme suit: premièrement, on définit µ m P tr par µ m argmax µP tr }γpµ; ¤q ¡ γ m¡1 pµ; ¤q} V pΩ tr q .

q m p¤q : r m p¤q r m px m q , (1.37)
et on calcule la dernière ligne de la matrice d'interpolation B donnée par

B mi : pq i px m qq, d1 ¤ i ¤ m.
(1.38)

La construction de l'approximation EIM ainsi décrite satisfait trois propriétés importantes:

(i) Les fonctions de base pq j q 1¤j¤M sont linéairement indépendantes.

(ii) La matrice d'interpolation B est triangulaire inférieure à diagonale unitaire, donc inversible.

(iii) Pour M grand, l'approximation γ M tend vers γ en norme infinie V pP tr ¢Ω tr q.

1.2.6 Brève bibliographie sur la RBM et ses applications L'idée d'utilisation de bases réduites a été introduite dans des travaux anciens [START_REF] Almroth | Automatic choice of global shape functions in structural analysis[END_REF][START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF]. Forgée sous sa forme actuelle, la RBM a été étudiée dans les travaux fondateurs [START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF][START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF] et deux monographies récentes présentent en détail ses différents aspects [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. La RBM a été appliquée à un large éventail de problèmes. Par exemple, elle a été utilisée pour la quantification d'incertitudes [START_REF] Barrett | On the reduced basis method[END_REF][START_REF] Chen | Sparse-grid, reduced-basis Bayesian inversion[END_REF] et pour l'échantillonnage par éléments finis avec comme but l'accélération des calculs spatiotemporels d'un problème d'élasticité [START_REF] Hoang | An efficient goaloriented sampling strategy using reduced basis method for parametrized elastodynamic problems[END_REF]. La méthode RBM a également été étendue aux problèmes de contrôle de flux optimaux pour les équations elliptiques noncoercives [START_REF] Negri | Reduced basis approximation of parametrized optimal flow control problems for the Stokes equations[END_REF], après avoir permis de générer des sorties en temps réel des équations de Navier-Stokes incompressibles avec des estimateurs d'erreur [START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF]. Un autre estimateur d'erreur a posteriori a été développé dans [START_REF] Casenave | Accurate and online-efficient evaluation of the a posteriori error bound in the reduced basis method[END_REF] en utilisant la méthode d'interpolation empirique (EIM) afin de surmonter les erreurs d'arrondi. La méthode EIM est utilisée dans plusieurs contextes et a donné lieu à plusieurs variantes [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Maday | A generalized empirical interpolation method: application of reduced basis techniques to data assimilation[END_REF]. Concernant l'application à la modélisation thermique, une autre application de la RBM a été la création de solutions fiables en temps réel pour la modélisation d'ailettes thermiques et autres structures complexes [START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF]. Par ailleurs, des travaux récents présentent une version adaptée à la solution des problèmes décrits par les EDP paraboliques [START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF]. Bien que dans un premier temps l'algorithme glouton fut utilisé seul pour construire l'espace réduit X N dans le cas parabolique, il est généralement couplé avec une POD pour améliorer l'efficacité de la méthode. Cette méthode couplée relativement récente a été baptisée POD-greedy dans la littérature [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF].

Objectifs de la thèse

Cette thèse a pour objectif de traiter trois problématiques principales.

1. Il est communément admis dans la littérature que la phase de calcul amont de la méthode des bases réduites, dite phase 'hors-ligne', peut être onéreuse, voire très onéreuse. Ce critère est rarement pris en compte dans les études académiques car la phase hors-ligne est réalisée une fois pour toutes. Cependant, le contexte industriel détaillé plus haut (et bien d'autres) impose des contraintes de ressources et de coût de calcul même pendant cette phase. Notre objectif sera donc de remédier à ce problème en concevant une méthode qui optimise la phase hors-ligne sans pour autant dégrader la qualité de l'approximation réduite finale. Cette problématique n'a été que très peu étudiée jusqu'à présent [START_REF] Daversin | Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF][START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF].

2. Les méthodes de réduction de modèles sont pensées pour réduire des équations variationnelles. L'extension de ces travaux aux inégalités variationnelles a fait l'objet de quelques travaux récents [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF][START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Haasdonk | A reduced basis method for parametrized variational inequalities[END_REF]. Une application d'intérêt majeur est la mécanique du contact. Néanmoins, les travaux sur ce sujet supposent des hypothèses simplificatrices que sont la linéarité des contraintes et la coïncidence des maillages. Par conséquent, notre second objectif sera d'établir une stratégie de réduction de modèles applicable pour une classe plus générale de problèmes de contact élastique.

3. Le troisième objectif de cette thèse est de valoriser les campagnes d'essai et leurs résultats expérimentaux pour des études par simulation numérique tout en réduisant les moyens et temps de calcul. Nous considérons plus particulièrement une stratégie d'assimilation de données couplée à la réduction de modèles [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF] que nous étendrons aux problèmes de thermique non-linéaire instationnaire.

État de l'art et principaux résultats

Nous présentons finalement un résumé des travaux effectués dans le cadre de cette thèse. Ces travaux sont présentés plus amplement dans les trois chapitres suivants. Dans le premier de ces chapitres, nous nous concentrons sur un problème de thermique non-linéaire instationnaire pour lequel nous concevons une méthode de réduction du temps de calcul hors-ligne afin de diminuer davantage les ressources requises pour la réduction de cette classe de problèmes. Dans le deuxième chapitre, nous étudions ensuite la réduction de modèles pour un problème de mécanique de contact. Finalement, le troisième chapitre est dédié à l'élaboration d'une stratégie d'assimilation de données pour les problèmes instationnaires, dont la classe de problèmes traitée au premier chapitre.

Réduction de modèle en thermique non-linéaire

Nous nous intéressons aux problèmes paraboliques non-linéaires pour lesquels une méthode RBM a été développée dans [START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. Ces travaux ne prennent pas en compte l'éventuel coût de la phase offline. Concernant l'optimisation des coûts hors-ligne de cette classe de problèmes, quelques travaux récents existent dans la littérature. L'idée d'un enrichissement progressif de l'approximation EIM et de l'espace réduit a été récemment proposée dans [START_REF] Daversin | Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF] pour les EDP non-linéaires stationnaires. La construction progressive de l'EIM et de la base réduite a été également abordée dans [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. Dans cette méthode, le critère d'enrichissement est commun à la fois à l'EIM et à la base réduite, et l'échantillon maximisant un estimateur d'erreur a posteriori est sélectionné pour enrichir les deux bases. Une autre méthode est développée dans [START_REF] Urban | Affine decompositions of parametric stochastic processes for application within reduced basis methods[END_REF]; on y propose une construction progressive de l'EIM à l'aide d'approximations basées sur une POD des trajectoires haute-fidélité. Nous étudions dans le chapitre 2 de nouveaux développements de la méthode des bases réduites et de la méthode d'interpolation empirique (RB-EIM) pour des problèmes paraboliques non-linéaires. Nous développons une nouvelle méthodologie: la méthode progressive RB-EIM (PREIM) pour les problèmes non-linéaires paraboliques. Ici, le but est de réduire le coût hors-ligne final tout en maintenant une bonne approximation RB dans la phase en-ligne. L'idée de base est un enrichissement progressif de l'approximation EIM et de l'espace RB, contrairement à l'approche standard où l'approximation EIM et l'espace RB sont construits séparément. PREIM utilise des calculs haute-fidélité chaque fois que ceux-ci sont disponibles et des calculs RB dans le cas contraire. Ce chapitre correspond à un article paru chez SIAM Journal on Scientific Computing (SISC).

Réduction de modèle pour les problèmes de contact

Dans la littérature, trois articles récents abordent la réduction de modèle pour les inéquations variationnelles. Le premier article [START_REF] Haasdonk | A reduced basis method for parametrized variational inequalities[END_REF] étend la méthode RBM standard aux inégalités variationnelles linéaires résolues par une formulation mixte. En ce qui concerne la construction des bases, les bases 'primale' (pour la solution primale) et duale (pour les multiplicateurs de Lagrange) sont directement composées d'échantillons bien choisis. Aucune phase de compression supplémentaire n'est considérée. Dans la méthode dite 'Projection-Based' (PB) [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] qui a été introduite spécifiquement pour résoudre les problèmes de contact avec des contraintes linéaires, les bases primale et duale sont construites différemment. Une base primale est obtenue par POD, comme c'est souvent le cas pour la réduction des égalités variationnelles. Vu que [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] se concentre sur des problèmes de dynamique avec de nombreux instants pour réaliser l'échantillonnage temporel, l'ensemble des échantillons des multiplicateurs de Lagrange s'avère souvent assez important. Sa compression devient alors une préoccupation majeure. L'idée est de construire la base duale réduite en appliquant l'algorithme de factorisation de matrices non-négatives (NMF) [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] à l'ensemble des échantillons de multiplicateurs de Lagrange. La NMF garantit la positivité des vecteurs de base ainsi qu'une dimension de base relativement faible, mais la base duale qui en résulte est beaucoup moins précise que la base primale. De plus, l'utilisateur ne spécifie pas une tolérance requise mais un nombre de vecteurs de base dominants à conserver. Enfin, [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF] étend un autre type de méthodes de réduction de modèles, dites méthodes d'hyper-réduction (HR), aux problèmes de contact avec des contraintes linéaires. L'extension proposée de la méthode HR aux problèmes de contact consiste à conserver peu de vecteurs de la base duale haute-fidélité, car le nombre de noeuds de contact est limité à un domaine d'intégration réduit (RID). Par conséquent, seuls les noeuds de contact du RID sont traités, mais avec une haute fidélité locale. Jusqu'à présent, tous les résultats existants se limitent aux contraintes linéaires.

Nous étudions dans le chapitre 3 de nouveaux développements de la RBM pour les inéquations variationnelles avec contraintes non-linéaires. Nous proposons une méthode de base réduite combinée à la méthode d'interpolation empirique pour traiter la contrainte non-linéaire. Dans ce contexte, une base réduite 'primale' est nécessaire pour la solution primale et une base réduite 'duale' est requise pour les multiplicateurs de Lagrange. Nous proposons de construire cette dernière en utilisant un algorithme hiérarchique qui conserve la non-négativité des vecteurs de la base duale. Cette stratégie de réduction est ensuite appliquée aux problèmes de contact élastique sans frottement pour les maillages non-coïncidents. Ce type de maillages induit une non-linéarité de discrétisation et son traitement est d'autant plus important que le cas de non-coïncidence est très majoritaire dans la pratique. Ces deux dernières conditions issues de la physique du problème n'ont pas encore été abordées, à notre connaissance, dans la littérature.

Assimilation de données

Une méthode d'assimilation de données a été récemment couplée à la réduction de modèle, à savoir la méthode dite 'Parameterized-Background Data-Weak' (PBDW) [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF]. Elle a pour but de traiter les divergences pouvant apparaître entre les prédictions d'une simulation numérique basée sur un modèle numérique et des mesures expérimentales. Cette problématique est fréquemment rencontrée dans la pratique, et il est judicieux de prendre en compte les données expérimentales afin d'améliorer les modèles mathématiques. La méthode PBDW est une formulation variationnelle pour les problèmes d'assimilation de données modélisés par des équations aux dérivées partielles. Plusieurs tests numériques indiquent une grande amélioration de la prédiction par le modèle mathématique grâce à la prise en compte des données expérimentales. Une étude approfondie des liens entre la méthode PBDW et d'autres techniques d'assimilation de données est proposée dans [START_REF] Taddei | An adaptive parametrized-background data-weak approach to variational data assimilation[END_REF][START_REF] Taddei | A localization strategy for data assimilation; application to state estimation and parameter estimation[END_REF]. D'autres travaux ont par la suite émergé dans ce même contexte. Dans [START_REF] Argaud | Sensor placement in nuclear reactors based on the generalized empirical interpolation method[END_REF], une technique basée sur la Generalized EIM (GEIM) est proposée pour le placement des capteurs de mesures expérimentales.

Le chapitre 4 a été réalisé à l'occasion d'une collaboration avec le Pr. Patera initiée par une mobilité d'un mois au Massachusetts Institute of Technology (MIT). Il est dédié à une extension de la méthode PBDW. Initialement introduite pour les problèmes elliptiques, ce chapitre élargit son cadre d'application aux problèmes paraboliques. Pour ce faire, nous exploitons la POD-greedy [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF] qui permet d'obtenir des espaces réduits convenables pour tout le transitoire temporel. Ensuite, nous présentons un nouvel algorithme pour la phase hors-ligne de la méthode PBDW. Cet algorithme présente deux intérêts majeurs. D'un côté, il permet d'obtenir des espaces réduits plus représentatifs de la solution recherchée. D'un autre côté, il conduit à une diminution significative du nombre de mesures expérimentales nécessaires lors de la phase en-ligne et en accélère davantage la résolution.

Développements informatiques

Cette thèse a donné lieu à plusieurs développements informatiques. Ceux-ci sont capitalisés dans deux outils phares d'EDF R&D. Une première partie est capitalisée dans la 'forge Pléiade', qui est un dépôt de codes informatiques permettant l'accès et la modification de codes sources à plusieurs acteurs impliqués dans un même projet. Quant à la seconde partie, elle est intégrée à la plateforme MAP (Materials Ageing Platform). Cette plateforme permet l'automatisation du code avec un processus de validation industriel plus avancé, une documentation et un accès à un public d'utilisateurs plus large. Ce processus inclut des tests réguliers de non-régression au fur et à mesure de l'évolution des versions des différents langages de programmation et outils impliqués dans l'élaboration du code. Ainsi, les codes informatiques réalisés pendant cette thèse se divisent comme suit: • Dans la forge Pléiade sont déposés les codes restants, à savoir les fichiers sources de réduction de modèles pour la mécanique de contact élastique nonfrottant et l'assimilation de données.

Introduction

The Reduced-Basis (RB) method devised in [START_REF] Machiels | Output bounds for reduced-basis approximations of symmetric positive definite eigenvalue problems[END_REF][START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF] (see also the recent textbooks [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]) is a computationally effective approach to approximate parametrized Partial Differential Equations (PDEs) encountered in many problems in science and engineering. For instance, the RB method is often used in real-time simulations, where a problem needs to be solved very quickly under limited computational resources, or in multi-query simulations, where a problem has to be solved repeatedly for a large number of parameter values. Let P denote the parameter set. The RB method is split into two stages: (i) an offline stage where a certain number of so-called High-Fidelity (HF) trajectories are computed for a training subset of parameters P tr P(typically a finite element space based on a fine mesh); (ii) an online stage for real-time or multi-query simulations where the parameter set P is explored more extensively. The output of the offline phase includes an approximation space of small dimension spanned by the so-called RB functions. The reduced space then replaces the much larger HF space in the online stage. The crucial point for the computational efficiency of the overall procedure is that computations in the HF space are allowed only in the offline stage.

In the present work, we are interested in nonlinear parabolic problems for which a RB method has been successfully developed in [START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF]. A key ingredient to treat the nonlinearity so that the online stage avoids HF computations is the Empirical Interpolation Method (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF]. The EIM provides an approximation of the nonlinear (or non-affine) terms in the PDE. This approximation is built using a greedy algorithm as the sum of M functions, where the dependence on the space variable is separated from the dependence on the parameter (and the time variable for parabolic problems). The integer M is called the rank of the EIM and controls the accuracy of the approximation. Although the EIM is performed during the offline stage of the RB method, its cost can become a critical issue since the EIM can require an important number of HF computations for an accurate approximation of the nonlinearity. The cost of the EIM typically scales with the size of the training set P tr .

The goal of the present work is to overcome this issue. To this purpose, we devise a new methodology, the Progressive RB-EIM (PREIM) method, which aims at reducing the computational cost of the offline stage while maintaining the accuracy of the RB approximation in the online stage. The key idea is a progressive enrichment of both the EIM approximation and the RB space, in contrast to the standard approach, where the EIM approximation and the RB space are built separately. In PREIM, the number of HF computations is at most M , and it is in general much lower than M in a time-dependent context where the greedy selection of the pair pµ, kq to build the EIM approximation (where µ is the parameter and k refers to the discrete time node) can lead to repeated values of µ for many different values of k. In other words, PREIM can select multiple space fields within the same HF trajectory to build the EIM space functions. In this context, only a modest number of HF trajectories needs to be computed, yielding significant computational savings with respect to the standard offline stage. PREIM is driven by convergence criteria on the quality of both the EIM and the RB approximation, as in the standard RB-EIM procedure. PREIM is devised in order to have a guaranteed termination, and in the worst-case scenario, the same number of HF trajectories is computed as in the standard RB-EIM algorithm, thus reaching the same level of accuracy for the representation of the nonlinearity and the construction of the RB functions (if this level of accuracy turns out to be insufficient, the parameter training set has to be enlarged as usual in the standard algorithm). In this worst-case scenario, the computational cost of PREIM may be slightly larger than that of the standard algorithm because of the way the intermediate calculations of trajectories are organized in PREIM. However, we expect that in many practical situations, e.g., when the computation of HF trajectories dominates the cost of the progressive construction of the EIM, PREIM can bring computational benefits with respect to the standard approach. These benefits, which are particularly sizeable whenever the nonlinearity can be represented by an EIM approximation of relatively modest rank, are illustrated in this work on three test cases, including one derived from a three-dimensional industrial prototype. Yet, the present study remains heuristic, and a theoretical analysis of the possible computational gains of PREIM can be pursued in future work.

The idea of a progressive enrichment of both the EIM approximation and the RB space has been recently proposed in [START_REF] Daversin | Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF] for stationary nonlinear PDEs, where it is called Simultaneous EIM/RB (SER). Thus, PREIM extends this idea to timedependent PDEs. In addition, there is an important difference in the greedy algorithms between SER and PREIM. Whereas SER uses only RB computations, PREIM uses HF computations whenever available, both for the greedy selection of the parameters and the time nodes, as well as for the space-dependent functions in the EIM approximation. These aspects are particularly relevant since they improve the accuracy of the EIM approximation. This is illustrated in our numerical experiments on nonlinear parabolic PDEs. The progressive construction of the EIM and the RB has been recently addressed within the Empirical Interpolation Operator Method in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF]. Therein, the enrichment criterion is common to both the EIM and the RB, and the snapshot maximizing an a posteriori error estimator is selected to enrich both bases. Instead, PREIM has dedicated criteria for the quality of the EIM approximation and for the RB approximation. Furthermore, PREIM systematically exploits the knowledge of the HF trajectories whenever available, and an update step is performed in order to confirm the current parameter selection. We also mention the Proper Orthogonal Empirical Interpolation Method from [START_REF] Urban | Affine decompositions of parametric stochastic processes for application within reduced basis methods[END_REF], where a progressive construction of the EIM approximation is devised using POD-based approximations of the HF trajectories.

The chapter is organized as follows. In Section 2.2, we introduce the model problem. In Section 2.3, we briefly recall the main ideas of the nonlinear RB method devised in [START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF], and in Section 2.4, we briefly recall the EIM procedure in the standard offline stage as devised in [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF]. The reader familiar with the material can jump directly to Section 2.5, where PREIM is introduced and discussed. Section 2.6 presents numerical results illustrating the performance of PREIM on nonlin-ear parabolic problems related to heat transfer including a three-dimensional valve prototpye for flow regulation. Finally, Section 2.7 provides a technical complement.

Model problem

In this section, we present a prototypical example of a nonlinear parabolic PDE. The methodology we propose is illustrated on this model problem but can be extended to other types of parabolic equations. We consider a spatial domain (open, bounded, connected subset) Ω R d , d ¥ 1, with a Lipschitz boundary, a finite time interval I r0, T s, with T ¡ 0, and a parameter set P R p , p ¥ 1, whose elements are generically denoted by µ P. Our goal is to solve the following nonlinear parabolic PDE for many values of the parameter µ P: find u µ : I ¢ Ω Ñ R such that

fu µ ft ¡ ∇ ¤ pκ 0 Γpµ, u µ qq∇u µ ¨ f, in I ¢ Ω, ¡ κ 0 Γpµ, u µ q ¨fu µ
fn φ e , on I ¢ fΩ, u µ pt 0, ¤q u 0 p¤q, in Ω, (2.1) where κ 0 ¡ 0 is a fixed positive real number, Γ : P ¢ R Ñ R is a given nonlinear function, f : I ¢ Ω Ñ R is the source term, φ e : I ¢ fΩ Ñ R is the time-dependent Neumann boundary condition on fΩ, and u 0 : Ω Ñ R is the initial condition. For simplicity, we assume without loss of generality that f , φ e , and u 0 are parameterindependent. We assume that f L 2 pI; L 2 pΩqq and φ e L 2 pI; L 2 pfΩqq (this means that f ptq L 2 pΩq and φ e ptq L 2 pfΩq for (almost every) t I), and we also assume that u 0 H 1 pΩq. We make the standard uniform ellipticity assumption β 1 ¤ κ 0 Γpµ, zq ¤ β 2 with 0 β 1 β 2 V, for all pµ, zq P ¢ R. With the above assumptions, it is reasonable to look for a weak solution u µ L 2 pI; Y q H 1 pI; Y I q.

Remark 2.1 (Initial condition). For parabolic PDEs, the initial condition is often taken to be in a larger space, e.g., u 0 L 2 pΩq. Our assumption that u 0 Y is motivated by the RB method, where basis functions in Y are sought as solution snapshots in time and for certain parameter values. In this context, we want to include the possibility to select the initial condition as a RB function.

Remark 2.2 (Heat transfer).

One important application we have in mind for (4.68) is heat transfer problems. In this context, the PDE can take the slightly more general form

αpu µ q fu µ ft ¡ ∇ ¤ pκ 0 Γpµ, u µ qq∇u µ ¨ f, in I ¢ Ω,
where αpu µ q stands for the mass density times the heat capacity. Moreover, the quantity pκ 0 Γpµ, u µ qq represents the thermal conductivity. Note also that φ e ¡ 0 means that the system is heated.

In practice, one way to solve (4.68) is to use a Y -conforming Finite Element Method [START_REF] Ern | Theory and practice of finite elements[END_REF] to discretize in space and a time-marching scheme to discretize in time.

The Finite Element Method is based on a finite element subspace X Y defined on a discrete nodal subset Ω tr Ω, where CardpΩ tr q N . To discretize in time, we consider an integer K ¥ 1, we let 0 t 0 ¤ ¤ ¤ t K T be pK 1q distinct time nodes over I, and we set K tr t1, ¤ ¤ ¤ , Ku, K tr t0u K tr , I tr tt k u kK tr , and ∆t k t k ¡ t k¡1 for all k K tr . As is customary with the RB method, we assume henceforth that the mesh-size and the time-steps are small enough so that the above space-time discretization method delivers HF approximate trajectories within the desired level of accuracy. These trajectories, which then replace the exact trajectories solving (4.68), are still denoted u µ for all µ P. Henceforth, we use the convention that the superscript k always indicates a time index; thus, we write u k µ p¤q u µ pt k , ¤q X, f k p¤q f pt k , ¤q L 2 pΩq, and φ k e p¤q φ e pt k , ¤q L 2 pfΩq.

Applying a semi-implicit Euler scheme, our goal is, given u 0

µ u 0 X, to find pu k µ q kK tr X K such that, for all k K tr , dv X, mpu k µ , vq ∆t k a 0 u k µ , v ¨ ∆t k n Γ µ, u k¡1 µ , v ¨ mpu k¡1 µ , vq ∆t k l k pvq, (2.2)
with the bilinear forms m :

Y ¢ Y Ñ R, a 0 : Y ¢ Y Ñ R and the linear forms l k : Y Ñ R such that mpv, wq » Ω vw, a 0 pv, wq κ 0 » Ω ∇v ¤∇w, l k pvq » Ω f k v » fΩ φ k e v, (2.3)
and the nonlinear form

n Γ : P ¢ Y ¢ Y Ñ R such that n Γ pµ, v, wq » Ω Γpµ, vq∇v ¤ ∇w, (2.4)
for all µ P and all v, w Y . In (2.2), the nonlinearity is treated explicitly, whereas the diffusive term is treated implicitly. This choice avoids dealing with a nonlinear solver at each time-step. The computation of derivatives of discrete operators within Newton's method is addressed, e.g., in [START_REF] Drohmann | Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation[END_REF].

The Reduced-Basis method

In this section, we briefly recall the Reduced-Basis (RB) method for the nonlinear problem (2.2) [START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Grepl | Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations[END_REF]. Let XN X be a so-called reduced subspace such that N dimp XN q 3 dimpXq N . Let pθ n q 1¤n¤N be a Y -orthonormal basis of XN . For all µ P and k K tr , the RB solution ûk µ XN that approximates the HF solution u k µ X is decomposed as

ûk µ N ņ1 ûk µ,n θ n , (2.5)
with real numbers ûk µ,n for all n t1, . . . , N u. Let us introduce the component vector ûk µ : pû k µ,n q 1¤n¤N R N , for all µ P and k K tr . Let û0 be the Y -orthogonal projection of the initial condition u 0 X onto XN with associated component vector û0 R N . Replacing u k µ X in the weak form (2.2) by the approximation ûk µ XN with associated component vector ûk µ R N , and using the test functions pθ p q 1¤p¤N , we obtain the following problem written in algebraic form: Given û0

µ û0 R N , find pû k µ q kK tr pR N q K such that, for all k K tr , pM ∆t k A 0 qû k µ ∆t k f k Mû k¡1 µ ¡ ∆t k gpû k¡1 µ q, (2.6)
with the matrices M, A 0 R N ¢N and the vectors

f k R N such that M ¡ mpθ n , θ p q © 1¤p,n¤N , A 0 ¡ a 0 pθ n , θ p q © 1¤p,n¤N , f k l k pθ p q ¨1¤p¤N , (2.7)
and the vector gpû k¡1 µ q R N such that

gpû k¡1 µ q ¢ N ņ1 ûk¡1 µ,n » Ω Γ ¢ µ, N ņI 1 ûk¡1 µ,n I θ n I ∇θ n ¤ ∇θ p 1¤p¤N . (2.8)
The difficulty is that the computation of gpû k¡1 µ q requires a parameter-dependent reconstruction using the RB functions in order to compute the integral over Ω.

To avoid this, we need to build an approximation γ M of the nonlinear function

γ : P ¢ K tr ¢ Ω Ñ R such that γpµ, k, xq : Γpµ, u k µ pxqq, (2.9)
in such a way that the dependence on x is separated from the dependence on pµ, kq. More precisely, for some integer M ¡ 0, we are looking for an (accurate) approximation γ M : P ¢ K tr ¢ Ω Ñ R of γ under the separated form

γ M pµ, k, xq : M j1 ϕ k µ,j q j pxq, (2.10)
where M is called the rank of the approximation and ϕ k µ,j are real numbers that we find by interpolation over a set of M points tx 1 , . . . , x M u in Ω tr by requiring that

γ M pµ, k, x i q γpµ, k, x i q Γpµ, u k µ px i qq, di t1, ¤ ¤ ¤ , M u.
(2.11)

The interpolation property (2.11) is achieved by setting

ϕ k µ,j pB ¡1 γ k µ q j , dj t1, ¤ ¤ ¤ , M u, where γ k µ : Γpµ, u k µ px i qq ¨1¤i¤M R M , (2.12)
and B pq j px i qq 1¤i,j¤M R M ¢M must be an invertible matrix. Therefore, (2.10) can be rewritten as follows:

γ M pµ, k, xq M j1
pB ¡1 γ k µ q j q j pxq.

(2.13)

The points px i q 1¤i¤M in Ω tr and the functions pq j q 1¤j¤M defined on Ω are determined by the EIM algorithm [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] which is further described in Section 2.4 below.

Let us now describe how we can use the EIM approximation (2.13) to allow for an offline/online decomposition of the computation of the vector gpû k¡1 µ q defined in (3.44). Under the (reasonable) assumptions ûk µ u k µ and Γpµ, ûk

µ pxqq Γ µ, u k µ pxq ¨, (2.14)
we obtain Γpµ, ûk

µ pxqq Γ µ, u k µ pxq ¨ γpµ, k, xq γ M pµ, k, xq (2.15)
M j1 pB ¡1 γ k µ q j q j pxq M j1 pB ¡1 γk µ q j q j pxq, with the vector γk µ : pΓpµ, ûk µ px i qqq 1¤i¤M R M . The problem (2.6) then becomes:

Given û0 µ û0 R N , find pû k µ q kK tr pR N q K such that, for all k K tr , pM ∆t k A 0 qû k µ ∆t k f k M ¡ ∆t k D k¡1 µ ¨û k¡1 µ , (2.16
)

with the matrix D k¡1 µ R N ¢N such that D k¡1 µ M j1 B ¡1 γk¡1 µ ¨jC j , (2.17)
where

C j ¢ » Ω q j ∇θ n ¤ ∇θ p 1¤p,n¤N R N ¢N , d1 ¤ j ¤ M.
(2.18)

The overall computational procedure can now be split into two stages:

(i) An offline stage where one precomputes on the one hand the RB functions pθ n q 1¤n¤N leading to the vectors û0 R N , pf k q kK tr pR N q K and the matrices M, A 0 R N ¢N , and on the other hand the EIM points px i q 1¤i¤M and the functions pq j q 1¤j¤M leading to the matrices B R M ¢M and C j R N ¢N , for all j t1, . . . , M u. The offline stage is discussed in more detail in Section 2.4.

(ii) An online stage to be performed each time one wishes to compute a new trajectory for a parameter µ P. All that remains to be performed is to compute the vector p γ k¡1 µ R M and the matrix D k¡1 µ R N ¢N and to solve the N -dimensional linear problem (2.16) for all k K tr . The online stage is summarized in Algorithm 2.1.

The standard offline stage

There are two tasks to be performed during the offline stage:

(T 1 ) Build the rank-M EIM approximation (2.10) of the nonlinear function γ defined by (2.9);

Algorithm 2.1 Online stage

Input : µ, pθ n q 1¤n¤N , û0 , pf k q kK tr , M, A 0 , px i q 1¤i¤M , pq j q 1¤j¤M , and pC j q 1¤j¤M 1: Set k 1 and û0 µ û0 2: while k K tr do 3:

Compute D k¡1 µ using (2.17 Set q m : r m {r m px m q and compute pB mi q 1¤i¤m by setting B mi : pq i px m qq 7:

γ 0 0 3: Search pµ m , k m q argmax pµ,kqP tr ¢K tr }Γpµ, u k µ p¤qq ¡ γ m¡1 pµ, k, ¤q} V pΩ tr q 4: Set r m p¤q : Γpµ m , u km µm p¤qq ¡ γ m¡1 pµ m , k m ,
Set m m 1

8:

Search pµ m , k m q argmax pµ,kqP tr ¢K tr }Γpµ, u k µ p¤qq ¡ γ m¡1 pµ, k, ¤q} V pΩ tr q 9:

Set r m p¤q : Γpµ m , u km µm p¤qq ¡ γ m¡1 pµ m , k m , ¤q and x m argmax xΩ tr |r m pxq| 10: end while 11: Set M : m ¡ 1 Output : px i q 1¤i¤M and pq j q 1¤j¤M (T 2 ) Explore the solution manifold in order to construct a linear subspace XN X of dimension N .

In the standard offline stage, these two tasks are performed independently. Let us first discuss Task (T 1 ), i.e., the construction of the rank-M EIM approximation. Recall from Section 2.3 that the goal is to find the interpolation points px i q 1¤i¤M in Ω tr Ω and the functions pq j q 1¤j¤M with q j : Ω Ñ R. The construction of the EIM approximation additionally uses a training set P tr P for the parameter values; in what follows, we denote by P the cardinality of P tr . For a real-valued function v defined on Ω tr , we define }v} V pΩ tr q : max xΩ tr |vpxq|. Given an iteration counter m ¥ 1 and a function γ m¡1 defined on P tr ¢ K tr ¢ Ω, with the convention that γ 0 0, an EIM iteration consists of the following steps. First, one defines pµ m , k m q P tr ¢ K tr by pµ m , k m q argmax pµ,kqP tr ¢K tr }Γpµ,

u k µ p¤qq ¡ γ m¡1 pµ, k, ¤q} V pΩ tr q , (2.19)
where we notice the use of the HF trajectories for all values of the parameter µ in the training set P tr . Once pµ m , k m q has been determined, one sets

r m p¤q : Γpµ m , u km µm p¤qq ¡ γ m¡1 pµ m , k m , ¤q, x m argmax xΩ tr |r m pxq|, (2.20)
and one checks whether |r m px m q| ¡ eim for some user-defined positive threshold eim . If this is the case, one sets

q m p¤q : r m p¤q r m px m q , (2.21)
and one computes the new row of the matrix B by setting B mi : pq i px m qq, for all 1 ¤ i ¤ m. The standard EIM procedure is presented in Algorithm 2.2.

Let us now briefly discuss Task (T 2 ) above, i.e., the construction of a set of RB functions with cardinality N . First, as usual in RB methods, the solution manifold is explored by considering a training set for the parameter values; for simplicity, we consider the same training set P tr as for the EIM approximation. This way, one only explores the collection of points tu k µ u pµ,kqP tr ¢K tr in the solution manifold. For this exploration to be informative, the training set P tr has to be chosen large enough. The exploration can be driven by means of an a posteriori error estimator (see, e.g., [START_REF] Rovas | Reduced-basis output bound methods for parabolic problems[END_REF]) which allows one to evaluate only N HF trajectories. However, in the present setting where HF trajectories are to be computed for all the parameters in P tr when constructing the EIM approximation, it is natural to exploit these computations by means of a Proper Orthogonal Decomposition (POD) [START_REF] Hinze | Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: error estimates and suboptimal control[END_REF][START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] to define the RB. This technique is often considered in the literature to build the RB in a time-dependent setting, see, e.g., [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. In practice, a POD of the whole collection of snapshots may be computationally demanding (or even unfeasible) when a very large number of functions is considered. Thus, we adopt a POD-based progressive construction of the reduced basis in the spirit of the PODgreedy algorithm from [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF]. Therein, one additional RB function is picked at a time, whereas here we can pick more than one function. The progressive construction of the RB is presented in Algorithm 2.3, where we have chosen an enumeration of the parameters in P tr from 1 to P . The initialization of Algorithm 2.3 is made by computing pθ n q 1¤n¤N 1 PODpS 1 , pod q for the trajectory S 1 associated with the parameter µ 1 . That is, we select the first N 1 POD modes out of the set S 1 with error threshold pod (for completeness, this procedure is briefly outlined in Section 2.7).

The next steps of the algorithm are performed in an iterative fashion. For each new trajectory, we first subtract its projection on the current RB, and then perform a POD on the projection and merge the result with the current RB. This specific part of the procedure, called UPDATE RB, is presented in Algorithm 2.4; this part of the procedure is presented separately since it will be re-used later on.

Algorithm 2.3 Progressive RB

Input : P tr , K tr , and pod ¡ 0 1: Compute S p ¨1¤p¤P pu k µp q kK tr ¨1¤p¤P P HF trajectories 2: Compute pθ n q 1¤n¤N 1 PODpS 1 , pod q 3: Set p 1 4: while p P do 5:

Set p p 1

6:

Compute pθ n q 1¤n¤N p UPDATE RB pθ n q 1¤n¤N p¡1 , S p , pod 7:

end while 8: Set N : N P 9: Compute û0 , pf k q kK tr , M, and A 0 10: Compute the matrices pC j q 1¤j¤M

Output : pθ n q 1¤n¤N , û0 , pf k q kK tr , M, A 0 , and pC j q 1¤j¤M

Remark 2.3 (Threshold pod ). For the initialization (line 2 of Algorithm 2.3), one can use a relative error threshold for pod (for instance, pod 1%). Instead, for the iterative loop (line 6 of Algorithm 2.3), the threshold pod can be set to the greatest singular value that has been truncated at the initialization step.

Remark 2.4 (Order of EIM and RB). Algorithms 2.2 and 2.3 can be performed in whatever order. If Algorithm 2.3 is performed first, the computation of the matrices pC j q 1¤j¤M is postponed to the end of Algorithm 2.2. Moreover, the HF trajectories pu k µ q pµ,kqP tr ¢K tr appearing in both algorithms are computed only once.

The Progressive RB-EIM method (PREIM)

In this section, we first present the main ideas of the PREIM algorithm. Then we describe one important building block called UPDATE EIM. Finally, using this building block together with the procedure UPDATE RB from Algorithm 2.4, we present the PREIM algorithm.

Main ideas

PREIM consists in a progressive construction of the EIM approximation and of the RB. The key idea is that, unlike the standard EIM for which HF trajectories are computed for all the parameter values in the training set P tr (Algorithm 2.2, line 1), PREIM works with an additional training subset P HF m P tr that is enriched progressively with the iteration index m of PREIM. The role of P HF m is to collect the parameter values for which a HF trajectory has already been computed. PREIM is designed such that CardpP HF m q ¤ m for all m t1, . . . , M u. This means that when the final rank-M EIM approximation has been computed, at most M HF trajectories have been evaluated, whence the computational gain with respect to the standard offline stage provided M 3 P . m . The RB trajectories can be modified at each iteration m of PREIM. This happens whenever a new value of the parameter is selected in the greedy stage of the EIM so that the approximation of the nonlinearity is modified. To reflect this dependency, we add a superscript m to the RB trajectories which are now denoted pû m,k µ q kK tr for all µ P tr zP HF (2.23)

The goal of every PREIM iteration is twofold:

(i) produce a set of RB functions pθ m n q 1¤n¤N m (the RB functions and their number depend on m); p φm q k µ,j qj pxq.

(2.24)

The notation γm rP HF m ,Xm,Qms in (2.24) indicates the data rP HF m , X m , Q m s that is used to build the approximation of the nonlinearity. More precisely, this construction uses the PREIM training set P HF m , the sequence of interpolation points X m : px i q 1¤i¤m

in Ω tr (with xm computed at iteration m), and the sequence of functions Q m :

pq j q 1¤j¤m defined on Ω (with qm computed at iteration m). The progressive construction of these three ingredients is described below. Then, considering the (invertible) lower-triangular matrix B R m¢m whose last row is calculated using Bmj qj px m q for all j t1, . . . , mu, we compute the real numbers p φm q k µ,j in (2.24) from the relations m j1 Bij p φm q k µ,j γm pµ, k, xi q, di t1, . . . , mu,

for all pµ, kq P ¢ K tr . All the real numbers p φm q k µ,j depend on m since the righthand side of (2.25) depends on m.

The procedure UPDATE EIM Algorithm 2.5 UPDATE EIM

Input : pθ n q 1¤n¤N m¡1 , P HF in , X m¡1 , Q m¡1 , and eim 1: Compute pū k µ q pµ,kqP tr ¢K tr using pθ n q 1¤n¤N m¡1 2: Search pµ m , k m q argmax pµ I ,k I qP tr ¢K tr }Γ µ I , ūk

I µ Ip¤q ¨¡ γm¡1 rP HF in ,X m¡1 ,Q m¡1 s pµ I , k I , ¤q} V pΩ tr q based on RB/HF 3: Define r r m p¤q Γ µ m , ūkm µm p¤q ¨¡ γm¡1 rP HF in ,X m¡1 ,Q m¡1 s pµ m , k m , ¤q. 4: if µ m P HF in then 5:
Compute S out pu k µm q kK tr and set P HF out P HF in tµ m u one HF trajectory

6:

Search pμ m , km q argmax pµ I ,k I qP HF out ¢K tr }Γ µ I , u k I µ Ip¤q ¨¡γ m¡1 rP HF in ,X m¡1 ,Q m¡1 s pµ I , k I , ¤q} V pΩ tr q 7: else 8:

Set S out r, P HF out P HF in , and pμ m , km q pµ m , k m q 9: end if 10: Define rm p¤q : Γ μm , u km μm p¤q ¨¡ γm¡1

rP HF in ,X m¡1 ,Q m¡1 s pμ m , km , ¤q 11: if }r m } V pΩ tr q eim then 12:

Set incr rk = FALSE Set incr rk = TRUE 17:

Define r m p¤q rm p¤q 18:

Set X out pX m¡1 , xm q and Q out pQ m¡1 , qm q with xm , qm as in Algorithm 2.2 (lines 6 and 9). 19: end if 20: Define δ eim m }r m } V pΩ tr q Output : incr rk, P HF out , X out , Q out , S out , and δ eim m An essential building block of PREIM is the procedure UPDATE EIM described in Algorithm 2.5. The input is the RB functions pθ n q 1¤n¤N m¡1 , the triple rP HF in , X m¡1 , Q m¡1 s describing the current approximation of the nonlinearity (the choice for the indices will be made clearer in the next section, and is not important at this stage), and the threshold eim . The output is the flag incr rk which indicates whether or not the rank of the EIM approximation has been increased, and if incr rk TRUE, the additional output is the triple rP HF out , X out , Q out s to devise the new EIM approximation, possibly a new HF trajectory S out , and a measure δ eim m on the EIM error. First (see line 2), one selects a new pair pµ m , k m q P tr ¢ K tr in a greedy fashion as follows:

pµ m , k m q argmax pµ I ,k I qP tr ¢K tr }Γ µ I , ūk I µ Ip¤q ¨¡ γm¡1 rP HF in ,X m¡1 ,Q m¡1 s pµ I , k I , ¤q} V pΩ tr q . (2.
pμ m , km q argmax pµ I ,k I qP HF out ¢K tr }Γ µ I , u k I µ Ip¤q ¨¡ γm¡1 rP HF in ,X m¡1 ,Q m¡1 s pµ I , k I , ¤q} V pΩ tr q . (2.27)
We notice that this re-selection criterion only handles HF trajectories since the parameter values are in P HF out . Moreover, (2.27) only requires to probe the values for µ m , since the values for the other parameters, which are in P HF in , have already been evaluated in (2.26). Finally, to prevent division by small quantities, the value of the residual }r m } V pΩ tr q is checked in line 11. If this value is too small, the pair pμ m , km q is rejected and the rank of the EIM approximation is not increased.

The PREIM algorithm

We are now ready to describe the PREIM procedure, see Algorithm 2.6. PREIM is an iterative method that builds progressively the RB and the EIM approximation. The iteration is controlled by three tolerances: pod ¡ 0 which is used in the progressive increment of the RB, eim ¡ 0 which is used to check the quality of the EIM approximation, and rb ¡ 0 which is used to check the quality of the RB. The termination criterion involves the quality of both the EIM and the RB approximations, see line 7. Note that this is the same criterion as in the standard RB-EIM approach.

Within each PREIM iteration, the two previously-described procedures UP-DATE EIM and UPDATE RB are called. First, one attempts to improve the EIM approximation (line 9). If this is successful (i.e., if incr rk TRUE), the RB is updated by using the possibly new HF trajectory S m (line 21). Otherwise (i.e., if incr rk FALSE), the RB is possibly updated (line 12) and a new improvement of the EIM is attempted (line 13). In general, the RB is improved because a new HF trajectory has been computed. Whenever this is not the case, a new HF trajectory is anyway computed in line 18 (cf. Remark 2.5 below). The choice of this new HF trajectory can be driven by a standard greedy algorithm based on the use of a classical a posteriori error estimator. More precisely, for a given reduced basis pθ n q 1¤n¤N and given sets of training points X and functions Q used for the current EIM approximation of the nonlinearity, the associated a posteriori error estimator for a given value of the parameter µ P is denoted by ∆ X ,Q pθnq 1¤n¤N pµq. Finally, we observe that the reduced matrices and vectors in line 22 of Algorithm 2.6 need to be updated since these quantities depend on the RB functions which can change at every iteration.

Remark 2.5 (Worst-case scenario). The worst-case scenario is that in which PREIM would compute as many trajectories as the standard EIM. In this situation, the RB space would be identical to that of the standard RB-EIM. Regarding the approximation of the nonlinearity, if PREIM is carried on until M M max : P ¢K, the resulting rank-M approximation would be exact for all the parameters in P tr . Hence, as RB , POD , and eim tend to zero, RB-EIM and PREIM produce the same approximations at termination (recall that termination is guaranteed for both algorithms).

Let us now discuss the initialization of PREIM. In line 2, one can choose an initial PREIM training set P HF 1 composed of a single parameter, as is often the case with greedy algorithms. Although the nonlinearity may not be well-described initially, one can expect that the description will improve progressively. Still, to allow for more robustness in the initialization, one can consider an initial set P HF 1 composed of several parameters. One can then compute the HF trajectories for all µ P HF 1 and compress them using the POD procedure with threshold pod (if P HF 1 contains more than one value, a progressive version is used). Finally, one selects

pµ 1 , k 1 q argmax pµ I ,k I qP HF 1 ¢K tr }Γ µ I , u k I µ Ip¤q ¨} V pΩ tr q ,
(2.28) one defines r 1 p¤q Γpµ 1 , u k 1 µ 1 p¤qq and computes X 1 px 1 q, Q 1 pq 1 q (as in the standard EIM procedure), and one sets δ eim 1 }r 1 } V pΩ tr q . Let us finally point out that a good initialization of PREIM can favor its early termination. For instance, one can try to select the first parameter as one for which the nonlinearity has a sizeable effect.

Remark 2.6 (PREIM-NR and U-SER variants). We can consider two variants in the procedure UPDATE EIM (Algorithm 2.5) and therefore in PREIM. A first variant consists in skipping the re-selection step in line 6 of Algorithm 2.5. This variant, which we call PREIM-NR (for 'no re-selection'), will be tested numerically in the next section in order to highlight the actual benefits brought by the re-selection.

A second variant is to replace ūm,k µ with ûm,k µ in lines 1 and 2 of Algorithm 2.5, and to skip the re-selection step in line 6. We call this variant U-SER since it can be viewed as an extension of SER [START_REF] Daversin | Simultaneous empirical interpolation and reduced basis method for non-linear problems[END_REF] to the unsteady setting. The crucial difference between PREIM-NR and U-SER is that U-SER uses RB trajectories to compute the space-dependent functions in the EIM approximation, whereas PREIM-NR uses HF trajectories.

Algorithm 2.6 Progressive RB-EIM (PREIM)

Input : P tr , K tr , Ω tr , pod ¡ 0, eim ¡ 0, and rb ¡ 0 1: Set m 1 2: Choose P HF 1 P tr of cardinality J and compute S 1 pu k µ q pµ,kqP HF 1 ¢K tr J ¥ 1 HF trajectories

3: Compute pθ 1 n q 1¤n¤N 1 PODpS 1 , pod q. 4: Compute û0 R N 1 , pf k q kK tr pR N 1 q K , M R N 1 ¢N 1 , and A 0 R N 1 ¢N 1 5: Compute pX 1 , Q 1 , δ eim 1 q INIT EIMpP HF 1 q and C 1 R N 1 ¢N 1 6: Compute δ rb 1 max µP tr ∆ 1 pµq 7: while (δ eim m ¡ eim or δ rb m ¡ rb ) do 8:
Set m m 1 and P HF in : P HF m¡1 9:

(incr rk,

P HF out , X out , Q out , S out δ eim m ) = UPDATE EIM (pθ m¡1 n q 1¤n¤N m¡1 , P HF in , X m¡1 , Q m¡1 , eim ) 10:
while incr rk = FALSE do

pθ m¡1 n q 1¤n¤N m¡1 UPDATE RB pθ m¡1 n q 1¤n¤N m¡1 , S out , pod 13: (incr rk, P HF out , X out , Q out , S out , δ eim m ) = UPDATE EIM (pθ m¡1 n q 1¤n¤N m¡1 , P HF in , X m¡1 , Q m¡1 , eim ) 14:
if incr rk = TRUE then 15:

Step to line 20 Compute pθ m n q 1¤n¤N m UPDATE RB pθ m¡1 n q 1¤n¤N m¡1 , S m , pod

22:

Update û0 R N m , pf k q kK tr pR N m q K , and the matrices M, A 0 , pC j q 1¤j¤m in R N m ¢N m

23:

Compute δ rb m max µP tr ∆ Xm,Qm pθ m n q 1¤n¤N m pµq 24: end while 25: Set M : m

Output : pθ n q 1¤n¤N M , û0 , pf k q kK tr , M, A 0 , X M , Q M , and pC j q 1¤j¤M

Numerical results

In this section, we illustrate the above developments on three test cases related to transient heat transfer problems. The first two test cases use the idealized 2D geometry of a perforated square plate; the first test case involves a nonlinearity on the solution, whereas the second test case considers a nonlinearity on its partial derivatives. The third test case is based on the three-dimensional geometry of an industrial valve prototype used for flow regulation in nuclear reactor operation, while we use the same type of nonlinearity as in the first test case. Our goal is to illustrate the computational performance of PREIM and to compare it to the standard EIM approach described in Section 2.4 and to the variants PREIM-NR and U-SER described in Remark 2.6. HF trajectories are computed using a Finite Element subspace X Y H 1 pΩq consisting of continuous, piecewise affine functions. The HF computations use the industrial software code_aster [START_REF] De | Finite element code aster , analysis of structures and thermomechanics for studies and research[END_REF] for the first test case, FreeFem++ [START_REF] Hecht | New developments in freefem++[END_REF] for the second test case, and a combination of the industrial software Salomé and FreeFem++ for the third test case. The reduced-order modeling algorithms have been developed in Python. In all the test cases, the dominant error component turns out to be the one resulting from the approximation of the nonlinearity, rather than the one resulting from the RB. For this reason, PREIM has been run using only the stopping criterion δ eim m ¡ eim in line 7 of Algorithm 2.6.

Test case (a): Nonlinearity on the solution

We consider a two-dimensional setting based on the perforated plate illustrated in

¨2©

, with u 0 293 K (20 o C) and u m 323 K (50 o C). We define κ 0 1.05 m 2 ¤K ¡2 ¤s ¡1 and φ e 3 K¤m¤s ¡1 (these units result from our normalization by the density times the heat capacity). For space discretization, we use a mesh containing N 1438 nodes (see Figure 2.1). Regarding time discretization, we consider the time interval I r0, 5s, the set of discrete times nodes K tr t1, ¤ ¤ ¤ , 50u, and a constant time step ∆t k 0.1 s for all k K tr . Finally, we consider the parameter interval P r1, 20s, the training set P tr t1, ¤ ¤ ¤ , 20u, and we use the larger set t0.25i | 0 ¤ i ¤ 80u to verify our numerical results. In Figure 2.2, we show the HF temperature profiles over the perforated plate at two different times and for two different parameter values. We can see that, as the simulation time increases, the temperature is, overall, higher for larger values of the parameter µ than for smaller values. Also, for larger values of µ, the temperature variation tends to be less uniform over the plate than for smaller values of µ. We now investigate PREIM, which we first run with thresholds pod 10 ¡3 and eim 5 ¤ 10 ¡2 . Table 2.2 shows the selected parameters and discrete time nodes at each stage of PREIM. We can make two important observations from this table. First, after 13 iterations, PREIM has only selected four different parameter values, and has therefore computed only four HF trajectories (the iterations for which a new parameter value is selected are indicated in gray in Table 2.2). In the other 9 out of the 13 iterations, a different time snapshot of an already existing HF trajectory has been selected. Second, by comparing the lines in Table 2.2 related to µ and μ, we can see that a parameter re-selection happened at iteration m 7. . Left: Errors for the standard RB-EIM procedure, PREIM, and U-SER. Right: Errors }r r m } V pΩ tr q and }r m } V pΩ tr q for PREIM.

The left panel of Figure 2.3 displays the error on the approximation of the nonlinear function Γ for the standard RB-EIM procedure and for PREIM as a function of the iteration number m (the additional curve concerning U-SER will be commented afterwards), i.e., we plot }r m } V pΩ tr q (line 3 of Algorithm 2.5) and }r r m } V pΩ tr q (line 10 of Algorithm 2.5) as a function of m, see (2.24). We can see that the quality of the approximation of the nonlinearity is almost the same for PREIM as for the standard RB-EIM procedure; yet, the former achieves this accuracy by computing 20% of the HF trajectories computed by the latter (4 instead of 20 HF trajectories). The right panel of Figure 2.3 shows the values of }r r m } V pΩ tr q and }r m } V pΩ tr q as a function of m.

The two quantities differ when the parameter µ m in line 2 of Algorithm 2.5 is not in the set P HF m¡1 so that }r r m } V pΩ tr q is computed using a RB approximation, whereas }r m } V pΩ tr q results from a HF trajectory. Discarding the initialization, this happens for m t6, 9, 10u. The fact that }r r m } V pΩ tr q and }r m } V pΩ tr q take rather close values indicates that the RB provides an accurate approximation of the HF trajectory. The error is generically denoted }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq . We observe an excellent agreement over the whole parameter range. In the right panel of Figure 2.4, we also consider the space-time errors on the trajectories produced using the approximation of the nonlinearity resulting from PREIM with the RB resulting from the standard algorithm. We do not observe any significant change with respect to the left panel, which indicates that the dominant error component is that associated with the approximation of the nonlinearity. We consider the tighter couple of thresholds pod 10 ¡5 and eim 5¤10 ¡3 in Figure 2.5. Here, we can observe some differences in the errors produced by the standard RB-EIM and PREIM procedures, although both errors remain comparable and reach similar maximum values over the parameter Figure 2.5 -Test case (a): RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq for pod 10 ¡5 and eim 5 ¤ 10 ¡3 . range. While the standard procedure is slightly more accurate for most parameter values, the conclusion is reversed for some other values. Moreover, the curves on the right panel of Figure 2.5 corroborate the fact that once again, the dominant error component is that associated with the approximation of the nonlinearity. Let us further explore the PREIM algorithm by comparing it to its variants U-SER and PREIM-NR introduced in Remark 2.6. Table 2.3 reports the selected parameters and time nodes in U-SER and PREIM-NR (compare with Table 2.2 for PREIM). Both U-SER and PREIM-NR need to compute five HF trajectories, which is only 25% of those needed with the standard RB-EIM procedure, but this is still one more HF trajectory than with PREIM. One difference between U-SER and PREIM-NR is that new parameters are selected earlier with U-SER. Interestingly, after 13 iterations, U-SER and PREIM-NR have selected the same five parameters.

m
Another interesting observation is that U-SER actually selects the same couple pµ, kq twice (this happens for m 2 and m 6); this can be interpreted by observing that owing to the improvement of the RB using HF trajectories between iterations m 2 and m 6, the algorithm detects the need to improve the approximation of the nonlinearity by using the same pair pµ, kq. The same observation can be made for PREIM-NR (this happens for m 4 and m 8). We emphasize that re-selecting the same pair pµ, kq is not possible within PREIM since the selection is based on HF trajectories. The left panel of Figure 2.3 displays the error on the approximation of the nonlinear function Γ obtained with U-SER and compares it to the error obtained with the standard RB-EIM and PREIM procedures that were already discussed. The U-SER error is evaluated as sup pµ,kqP tr ¢K tr }Γpµ, ûk µ p¤qq¡γ m rP HF m ,Xm,Qms pµ, k, ¤q} V pΩ tr q .

We observe that the approximation of the nonlinearity is somewhat less sharp with U-SER than with PREIM. Figure 2.6 reports the space-time errors (measured using the 2 -norm in time and the H 1 -norm in space) on the trajectories produced by PREIM and U-SER for the whole parameter range. We observe that the U-SER error is always larger, sometimes up to a factor of five, but for the larger parameter values which produce the larger errors, the quality of the results produced by PREIM and U-SER remains comparable. Finally, we provide an assessment of the runtimes in Table 2. [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF]. We can see that for the standard RB-EIM procedure, the computation of the HF trajectories dominates the cost of the offline phase. For both PREIM and U-SER, the cost of these HF computations is substantially reduced. At the same time, the cost of the greedy algorithm (which includes the construction of the EIM and of the RB) is increased by 50% with respect to the standard RB-EIM procedure. However, the impact on the total runtime is marginal. We consider the nonlinear parabolic problem (4.68) with the nonlinear function Γpµ, zq : sin ¡ ωµ p fu fx q 2 p fu fy q 2 ¨©2 , where ω 6.25 ¤ 10 ¡3 . We define u 0 293 K (20 o C), κ 0 1 m 2 ¤K ¡2 ¤s ¡1 and φ e 3 K¤m¤s ¡1 (these units result from our nor-malization by the density times the heat capacity). For the space discretization, we use a mesh containing N 1429 nodes. Regarding time discretization, we consider the time interval I r0, 2.5s, the set of discrete times nodes K tr t1, . . . , 50u, and a constant time step ∆t k 0.05 s for all k K tr . Finally, we consider the parameter interval P r1, 40s and the training set P tr t1, . . . , 40u. In Figure 2.7, we show the temperature isovalues over the perforated plate at two different times for µ 1. We can observe different boundary layers depending on the time (the same observation can be made by varying the parameter value). During the standard offline stage, we perform P 40 HF computations. Knowing that K 50, the set S (Algorithm 2.2, line 1) contains 2040 fields, each consisting of N 1429 nodal values. Applying Algorithm 2.3 based on the H 1 -norm, a truncation threshold pod 5 ¤ 10 ¡2 , and parameters enumerated with increasing values, we obtain N 15 RB functions. Table 2.5 shows the dimension of the RB space as a function of the enumeration index p. Table 2.6 shows the evolution of the error on the nonlinearity within the standard EIM for selected values of the rank of the EIM approximation. The fact that the nonlinearity depends on the partial derivatives of the solution challenges the EIM; indeed, the error decay is not as fast as in the previous test case. This observation is corroborated by the fact that the functions pq j q 1¤j¤M all look quite different (not shown for brevity).

Test case (b): Nonlinearity on the partial derivatives

p
We now investigate the performance of PREIM, which we run with thresholds pod 5 ¤ 10 ¡2 and either eim 10 ¡1 or eim 10 ¡3 . Table 2.7 shows the selected parameters and time nodes at each iteration. For eim 10 ¡1 , PREIM performs 9 iterations, and three parameters are selected for HF computations, whereas for eim 10 ¡3 , PREIM performs 11 further iterations and six more HF computations to reach the requested threshold. Moreover, the evolution of the size of the reduced basis within PREIM is shown in Table 2.8 for selected values of the rank of the EIM approximation. As can be noticed, the approximation of the nonlinearity requires more computational effort than that of the solution manifold.

Figure 2.8 shows the decrease of the EIM approximation error on the nonlinearity for PREIM with pod 5 ¤ 10 ¡2 and eim 10 ¡3 . We observe that each time a new HF trajectory is computed, i.e., whenever the quantities }r r m } V pΩ tr q and }r m } V pΩ tr q differ, the difference is actually rather small, thereby confirming the already accurate approximation of the nonlinearity by the RB solutions in PREIM. The left panel of Figure 2.9 illustrates the space-time errors (measured in the 2 pI tr ; H 1 pΩ tr qq-norm) on the trajectories produced by the standard RB-EIM and the PREIM procedures for the whole parameter range. We observe that for lower parameter values, PREIM delivers somewhat less accurate results, whereas the conclusion is reversed for higher parameter values. Altogether, both errors stay within comparable upper bounds. The right panel of Figure 2.9 shows that the error component associated with the approximation of the nonlinearity is still the dominant one, except for the parameter range r1, 5s, where the RB from the standard algorithm improves the error. Incidentally, we observe that these smaller values of the parameter were not selected within PREIM for approximating the nonlinearity. Finally, Figure 2.10 shows the same results for the tighter thresholds pod 5 ¤ 10 ¡2 and eim 10 ¡4 . Here, 14 HF computations and 100 PREIM iterations were needed. We can see that the PREIM error closely matches that of the standard RB-EIM procedure.

Test case (c): 3D industrial valve prototype

Here, we present a three-dimensional test case whose geometry is based on a flow regulation valve used in nuclear reactor operation. We consider the nonlinear parabolic problem (4.68) with the nonlinear function Γpµ, vq : sin }r m } V pΩ tr q 1.0 5.5E¡1 2.7E¡1 5.2E¡2 1.1E¡2 7.5E¡3 1.5E¡3 We now investigate PREIM, which we first run with thresholds pod 10 ¡3 and eim 5 ¤ 10 ¡3 . Table 2.10 shows the selected parameters and discrete time nodes at each stage of PREIM. Out of ten iterations, we can see that a parameter re-selection happened at four iterations. For some of the remaining iterations, the selected time nodes have been changed at the re-selction step. Figure 2.13 -Test case (c): Left: EIM approximation errors }r r m } V pΩ tr q and }r m } V pΩ tr q as a function of m for PREIM with pod 10 ¡4 and eim 5¤10 ¡3 . Right: RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq {}u µ } 2 pI tr ;H 1 pΩ tr qq for pod 10 ¡3 and eim 5 ¤ 10 ¡3 .

The left panel of Figure 2.13 shows the decrease of the EIM approximation error on the nonlinearity within PREIM for pod 10 ¡4 and eim 5 ¤ 10 ¡3 . Overall, the difference between }r r m } V pΩ tr q and }r m } V pΩ tr q is rather small except for m 9 where RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq {}u µ } 2 pI tr ;H 1 pΩ tr qq . Left: pod 10 ¡4 and eim 5 ¤ 10 ¡1 . Right: pod 10 ¡4 and eim 5 ¤ 10 ¡3 . values. Altogether, both errors stay within comparable upper bounds. Figure 2.14 displays the space-time errors for the tighter tolerance pod 5 ¤ 10 ¡4 that delivers more accurate reluts as expected. Moreover, tightening eim makes both the RB-EIM and PREIM errors overlap. Thus, one can infer that the dominant error in this test case is rather the RB error; whence the numerous HF computations that need to be performed, as seen in Table 2 Finally, we provide an assessment of the runtimes in Table 2.11. One can notice that the greedy procedure accounts for a slightly greater percentage of the offline stage in PREIM compared to the standard RB-EIM. This is mainly due to the additional intermediate calculations in PREIM. However, as previously shown, the dominant part of the offline stage are the HF computations; this illustrates again the relevance of using PREIM.

Technical complement : Proper Orthogonal Decomposition

The goal of this complement is to briefly describe the procedure associated with the notation pθ 1 , . . . , θ N q PODpS, pod q,

(2.29)
which is used in Algorithms 2.3, 2.4, and 2.6, where S pv 1 , . . . , v R q is composed of R ¥ 1 functions in the space X and pod is a user-prescribed tolerance. For simplicity, we adopt an algebraic description, and we refer the reader to [START_REF] Gubisch | Chapter 1: Proper orthogonal decomposition for linear-quadratic optimal control[END_REF] for further insight. Let p 1 , . . . , N q be a basis of X, where dimpXq N . For a function w X, we denote by w : pw j q 1¤j¤N its coordinate vector in R N , so that w °N j1 w j j . The algebraic counterpart of (2.29) is that we are given R vectors forming the rectangular matrix S : pv 1 , . . . , v R q R N ¢R , and we are looking for N vectors forming the rectangular matrix Θ : pθ 1 , . . . , θ N q R N ¢N . The vectors pθ 1 , . . . , θ N q are to be orthonormal with respect to the Gram matrix of the inner product in X. In the present setting, we consider the Gram matrix

C N R N ¢N such that C N ¡ mp n , p q ηa 0 p n , p q © 1¤p,n¤N , (2.30)
where η ¡ 0 is a user-prescribed weight and the bilinear forms m and a 0 are defined in (2.3). Thus, we want to have θ T n C N θ p δ n,p , the Kronecker delta, for all n, p t1, . . . , N u.

Let us set T : pC N q 1 2 S R N ¢R and consider the integer D minpN , Rq (in most situations, we have D R and D 3 N ). Computing the Singular Value Decomposition [START_REF] Noble | Applied Linear Algebra[END_REF] of the matrix T, we obtain the real numbers σ 1 ¥ σ 2 ¥ ¤ ¤ ¤ ¥ σ D ¥ 0, the orthonormal family of column vectors pξ n q 1¤n¤D pR N q D (so that ξ T n ξ p δ p,n ) and the orthonormal family of column vectors p ψn q 1¤n¤D pR R q D (so that ψT n ψp δ p,n ), and we have

T D ņ1 σ n ξ n ψT n . (2.31)
From (2.31), it follows that T ψn σ n ξ n and T T ξ n σ n ψn for all n t1, . . . , Du.

The vectors we are looking for are then given by θ n : pC N q ¡ 1 2 ξ n for all n t1, . . . , N u with N : maxt1 ¤ n ¤ D | σ n ¥ POD u. It is well-known that the N -dimensional space spanned by the vectors pθ n q 1¤n¤N minimizes the quan-

tity °R r1 inf zZ N pv r ¡ zq T C N pv r ¡ zq among all the N -dimensional subspaces Z N of R N . Moreover, we have }v ¡ Π Z N v} X ¤ σ N 1 }v} X , for all v S.
In practice, when D R, we can avoid the computation of the matrix pC N q 1 2 and of its inverse by considering the matrix of smaller dimension T T T S T C N S R R¢R . Solving for the eigenvalues of T T T, we obtain the vectors ψn with associated eigenvalues σ 2 n since we have T T T ψn σ n T T ξ n σ 2 n ψn . Then, the vectors pθ n q 1¤n¤N are obtained as

θ n C N ¨¡ 1 2 ξ n 1 σn C N ¨¡ 1 2 T ψn 1
σn S ψn .

CHAPTER 3

A REDUCED-BASIS METHOD FOR PARAMETRIZED VARIATIONAL INEQUALITIES WITH NONLINEAR CONSTRAINTS

Abstract

We investigate new developments of the Reduced-Basis (RB) method for parametrized optimization problems with nonlinear constraints. In this chapter, we propose a reduced-basis scheme in a saddle-point form combined with the Empirical Interpolation Method to deal with the nonlinear constraint. In this setting, a 'primal' reduced-basis is needed for the primal solution and a 'dual' one is needed for Lagrange multipliers. We suggest to construct the latter using a 'cone-projected' hierarchical algorithm that conserves the non-negativity of the dual basis vectors. The reduction strategy is applied to elastic frictionless contact problems including the possibility of using non-matching meshes. We study test cases that are inspired from existing work on finite elements for contact mechanics. The numerical examples confirm the efficiency of the reduction strategy.

Introduction

Constrained optimization problems are of great importance in numerous engineering applications. Owing to the nonlinear nature of some constraints, the algorithms designed for solving nonlinearly-constrained optimization problems often suffer from slow convergence; thereby entailing subsequent computational costs. Besides, the Reduced-Basis (RB) method [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF] is a computationally effective approach to approximate the solutions of parametrized Partial Differential Equations (PDEs) encountered in many problems in science and engineering. In particular, it is highly beneficial in nonlinear settings that tend to substantially increase computational complexity. Although significant progress has been achieved in this field [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF][START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Haasdonk | A reduced basis method for parametrized variational inequalities[END_REF], the model reduction of parametrized optimization problems involving nonlinear constraints remains in need of further advances.

Here, the problem of interest is a parametrized optimization problem with nonlinear constraints which is formulated as a saddle-point problem and numerically solved using Lagrange multipliers. We consider a situation where this problem must be solved in a multi-query or real-time context, i.e., the problem has to be solved repeatedly for a large number of parameter values or it needs to be solved very quickly under limited computational resources. For standard PDEs in variational form, RB methods provide efficient tools for complexity reduction. More precisely, instead of the High-Fidelity (HF) problem, which is typically infinite-dimensional or rather high-dimensional after a finite element discretization, a low-dimensional model is generated. This low-dimensional problem can then be solved significantly faster for a wide range of parameters. Lately, the reduced modeling of variational inequalities has gained growing interest. In the literature, three recent papers address somewhat related problems. The first paper [START_REF] Haasdonk | A reduced basis method for parametrized variational inequalities[END_REF] extends the standard RB method to linear variational inequalities solved through a mixed formulation. Regarding the construction of the bases, the 'primal' basis (for the primal solution) and the 'dual' one (for the Lagrange multipliers) are directly composed of well-chosen snapshots. No additional compression phase is considered. In the so-called Projection-Based method of [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF], which has been specifically introduced to address time-dependent contact problems with linear constraints, the primal and the dual bases are built differently. An efficient primal RB is obtained using the POD. Since [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] focuses on dynamic problems with numerous time instants, the set of Lagrange multiplier snapshots can rapidly become sizable and its compression is therefore an important task. Therein, a dual basis is built by applying the Non-negative Matrix Factorization (NMF) algorithm [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] to the set of Lagrange multiplier snapshots. The NMF guarantees non-negative basis vectors and a limited RB dimension, but the resulting dual RB can be (far) less accurate than the primal RB. Another concern is that the user does not specify a required error tolerance as an input but a number of dominant basis vectors to retain. Finally, the work in [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF] extends another type of model-order reduction methods called Hyper-Reduction (HR) to contact problems with linear constraints. The proposed extension of the HR method consists in conserving a few vectors of the High-Fidelity (HF) dual basis because the number of contact nodes is limited to a Reduced Integration Domain (RID). Hence, only the contact nodes in the RID are treated but with a local high fidelity. So far, all the existing results are restricted to linear constraints.

In this paper, we propose to extend constrained model reduction to the framework of nonlinear constraints. Motivated by industrial applications in contact mechanics, we address a nonlinear type of constraints that can be written in a quasilinear form so as to use a fixed-point iterative scheme. We express the problem of interest in a saddle-point form, and we apply the Empirical Interpolation Method (EIM) [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Maday | A general multipurpose interpolation procedure: the magic points[END_REF] to allow for an offline/online decomposition of the nonlinear con-straints. Regarding basis construction, a 'primal' reduced basis is needed for the primal solution and a 'dual' one is needed for the Lagrange multipliers. Whereas the primal basis is constructed using a standard POD, we introduce a 'cone-projected' algorithm that builds nested dual bases while preserving the non-negativity of the basis vectors. The forthcoming analysis is meant to address static parametrized optimization problems, the extension to a time-dependent setting being straightforward once a time-discretization scheme has been chosen. An important application we have in mind is frictionless contact in a generic framework, namely without the small displacement assumption and for various types of contact constraints, including the intricate case of non-matching meshes. Both of the previous features lead to nonlinear operators in the definition of the constraint. On the one hand, the small displacement hypothesis allows one to consider the same normal vector on both contact boundaries. This work aims at going beyond this assumption. On the other hand, assuming that the meshes match eludes the nonlinearity induced by the spatial discretization of the constraint. Unfortunately, the latter assumption is not realistic in many engineering scenarios. Therefore, we also aim at addressing the case of non-matching meshes.

This chapter is organized as follows. In Section 3.2, we introduce the model problem. In Section 3.3, we consider elastic contact problems. Since we do not consider the simplifying hypotheses discussed above, we describe the way we derive the nonlinear non-interpenetration condition in some detail. In Section 3.4, we return to the general setting and we apply the RB method to derive a reduced resolution scheme. In Section 3.5, we discuss the offline stage in some detail, we present the EIM procedure for the nonlinear constraint, and we describe the construction of the primal and dual RB spaces. In Section 3.6, we present numerical results illustrating the performance of the method in the framework of elastic frictionless contact. We consider the contact problem between two spheres introduced by Hertz [START_REF] Hertz | Über die berührung fester elastischer körper[END_REF] with a parametrization of the radius of one of the spheres. Finally, Section 3.7 collects some technical results.

Model problem

Let V be a separable Hilbert space composed of functions defined on a spatial domain (open, bounded, connected subset) Ω R d , d ¥ 1, with a Lipschitz boundary fΩ.

Let Ω denote the closure of Ω and let P denote a parameter set. We define a continuous, symmetric and coercive bilinear form a : P ¢V ¢V Ñ R (the attributes of a are with respect to its second and third arguments), and a continuous linear form f : P ¢ V Ñ R (the attributes of f are with respect to its second argument). We also define a nonlinear continuous map k : P ¢ V ¢ V Ñ L 2 pΓ c 1 q and a nonlinear continuous map g : P¢V Ñ L 2 pΓ c 1 q, for a subset Γ c 1 fΩ. For simplicity, we consider at this stage that the domain Ω and the subset Γ c 1 are parameter-independent. A more general setting with parameter-dependent Ωpµq and Γ c 1 pµq will be considered from Section 3.3.2 onwards. For all µ P, we want to solve the following nonlinear minimization problem:

Find upµq V such that (3.1)

The semicolons in (3.1) indicate that a is linear with respect to its second and third arguments and that k is linear with respect to its third argument. Remark 3.1 (Nonlinear constraint). The nonlinear constraint in (3.1) can be formulated more compactly as ζpµ, upµqq ¤ 0 for the nonlinear continuous map ζpµ, vq : P ¢ V Ñ L 2 pΓ c 1 q defined as ζpµ, vq : kpµ, v; vq ¡ gpµ, vq.

(3.

2)

The

¢ V ¢ V Ñ L 2 pΓ c 1 q.
In the present setting, we make three assumptions. First, we assume that the inequality constraint in (3.1) is quasi-linear, i.e., that k is linear with respect to its third argument. This assumption will be exploited below in setting up an iterative solver for the discrete version of (3.1). Second, we assume that g satisfies gpµ, 0q ¥ 0.

Third, we assume that the problem (3.1) is well-posed. Note that the functional minimized in (3.1) is strongly convex and continuous. Moreover, the set of admissible states

K tv V | kpµ, v; vq ¤ gpµ, vqu (3.3)
is non-empty since 0 K because gpµ, 0q ¥ 0, and the set is closed owing to the continuity of k and g. Therefore, the existence of a minimizer is guaranteed. Our third assumption then means that we assume the uniqueness of the searched minimizer in K. In fact, the above setting is motivated by contact constrained problems that will be described in more detail in Section 3.3 below.

Let W be a non-empty closed convex cone composed of functions defined on the subset Γ c 1 fΩ. We assume that W : L 2 pΓ c 1 , R q, with R : r0, Vq. Using the test space W, the variational formulation of the inequality constraint in (3.1) reads as follows: Lpµqpv, ηq,

where the Lagrangian Lpµq : V ¢ W Ñ R is defined as

Lpµqpv, ηq : 1 2 apµ; v, vq ¡ f pµ; vq £ » Γ c 1 kpµ, v; vqη ¡ » Γ c 1 gpµ, vqη , (3.6)
and upµq and λpµq are respectively called the primal and the dual solutions of the saddle-point problem (3.5).

In practice, one uses a conforming Finite Element Method (FEM) [START_REF] Ern | Theory and practice of finite elements[END_REF] to discretize (3.5) in space. The FEM is based on a finite element subspace V N : spantφ 1 , . . . , φ N u V defined using a discrete nodal subset Ω tr Ω, where CardpΩ tr q N . Besides, one introduces the subset W R : span tψ 1 , . . . , ψ R u W defined using a discrete nodal subset Γ c,tr 1 Γ c 1 , where CardpΓ c,tr 1 q R. The notation span means that linear combinations are restricted to non-negative coefficients.

The discrete sadlle-point problem reads:

Find pu N pµq, λ R pµqq V N ¢ W R such that pu N pµq, λ R pµqq arg minmax vV N ,ηW R Lpµqpv, ηq, (3.7)
with the Lagrangian defined in (3.6). Note that the discrete inequality constraint reads

» Γ c 1 kpµ, u N pµq; u N pµqqψ i ¤ » Γ c 1 gpµ, u N pµqqψ i , di t1, . . . , Ru. (3.8)
As is customary with the RB method, we assume henceforth that the mesh-size is small enough so that the above space discretization method delivers HF approximate primal and dual solutions within the desired level of accuracy. Introducing the component vectors upµq : pu n pµqq 1¤n¤N R N and λpµq : pλ n pµqq 1¤n¤R R R of u N pµq and λ N pµq, the algebraic formulation of (3.7) reads: Find pupµq, λpµqq R N ¢ R R satisfying pupµq, λpµqq arg minmax

vR N ,ηR R 1 2 v T Apµqv ¡ v T f pµq η T Kpµ, vqv ¡ gpµ, vq ¨, (3.9)
with the matrices Apµq R N ¢N and Kpµ, wq R R¢N such that Apµq ij apµ; φ j , φ i q, Kpµ,

wq ij » Γ c 1 kpµ, w; φ j qψ i , (3.10)
and the vectors f pµq R N and gpµ, wq R R such that f pµq j f pµ; φ j q, gpµ, wq j

» Γ c 1 gpµ, wqψ j . (3.11)
In the sequel, we will solve (3.9) using an iterative algorithm, where the terms Kpµ, vq and gpµ, vq are treated explicitly. This amounts to a so-called 'secant method' or Kačanov iterative method [START_REF] Fučík | Kačanov-Galerkin method and its application[END_REF]. The Kačanov iteration consists in solving the following problem : For all k ¥ 1, pu k pµq, λ k pµqq arg minmax

vR N ,ηR R 1 2 v T Apµqv ¡ v T f pµq η T Kpµ, u k¡1 pµqqv ¡ gpµ, u k¡1 pµqq ¨. (3.12)
For a user-defined tolerane ka ¡ 0, the stopping criterion of the Kačanov iterative algorithm reads

}u k pµq ¡ u k¡1 pµq} R N }u k¡1 pµq} R N ¤ ka . (3.13)
Depending on the problem and output of interest, an additional check on the dual increment }λ k pµq ¡ λ k¡1 pµq} R R {}λ k¡1 pµq} R R can be performed. A brief comparison with the Newton method is presented in Section 3.7.3. The advantage of the Kačanov iterative method is its simplicity. Indeed, unlike the standard Newton method, Kačanov iterations do not require any computation of Jacobian preconditioners, thereby achieving significant computational savings when solving (3.9). On the other hand, if the Newton method converges, it is (much) faster than the Kačanov iteration. In Section 3.4 below, the reduced problem will be solved using the Kačanov iteration as well. Therein, we will shortly discuss the influence of the solver on the reduction scheme.

Prototypical example: elastic contact

The model reduction of mechanical problems involving contact remains an important issue in computational solid mechanics. In this section, we consider the case of frictionless contact for linear elasticity. An important difference with regard to the previous section is that the domain over which the problem is posed is now parameter-dependent.

Linear elasticity

For all µ P, the domain Ωpµq R d , d t2, 3u represents the initial configuration of a deformable medium initially at equilibrium and to which an external load pµq : Ωpµq Ñ R d is applied. The standard linear elasticity problem consists in finding the displacement field upµq : Ωpµq Ñ R d induced by the externally applied force field pµq once the system has reached equilibrium. Let σpupµqq : Ωpµq Ñ R d¢d be the stress tensor in the medium. The equilibrium conditions can be expressed as ∇ ¤ σpupµqq pµq, in Ωpµq.

(3.14)

We define the functional space Vpµq such that Vpµq : H 1 pΩpµq; R d q.

(3.15)

For all v Vpµq, let ε : Ωpµq Ñ R d¢d be the linearized strain tensor defined as εpvq : 1

2 p∇v ∇v T q.

(3.16)

In the framework of linear isotropic elasticity, the stress tensor is related to the linearized strain tensor by the formula

Non-interpenetration condition

We intend to model the non-interpenetration condition in a general framework. We use the previously introduced notation and only define the new quantities dedicated to this particular formulation. We consider that the domain Ωpµq can be partitionned as

Ωpµq Ω 1 pµq Ω 2 pµq, where Ω 1 pµq and Ω 2 pµq represent the initial configuration of the two disjoint deformable media. For all µ P, let Γ c 1 pµq and Γ c 2 pµq be the potential contact boundaries of Ω 1 pµq and Ω 2 pµq, such that Γ c pµq : Γ c 1 pµq Γ c 2 pµq. For all v Vpµq and all i t1, 2u, we introduce the following functions v i : Ω i pµq Ñ R d such that

v i : v |Ω i pµq . (3.20)
In order to formulate the contact conditions in a general setting, we need to introduce some auxiliary geometric mappings. An illustration of the various geometric mappings is given in Figure 3.1. For all v V, all µ P and all i t1, 2u, we define the geometric mappings

ψ i pµ, v i q : Γ c i pµq Ñ Υ c i pµ, v i q z Þ Ñ z v i pzq, (3.21) where Υ c i pµ, v i q : ψ i pµ, v i q Γ c i pµq ¨.
In what follows, we assume implicitly that v i is injective so that ψ i pµ, v i q is injective as well. Therefore, pψ i pµ,

v i qq ¡1 : Υ c i pµ, v i q Ñ Γ c
i pµq is well defined. This assumption is natural in the context of solid mechanics.

Under a local convexity assumption on Υ c 2 pµ, v 2 q, the contact mapping

ϑpµ, vq : Υ c 1 pµ, v 1 q Ñ Υ c 2 pµ, v 2 q z 1 Þ Ñ argmin z 2 Υ c 2 pµ,v 2 q }z 1 ¡ z 2 }, (3.22)
is well defined. The contact mapping ϑpµ, vq can be physically interpreted as the function relating every point on Υ c To be more precise, for all z Υ c 1 pµ, v 1 q, ϑpµ, vqpzq is the orthogonal projection of

z onto Υ c
2 pµ, v 2 q. The contact mapping ϑpµ, vq depends on the displacement field v and is therefore unknown a priori for the solution v upµq. For all v V and all µ P, we define the mapping

5 ρpµ, vq : Γ c 1 pµq Ñ Γ c 2 pµq ρpµ, vq : pψ 2 pµ, v 2 qq ¡1 ¥ ϑpµ, vq ¥ ψ 1 pµ, v 1 q, (3.23)
and the vector field of outward normals on Υ c

2 pµ, v 2 q n 2 pµ, v 2 q : Υ c 2 pµ, v 2 q Ñ R d . (3.24)
It is also convenient to introduce the vector field

r n 2 pµ, vq : Γ c 1 pµq Ñ R d r n 2 pµ, vq : n 2 pµ, v 2 q ¥ ϑpµ, vq ¥ ψ 1 pµ, v 1 q, (3.25)
which corresponds to the outward normal on Υ c 2 pµ, v 2 q but defined at the corresponding point in Γ c

1 pµq through the mapping ϑpµ, vq ¥ ψ 1 pµ, v 1 q. For an admissible solution upµq pu 1 pµq, u 2 pµqq Vpµq, the noninterpenetration condition reads: For all z Γ c 1 pµq, pu 1 pµqpzq ¡ pu 2 pµq ¥ ρpµ, upµqqq pzqq ¤r n 2 pµ, upµqqpzq ¥ ρpµ, upµqqpzq ¡ z ¨¤ r n 2 pµ, upµqqpzq.

(

The reduced-basis model

In this section, we return to the general setting of Section 3.2 and we derive a general RB formulation for the nonlinear minimization problem (3.1), and more precisely its algebraic saddle-point formulation (3.9). Yet, following Section 3.3, we now consider the more general setting of a parameter-dependent domain Ωpµq.

Reference domain

In view of model reduction, we assume that there exists a bi-Lipschitz diffeomorphism called geometric mapping hpµq and defined on a parameter-independent reference domain q Ω such that

hpµq : q Ω Ñ Ωpµq x Þ Ñ I i1 h i pµ, xq1 q Ω i pxq, (3.33)
where t q Ω i u I i1 is a partition of q Ω. Using this geometric mapping, we introduce the reference Hilbert space q V : H 1 p q Ω; R d q composed of functions defined on q Ω so that

Vpµq q V hpµq ¡1 tv ¥ hpµq ¡1 | v q Vu. (3.34)
In what follows, for all i t1, . . . , Iu, we set Ω i pµq hpµqp q Ω i q,

(3.35) so that q Ω i is parameter-independent.

We assume for simplicity that I 2 in (3.33), which corresponds to the situation where there are two disjoint solids Ω 1 pµq and Ω 2 pµq that can come into contact. We then define the reference contact boundaries such that Γ c i pµq hpµqp q Γ c i q, di t1, 2u. (3.36) We also define Wpµq : L 2 pΓ c 1 pµq; R q and | W :

L 2 p q Γ c 1 ; R q such that Wpµq | W h ¡1 pµq |Γ c 1 pµq tq η ¥ h ¡1 pµq |Γ c 1 pµq | q η |
Wu.

(3.37)

Reduced basis spaces

Let V N pµq and W R pµq be the FEM discretizations of Vpµq and Wpµq, respectively.

The reference counterparts of these spaces are q V N and | W R . Towards an accurate approximation of the solution manifold, we consider a so-called 'primal' RB subspace p V N pµq and a so-called 'dual' RB subcone

x W R pµq such that p V N pµq V N pµq Vpµq, and x W R pµq W R pµq Wpµq. (3.38)
The dimensions of these spaces are such that N dimp p V N pµqq 3 dimpV N pµqq N and R dimp x W R pµqq 3 dimpW R pµqq R. The reference counterparts of the parameter-dependent RB spaces in (3.38) are spaces q V N and | W R that satisfy

q V N q V N q
V, and

| W R | W R | W. (3.39)
Let p q θ n q 1¤n¤N be an orthonormal basis of q V N and let p q ξ n q 1¤n¤R be generating vectors of the cone | W R , i.e., | W R span t q ξ 1 , . . . , q ξ R u t °R n1 ω n q ξ n | ω n ¥ 0u. For all µ P, the primal RB solution p upµq p V N pµq and the dual RB solution (Lagrange multipliers) p

λpµq x W R pµq that approximate the HF solution pu N pµq, λ R pµqq V N pµq ¢ W R pµq are decomposed as p upµq N ņ1 p u n pµq q θ n ¥ hpµq ¡1 , p λpµq R ņ1 p λ n pµq q ξ n ¥ hpµq ¡1 , ( 3
θ n ¥ hpµq ¡1 , q θ p ¥ hpµq ¡1 q, (3.42a) p Kpµ, p vq pn » Γ c 1 pµq k ¢ µ, N i1 p v i q θ i ¥ hpµq ¡1 ; q θ n ¥ hpµq ¡1 q ξ p ¥ hpµq ¡1 , (3.42b)
and the vectors p f pµq R N and p gpµ, p vq R R such that p f pµq p f pµ; q θ p ¥ hpµq ¡1 q,

(3.43a) p gpµ, p vq p » Γ c 1 pµq g £ µ, N i1 p v i q θ i ¥ hpµq ¡1 q ξ p ¥ hpµq ¡1 .
(3.43b)

Separation of the elastic energy

We assume the existence of two integers J a and J f and of matrices p A j R N ¢N , with 1 ¤ j ¤ J a , and vectors p f j R N , with 1 ¤ j ¤ J f , such that p A j,pn a j p q θ n , q θ p q, and p f j,pn f j p q θ p q, (3.44) defined using the reference basis functions such that the matrix p Apµq defined in (3.42a) and the vector p f pµq defined in (3.43a) can be affinely decomposed under the form

¡ p Apµq © np J a j1 α a j pµq p
A j,np , and

¡ p f pµq © p J f j1 α f j pµq p f j,p , d1 ¤ n, p ¤ N, (3.45)
where the dependencies on µ and n, p are separated. The key point is that the quantities in (3.44) are offline-computable. During the online stage, all that remains to be performed is the assembly of p Apµq and p f pµq using (3.45) for each new parameter value µ P.

In order to clarify how the separated representations in (3.45) are derived, we consider the elastic problem defined in Section 3.3 with homogenous load function pµq and a test function w. We have

f pµ; wq » Ωpµq pµqwpxq dx » q Ω pµqwphpµqpq xqqdet pJacphpµqqpq xqq dq x, (3.46)
where the notation det pJacphpµqqq refers to the determinant of the Jacobian matrix of the geometric mapping hpµq. Using the definition of hpµq in (3.33), we obtain

f pµ; wq 2 j1 » q Ω j h j pµq pµq q wpq xqdet Jacp1 q Ω j qpxq ¨dq x 2 j1
h j pµq pµq » q Ω j q wpq xq dq x.

(3.47)

Consequently, J f 2, and in (3.45), we have α f j pµq h j pµq pµq and p f j,p ³ q Ω j q θ p pq xq dq x, for all j t1, 2u.

Separation of the constraint

The remaining bottleneck is the computation of p Kpµ, p vpµqq and p gpµ, p vpµqq in (3.42b) and (3.43b) which requires parameter-dependent reconstructions using the FEM basis functions in order to compute the integrals over Γ c 1 pµq.

The key idea is to introduce a further approximation. Namely, we search for approximations κ M k and γ M g of the nonlinear maps κ : P ¢ t1, . . . , N u ¢ q Γ c 1 Ñ R and γ : P ¢ q Γ c 1 Ñ R defined such that κpµ, n, q xq : kpµ, upµq; φ n qphpµqpq xqq, and γpµ, q xq : gpµ, upµqqphpµqpq xqq. (3.48) Our goal in building the approximations κ M k and γ M g is to separate the dependence on µ from the dependence on the other variables. More precisely, for some integers

M k , M g ¥ 1, we look for (accurate) approximations κ M k : P ¢ t1, . . . , N u ¢ q Γ c 1 Ñ R of κ and γ M g : P ¢ q Γ c 1 Ñ R of γ in the separated form κ M k pµ, n, q xq : M k j1 ϕ κ j pµqq κ j pn, q xq, γ M g pµ, q
xq :

M g j1 ϕ γ j pµqq γ j pq xq, (3.49)
where M k (resp. M g ) is called the rank of the approximation and ϕ κ j pµq (resp. ϕ γ j pµq) are real numbers that are found by interpolation. For κ M k , we interpolate over a set of M k pairs tpn κ

1 , q x κ 1 q, . . . , pn κ M k , q x κ M k qu in t1, . . . , N u ¢ q Γ c 1 , whereas for γ M g , we interpolate over a set of M g points tq x γ 1 , . . . , q x γ M g u in q Γ c 1 .
The interpolation is performed using the EIM [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] and leads to the quantities p κpµ,

p vq R M k , B κ R M k ¢M k , p γpµ, p
vq R M g and B γ R M g ¢M g defined as follows:

p vR N , p ηR R 3 1 2 p v T p Apµqp v ¡ p v T p f pµq p η T D κ pµ, p vqp v ¡ D γ pp vqp γpµ, p vq ¨A, (3.51)
with the matrices

D κ pµ, p vq M k j1 C κ j ppB κ q ¡1 p κpµ; p vqq j , and D γ C γ pB γ q ¡1 , (3.52)
where C κ j R R¢N and C γ R R¢M g are given by

C κ j,pn ¢ N i1 » q Γ c 1 q
θ n,i q κ j pi, ¤q q ξ p pn , and

C γ pj ¢ » q Γ c 1 q γ j q ξ p pj . (3.53)
The overall computational procedure can now be split into two stages:

(i) An offline stage where one precomputes on the one hand the RB subspace p V N and the RB subcone x W R leading to the vectors t p f r u 1¤j¤J a in R N and the matrices t p A r u 1¤j¤J f in R N ¢N , and on the other hand the EIM pairs tpn κ i , q

x κ i qu 1¤i¤M k , the EIM points tq x γ i u 1¤i¤M g , the EIM functions tq κ j u 1¤j¤M k , and the EIM functions tq γ j u 1¤j¤M g , leading to the matrices

B κ R M k ¢M k , B γ R M g ¢M g , tC κ j u 1¤j¤M k in R R¢N and C γ R R¢M g .
The offline stage is discussed in more detail in Section 3.5.

(ii) An online stage to be performed each time one wishes to compute a new solution for a parameter µ P. All that remains to be performed is to assemble the vector p f pµq R N and the matrix p Apµq R N ¢N using (3.45), to compute the vectors p κpµ, p vq R M k and p γpµ, p vq R M g defined in (3.50), to assemble the matrix D κ pµ, p vq defined in (3.52), and to solve the reduced saddle-point problem (3.51). The online stage is summarized in Algorithm 3.1. upµq and p λpµq Remark 3.4 (EIM matrix). The computations in Algorithm 3.1 only require the knowledge of the matrix pB κ q ¡1 . In order to optimize the computational costs, B κ is inverted prior to the algorithm, i.e. during the offline stage. The matrix B γ is also inverted when computing the matrix D γ (see (3.52)), which is an input of the offline stage.

Algorithm 3.1 Online stage

Input : µ, t p f j u 1¤j¤J f , t p A j u 1¤j¤J a , tpn κ i , q x κ i qu 1¤i¤M k , tq x γ i u 1¤i¤M g , tq κ j u 1¤j¤M k , tq γ j u 1¤j¤M g ,
Remark 3.5 (EIMs on k and g). Owing to the quasi-linear structure of the inequality constraint, the reduced problem (3.51) is solved using the Kačanov algorithm. At first glance, the influence of this resolution choice is that we have to perform the EIM twice since the maps k and g are separated one at a time. Were we to use a standard Newton method by considering the one-term constraint ζpµ, upµqq ¤ 0 (see (3.2)), we would only perform a single EIM. However, an additional EIM would be needed in the Newton method in order to compute the Jacobian preconditioning matrix. Thus, both options (Kačanov or Newton) lead to two distinct EIMs and the storage cost is essentially the same.

The offline stage

There are two tasks to be performed during the offline stage:

(T 1 ) Build the rank-M k and the rank-M g EIM approximations in (3.49);

(T 2 ) Explore the solution manifold in order to construct the linear subspace q

V N q V N of dimension N and the subcone | W R | W R of dimension R.
Tasks (T 1 ) and (T 2 ) can be performed independently and in whatever order. Since Task (T 1 ) can be considered to be standard, we only discuss Task (T 2 ), i.e., the construction of the sets of primal and dual RB functions with cardinalities N and R respectively. First, as usual in RB methods, the solution manifold is explored by considering a training set for the parameter values. For simplicity, one can consider the same training set P tr as for the EIM approximations. This way, one only explores the collection of points S pri tupµqu µP tr and S du tλpµqu µP tr respectively in the primal and dual solution manifolds. For this exploration to be informative, the training set P tr has to be chosen large enough. In the present setting where HF solutions are to be computed for all the parameters in P tr when constructing the EIM approximations, it is natural to compress these computations by means of a Proper Orthogonal Decomposition (POD) [START_REF] Hinze | Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: error estimates and suboptimal control[END_REF][START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF] to define the primal RB subspace q V N . This technique is often considered in the literature to build the RB, see, e.g., [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations[END_REF]. Bearing in mind that the dual RB cone | W R is meant to represent the set of Lagrange multipliers, its spanning vectors should all have nonnegative components. Consequently, the POD is not appropriate to build the dual RB cone. If the training set has a reasonable size, one could keep all the snapshots, especially if they have been computed via a posteriori error estimation. In [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF], it is suggested to use the Non-negative Matrix Factorization (NMF) algorithm [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] whenever the number of training snapshots is relatively large, for instance in the case of a time-dependent problem. For a set of positive snapshots S du and an integer R, the procedure NMFpS du , Rq returns R vectors pw 1 , . . . , w R q with nonnegative components (cf. Section 3.7.2). Nonetheless, the resulting dual RB can be less accurate than the primal RB and the user does not specify a required error tolerance as an input but only the cardinality of the dual RB. In practice, it is often difficult to anticipate the approximation capacity of the dual RB from its cardinality.

Here, we suggest to build a dual hierarchical basis from the Lagrange multiplier snapshots computed offline. In the spirit of weak greedy algorithms, the idea is to order the snapshots depending on their relevance to represent the entire set of snapshots. For all µ P tr , we define the reference Lagrange multiplier snapshot q λpµ; ¤q : λpµ; ¤q ¥ hpµq : q Γ c 1 Ñ R.

(3.54) First, we choose µ 1 P tr such that

µ 1 argmax µP tr } q λpµ; ¤q} V p q Γ c,tr 1 q , (3.55)
where the discrete subset q Γ c,tr

1 q
Γ c 1 is introduced to compute the maximizer in space.

Afterwards, at each iteration n ¥ 2, we define the convex cone q K n span tλpµ 1 , ¤q¥ hpµq, . . . , λpµ n¡1 , ¤q ¥ hpµqu and select a new parameter value µ n P tr using the criterion µ n argmax µP tr } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr 1 q , (3.56)

where Π q K n¡1 is the orthogonal projector onto the convex cone q K n¡1 . At each iteration, we check whether or not the stopping criterion max µP tr } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr

1 q ¤ du , (3.57)
is fulfilled. The steps of the 'cone-projected' weak greedy algorithm are summarized in Algorithm 3.2.

Algorithm 3.2 Cone-projected weak greedy algorithm

Input : P tr , q Γ c,tr

1
and du ¡ 0 1: Compute S du t q λpµ; ¤qu µP tr # HF solutions

2: Set q K 0 t0u 3: Set n 1 4: Set r 1 max µP tr } q λpµ; ¤q} V p q Γ c,tr 1 q 5: while (r n ¡ du ) do 6:
Search µ n argmax

µP tr } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr 1 q 7:
Set q K n : span t q λpµ 1 ; ¤q, . . . , q λpµ n ; ¤qu 8:

Set n n 1

9:

Set r n : max µP tr } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr 1 q 10: end while

11: Set R : n ¡ 1 Output : | W R : q K R .
Remark 3.6 (Cone projections). The projection onto the cone q K n¡1 in line 6 of Algorithm 3.2 is not trivial. We use the off-the-shelf solver from [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF]. Remark 3.7 (Relative error). Algorithm 3.2 can be run using a relative error criterion instead of an absolute one. Towards this end, one replaces the absolute error } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr 1 q in lines 6 and 9 by the relative error } q λpµ; ¤q ¡ Π q K n¡1 p q λpµ; ¤qq} V p q Γ c,tr

1 q {} q λpµ; ¤q} V p q Γ c,tr 1 q .
(3.58) Remark 3.8 (Elementary compression). For the reduced problem, if we choose to conserve all the snapshots of the Lagrange multipliers, we can still check for some computational savings. This can be achieved by suppressing the constraints that are never saturated for any of the parameters in the training set P tr but were initially introduced in the HF model. In practice, we reduce the dimensions of the matrix Kpµ, upµqq and the vector gpµ, upµqq by removing the lines and columns of Kpµ, upµqq and the components of gpµ, upµqq that always vanish no matter the value of the parameter µ P tr .

Numerical results

In this section, we illustrate the above developments by a numerical example related to elastic contact in a two-dimensional framework. The investigated test case is the seminal contact problem between two spheres introduced by Hertz [START_REF] Hertz | Über die berührung fester elastischer körper[END_REF]. The HF computations use the software Freefem++ [START_REF] Hecht | New developments in freefem++[END_REF] in order to generate the mesh and the matrices and vectors of the saddle-point problem. The HF solution of the saddlepoint problem is then retrieved using Python and its convex optimization package cvxopt [START_REF] Andersen | Cvxopt: A python package for convex optimization[END_REF]. The RB algorithms are developed in Python. Our goal is to illustrate the computational performance of the method. We choose a Galerkin discretization for the primal solution and a collocation method for the contact nodes. We introduce a reference mesh grid q Ω tr q

Ω with Cardp q Ω tr q N and the corresponding finite element test functions q φ j,x p q φ j , 0q and q φ j,y p0, q φ j q for d 2 (the extension being straightforward for d 3). Also, we introduce a set of ordered sampling points q Γ c,tr

1 : tq z 1 , ¤ ¤ ¤ , q z R u q Γ c 1 q
Ω tr with Cardp q Γ c,tr 1 q R. By means of collocation over this set of points, the expressions of the matrices C κ j and C γ in (3.52) and (3.53) are recast as

6 9 9 8 9 9 7 C κ j x W R Q κ j , where C κ j ¢ N i1 » Γ c 1 θ n,i q κ j pi, q z p q 1¤p¤R,1¤n¤N R R¢N , C γ x W R Q γ , where Q γ q γ j pq z p q ¨1¤p¤R,1¤j¤M g R R¢N , (3.59) where x
W R R ¢ R is the matrix of dual snapshots. Consider the parameter set P r0.9, 1.12s, the discrete training set P tr t0.905 0.01i| 0 ¤ i ¤ 22u and a two-dimensional setting based on two half-disks with an upper constant radius R up 1m, a lower parametric radius R low µ and an initial gap between the centers of the disks γ 0 0.1m. The left panel of Figure 3.2 shows the reference domain q Ω Ωpµ 1q. HF solutions are computed using a finite element subspace defined on a mesh of Ω tr with 675 nodes consisting of continuous, piecewise affine functions and the potential contact zone contains 51 nodes on each disk. Consequently, our problem has N 1350 degrees of freedom and R 51 Lagrange multipliers. The materials of both hald-disks are identical and the HF computations are run with a Young modulus E 15Pa and a Poisson coefficient ν 0.35. Regarding boundary conditions, we consider a homogeneous Dirichlet condition u x 0 and u y 0 on the lower horizontal edge, a homogeneous Dirichlet condition u x 0 on the upper horizontal edge, and an imposed displacement u y ¡0.4m on the upper horizontal edge. The right panel of Figure 3.2 displays the HF displacement field for the parameter value µ 1.12m, whereas Figure 3.3 displays the normal contact stress as a function of the reference configuration's abscissas for the parameter values µ P tr . The normal contact is zero on the nodes where the contact between the two half-disks is not established at equilibrium.

During the offline stage, we perform P 23 HF computations. Applying the POD to S pri based on the energy norm and an absolute truncation threshold pod 10 ¡3 , the primal space q V N is composed of 11 RB functions. Table 3.1 shows the size of the reduced basis as a function of the tolerance pod . As can be seen, the number of offline computations is equal to the dimension of the primal space for pod 4.10 ¡6 . The left panel of Figure 3.4 illustrates the decrease of the singular values associated with the POD modes. The decrease is not as sharp as is often the case for variational equalities. Moreover, the higher the rank of the singular value, the milder the decrease of the error. In order to build the dual reduced basis, we test both the NMF suggested in [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] and the cone-projected algorithm (see Algorithm 3.2). Table 3.2 shows the dimension R of the dual space as a function of the truncation threshold du for both algorithms. The cone-projected greedy algorithm achieves the same accuracies with less basis function than the NMF. Notice that, at the first iterations of the procedure, the NMF uses at least twice as much functions as the cone-projected greedy algorithm. The right panel of Figure 3.4 displays the approximation error for the dual | W R space as its dimension R increases. This figure clearly shows that the cone-projected greedy algorithm outperforms the NMF. Note that it is actually pointless to perform an NMF with an input number of modes equal to the total number of modes. The same reasoning applies to the cone-projected greedy algorithm as well. However, as the cone-projected greedy algorithm is (in principle) steered by an accuracy threshold rather than a number of modes, the case in which all the modes are retained can still be justified. Hence, the right panel of figure 3.4 is meant to compare the performances of both algorithms at convergence. We now perform the EIM twice so as to allow for an offline/online decomposition V N . Right: Approximation error for the dual space | W R .

of both terms in the inequality constraint. The convergence of the EIM is reported in Figure 3.5. The approximation error decreases clearly faster for the nonlinear gap function gpµ, upµqq than for the nonlinear contact operator Kpµ, upµqq. This observation is not counter-intuitive since the contact map is a trivariate function, whereas the gap map is a bivariate function. Therefore, the EIM on Kpµ; uq needs a higher rank M k to be accurate.

Let us now investigate the online stage for the entries du 10 ¡4 , k eim 10 ¡2 , for the cone-projected greedy algorithm is always below that of the NMF. Moreover, the right panel of Figure 3.6 shows the relative H 1 -error error on the displacement field defined as follows:

e displ pµq : }p upµq ¡ u N pµq} H 1 pΩ tr q }u N pµq} H 1 pΩ tr q . (3.61)
The conlusion is similar as for e ener pµq. for larger parameter values, i.e. µ t1.01, . . . , 1.12u. The reason for this is that the spatial discretization becomes coarser with the increase of the parameter value. Notice that the interpenetration curves for the NMF and the cone-projected greedy almost overlap.

As a second test, we run both methods with the looser tolerances du 1 ¤ 10 ¡2 , k eim 2 ¤ 10 ¡1 , g eim 10 ¡2 and pod 10 ¡2 . Figure 3.8 displays the errors e ener pµq and e displ pµq using only R 4 basis vectors for the cone-projected greedy and R 10 basis vectors for the NMF. The overall error is larger than in the previous cases since all the tolerances have been loosened. In order to get a clearer insight on the impact of the dual space | W R , we display the same plots as in the first simulation with the larger truncation threshold du 5 ¤ 10 ¡3 . The other tolerances are k eim 10 ¡2 , g eim 10 ¡4 and pod 10 ¡4 . In this configuration, the NMF conserves R 18 dual basis vectors, whereas the coneprojected greedy only conserves R 5 dual basis vectors. The minimum energy Figure 3.9 -Error quantification with R 18 for the NMF and R 5 for the cone-projected greedy. Left: Error on the minimum energy e ener pµq. Right: Relative H 1 -error for the displacement field e displ pµq.

error e ener pµq and the relative H 1 -error e displ pµq are plotted in Figure 3.9. In spite of the substantial difference between the sizes of the NMF space and the coneprojected greedy space, Figure 3.9 shows that the cone-projected greedy algorithm still delivers accurate approximations and, in average, produces smaller errors. For the comparison between the two algorithms to be fairer, we keep R 5 basis vectors for the NMF and display the error indicators e ener pµq and e displ pµq in Figure 3.10.

In this situation, the error for the cone-projected greedy algorithm is always below

Technical complements

This section gathers three technical complements that give further insight on some of the aspects discussed within this chapter.

Kačanov vs. Newton

Suppose that we need to solve numerically the nonlinear problem: Find u R N such that Fpuq 0, where F : R N Ñ R N . Given an initialization u 0 R N , the standard Newton iteration reads

DFpu k qpu k 1 ¡ u k q ¡Fpu k q, dk ¥ 0. (3.63)
If the nonlinear operator F is of the form Fpvq Apvqv ¡ gpvq, with A : R N Ñ R N ¢N and g : R N Ñ R N , the iteration (3.63) becomes

DApu k q ¤ u k Apu k q ¡ Dgpu k q ¨pu k 1 ¡ u k q ¡Apu k qu k gpu k q, (3.64) leading to DApu k q ¤ u k ¡ Dgpu k q ¨pu k 1 ¡ u k q Apu k qu k 1 gpu k q.
(3.65)

Provided that DApvq ¤ v ¡ Dgpvq 0, (3.65) leads to the Kačanov iteration Apu k qu k 1 gpu k q.

(3.66)

Non-negative Matrix Factorization (NMF)

Let us give some details on the NMF for completeness. The following results can be found in the standard literature on clustering [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. The goal is to briefly describe the procedure associated with the notation pw 1 , . . . , w R q NMFpT, Rq,

which is mentioned in Section 3.5, where we are given P vectors pt 1 , ¤ ¤ ¤ , t P q forming the rectangular matrix T R R¢P whose entries are all non-negative, and we are looking for R positive vectors forming the rectangular matrix W : pw 1 , . . . , w R q R R¢R . We define the error function that quantifies the quality of the approximation of the matrix A by a matrix B as the Frobenius norm of the difference between A and B:

}A ¡ B} 2 : ij pA ij ¡ B ij q 2 , (3.68)
which clearly vanishes if and only if A B. The NMF optimization problem then reads:

pW, Hq argmin

WR R¢R HR R¢P }T ¡ W H} 2 , (3.69)
where W R R¢R and H R R¢P . The function }T ¡ WH} is not convex in both variables W and H together. Thus, only the recovery of local minima is feasible.

Regarding the search algorithm, the decrease of the Frobenius norm is proven in [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] for the following iterative update rules

H ij Ð H ij pW T Tq ij pW T WHq ij , W ij Ð W ij pTH T q ij pWHH T q ij . (3.70)
The first update rule is equivalent to a gradient descent algorithm

H ij Ð H ij η ij pW T Tq ij ¡ pW T WHq ij ¨, (3.71)
with η ij H ij {pWHH T q ij . The reasoning for W ij is similar. The motivation for (3.70) is that if the pair pW, Hq yields an exact reconstruction, i.e. T WH, then pW, Hq is a fixed-point of the algorithm. Finally, note that, in contrast to the POD, the integer R is a required input for the NMF. Moreover, the uniqueness of the NMF is not guaranteed. In fact, any positive matrix D R R¢R satisfies T WDD ¡1 H, thereby leading to another NMF decomposition.

Proof of Lemma 3.1

Proof. We first prove the converse implication, i.e. (3.30) ùñ (3.26). It is readily verified that

(3.30) ùñ Υ c 1 pµ, u 1 pµqq e Ωpµ,upµqq Ω 2 pµ, u 2 pµqq Υ c 2 pµ, u 2 pµqq, (3.72)
where the union in (3.72) is disjoint. For all z Γ c 1 pµq, we have ψ 1 pµ, u 1 pµqqpzq z u 1 pµqpzq Υ c 1 pµ, u 1 pµqq. Thus, the inclusion (3.72) implies that, for all z Γ c 1 pµq :

ψ 1 pµ, u 1 pµqqpzq Ω 2 pµ, u 2 pµqq or ψ 1 pµ, u 1 pµqqpzq Υ c 2 pµ, u 2 pµqq. (3.73)
By definition of ϑpµ, upµqq, we have

ψ 1 pµ, u 1 pµqqpzq Ω 2 pµ, u 2 pµqq ðñ pψ 1 pµ, u 1 pµqqpzq ¡ pϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqqq pzq ¤ r n 2 pµ, upµqqpzqq ¡ 0,
and

ψ 1 pµ, u 1 pµqqpzq Υ c 2 pµ, u 2 pµqq ðñ pψ 1 pµ, u 1 pµqq ¡ ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqqq pzq 0 ðñ ppψ 1 pµ, u 1 pµqq ¡ ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqqq pzq ¤ r n 2 pµ, upµqqpzqq 0.
The latter equivalence follows from the colinearity of the two vectors involved in the scalar product. Thus, (3.30) ùñ ppψ 1 pµ, u 1 pµqq ¡ ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqqq pzq ¤ r n 2 pµ, upµqqpzqq ¥ 0.

(3.74)

Moreover, by definition of ψ 2 pµ, u 2 pµqq, we have

ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqq pψ 2 pµ, u 2 pµqqq ¡1 ¥ ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqq u 2 pµq ¥ pψ 2 pµ, u 2 pµqqq ¡1 ¥ ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqq, leading to ϑpµ, upµqq ¥ ψ 1 pµ, u 1 pµqq ρpµ, upµqq u 2 pµq ¥ ρpµ, upµqq.
A straightforward replacement in (3.74) yields that, for all z Γ c 1 pµq, pψ 1 pµ, u 1 pµqq ¡ u 2 pµq ¥ ρpµ, upµqq ¡ ρpµ, upµqqp zq ¤ r n 2 pµ, upµqqpzq ¥ 0. (

Using the definition of ψ 1 pµ, u 1 pµqq, (3.75) yields, for all z Γ c 1 pµq, pu 1 pµq ¡ u 2 pµq ¥ ρpµ, upµqqq pzq ¤ n 2 pµ, upµqqpzq ¥ pρpµ, upµqqpzq ¡ zq ¤ r n 2 pµ, upµqqpzq, thereby showing that (3.30) ùñ (3.26).

Let us now prove the direct implication, i.e. (3.26) ùñ (3.30). Let z Ω 1 pµ, u 1 pµqq Ω 2 pµ, u 2 pµqq. Since in particular, z Ω 1 pµ, u 1 pµqq, we consider two cases:

• z is a boundary point of the first solid, i.e. z Υ c 1 pµ, u 1 pµqq; • z is an interior point of the first solid, i.e. z Ω 1 pµ, u 1 pµqq.

If z is a boundary point, we have z Υ c 1 pµ, u 1 pµqq ùñ z ¡ ϑpµ, upµqqpzq ¤ r n 2 pµ, upµqq ¥ pψ 1 pµ, u 1 pµqqq ¡1 pzq ¨¥ 0 ðñ z ϑpµ, upµqqpzq or z ¡ ϑpµ, upµqqpzq ¤ r n 2 pµ, upµqq ¥ pψ 1 pµ, u 1 pµqqq ¡1 pzq ¨¡ 0 ùñ z Υ c 1 pµ, u 1 pµqq Υ c 2 pµ, u 2 pµqq,
where we have used that z ¡ ϑpµ, upµqqpzq ¤ r n 2 pµ, upµqq ¥ pψ 1 pµ, u 1 pµqqq ¡1 pzq ¨¥ 0 ùñ z Ω 2 pµ, u 2 pµqq,

Introduction

State estimation is a task in which the quantity of interest is the 'true' state u true of a physical system over a space or space-time domain of interest. However, numerical prediction based on a given mathematical model may be deficient due to limitations imposed by available knowledge. In other words, the mathematical model can only take anticipated or parametric uncertainty into account. A more accurate prediction requires the incorporation of experimental observations in order to accommodate unanticipated or nonparametric uncertainty. The Parameterized-Background Data-Weak (PBDW) formulation for variational data assimilation is a data-driven reduced order modeling approach that was initially devised in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF] so as to merge prediction by model with prediction by data. The PBDW approach has been developed in order to estimate the true state u true for several configurations of a physical system. Supposing that the true state u true depends on some unknown parameter ω in an unknown parameter set Θ that represents the unanticipated uncertainty, the goal is to account for the dependency of the true state u true pωq on uncertain parameters by means of the sole knowledge of data. In this chapter, whenever the context is unambiguous, the parameter ω is dropped.

The formulation combines a so-called 'best-knowledge' (bk) model represented by a parametrized partial differential equation (PDE) and experimentally observable measurements. The use of data in the PBDW approach is fundamental not only to reconstruct the quantities of interest, but also to correct the possible bias in the mathematical bk model. The PBDW approach provides the following attractive features:

• The PBDW variational formulation simplifies the construction of a priori error estimates which can guide the optimal choice of the experimental observations.

• The PBDW formulation uses a bk model that accommodates anticipated uncertainty associated with the parameters of the model in a computationally convenient way. This bk model is typically built using model-order reduction techniques.

Note that the PBDW formulation does not explicitly include the equations of the bk model, but only a finite collection of solutions to the bk model. Thus, another important feature of the PBDW approach is its non-intrusiveness. In fact, once the subspace Z N has been generated, we no longer need the bk model. The PBDW approach was devised in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF] for steady problems. In this chapter, we propose an extension of the PBDW approach to time-dependent state estimation. Two main contributions to the standard PBDW approach are presented:

• We build appropriate background spaces for the time-dependent setting using the POD-greedy algorithm [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF].

• We propose a modified offline stage so as to alleviate its computational cost which can be sizeable in a time-dependent setting. The new offline stage allows for a better computational efficiency owing to a smaller online system. Moreover, it achieves substantial cost savings associated with data collection since it diminishes the number of observation sensors needed online. Note that this modified offline stage can also be applied to a steady framework.

This chapter is organized as follows. In Section 4.2, we recall the standard PBDW approach for steady problems as introduced in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF]. In Section 4.3, we extend the PBDW approach to the time-dependent framework. In Section 4.4, we discuss the offline stage. Finally, in Section 4.5, we illustrate our method by numerical results.

Parametrized-Background Data-Weak (PBDW) approach

In this section, we first introduce the notation that will be used throughout the chapter. Here, we focus on a time-independent setting. We consider a spatial domain (open, bounded, connected subset) Ω R d , d ¥ 1, with a Lipschitz boundary. We introduce a Hilbert space U composed of functions defined over Ω. The space U is endowed with an inner product p¤, ¤q and we denote by } ¤ } the induced norm; U consists of functions tw : Ω Ñ R | }w} Vu. To fix the ideas, we assume that H 1 0 pΩq U H 1 pΩq, and we denote the dual space of U by U I . The Riesz operator R U : U I Ñ U satisfies, for each U I , pR U p q, vq pvq, dv U. For any closed subspace Q U, the orthogonal complement of Q is defined as Q u : tw U | pw, vq 0, dv Qu.

(4.2)
Finally, we introduce a parameter set P R p , p ¥ 1, whose elements are generically denoted by µ P.

Best-knowledge (bk) model

The first source of information we shall afford ourselves in the PBDW approach is a so-called 'best-knowledge' (bk) mathematical model in the form of a parameterized PDE posed over the domain Ω (or more generally, over a domain Ω bk such that Ω Ω bk ). Given a parameter value µ in the parameter set P, we denote the solution to the bk parameterized PDE as u bk pµq U. Then, we introduce the manifold associated with the solutions of the bk model

M bk : tu bk pµq | µ Pu U. (4.3)
In ideal situations, the true solution u true is well approximated by the bk manifold, i.e., the model error bk mod pu true q : inf zM bk }u true ¡ z},

is very small.

We introduce nested background subspaces Z 1 . . . Z N . . . U that are generated to approximate the bk manifold M bk to a certain accuracy. These subspaces can be built using various model-order reduction techniques, for instance, the RBM described in the previous chapters. Note that the indices of the subspaces conventionally indicate their dimensions. To measure how well the true solution is approximated by the background space Z N , we define the quantity bk N pu true q : inf zZ N }u true ¡ z}.

(4.5)

When N is large enough, we have bk N pu true q bk mod pu true q.

(4.6)

Moreover, we introduce the reduction error bk red,N : sup

uM bk inf zZ N }u ¡ z}, (4.7)
which encodes the loss of accuracy caused by solving the bk model in the Ndimensional background space Z N . Figure 4.1 illustrates both the model and reduction errors, where Π Z N pu true q and Π M bk pu true q are closest points to u true in Z N and M bk , respectively. Note that Π Z N is the U-orthogonal projection onto Z N . The background space Z N can be interpreted as a prior space that approximates the bk manifold which we hope approximates well the true state u true . As previously alluded to, u true rarely lies in M bk in realistic engineering study cases.

Unlimited-observations statement

Let us first describe an ideal situation. The unlimited-observations PBDW statement reads: find pu ¦

N , z ¦ N , η ¦ N q U ¢ Z N ¢ U such that pu ¦ N , z ¦ N , η ¦ N q arginf u N U z N Z N η N U }η N } 2 , (4.8)
subject to pu N , vq pη N , vq pz N , vq, dv U, (4.9a)

pu N , φq pu true , φq, dφ U.

(4.9b)

A direct consequence of (4.9a) is that

u ¦ N z ¦ N η ¦ N , (4.10)
whereas (4.9b) implies that u ¦ N u true . We will be using the following terminology:

• State estimate: The goal of the PBDW statement (4.8)-(4.9) being to estimate the true state u true , the first component of its solution, u ¦ N , is called the 'state estimate'. In the present ideal situation of unlimited observations, the state estimate coincides with the true state. Since the bk mathematical model of a physical system is often deficient, one cannot realistically assume that the state estimate u ¦ N of u true lies completely in the bk manifold (or in the background space Z N ). Therefore, the update estimate η ¦ N is meant to cure the deficiency of the bk model by capturing unanticipated uncertainty. In other words, the key idea of the PBDW statement (4.8) is to search for the smallest correction to the bk manifold. The following result is proved in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF].

¦ N u true , z ¦ N Π Z N pu true q, η ¦ N Π Z u N pu true q.
(4.11)

Proof. We have already seen that u ¦ N u true . Next, we deduce from (4.9a) that

u true z ¦ N η ¦ N . Since (4.8) is a minimization of }η ¦ N }, it follows that z ¦ N Π Z N pu true q. Thus, η ¦ N Π Z u N pu true q.
The Euler-Lagrange saddle-point problem associated with the PBDW statement (4.8)-(4.9) reads: find pz ¦

N , η ¦ N q Z N ¢ U such that 5 pη ¦ N , qq pz ¦ N , qq pu true , qq, dq U, pη ¦ N , pq 0, dp Z N , (4.12)
and set

u ¦ N z ¦ N η ¦ N . (4.13)
As mentioned earlier, the saddle-point problem (4.12) is purely geometric and does not include any explicit reference to the bk model. The unique link to the bk model is through the background space Z N . Therefore, the PBDW approach is applicable to a wide class of engineering problems. Moreover, the non-intrusiveness of (4.12) simplifies its implementation.

Observable space

The evaluation of the right-hand side pu true , qq in (4.12) requires the full knowledge of the true state u true which is unrealistic. In practice, one can only afford a limited number of experimental observations of the true state u true . In the present setting, the experimental observations are interpreted as the application of prescribed observation functionals obs m U I for all m t1, . . . , M u such that the m-th experimental observation is given by obs m pu true q R, dm t1, . . . , M u. where x c m is the center of the Gaussian that reflects the location of the sensor, and

r m 3 |R m | 1
d is the standard deviation of the Gaussian that reflects the filter width of the sensor.

Generally, we introduce the M -dimensional experimentally observable space

U M : Spantq m u 1¤m¤M U, (4.17)
where q m : R U p obs m q is the Riesz representation of obs m U I , for all m t1, . . . , M u.

The experimental observations of the true state satisfy u true , q m ¨ obs m pu true q, dm t1, . . . , M u. Hence, for all q U M such that

q M m1 α m q m , with pα m q 1¤m¤M R M , (4.19)
the inner product pu true , qq can be deduced from the experimental observations as a linear combination of the M available observations:

u true , q ¨ M m1 α m u true , q m ¨ M m1
α m obs m pu true q.

(4.20)

Limited-observations statement

Let us now describe the PBDW statement in the case of limited observations. Henceforth, we make the crucial assumption that

Z N U u M t0u, (4.21)
which is meant to ensure the well-posedness of the PBDW statement with limited observations (cf. Proposition 4.3 below). This assumption can be viewed as a requirement to have enough sensors (note that Z N U u t0u). The limitedobservations PBDW statement reads: find pu ¦

N,M , z ¦ N,M , η ¦ N,M q U ¢ Z N ¢ U such that pu ¦ N,M , z ¦ N,M , η ¦ N,M q arginf u N,M U z N,M Z N η N,M U }η N,M } 2 ,
¦ N,M z ¦ N,M η ¦ N,M Z N pZ u N U M q.
Then, (4.25a) yields pη ¦ N,M , qq pu true , qq, dq Z u N U M .

We conclude that η ¦ N,M Π Z u N U M pu true q.

The decomposition of the state estimate u ¦ N,M is illustrated in Figure 4.2. Remark 4.1 (Perfect background space). The choice of the background space Z N and of the observable space U M may lead to several specific configurations. In particular, the background space Z N is said to be perfect if the reduction error (cf.

Algebraic formulation

We now present the algebraic formulation of the limited-observations PBDW statement. We first introduce an N -dimensional approximation space U N of the infinitedimensional space U as well as discrete approximation spaces Z N U N and U M U N of the subspaces Z N and U M , respectively. These spaces are built using finite elements [START_REF] Ern | Theory and practice of finite elements[END_REF]. We assume that the size of the mesh is small enough so that the N -dimensional space discretization delivers High-Fidelity (HF) approximations within the requested level of accuracy. To alleviate the notation, we have dropped the superscript N ; hence, the discrete FEM spaces are denoted Z N and U M instead of Z N N and U N M , but we still keep the notation U N for the FEM-discretization space.

Then, we introduce a basis for the background space Z N : Spantζ n u 1¤n¤N . The update space is spanned by the Riesz representations of the observation functionals in U N , i.e., U M : Spantq m u 1¤m¤M , where q m U N , for all m t1, . . . , M u.

The high-fidelity (HF) discretization of the saddle-point problem (4.25) is: Find

pz ¦ N,M , η ¦ N,M q Z N ¢ U M such that pη ¦ N,M , qq pz ¦ N,M , qq pu true , qq, dq U M , (4.27a) pη ¦ N,M , pq 0, dp Z N . (4.27b)
The solution of (4.27) is then searched under the form

z ¦ N,M N ņ1 z n ζ n , and η ¦ N,M M m1 η m q m , (4.28)
and we introduce the component vectors z ¦

N,M : pz n q 1¤n¤N R N and η ¦ N,M : pη m q 1¤m¤M R M .

(4.29)

We also introduce the basis matrices Z N R N ¢N and U M R N ¢M whose column vectors are the components of the functions tζ n u 1¤n¤N and tq m u 1¤m¤M respectively in the basis of U N . In algebraic form, the FEM-discretized saddle-point problem (4.27) reads: find pz ¦

N,M , η ¦ N,M q R N ¢ R M such that ¢ A B B T 0 ¢ η ¦ N,M z ¦ N,M ¢ obs M 0 , (4.30)
with the matrices is equal to the cosine of the angle between the linear subspaces Z N and U u M (cf. Figure 4.2). Furthermore, it is readily verified that

A ¡ pq m I, q m q © 1¤m,m I ¤M R M ¢M , B ¡ pζ n , q m q © 1¤m¤M,1¤n¤N
β N,M 0 ðñ Z N U u M $ t0u, (4.35a) β N,M 1 ðñ Z N U M . (4.35b)
The case (4.35b) can hardly occur in practice with a reasonable (not too high) number of observation sensors. Loosely speaking, a sensor is localized in space. Thus, it concerns only a limited number of degrees of freedom.

Remark 4.3 (Insufficient observations

). If M N , then (4.30) is necessarily illposed. Indeed, we have M N ðñ dimpU M q dimpZ N q.

Moreover, we have

dimpZ N q ¡ codimpU u M q ¤ dimpZ N U u M q ðñ dimpZ N q ¡ dimpU M q ¤ dimpZ N U u M q ùñ 1 ¤ dimpZ N U u M q ðñ Z N U u M $ t0u ðñ β N,M 0,
where the last equivalence follows from Proposition 4.3.

In practice, the matrices A and B are computed using the algebraic formulas

A U M MU M , and B Z N MU M , (4.36)
where M is the Gram matrix of the inner product in U. Thus, solving (4.30) allows for a straightforward reconstruction of the components of the state estimate in the basis of U N as follows:

u ¦ N,M Z N z ¦ N,M U M η ¦ N,M .
(4.37)

Time-dependent PBDW

Consider a finite time interval I r0, T s, with T ¡ 0. To discretize in time, we consider an integer K ¥ 1, we define 0 t 0 ¤ ¤ ¤ t K T as pK 1q distinct time nodes over I, and we set K tr t1, . . . , Ku, K tr t0u K tr and I tr tt k u kK tr .

Remark 4.5 (Initial condition). In the present setting, we choose not to solve the PBDW statement for the initial time node k 0. It is straightforward to consider a setting where the initial time node is also included.

Unlimited-observations statement

In this ideal setting, we assume that u true C 0 pI; Uq. The time-dependent unlimitedobservations PBDW statement reads: for each k K tr , find pu

k,¦ N , z k,¦ N , η k,¦ N q U ¢ Z N ¢ U such that pu k,¦ N , z k,¦ N , η k,¦ N q arginf u N U z N Z N η N U }η N } 2 , (4.43)
subject to pu N , vq pη N , vq pz N , vq, dv U,

pu N , φq pu k,true , φq, dφ U.

where u k,true : u true pt k , ¤q. For each k K tr , the solution of (4.43)-(4.44) is given by

u k,¦ N u k,true , z k,¦ N Π Z N pu k,true q, η k,¦ N Π Z u N pu k,true q. ( 4
, η k,¦ N q Z N ¢ U such that pη k,¦ N , qq pz k,¦ N , qq pu k,true , qq, dq U, (4.46a) pη k,¦ N , pq 0, dp Z N . (4.46b)
The unlimited-observations state estimate is then

u k,¦ N z k,¦ N η k,¦ N , dk K tr .
(4.47)

Limited-observations statement

We now weaken the regularity assumption on the true state and only assume that u true L 1 pI; Uq. We introduce the time-integration intervals

I k rt k ¡ δt k , t k δt k s, dk K tr , (4.48)
where δt k ¡ 0 is a parameter related to the precision of the sensor (ideally, δt k minpt k 1 ¡ t k , t k ¡ t k¡1 q with obvious adaptation if k=K). Then, for any function v L 1 pI; Uq, we define the time-averaged snapshots

v k pxq : 1 |I k | » I k vpt, xq dt U, dk K tr . (4.49)
As in the steady case, we consider observation functionals that render the behaviour of given sensors. We use the same observation functionals as in the time-independent setting, but we let them act on the time-averaged snapshots of the true solution, i.e., we consider k,obs m pu true q : obs m pu k,true q, dm t1, . . . , M u, dk K tr .

(

pu true q 1 2πr 2 m » Rm u k,true pxq exp ¢ ¡px ¡ x c m q 2 2r 2 m dxdt, 1 |I k | 1 2πr 2 m » I k » Rm u true pxq exp ¢ ¡px ¡ x c m q 2 2r 2 m dxdt. (4.52)
Generally, we introduce the time-independent observable space U M U such that U M Spantq 1 , . . . , q M u. The observation functionals in U I are then defined as k,obs m pu true q pu k,true , q m q, dm t1, . . . , M u, dk K tr .

(4.54)

Note that, for fixed sensor locations, the computational effort to compute the Riesz representations of the observation functionals is time-independent and is incurred only once so that the experimental observations of the true state satisfy k,obs m

pu true q u k,true , q m ¨ 1 |I k | » I k
obs m pu true pt, ¤qqdt, dm t1, . . . , M u, dk K tr .

(4.55)

Hence, for all q U M such that,

q M m1 α m q m , (4.56)
the inner product u k,true , q ¨is deduced from the experimental observations as follows:

u k,true , q ¨ 1 |I k | » I k M m1 α m u true pt, ¤q, q m ¨dt 1 |I k | M m1 α m » I k
obs m pu true pt, ¤qqdt.

(4.57) We are now ready to write the limited-observations PBDW statement: for each pu N,M , φq pu k,true , φq, dφ U M .

k K tr , find pu k,¦ N,M , z k,¦ N,M , η k,¦ N,M q U ¢ Z N ¢ U such that pu k,¦ N,M , z k,¦ N,M , η k,¦ N,M q arginf u N,M U z N,M Z N η N,M U }η N,M } 2 , ( 4
(4.59b)

The limited-observations saddle-point problem associated with (4.58) reads: for

each k K tr , find pz k,¦ N,M , η k,¦ N,M q Z N ¢ U M such that pη k,¦ N,M , qq pz k,¦ N,M , qq pu k,true , qq, dq U M , (4.60a
) pη k,¦ N,M , pq 0, dp Z N , (4.60b)
and the limited-observations state estimate is The assumption (4.62) is typically reasonable for a sensor of small precision δt k .

u k,¦ N,M z k,¦ N,M η k,¦ N,M , dk K tr . ( 4
In algebraic form, the limited-observations PBDW statement reads: for each

k K tr , find pz k,¦ , η k,¦ q R N ¢ R M such that ¢ A B B T 0 ¢ η k,¦ z k,¦ ¢ k,obs 0 , (4.63)
with the matrices

A ¡ pq m I, q m q © 1¤m,m I ¤M R M ¢M , B ¡ pζ n , q m q © 1¤m¤M,1¤n¤N R M ¢N , (4.64)
and the vector of observations k,obs obs m pu k,true q ¨1¤m¤M R M .

(4.65)

Similarly to the steady PBDW linear system (4.30), we solve (4.63) through an offline/online decomposed computational procedure whenever several realizations u true pωq of the true state are to be considered.

Offline stage

In this section, we discuss the offline stage. Our main goal is to address the construction of the background space Z N .

Background space construction via POD-greedy

Suppose that we have computed a set of high-fidelity (HF) trajectories S pS k q kK tr pu k pµqq µP tr ¨kK tr ,

where u k pµq : upµqpt k , ¤q, for all k K tr . If we were to consider the PBDW statement (4.58)-(4.59) for each k K tr as an independent steady PBDW statement, we would be working with the time-dependent background spaces

Z k N k PODpS k , pod q, dk K tr , (4.67)
where the procedure POD is defined in Chapter 1. However, this strategy is not convenient since the sizes N k of the background spaces Z k N k would depend on k. Since the observable space U M is fixed, the same non-homogeneity between time nodes would also arise in the stability constant β N k ,M . Thus, we propose to apply a POD-greedy algorithm [START_REF] Haasdonk | Convergence rates of the POD-greedy method[END_REF] in order to build a time-independent background space Z N that will be used for all k K tr . The advantage is that the PBDW matrices A and B and the stability constant β N,M remain unchanged regardless of the discrete time node. The offline stage using the POD-greedy algorithm is summarized in Algorithm 4.1.

Algorithm 4.1 Offline stage via POD-greedy

Input : S and pod .

Q init : an initial set of Riesz representations for the observations.

1: Compute Z N : POD-greedypS, pod q. 2: Set U M : SpantQ init u.

3: Compute the matrices A and B using Z N and U M .

Output : Z N , U M , A and B.

Background space construction via state estimation

We now devise a new algorithm in the context of time-dependent PBDW to perform the offline stage. Here, the construction of the background space Z N , the choice of the observation space U M and the PBDW matrices are modified. The key idea of the new procedure is to precompute the PBDW state estimates of the parameters in the training set P tr during the offline stage. The background space is then deduced from these PBDW state estimates. The benefit is that the newly created background space incorporates data-based knowledge. The modified offline stage of the PBDW for time-dependent problems is described in Algorithm 4.2. Within the modified offline stage, we use the so-called 'Greedy stability maximization' (S-Greedy) algorithm (considered in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF]) in line 6 in order to identify the least stable mode and then take the best measurement. The algorithm uses an input space Z N that results from a POD-greedy procedure so that Z 1 contains the dominant mode, and so forth. The S-greedy algorithm selects the observations progressively. Thus, the enrichement of the observable space U M stops once the minimum stability β min is reached. the modified offline stage in the proposed algorithm offers four major advantages:

• Improved background space: Since the background space Z N is built using both the bk model and the observations, it is expected to have better approximation capacities of the true state.

• Reduced number of online observations: In line 6 of Algorithm 4.2, we select each new data point so as to maximize the stability constant β N,M . Thus, the observations that will be used during the online stage are mainly needed only for stability and not for accuracy.

• Reduced dimension of the online PBDW statement: Since the number of observations is significantly reduced, the modified PBDW matrices are of smaller size compared to the matrices of the standard PBDW. Thus, using the modified offline algorithm, the online PBDW formulation is solved faster.

• Reduced storage cost: Owing to the reduced number of online measurements, the dimensions of the observable space U M and of the matrices A and B are smaller, whence the storage gain.

Regarding computational efficiency, the modified procedure consists of more steps than in the standard PBDW. However, all the additional steps of the algorithm are Q init : an initial set of Riesz representations of the observations.

1: Choose a random q 1 Q init . 2: Set U 1 : Spantq 1 u. Compute the least stable mode and the associated supremizer

w inf arginf wZ N sup vU m¡1
pw, vq }w} }v} , and v sup Π U m¡1 pw inf q.

7:

Identify the least well-approximated vector q m argsup qQ init |pq, w inf ¡ v sup q|.

8:

Set U m : SpantU m¡1 , q m u.

9:

Compute the stability constant β N,m .

10:

m m 1.

11: end while 12: M : m.

Output : U M . performed offline. As for all reduced-order modeling techniques, the goal of the algorithm is to further improve the online efficiency. Hence, the computational savings brought by the new PBDW formulation come, in our opinion, at a reasonable offline price. Indeed, the resolution of the (online) standard PBDW statement for each parameter µ P tr has a reduced computational cost. The only relevant additional computational cost incurred offline is related to the second POD-greedy (cf. line 5 of Algorithm 4.2). We believe this computational effort remains acceptable.

Remark 4.8 (Least stable mode). Line 6 of Algorithm 4.3 may return several infimums. Among these infimums, we select a function whose norm in U is maximal.

Remark 4.9 (Steady setting). In a time-dependent framework, the computational savings induced by the modified offline stage are substantial; in particular because of the influence of the time steps. However, Algorithm 4.3 can be applied in the steady setting as well.

Numerical results

In this section, we illustrate the above developments on test cases related to the heat equation. The goal is to illustrate the computational performance of our algorithms. In all our test cases, we consider a two-dimensional setting based on the tions use the software FreeFem++ [START_REF] Hecht | New developments in freefem++[END_REF], whereas the reduced-order modeling and the PBDW-related algorithms have been developed in Python.

Physical model problem

We apply the above methodology to the following parabolic PDE: For many values of the parameter µ P, find upµq : I ¢ Ω Ñ R such that

Synthetic data generation

We synthesize the data by first synthesizing a true solution and then applying to it the linear functionals by means of their Riesz representations in the observable space U M . In order to synthetize the true solution, we consider a 'true model' based on the bi-material plate (cf. Figure 4.4) where we choose a fixed internal diffusivity

D

Test configurations

In order to investigate the PBDW formulation, we perform test cases on two distinct configurations:

1. Perfect model: The bk model is said to be perfect when bk mod pu true pωqq 0, for every ω Θ (see (4.4)) (we recall that ω represents the unanticipated uncertainty). In this situation, u true pωq M bk for all ω Θ. Although the model is perfect, some discrepancies between the HF solutions and the measurements might arise from model-order reduction since M bk $ Z N (cf. Note that this scenario seldom occurs in engineering situations. This test configuration is meant to assess the accuracy of the PBDW formulation when the observable space U M scarcely has additional information compared to Z N .

Imperfect model:

The bk model is said to be imperfect when the modeling error does not vanish. In this situation, there exists at least one (and in general many) ω Θ such that bk mod pu true pωqq $ 0, i.e., u true pωq M bk . Consider for instance the plates in Figures 4.3 and 4.4. If the true solution is generated synthetically using the bimaterial plate, an example of an imperfect bk model can be the one for which we solve the same PDE that has generated the true states without accounting for the difference in diffusivity between the subdomains of the plate.

Background space construction via POD-greedy

In this section, four test cases are considered to study the PBDW approach.

• Test case (a): Linear perfect.

• Test case (b): Linear imperfect.

• Test case (c): Nonlinear perfect.

• Test case (d): Nonlinear imperfect. profiles for the model problem (4.68)¡(4.69) over the homogeneous plate at the end of the simulation, i.e., for t K 10s and for two parameter values. We recall that these solutions will be used as true solutions for the perfect linear case. As the time evolves, the energy related to the flux φ e propagates through the plate which is progressively heated. Moreover, the overall temperature is higher for smaller values of the parameter µ than for larger values. As physically expected, the thermal diffusion over the plate is stronger for larger values of µ than for smaller values.

Test case (a): Linear perfect model We consider the case of a perfect bk model for which the diffusivity is uniform over the entire domain Ω. Thus, the true solutions correspond to the HF computations of the bk model. The resulting trajectories are reduced using the POD-greedy algorithm. For instance, for a tolerance value pod 10 ¡2 , the background space Z N is composed of N 5 modes. Regarding observations, the initial set Q init is obtained using M CardpQ init q 121 sensors that are uniformly placed over the plate (see Figure 4.3). Using both the background space Z N and the observable space U M pQ init q, we build the offline matrices A and in diffusivity between Ω int and Ω ext leads as expected to a kink in the temperature isolines. When µ 1, the thermal diffusion is stronger in the inner plate corresponding to Ω int , whereas for µ ¡ 1, the thermal diffusion is weaker in the inner plate.

Using the HF trajectories produced by the bk model, we generate a background space Z N by means of a POD-greedy algorithm. We use M 121 observations to build the observable space U M . The relative H 1 -errors e k pµq defined in (4.76) are shown in Figure 4.8 as a function of the value of the parameter µ. For instance, for a tolerance value pod 10 ¡3 , Z N is spanned by N 7 vectors. Notice that the error vanishes for µ 0.5 since this configuration is equivalent to a perfect bk model. However, the bottom panels of Figure 4.8 show a gradual error increase with the dimension N of the bk space. This tendency was already observed for the linear perfect test case, although in smaller proportions. As before, the stability constant β N,M is degraded when increasing the dimension N of the background space Z N . Moreover, the enrichment of Z N does not add relevant modes anymore (in terms of associated singular values). For the sake of comparison, we enrich the observable space U M such that M 676 and plot the relative H 1 -errors e k pµq for the same values of pod in Figure 4.9. Our interpretation is confirmed since the stability issues do not arise anymore for pod 10 ¡4 . Owing to the increase of M , the stability decrease with respect to N is somewhat compensated. Finally, the bottom-right panel of Figure 4.9 shows the beginning of an error increase. Using the same reasoning as above, we conclude that more observations are needed for pod 10 ¡6 .

Let us now investigate the relative H 1 -errors as a function of the dimension M of the observable space U M . Figure 4.10 shows that the larger the set of observations, the smaller the error. Finally, we visualize the stability constant β N,M as a function of the number of observations M in Figure 4.11. The left panel of the figure shows a single curve for clarity, whereas the right panel includes curves for several values of the tolerance pod (note that the two panels do not use the same rule). As expected, for a constant value of N , the more the observations, the better the stability. For a number of observations M 3000, the PBDW formulation is perfectly stable (or

Test case (c): Nonlinear perfect model

We consider the case with a perfect bk model. Thus, the true solutions correspond to the HF computations of the bk model (cf. left panel of Figure 4.12). The resulting trajectories are reduced using the POD-greedy algorithm. For instance, for a tolerance value pod 10 ¡2 , the background space Z N consists of N 3 modes. Regarding observations, the initial set Q init is obtained using M CardpQ init q 121 sensors that are uniformly placed over the plate (see Figure 4.3). During the online stage, we estimate the state u ¦ N,M for every parameter µ in the training set P tr . In Figure 4.13, we display the state estimation relative H 1 -error e k pµq defined in (4.76) as a function of the value of the parameter µ for several values of pod . In contrast to the linear case, the error always decreases for smaller tolerances pod , i.e., with the dimension N of the background space Z N . However, we expect that, for some very small tolerance value (e.g. pod such that N ¡ M ), the stability issues mentioned above would arise again. observations to build the observable space U M . For pod 10 ¡3 , Z N is spanned by N 5 vectors. Notice that the error vanishes for µ 0.25 since this configuration is equivalent to a perfect bk model. We notice that the relative H 1 -error e k pµq increases because the stability constant decreases. Figure 4.15 visualizes the relative H 1 -error e k pµq for a higher number of observations M 676. We observe that augmenting the dimension of the observable space U M cures the stability issues. Also, the errors are lower owing to the higher number of observations. Finally, Figure 4.16 shows the stability constant β N,M as a function of the number of observations M . The behaviour is quite similar to the linear case. Hence, the nonlinear character of the problem does not influence the overall features of the PBDW statement. This observation corroborates the independence with regard to the bk model. N,M has a larger norm. This is due to the offline inclusion of observations in the new background space Z N through offline state estimation. Therefore, we deduce that the modified offline algorithm achieves the expected objective. CHAPTER 5

CONCLUSIONS AND PERSPECTIVES

In this thesis, we have devised three new methodologies in the Reduced Basis (RB) context. First, we have introduced a Progresssive RB Empirical Interpolation Method (PREIM) that allows one to diminish the offline costs incurred in the nonlinear RB method applied to unsteady nonlinear PDEs. The main reason for this computational gain is that the computation of high-fidelity trajectories is the dominant part of the offline cost. Numerical tests on both two-dimensional and industrial three-dimensional heat transfer problems with nonlinear thermal conductivities have illustrated the computational efficiency and the accuracy of the proposed algorithm.

In the present study, the computational benefits of PREIM have been evaluated by comparing the results to those of the reference method, i.e. the standard RB-EIM and also to the Simultaneous EIM-RB method (SER) which is the closest progressive method to PREIM available in the literature. Comparisons with other progressive RB-EIM methods, such as those mentioned in the introductive section of Chapter 2, and a more theoretical study of PREIM can be considered for future work. Another relevant perspective is the application of PREIM to a more systematic study of 3D flow regulation in industrial applications related to nuclear reactor operation. Second, we have presented a RB scheme for parametrized variational inequalities, more particularly in the framework of contact mechanics. We have addressed the conceptional issue raised by nonlinear inequality constraints in the RB scheme when considering a general setting where neither small displacements nor matching meshes are assumed, as was the case so far in the literature. The present RB scheme allows one to reduce the dual basis of Lagrange multipliers while maintaining their positivity by means of a greedy algorithm based on the projection onto a convex cone. A systematic comparison of the present method to the few methods so far available in the literature such as the Non-negative Matrix Factorization (NMF) and Hyper-Reduction (HR) is a relevant perspective. The method that we have proposed can be used under a local convexity assumption on the solids that come into contact. A challenging research perspective is the extension of our developments to concave solids and to multi-body contact problems. Moreover, we have only considered equilibrium problems. Thus, the treatment of dynamic contact problems remains on the agenda.

Third, we have addressed data assimilation in the context of RB schemes. Our first contribution is the extension of the Parametrized-Background Data-Weak (PBDW) approach to the time-dependent setting by means of a POD-greedy reduced basis construction. Since the construction of the basis is performed offline, the algorithm renders the time dependence of the problem we are addressing while the time stepping scheme remains unchanged. An interesting research direction is to devise an online time-stepping PBDW scheme so as to take advantage of the potential interactions between the time steps. As a second contribution, we have devised a modified offline algorithm that exploits offline state estimates in order to diminish both the dimension of the online PBDW statement and the number of required sensors collecting data. The idea is to exploit in situ observations in order to update the best-knowledge model, thereby improving the approximation capacity of the background space. The proposed algorithm is sequential, i.e., we first build the background space before chosing the observable space. An interesting perspective is to study a progressive construction of both spaces. As regards the quality of the spaces built offline, the PBDW suggests separate criteria for the background space and the observable space. Hence, another relevant research direction is to study a joint stopping criterion for the construction of both spaces.

Figure 1 . 1 -

 11 Figure 1.1 -Prototype du robinet RAMA instrumenté. À gauche: image réelle du robinet. À droite: Simulation numérique.

Figure 1 . 2 -

 12 Figure 1.2 -Exemple de contact entre six goujons et une pièce du robinet RAMA (Remerciement à VELAN SAS pour la fourniture des plans nécessaires pour le maillage).

( 1 . 35 )

 135 Une fois le paramètre µ m déterminé, on pose r m p¤q : γpµ m ; ¤q ¡ γ m¡1 pµ m ; ¤q, et x m argmax xΩ tr |r m pxq|.

( 1 . 36 )

 136 Est vérifiée ensuite la satisfaction ou pas du critère |r m px m q| eim pour une tolérance eim ¡ 0 définie par l'utilisateur. Si ce critère est satisfait, on définit la nouvelle fonction EIM par

•

  Dans MAP sont intégrés les codes de réduction de modèle pour la thermique linéaire pour des problèmes d'échange thermique avec conditions aux limites de Robin. De plus, la partie hors-ligne d'un problème de thermique non-linéaire avec une non-linéarité qui porte sur la conductivité thermique y est finalisée et la partie en-ligne est bien avancée. L'intégration a été initiée pour une paramétrisation qui porte sur la conductivité thermique pour un clapet de régulation de débit (cf Figure1.4). Ensuite, elle a été élargie à une paramétrisation qui porte sur le coefficient d'échange thermique aussi bien pour le clapet de régulation de débit que pour le robinet RAMA (cf Figure1.1) dans le cadre du stage de fin d'études d'Alaeddine Jlaiel, encadré par l'auteur du manuscrit et Sébastien Meunier.

Figure 1 . 4 -

 14 Figure 1.4 -Demi-volume du clapet de régulation de débit.

) and ûk¡1 µ 4 : 5 :k k 1 6 2 :

 4512 Solve the reduced system (2.16) to obtain ûk µ Set : end while Output : pû k µ q kK tr Algorithm 2.2 Standard EIM Input : P tr , K tr , Ω tr , and eim ¡ 0 1: Compute S pu k µ q pµ,kqP tr ¢K tr P HF trajectories Set m 1 and

¤q and x m argmax xΩ tr |r m pxq| 5 :

 5 while (|r m px m q| ¡ eim q do 6:

Chapter 2 .

 2 PREIM for nonlinear parabolic problems Algorithm 2.4 UPDATE RB Input : Θ pθ n q 1¤n¤N , S, and pod ¡ 0 1: if S r then 2: end if Output : Θ At the iteration m ¥ 1 of PREIM, the trajectories for all µ P HF m are HF trajectories, whereas they are approximated by RB trajectories for all µ P tr zP HF

  function γm pµ, k, xq : Γpµ, ūm,k µ pxqq.

(

  ii) produce a rank-m approximation of the nonlinear function γm defined by (2.23) in the form γm rP HF m ,Xm,Qms pµ, k, xq : m j1

13 :

 13 Define r m p¤q r r m p¤q discard the EIM selection 14:Set X out X m¡1 and Q out Q m¡1 15: else 16:

Figure 2 . 1 -

 21 Figure 2.1 -Test cases (a) and (b): The computational domain is a perforated plate.

Figure 2 .

 2 Figure 2.1 with Ω p¡2, 2q 2 zr¡1, 1s 2 R 2 . We consider the nonlinear parabolic problem (4.68) with the nonlinear function Γpµ, vq : sin

Figure 2 . 2 -

 22 Figure 2.2 -Test case (a): HF solutions for the parameter values µ 1 (left) and µ 18 (right) at t 2 s (top) and t 5 s (bottom).

  m } V pΩ tr q 2.0 8.1E¡1 1.1E¡1 5.2E¡3 2.6E¡3 1.1E¡3 1.6E¡4

Figure 2 . 3 -

 23 Figure 2.3 -Test case (a): EIM approximation error as a function of m for pod 10 ¡3 and eim 5 ¤ 10 ¡2. Left: Errors for the standard RB-EIM procedure, PREIM, and U-SER. Right: Errors }r r m } V pΩ tr q and }r m } V pΩ tr q for PREIM.

Figure 2 . 4 -

 24 Figure 2.4 -Test case (a): RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq for pod 10 ¡3 and eim 5 ¤ 10 ¡2 .

Figure 2 . 6 -

 26 Figure 2.6 -Test case (a): RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq for pod 10 ¡3 and eim 5 ¤ 10 ¡2 .

Figure 2 . 7 -

 27 Figure 2.7 -Test case (b): HF solutions for the parameter values µ 1 at t 1 s (left, values from 20.2 to 22.1) and at t 2.5 s (right, values from 34.5 to 37.3).

Figure 2 . 8 -

 28 Figure 2.8 -Test case (b):EIM approximation errors }r r m } V pΩ tr q and }r m } V pΩ tr q as a function of m for PREIM with pod 5 ¤ 10 ¡2 and eim 10 ¡3 .

Figure 2 . 9 -

 29 Figure 2.9 -Test case (b): RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq for pod 5 ¤ 10 ¡2 and eim 10 ¡3 .

Figure 2 . 10 -

 210 Figure 2.10 -Test case (b): RB approximation error }u µ ¡ ûµ } 2 pI tr ;H 1 pΩ tr qq for pod 2.5 ¤ 10 ¡2 and eim 10 ¡4 .

Figure 2 . 11 -

 211 Figure 2.11 -Test case (c): half-section of an industrial flow-regulation valve (Acknowledgment to VELAN SAS for providing the design data necessary for the mesh).

Figure 2 . 14 -

 214 Figure 2.14 -Test case (c):

  , vq ¡ f pµ; vq kpµ, upµq; upµqq ¤ gpµ, upµqq a.e. on Γ c 1 .

  Lagrangian formulation, the optimization problem (3.1) is rewritten as a saddle-point problem. Specifically, (3.1) can be recast as: Find pupµq, λpµqq V ¢W such that pupµq, λpµqq arg minmax vV,ηW

Figure 3 . 1 -

 31 Figure 3.1 -Generic two-body contact problem. For simplicity of the scheme, the entire boundary is taken to be the potential contact boundary.

Figure 3 . 2 -

 32 Figure 3.2 -Left: Reference domain q Ω and mesh with N 1350. Right: HF displacement field u N pµq for µ 1.12m.

Figure 3 . 3 -

 33 Figure 3.3 -Normal contact stress components for the HF solutions. Vanishing values correspond to nodes where the contact between the two disks is not established at equilibrium.

pod 10 ¡2

 10 10 ¡3 10 ¡4 10 ¡5 4.10 ¡6

Figure 3 . 4 -

 34 Figure 3.4 -Offline basis construction. Left: Singular values resulting from the POD for the primal space qV N . Right: Approximation error for the dual space | W R .

Figure 3 . 5 -

 35 Figure 3.5 -EIM error as a function of the rank M of the EIM approximation. Left: for the nonlinear gap map gpµ, upµqq. Right: for the nonlinear contact map Kpµ, upµqq.

Finally, Figure 3 . 7 Γ c 1 min 0 ,Figure 3 . 7 -

 371037 Figure 3.7 -Estimation of the interpenetration e inter pµq.

Figure 3 . 8 -

 38 Figure 3.8 -Error quantification with R 10 for the NMF and R 4 for the cone-projected greedy. Left: Error on the minimum energy e ener pµq. Right: Relative H 1 -error for the displacement field e displ pµq.

Figure 3 . 10 -

 310 Figure 3.10 -Error quantification with R 5 for both algorithms. Left: Error on the minimum energy e ener pµq. Right: Relative H 1 -error for the displacement field e displ pµq.

Figure 3 . 11 -

 311 Figure 3.11 -Estimation of the interpenetration e inter pµq. Left: R 18 for the NMF and R 5 for the cone-projected greedy. Right: R 5 for both algorithms.

Figure 4 . 1 -

 41 Figure 4.1 -Model and reduction errors.

•

  Deduced background estimate: The first contribution z ¦ N in (4.10) lies in the background space Z N and is deduced from the PBDW statement, which takes the observations into account. Hence, z ¦ N is called the 'deduced background estimate'. • Update estimate: The second contribution η ¦ N in (4.10) is brought by the inclusion of the observations in the PBDW statement. The observations supplement the bk model. Thus, η ¦ N is called the 'update estimate'. The deduced background estimate z ¦ N can only represent anticipated uncertainty.

Proposition 4 . 1 (

 41 Unlimited observations). The solution of (4.8)-(4.9) is given by u

( 4 . 14 )

 414 One can consider any observation functional that renders the behaviour of some physical sensor. In the case of sensors measuring the state locally over user-defined subsets R m Ω, where m t1, . . . , M u, one possibility is to model each sensor through uniform local integration obs m pvq 1 |R m | » Rm vpxqdx, dv U.

(4. 15 )

 15 Another plausible option is, as introduced in[START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF], to consider obs

24 )

 24 pu N,M , vq pη N,M , vq pz N,M , vq, dv U, (4.23a) pu N,M , φq pu true , φq, dφ U M . (4.23b) As above, (4.23a) implies that the limited-observations state estimate u ¦ One can show (e.g., by introducing the Lagrangian) that the limited-observations problem (4.22)-(4.23) is equivalent to the limited-observations saddle-point problem: find pz ¦ N,M , η ¦ N,M q Z N ¢ U M such that pη ¦ N,M , qq pz ¦ N,M , qq pu true , qq, dq U M , (4.25a) pη ¦ N,M , pq 0, dp Z N , (4.25b) and define u ¦ N,M according to (4.24). We will see in Proposition 4.3 below that the linear system (4.25) is well posed under the assumption (4.21).

Figure 4 . 2 -

 42 Figure 4.2 -PBDW state estimation (courtesy of A. T Patera).

Proposition 4 . 2 (

 42 Update estimate). The update estimate is given byη ¦ N,M Π Z u N U M pu true q.

(4. 26 )

 26 Proof. By definition of the saddle-point problem (4.25), η ¦ N,M U M . From (4.25b), we infer that η ¦ N,M Z u N . Hence, η ¦ N,M Z u N U M . Thus, we have u

Figure 4 .

 4 Figure 4.1) vanishes, i.e., u true Z N . In this case, the pair pu true , 0q Z N ¢ U M is the unique solution to (4.25). Hence, the state estimate u ¦ N,M also belongs to Z N

Remark 4 . 7 (

 47 PBDW matrices). Notice that the PBDW matrices A and B are time-independent; only the right-hand side in (4.63) depends on k.

7 :

 7 Compute the matrices A and B using Z N and U M . Output : Z N , U M , A and B.

Algorithm 4 . 3 S

 43 -Greedy: Stability-maximization algorithm Input : N , Z N and β min p0, 1s.

Figure 4 . 3 -

 43 Figure 4.3 -Computational domain and mesh with N 6561. The little black squares are the observation subsets tR m u 121 m1 .

Figure 4 . 4 -

 44 Figure 4.4 -Computational domain and mesh for a bi-material plate with N 6561. The little black squares are the observation subsets tR m u 121 m1 .

Figure 4 .

 4 1).

4. 5 . 3 . 1 Figure 4 . 5 -

 53145 Figure 4.5 -Test cases (a) and (b) : HF solutions for the bk model with Neumann boundary conditions. Left: µ 1 in D uni (values from 17.3 o C to 118.3 o C). Right: µ 20 in D uni (values from 30.7 o C to 43.3 o C).

B.Figure 4 . 6 - 4 .

 464 Figure 4.6 -Test case (a): Relative H 1 -error e k pµq for some time nodes k K tr . Top left: pod 10 ¡2 (N 5). Top right: pod 10 ¡3 (N 7). Bottom left: pod 10 ¡4 (N 10). Bottom right: pod 10 ¡6 (N 15).

Figure 4 . 7 -

 47 Figure 4.7 -Test case (b) : Synthetic true solutions with Neumann boundary conditions. Left: µ 1 in D syn pµq (values from 17.4 o C to 117.9 o C). Right: µ 20 in D syn pµq (values from 29.8 o C to 44.3 o C).

Figure 4 . 8 -

 48 Figure 4.8 -Test case (b) : Relative H 1 -error e k pµq for some time nodes k K tr and M 121. Top left: pod 10 ¡2 (N 5). Top right: pod 10 ¡3 (N 7). Bottom left: pod 10 ¡4 (N 10). Bottom right: pod 10 ¡6 (N 15).

Figure 4 . 9 -

 49 Figure 4.9 -Test case (b) : Relative H 1 -error e k pµq for some time nodes k K tr and M 676. Top left: pod 10 ¡2 (N 5). Top right: pod 10 ¡3 (N 7). Bottom left: pod 10 ¡4 (N 10). Bottom right: pod 10 ¡6 (N 15).

Figure 4 . 10 -

 410 Figure 4.10 -Test case (b) : Relative H 1 -error e k pµq as a function of the number M of observations for t K 10.

Figure 4 . 11 -

 411 Figure 4.11 -Test case (b) : Stability constant β N,M as a function of M . On the right panel, the values of N are respectively 3, 5, 7, 10, 13, 15 for the values of pod in decreasing order.

4. 5 . 3 . 2

 532 Nonlinear caseHere, we consider the PDE (4.68)¡(4.72) with u r 303.15K, σ 5.67 ¢ 10 ¡8 W.m ¡2 .K ¡4 and ε 3.10 ¡3 . Except for the parameter interval P r0.1, 2s, the set P tr t0.1i, 1 ¤ i ¤ 20u and the time step ∆t k 0.1, all the other numerical data remain the same as for the linear test case from the previous section.

Figure 4 .

 4 Figure 4.12 -Test cases (c) and (d) : Left: HF solution for the bk model (values from 17.80 o C to 18.25 o C). Right: Synthetic true solution using a bi-material plate (values from 17.90 o C to 18.23 o C).

4. 5 . 3 . 4 Figure 4 .

 5344 Figure 4.14 shows the relative H 1 -error e k pµq defined in (4.76) using M 121

Figure 4 . 13 -

 413 Figure 4.13 -Test case (c) : Relative H 1 -error e k pµq for some time nodes k K tr and M 121. Top left: pod 10 ¡2 (N 3). Top right: pod 10 ¡3 (N 5). Bottom left: pod 10 ¡4 (N 7). Bottom right: pod 5.10 ¡6 (N 11).

Figure 4 . 14 -

 414 Figure 4.14 -Test case (d) : Relative H 1 -error e k pµq for some time nodes k K tr and M 121. Top left: pod 10 ¡2 (N 3). Top right: pod 10 ¡3 (N 5). Bottom left: pod 10 ¡4 (N 7). Bottom right: pod 10 ¡6 (N 11).

Figure 4 . 15 -

 415 Figure 4.15 -Test case (d) : Relative H 1 -error e k pµq for some time nodes k K tr and M 676. Top left: pod 10 ¡2 (N 3). Top right: pod 10 ¡3 (N 5). Bottom left: pod 10 ¡4 (N 7). Bottom right: pod 5.10 ¡6 (N 11).

Figure 4 . 16 -

 416 Figure 4.16 -Test case (d) : Stability constant β N,M . On the right panel, the values of N are respectively 2, 3, 5, 7, 11 for the values of pod in decreasing order.

Figure 4 . 18 -

 418 Figure 4.18 -Test case (e) : Relative H 1 -error e k pµq for the state estimate as a function of the time nodes. The various curves correspond to the different values of µ.

Figure 4 . 19 -

 419 Figure 4.19 -Test case (e) : Absolute H 1 -norms of the contributions z ¦ N init ,M init and

Figure 4 .

 4 Figure 4.20 -Test case (e) : Relative H 1 -error e k pµq for the state estimate as a function of the time nodes during the online stage. The various curves correspond to the different values of µ. Left: for all µ P tr . Right: for all µ P verif .

Figure 4 . 21 -

 421 Figure 4.21 -Test case (e) : Absolute H 1 -norms of the contributions z ¦ N,M and η ¦ N,M as a function of the time nodes. The various curves correspond to the different values of µ.

  . . . , N u avec N : maxt1 ¤ n ¤ D | σ n ¥ POD u. Il est intéressant de remarquer que, parmi tous les sous-espaces Z N de dimension N dans R N , l'espace de dimension N engendré par les vecteurs pθ n q 1¤n¤N est celui qui minimise la quantité En notant Z N vecttθ 1 , . . . , θ N u, nous avons de manière équivalente, que Z N De plus, nous avons }v ¡ Π Z N v} X ¤ σ N 1 }v} X , pour tout v S, où Π Z N désigne le projecteur X-orthogonal sur Z N . Dans la pratique, il est possible d'éviter le calcul de la matrice pC N q 1

	minimise la quantité
	R ŗ1	inf zZ N }v r ¡ z} 2 X	(1.19)

.17) À partir de

(2.31)

, nous déduisons que T ψn σ n ξ n et T T ξ n σ n ψn pour tout n t1, . . . , Du. Enfin, les vecteurs recherchés sont donnés par θ n : pC

N q ¡ 1 tout n t1, R ŗ1 inf zZ N pv r ¡ zq T C N pv r ¡ zq (1.18)

2 et de son inverse en considérant la matrice de plus petite dimension T T T S T C N S R R¢R , dite matrice d'auto-corrélation. Comme T T T ψn σ n T T ξ n σ 2 n ψn , la résolution du problème aux valeurs propres associé à T T T donne les vecteurs ψn et leurs valeurs propres associées σ 2 n . Ainsi, les vecteurs

  est une méthode itérative dans laquelle, à chaque itération, une nouvelle fonction de base réduite est ajoutée et la précision globale de la base est en général améliorée. Pour cela, il faut calculer une solution du problème haute-fidélité par itération et un total de N solutions haute-N dans le problème haute-fidélité par l'espace réduit X N . Nous détaillerons le calcul de ∆ N : P Ñ R en Section 1.2.4 ci-dessous. L'estimation d'erreur permet non seulement de quantifier et de contrôler l'erreur commise par la RBM, mais aussi d'explorer l'espace des paramètres de manière optimale. En général, l'estimateur d'erreur ∆ N vérifie l'inégalité }u N pµq ¡ u N pµq} X ¤ ∆ N pµq, dµ P. P ne pouvant en général pas être exploré de manière exhaustive, on explore uniquement un sous-ensemble P tr P qu'on appelle ensemble d'apprentissage.L'initialisation de la base réduite se fait par un premier calcul haute-fidélité pour un paramètre quelconque µ 1 permettant de définir X 1 vecttu N pµ 1 qu. Ensuite, pour chacune des itérations n ¥ 1 de cet algorithme, étant donnée une base tu

	(1.22)
	L'ensemble des solutions réduites pour µ

fidélité pour générer l'espace réduit de dimension N . On notera au passage que la POD nécessite le calcul de R ¥ N solutions haute-fidélité. L'algorithme glouton nécessite un estimateur d'erreur ∆ N pµq qui prédit l'erreur due à la réduction de dimension du problème, c'est-à-dire qui fournit une estimation de l'erreur induite en remplaçant X N pµ 1 q, . . . , u N pµ n qu, la prochaine fonction de base à sélectionner est celle pour laquelle nous approchons le moins bien la sortie du modèle, i.e. celle qui maximise l'estimateur sur l'erreur de réduction par rapport à l'ensemble des paramètres d'apprentissage P tr . Ainsi, l'algorithme glouton sélectionne un paramètre

  HF in . Therefore, the selection criterion (2.26) exploits the knowledge of the HF trajectory for all the parameter values in P HF in , and otherwise uses a RB trajectory. This is an important difference with respect to the standard offline stage. There are now two possibilities: (i) either µ m is already in P HF in ; then, no new HF trajectory is computed and we set P HF

	26)
	In (2.26), ūk I

µ I is defined as in (2.22) using the set P out : P HF in (line 8); (ii) or µ m is not in P HF in ; then we compute a new HF trajectory for the parameter µ m and we set P HF out : P HF in tµ m u (line 5). Our numerical experiments reported in Section 2.6 below will show that at many iterations of PREIM, the pair pµ m , k m q selected in (2.26) differs from the previously selected pair by the time index and not by the parameter value; this means that for many PREIM iterations, no additional HF computation is performed. In case of nonuniqueness of the maximizer in (2.26), one selects, if possible, a trajectory for which the parameter is not already in the set P HF in in order to trigger a computation of a new HF trajectory.

An additional feature of PREIM is that, whenever a new HF trajectory is actually computed, one can either confirm or update the selected pair pµ m , k m q using the following HF-based re-selection criterion (see line 6):

Table 2 .

 2 1 -Test case (a): Evolution of the standard EIM error. m is the rank of the EIM approximation.During the standard offline stage, we perform P 20 HF computations. Knowing that K 50, the set S (Algorithm 2.2, line 1) contains 1020 fields, each consisting of N 1438 nodal values. Applying the POD in a progressive manner (see Algorithm 2.3 with the parameters enumerated using increasing values) based on the H 1 -norm and a truncation threshold pod 10 ¡3 , we obtain N 6 RB functions. Afterwards, we perform the standard EIM algorithm whose convergence is reported in Table2.1 for selected values of the rank of the EIM approximation. For eim 5 ¤ 10 ¡2 , the final rank of the EIM approximation is M 8, whereas for eim 5 ¤ 10 ¡3 , the final rank of the EIM approximation is M 15.

	m	1 2 3 4 5 6 7 8 9 10 11 12 13
		μ 1 20 20 20 20 20 20 20 16 20 20 18 20
	PREIM	µ 1 20 20 20 20 18 20 20 16 20 20 18 20
		k 50 45 48 50 43 42 39 46 50 49 33 50 47
	Table 2.2 -Test case (a): Selected parameters and time nodes in PREIM. The
	gray cells correspond to a new parameter selection and, therefore, to a new HF
	computation.

Table 2 .

 2 3 -Test case (a): Selected parameters and time nodes in U-SER and PREIM-NR. The gray cells correspond to a new parameter selection and, therefore, to a new HF computation.

		1 2 3 4 5 6 7 8 9 10 11 12 13
	U-SER	µ 1 20 20 20 16 20 19 20 20 19 17 20 19 k 50 49 50 46 42 49 44 39 50 49 48 47 50
	PREIM-NR	µ 1 20 20 20 20 16 20 20 20 20 20 17 19 k 50 47 50 46 42 49 48 46 39 50 45 50 50

Table 2 .

 2 ¡2 .

		RB-EIM PREIM U-SER
	HF computations	99% 20.0% 25.0%
	greedy runtime	1%	1.5% 2.3%
	Total runtime	100% 21.5% 26.3%

4 -Test case (a): Runtime measurements.

Table 2 .

 2 5 -Test case (b): Size of the reduced basis in the standard algorithm with pod 5 ¤ 10 ¡2 .

	1 2 3 8 20 23 24 26 32 33 36 37 39 40
	RB size 3 4 5 6 7 8 9 10 11 12 13 14 15 15
	m }r m } V pΩ tr q 1.6 1.3 9.7E¡1 4.7E¡1 1.7E¡1 1.2E¡1 8.0E¡2 9.1E¡3 4.6E¡3 9.4E¡4 2 10 13 20 30 36 37 79 96 144

Table 2 .

 2

6 -Test case (b): Evolution of the standard EIM error. m is the rank of the EIM approximation and }r m } V pΩ tr q is the residual norm in (2.20).

Table 2 .

 2 7 -Test case (b): Selected parameters and time nodes in PREIM for eim 10 ¡1 (top) and eim 10 ¡3 (bottom). The gray cells correspond to a new parameter selection and, therefore, to a new HF computation.

	m	1 2 9 17 18 20
	RB size 5 6 6 7 9 9

Table 2 .

 2 8 -Test case (b): Size of RB generated within PREIM for pod 5 ¤ 10 ¡2 ; for eim 10 ¡1 , one stops at m 9, and for eim 10 ¡3 , one stops at m 20.

Table 2 . 9

 29 Applying the POD in a progressive manner based on the H 1 -norm and a relative truncation threshold pod 10 ¡3 defined as suggested in Remark 2.3, we obtain N 4 RB functions. Figure2.12 shows three POD modes. Afterwards, we perform the standard EIM algorithm whose convergence is reported in Table2.9.

	m	1 2 3 4 5 6 7 8 9 10
		μ 20 20 10 20 20 20 8 20 20 13
	PREIM	k 30 14 30 21 5 26 22 3 9 30
		µ 20 11 10 20 19 9 8 12 13 7
		k 31 31 31 21 31 31 31 31 31 31

-Test case (c): Evolution of the standard EIM error. m is the rank of the EIM approximation and }r m } V pΩ tr q is the residual norm in

(2.20)

.

stage, we perform P 20 HF computations. Knowing that K 31, the set S (Algorithm 2.2, line 1) contains 620 fields, each consisting of N 46, 018 nodal values.

Table 2 .

 2 10 -Test case (c): Selected parameters and time nodes in PREIM for pod 10 ¡3 and eim 5 ¤ 10 ¡3 . The gray cells correspond to a new parameter selection.

Table 2 .

 2 10 shows that the first selection has been discarded. The right panel of Figure 2.13 illustrates the space-time errors (measured in the relative 2 pI tr ; H 1 pΩ tr qq-

norm) on the trajectories produced by the standard RB-EIM and the PREIM procedures for the whole parameter range; here the tolerances are set to pod 10 ¡3 and eim 5 ¤ 10 ¡3 . We observe that for higher parameter values, PREIM delivers somewhat more accurate results, whereas the conclusion is reversed for lower parameter

Table 2 .

 2 .10. Still, for such tight tolerances and with quasi-identical output errors (cf. right panel of Figure 2.14), PREIM makes less than half of the HF computations incurred in the standard RB-EIM. 11 -Test case (c): Runtime measurements with pod 10 ¡4 . Left: eim 5 ¤ 10 ¡1 . Right: eim 5 ¤ 10 ¡2 .

	eim 5.10 ¡1 RB-EIM PREIM	eim 5.10 ¡2 RB-EIM PREIM
	HF computations	99.8% 10.0%	HF computations	99.6% 20.0%
	greedy runtime	0.2%	0.6%	greedy runtime	0.4%	2.4%
	Total runtime	100% 10.6%	Total runtime	100% 22.4%

  is the Young modulus, ν is the Poisson coefficient and I is the identity tensor in R d¢d . For simplicity, we have supposed that E and ν are parameterindependent. At this stage, we define the bilinear form a : P ¢ Vpµq ¢ Vpµq Ñ R

	where E introduced in (3.1) as	»
		apµ; v, wq	σpvq : εpwq,	(3.18)
		Ωpµq
	and the linear form f : P ¢ Vpµq Ñ R introduced in (3.1) as
		f pµ; wq	»	pµqw.	(3.19)
			Ωpµq
	σpvq	Eν p1 νqp1 2νq trpεpvqqI	E p1 νq εpvq,	(3.17)

  The constraint(3.26) is equivalent to the physical non-interpenetration conditionΩ 1 pµ, u 1 pµqq Ω 2 pµ, u 2 pµqq Υ c pµ, upµqq, pµ, u i pµqq : pId u i pµqqpΩ i pµqq, di t1, 2u,(3.31)andΥ c pµ, upµqq : Υ c 1 pµ, u 1 pµqq Υ c 2 pµ, u 2 pµqq,

	for all z Γ c 1 pµq. Hence, (3.26) can be recast as
	kpµ, upµq; upµqq ¤ gpµ, upµqq,	(3.29)
	leading to the same inequality constraint as in (3.1).
	Lemma 3.1. (3.30)
	where
	Ω i
		(3.28)

.26) At this stage, we can define the displacement map k and the gap map g in (3.1) as kpµ, w; vqpzq ¡ pv 1 pzq ¡ pv 2 ¥ ρpµ, wqq pzqq ¤ r n 2 pµ, wqpzq, (3.27) and gpµ, wqpzq ¡ pρpµ, wqpzq ¡ zq ¤ r n 2 pµ, wqpzq, (3.32) Proof. The proof is postponed to Section 3.7.3. Remark 3.2 (Physical interpretation). Incidentally, (3.30) means that the intersection of the two deformed solids Ω 1 pµ, u 1 pµqq and Ω 2 pµ, u 2 pµqq is necessarily a subset of their contact boundaries. Remark 3.3 (Symmetry). Note that the indices 1 and 2 play symmetric roles in (3.30).

  .40) with real numbers p u n pµq for all n t1, . . . , N u and non-negative real numbers p λ n pµq for all n t1, . . . , Ru. Let us introduce the component vectors p upµq : pp u n pµqq 1¤n¤N R N and p λpµq : p p λ n pµqq 1¤n¤R R R , for all µ P. The RB formulation of (3.9) reads: Find pp upµq, p λpµqq R N ¢ R R such that Apµq R N ¢N and p

		Kpµ, p vq R R¢N such that
	pp upµq, p λpµqq arg minmax p vR N , p ηR R	1 2

p v T p Apµqp v ¡ p v T p f pµq p η T p Kpµ, p vqp v ¡p gpµ, p vq ¨, (3.41) with the matrices p p Apµq pn apµ; q

  B κ , tC κ j u 1¤j¤M k and D γ . 1: Assemble the vector p f pµq and the matrix p Solve the reduced saddle-point problem (3.51) to obtain p upµq and p λpµq Output : p

	Apµq using (3.45)
	2: Compute p κpµ, p vq and p γpµ, p vq using (3.50)

3: Compute D κ pµq using p κpµ, p vq and (3.52) 4:

Table 3 .

 3 1 -Primal basis dimension N as a function of the truncation threshold pod . 10 ¡2 10 ¡2 5 ¤ 10 ¡3 10 ¡3

	NMF	R	4	10	19	20
	Cone-projected greedy R	2	5	6	13

du 5 ¤

Table 3 .

 3 2 -Dual basis dimension R as a function of the truncation threshold du .

  ) which cannot hold true since z Ω 2 pµ, u 2 pµqq by assumption. Actually, the implication (3.76) means thathz b Υ c 2 pµ, u 2 pµqq : 5 z ¡ z b ¤ r n 2 pµ, upµqq ¥ pψ 1 pµ, u 1 pµqqq ¡1 pzq pµ, u 1 pµqq Ω 2 pµ, u 2 pµqq ùñ hz b Υ c 1 pµ, u 1 pµqq : z b Ω 2 pµ, u 2 pµqq. (3.77)The implication (3.77) means that a non-empty intersection of both solids necessarily induces that at least a boundary point of each solid is in the interior of the other.From (3.77), we infer that the constraint (3.26) is violated. Therefore, Ω 1 pµ, u 1 pµqq Ω 2 pµ, u 2 pµqq r, which concludes the proof.

	¨	0
	z b ϑpµ, upµqqpzqq,
	which is in contradiction with the definition of the contact mapping ϑpµ, upµqq. If	z
	is an interior point of the first solid, we have
	z Ω 1

  R M ¢N , (4.31) Proof. The system (4.30) is a saddle-point problem with a symmetric positive definite matrix A. Therefore, it has a unique solution if and only if the matrix B is injective, i.e., if and only if kerpBq t0u. Using the definition of B in (4.31), we have kerpBq $ t0u ðñ hz Z N zt0u : dq U M : pz, qq 0, ðñ Z N U u M $ t0u. t0u and this statement is equivalent to β N,M ¡ 0. Finally, we readily verify that β N,M ¤ 1 owing to the Cauchy-Schwarz inequality. Remark 4.2 (Stability constant). In terms of geometry, the stability constant β N,M

					(4.34)
	Thus, (4.30) is well-posed if and only if Z N U u
	and the vector of observations
		obs M obs m pu true q ¨1¤m¤M R M .	(4.32)
	Proposition 4.3 (Well-posedness). The PBDW statement (4.30) has a unique solu-tion if and only if Z N U u M t0u. Equivalently, under this assumption, the stability
	constant	β N,M : inf wZ N	sup vU M	pw, vq }w} }v} p0, 1s.	(4.33)

M

  , vq pη N,M , vq pz N,M , vq, dv U,

	.58)
	subject to
	pu N,M (4.59a)

  The procedure S-greedy is described in Algorithm 4.3 below. Altogether, Algorithm 4.2 Modified offline stage of the time-dependent PBDW Input : P tr , K tr , S, pod , init pod and β min . Q init : an initial set of Riesz representations of observations. Compute the matrices A init and B init using Z init N init and U init M init . Estimate the state u k,¦ pµq for all pµ, kq P tr ¢ K tr .

1: Compute Z init N init : POD-greedypS, init pod q. 2: Set U init M init : SpantQ init u.

3:

4: 5: Compute Z N : POD-greedyptu k,¦ pµqu µP tr ,kK tr , pod q. 6: Compute U M : S-GreedypP tr , K tr , N , tZ n u N n1 , β min , Q init ).

3 :

 3 Compute the stability constant β 1,1 using Z 1 and U 1 .

4: Set m : 2. 5: while β N,m¡1 β min or m M do 6:

  Linear heat equation: We apply a homogeneous Neumann boundary condition on fΩ 0 and a non-homogeneous Neumann boundary condition on fΩ n , Note that fΩ 0 consists of the upper and right sides of the plate and fΩ n consists of its lower and left sides, so that fΩ fΩ 0 fΩ n . Dpµq fu fn σεpu 4 ¡ u 4 r q, on I ¢ fΩ,(4.72)where u r 303.15K (30 o C) is an enclosure temperature, σ 5.67 ¢ 10 ¡8 W.m ¡2 .K ¡4 is the Stefan-Boltzmann constant and ε 3.10 ¡3 is the emis-

	6 9 9 9 8 9 9 9 7	fupµq ft ¡ ∇ ¤ pDpµq∇upµqq 0, in I ¢ Ω, Boundary conditions, on I ¢ fΩ, upµqpt 0, ¤q u 0 , in Ω,	(4.68)
	where u 0 293.15K (20 o C). We will supplement (4.68) with two types of boundary
	conditions:
	1. i.e.,
		6 9 8 9 7	¡Dpµq fupµq fn 0, fn φ e , on I ¢ fΩ n , on I ¢ fΩ 0 , ¡Dpµq fupµq	(4.69)
	with φ e 3K¤m¤s ¡1 and
		fΩ 0 p¡2, 2q ¢ t2u t2u ¢ p¡2, 2q, fΩ n p¡2, 2q ¢ t¡2u t¡2u ¢ p¡2, 2q.	(4.70) (4.71)
	Thus, the resulting problem (4.68)¡(4.69) is linear. sivity. The Stefan-Boltzmann boundary condition is nonlinear and so is the
	resulting problem (4.68)¡(4.72).
	In what follows, the background spaces Z

2. Nonlinear heat equation: We apply Stefan-Boltzmann boundary conditions on fΩ, i.e., ¡ N will be generated by solving either the linear PDE (4.68)¡(4.69) or the nonlinear PDE (4.68)¡(4.72) with a uniform diffusivity function Dpµq such that Dpµqpxq D uni pµqpxq : µ1 Ω pxq, dx Ω. (4.73)

  int 1 and define, for each µ P, the diffusivity function Dpµq as Dpµqpxq D syn pµqpxq : µD int 1 Ωext pxq D int 1 Ω int pxq, dx Ω, int Ω ext and Ω int Ω ext r. The synthetic true solutions are then defined as the solutions of (4.68) for all µ P, with either the linear boundary condition (4.69) or the nonlinear boundary condition (4.72).

		(4.74)
	where
	Ω int p¡1, 1q 2 , and Ω ext p¡2, 2q 2 zp¡1, 1q 2 ,	(4.75)
	so that Ω Ω

Cette thèse se consacre à une classe de problèmes thermo-mécaniques, qui est en particulier celle rencontrée dans les études de robinetterie au sein d'EDF R&D. Dans cette section, nous présentons brièvement le cadre de ces études.1.1.1 Études de robinetterie et enjeux industrielsEn tant qu'équipements importants pour la sûreté et la disponibilité des centrales nucléaires, certains matériels de robinetterie doivent être qualifiés aux conditions normales, incidentelles et accidentelles de fonctionnement. La qualification a pour objectif de vérifier la capacité du matériel à assurer sa fonction dans des conditions de fonctionnement et pour une durée de vie spécifiées.Dans le cadre d'un partenariat entre le département Mécanique des Matériaux et Composants (MMC) d'EDF R&D et un fabricant de robinets, un essai sur un robinet à maintenance allégée (RAMA) fortement instrumenté a été mis en place sur la boucle de cyclage thermique (CYTHER) (cf. Figure1.1). L'expérience consiste en une succession de chocs thermiques alternés par passage d'un fluide sous pression. Il résulte de cet essai une quantité importante de données expérimentales, telles que des températures mesurées par des thermocouples, des efforts de serrage dans des goujons, des déformations résiduelles et bien d'autres mesures. Une étude récente

1.1. Contexte industriel

ξ n pour

pµ, v 1 q to a potential contact point on Υ c

pµ, v 2 q.

The ideas introduced in this chapter have been devised in collaboration with A. T. Patera during two two-week visits of the author of the manuscript to the Massachusetts Institute of Technology.

Offline/online procedure: Since several realizations u true pωq of the true state are considered, an offline/online procedure can be employed. During the offline stage, one precomputes the RB functions pζ n q 1¤n¤N and the Riesz respresenters pq m q 1¤m¤M leading to the matrices A R M ¢M and B R N ¢M once and for all.

Then, during the online stage, for each new set of observations corresponding to a new realization of the true state u true pωq, all that remains to be performed is to form the vector of observations obs M pωq and to retrieve the deduced background estimate z ¦ N,M pωq and the update estimate η ¦ N,M pωq by solving the pN M q-dimensional linear problem (4.30). The PBDW state estimate u ¦ N,M pωq is then computed using (4.37).

A priori error analysis

Here, we state an important result proved in [START_REF] Maday | A parameterized-background data-weak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF] and related to the a priori error analysis of the PBDW statement.

Hence, the error on the update estimate depends on the quality of approximation of the true state u true in Z N , i.e., inf zZ N }u true ¡ z}. The same comment is valid for the error on the deduced background estimate and the state estimate. Moreover, since

which shows that the three errors from Proposition 4.4 also depend on how well the subspace Z u N U M cures the lack of information in the background space Z N .

Remark 4.4 (Choice of spaces). The subspaces Z N and U M must be chosen carefully. In fact, we want a small angle between the spaces Z N and U M in order to increase the stability constant, but we need some overlap between the spaces Z u N and U M to improve the approximation capacity of Z u N U M .

Background space construction via state estimation

We now illustrate the performances of Algorithm 4.2 for the following linear imperfect case:

• Test case (e): We consider a simulation duration T 4s and a time step ∆t 0.1s. Test truths are synthetized with an internal diffusivity D int 0.2.

As opposed to the previous section, we choose a non-parametric bk model based on an HF computation for µ 0.5. The resulting unique trajectory is then reduced using a POD algorithm, which is equivalent to a POD-greedy for a single trajectory (cf. line 1 of Algorithm 4.2). For a tolerance value init pod 10 ¡2 , we obtain a background space Z init N init composed of N init 34 modes. Figure 4.17 shows the singular values that are retained. As regards observations, the initial set Q init consists of N init ,M init . One can notice that the latter is non-negligible compared to the former. Once the first part of the modified offline stage has been performed, we use the resulting state estimates in order to build the modified background space (cf. line 5 of Algorithm 4.2). For a tolerance pod 5.10 ¡2 , the POD-greedy algorithm selects four modes. Then, we build the observable space U M using M 121 uniformly distributed sensors (the optimal choice can be made using the S-Greedy algorithm, see Algorithm 4.3). Figure 4.20 displays the errors for the verification set P verif t0, . . . , 19u. The state estimation relative H 1 -error e k pµq remains comparable to that of the five parameters used for the offline construction. Regarding the online observations, we highlight that the online results are achieved using only M 8%M init . Finally, Figure 4.21 shows the absolute H 1 -norms of