Mireille Bossy

Olivier Glass

Philippe Moireau

Jacques Sainte-Marie

Keywords: 4, 2

Working on this thesis was equivalent to living in a roller coaster during almost 4 years. Exciting times, I must admit, in which I visited amazing places as Svalbard or Hong Kong while enjoying my work, and I also had a life in Paris I've never dreamed of. Even so, in an almost unnoticeable way, I began to forget why I chose to be a mathematician, since that willingness to learn slowly transformed into a distaste for what I was doing. At some point, near to what was supposed to be the end, I was only worried I wouldn't finish this manuscript in time, but then I simply faced reality and said : if Brexit was postponed several times, why can't I do the same ? Still, sometimes it seemed closer than my long-awaited thesis defense, even Notre-Dame got burned down in between, and I was also afraid that GRRM, a professional procrastinator, could publish his last book at any moment. By the end this PhD was a mix between uneasiness, impostor syndrome, total relief and more than a bit of hair turning white, but still, we made it.

You may be wondering why I took the liberty of writing all this. It's just that we inevitable soften our memories, and among many other things, writing is our small rebellion against time. In any case, it's time to come back to the usual purpose of these lines. I would like to thank all of whom were part of this process, because I couldn't have finished this thesis alone. Sure enough, an enormous cliché.

I would also like to extend my gratitude to all members of CEREMADE. Even if the architecture of the university building wasn't appealing, as military complexes often are, I really enjoyed my stay there. Many thanks to my office mates Nadia, Fang, Long and the aircon machine, our invaluable summer companion. Working together was like having a support group, and I definitely learned a lot from our interesting conversations about each other's respective cultures ! Thanks also to Jorge and the rest of the ephemeral pisco sour gang. Speaking of offices, being four floors below the rest of the lab was an excellent way to have a healthier lifestyle, based on climbing stairs, drinking coffee afterwards and hopefully, turned it into theorems; which also reminds me it wasn't until last year that I realized coffee was free for us. Then I guess it's not a surprise to anyone that I wasn't aware of the usual -and most of the time, nightmarish-paperwork, so I gratefully acknowledge the assistance of all administrative people, especially Cesar Faivre and Isabelle Bellier, Linda Mammoudi from the École doctorale, as well as Vincent Rivoirard and Jacques Féjoz on the directive side. ** * Finalmente, en la vorágine de ideas, fórmulas y técnicismos propios de este doctorado, es fácil creer que sus resultados son un fruto más de la "meritocracia" y omitir por completo el rol que jugaron mi familia y amigos en toda esta historia. Me gustaría agradecerles, aunque sea brevemente, en las próximas líneas.

A Milena, por darme esa claridad tan propia de ella que siempre termina poniéndome en perspectiva, convencerme de aceptar una beca inesperada y ser parte fundamental de estos años intensos y viajeros. Que nuestro actual paso por Francia y el futuro nos deparen más cosas buenas! A mi mamá y hermanos, por su generosidad al apoyar incondicionalmente cada una de mis decisiones, muchas veces precipitadas, confiando que estaba en lo correcto y aún sin entender porqué decidí dedicarme a esto. Gracias a los amigos matemáticos que decidieron repartirse por el mundo y de paso compartir un nuevo capítulo de esta estafa piramidal: Mónica y sus precisos consejos; Sandra y Jean-marc, por adoptarme cual gato en su hogar; Johan, Dani, señora Nicole y Abelino; Pedro y Pita, por las escapadas al château; y la capitana Anne, siempre más visionaria que todos nosotros. También a quienes de, a pesar de lo poco que nos vemos, nuestra amistad sigue intacta: Juan Carlos, Felipe, Karla, Gabriel, Pepo, por mencionar sólo a algunos. Y a los amigos con los cuales he compartido la experiencia de vivir acá, especialmente Víctor y Nico, Carlos y Garrido tiempo después. Las anécdotas del programa capital humano avanzado en París son interminables.

Les algorithmes de parallélisation en temps peuvent être très utiles lorsqu'il s'agit de intervalles de temps plus longs, comme c'est le cas des méthodes d'AD séquentielles, où les informations peuvent arriver sans interruption. Une question naturelle se pose alors : pouvons-nous combiner les deux procédures ? Nous commençons à répondre à cette question au Chapitre 2, en étudiant l'observateur de Luenberger et son couplage avec l'algorithme Pararéel.

Optimisation d'une bathymétrie

Bien que la bathymétrie puisse être mal connue dans de nombreuses situations, les modèles de propagation des ondes dépendent fortement de ce paramètre pour capturer le comportement de l'écoulement, ce qui souligne l'importance d'étudier les problèmes inverses concernant sa reconstruction à partir des données observées en surface libre.

Ces types de problèmes sont habituellement résolus en discrétisant simplement les équations qui les régissent ou à l'aide de méthodes d'AD séquentielles. Une alternative consiste à considérer la bathymétrie comme variable de contrôle d'un problème d'optimisation sous contrainte d'EDP, une approche utilisée en ingénierie côtière en raison des contraintes mécaniques associées aux structures des bâtiments et à leur interaction avec les vagues de mer. Cependant, sa résolution repose principalement sur des analyses de sensibilité, des méthodes numériques ou des simplifications du modèle qui conduisent à des solutions explicites et ensuite, les questions concernant un espace de contrôle approprié, la continuité de la fonction contrôle-état, la régularité des solutions et le caractère bien posé du problème ne sont en général pas abordées. Récemment, Dalphin et Barros [START_REF] Dalphin | Shape optimization of a moving bottom underwater generating solitary waves ruled by a forced KdV equation[END_REF] ont réalisé cette analyse théorique pour modéliser un générateur de vagues.

Au Chapitre 3 nous essayons de répondre aux questions ci-dessus lorsque nous envisageons une reformulation de l'équation de Helmholtz pour la modélisation de la propagation des ondes et une fonctionelle générale de coût qui peut être identifiée, par exemple, avec le décalage entre la solution d'onde prévue et celle observée.

Analyse mathématique de la méthode de l'élément de pale

Hydrotube Énergie, entreprise dédiée à la conception nautique et aux systèmes électroniques embarqués, a testé en 2015 le prototype d'une turbine hydraulique flottante sur la Garonne, à Bordeaux. Mais l'évolution vers une production à l'échelle industrielle nécessite le développement d'un logiciel numérique pour simuler un dispositif optimisé, qui est le but de son partenariat avec l'équipe ANGE (INRIA), un groupe de recherche en modélisation, analyse et simulation des écoulements géophysiques.

Parmi les problèmes à résoudre, nous étudions l'optimisation de l'efficacité de la turbine via la méthode de l'élément de pale. Présenté par Glauert [START_REF] Glauert | The Elements of Aerofoil and Airscrew Theory[END_REF], cette méthode est largement utilisée pour déterminer l'efficacité et donc les paramètres de conception d'une pale, en fonction de ses caractéristiques géométriques et mécaniques et du courant auquel elle est exposée. Elle résulte de la combinaison de deux méthodes indépendantes qui traitent le système fluide/turbine d'un point de vue macroscopique et local, ce qui donne un ensemble d'équations qui sont résolues en appliquant un algorithme itératif. Dans certains cas, un facteur de correction est nécessaire pour assurer l'existence des solutions du modèle.

Bien qu'ancienne, cette méthode est encore utilisée en raison de sa relative simplicité par rapport à la complexité du phénomène hydrodynamique développé dans le système fluide/turbine. Néanmoins, l'existence de solutions et la convergence de l'algorithme n'ont jamais été analysées d'un point de vue mathématique, ce qui est le but du Chapitre 4.

Contributions de cette thèse

Ce travail est divisé en trois chapitres indépendants portant sur les sujets susmentionnés, dans lequel j'ai eu le plaisir de collaborer avec Julien Salomon (INRIA Paris), Felix Kwok (Hong Kong Baptist University), Pierre-Henri Cocquet (Université de Pau et des Pays de l'Adour) et Jérémy Ledoux (Hydrotube Énergie). Dans ce qui suit, nous résumons nos principales contributions à chacun d'eux.

Parallélisation en temps de l'assimilation de données

Nous commençons au Chapitre 2 en proposant une procédure pour coupler des méthodes d'assimilation de données séquentielles avec des algorithmes de parallélisation en temps, qui consiste à diviser l'intervalle de temps non borné en sous-intervalles de même longueur (fenêtres) et à appliquer ensuite, selon un ordre séquentiel, le solveur parallèle en temps sur chacun d'eux. En considérant l'observateur de Luenberger comme méthode d'assimilation, nous fournissons un critère d'arrêt qui préserve son taux exponentiel de convergence, ce qui donne une estimation a posteriori de la précision du solveur. Afin d'aller plus loin, nous avons étudié plus précisément le cas d'une parallélisation en temps par l'algorithme Pararéel comme solveur parallèle en temps. Cela nous permet de concevoir un algorithme alternatif qui fournit une estimation a priori du nombre d'itérations nécessaires sur chaque fenêtre, ce qui nous permet également d'étudier l'efficacité théorique de l'ensemble de la procédure. Ces résultats sont basés sur une nouvelle estimation de convergence que nous dérivons pour Pararéel lorsque le solveur grossier est contractif.

Optimisation d'une bathymétrie

Nous passons à l'étude au Chapitre 3 de la détermination d'une bathymétrie à partir d'un problème d'optimisation, où une reformulation de l'équation de Helmholtz agit comme une contrainte. Même si cette équation est limitée à la description d'ondes de faible amplitude, elle est souvent utilisée en ingénierie en raison de sa simplicité, qui conduit à des solutions explicites lorsqu'une bathymétrie plate est considérée. En levant cette hypothèse, on obtient une formulation différente dans laquelle cette variable joue le rôle d'un diffuseur.

Sous des hypothèses appropriées sur la fonctionnelle de coût et l'ensemble admissible des bathymétries, nous sommes en mesure de prouver la continuité de la fonction contrôle-état et l'existence d'une solution optimale, en plus de la continuité et et le caractère borné de la vague résultante. Le problème de l'optimisation discrète est également abordé, en étudiant la convergence vers la solution optimale discrète ainsi que la convergence d'une approximation par éléments finis.

Analyse mathématique de la méthode de l'élément de pale

Nous abordons enfin au Chapitre 4 l'existence de solutions et la convergence des procédures de résolution pour la méthode de l'élément de pale. Le point clé de notre travail consiste à montrer que la décomposition proposée par Glauert peut être utilisée pour reformuler son ensemble original d'équations en une seule expression contenant deux termes bien distincts : un terme universel, indépendant de la turbine considérée et associé aux aspects macroscopiques du modèle ; et un terme expérimental, concernant les caractéristiques de la turbine et associé à la partie locale du modèle.

L'avantage de notre approche est qu'elle identifie explicitement les hypothèses sur les paramètres de la turbine qui garantissent l'existence d'une solution. De plus, elle nous aide aussi à présenter des critères de convergence pour différents algorithmes de résolution. Le modèle de Glauert est également utilisé pour optimiser la géométrie des pales, en ce sens qu'il maximise le rendement de la turbine. Nous rappelons les détails de la procédure de conception habituelle et discutons brièvement du cas où une correction du modèle est introduite.

État de l'art

Cette section est consacrée à un aperçu des différents sujets abordés dans cette thèse. Nous commençons par deux sujets étudiés dans le Chapitre 2 : les méthodes d'assimilation des données, en continu et en discret suivi d'un résumé des méthodes de parallélisation en espace et en temps, en particulier l'algorithme Pararéel. Nous passons ensuite à la discussion des éléments clés du Chapitre 3, l'estimation de la bathymétrie et modélisation des ondes, en mettant l'accent sur les principales hypothèses nécessaires pour dériver différents modèles de propagation des ondes. Enfin, le Chapitre 4 donne une brève description de la méthode de l'élément de pale.

Assimilation des données (AD)

Les modèles mathématiques sont largement utilisés pour décrire des systèmes complexes. Ils reposent sur des approximations et des simplifications d'un phénomène réel, ce qui définit en dernière instance son champ d'application. Après une validation appropriée, le but est de les utiliser pour décrire ou même prévoir une situation réelle, comme la surveillance sismique, la prévision de la grippe saisonnière ou l'estimation de l'état de charge d'une batterie [START_REF] Hu | Estimation of state of charge of a lithium-ion battery pack for electric vehicles using an adaptive Luenberger observer[END_REF]. Pour ce faire, nous devons intégrer des données réelles dans notre cadre. Les différentes techniques qui combinent les modèles mathématiques avec les observations disponibles pour améliorer la connaissance d'un système sont connues sous le nom d'Assimilation des données (AD). Parmi ces approches, nous rappelons ici les méthodes séquentielles et variationnelles.

Méthodes séquentielles

Supposons que pour une raison quelconque, par exemple des contraintes physiques ou budgétaires, nous n'ayons accès qu'à une information partielle sur l'état du système. Par exemple, dans certains cas, la condition initiale n'est pas connue avec précision, comme c'est le cas en science du climat [START_REF] Werndl | Initial-condition dependence and initial-condition uncertainty in climate science[END_REF]. Nous pouvons traiter ce manque d'information et cette incertitude en construisant un nouveau système qui utilise les observations disponibles pour se rapprocher de l'état réel. Dans un contexte déterministe, ce dispositif est appelé un observateur. +∞) x(0) = x 0 , y(t) = Cx(t), [START_REF] Ackermann | Der entwurf linearer regelungssysteme im zustandsraum[END_REF] où x ∈ R m est le vecteur d'état réel, y ∈ R q représente les observations (avec q < m, m, q ∈ N * ), u est une entrée et x(0) = x 0 une condition initiale inconnue. Les matrices A, B et C sont connues et leurs dimensions sont cohérentes. L'observateur de Luenberger [START_REF] Luenberger | Introduction to Dynamic Systems: Theory, Models, and Applications[END_REF] imite le modèle précédent, mais il inclut un terme supplémentaire dans la dynamique qui mesure l'écart entre les observations et ses propres prédictions. Il produit une estimation d'état x satisfaisant +∞) x(0) = x0 , ŷ(t) = C x(t), avec x0 une condition initiale arbitraire. Tant que le modèle original [START_REF] Ackermann | Der entwurf linearer regelungssysteme im zustandsraum[END_REF] est observable (c-à-d que l'état initial x(0) peut être déterminé uniquement à partir des observations dans [0, T ], pour tout T ), l'erreur résultante x(t) -x(t) peut être mise à zéro à un taux exponentiel en choisissant correctement la matrice de gain L, puis l'état réel est récupéré asymptotiquement.

Nous supposons que le système est régi par

       ẋ(t) = Ax(t) + Bu(t) t ∈ [0,
       ẋ(t) = Ax(t) + Bu(t) + L [y(t) -ŷ(t)] t ∈ [0,
Une autre alternative est le filtre de Kalman [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], qui prend en compte les erreurs de mesure et modélise les incertitudes représentées par des bruits blancs gaussiens (à la fois stationnaires et mutuellement non corrélés), afin de construire une estimation d'état qui minimise l'erreur quadratique moyenne. Notons que des extensions au cas non linéaire ont été développées, e.g. le observateur non linéaire de Luenberger [START_REF] Afri | State and parameter estimation: A nonlinear Luenberger observer approach[END_REF] et le filtre de Kalman étendu [START_REF] Ide | Unified notation for data assimilation : Operational, sequential and variational[END_REF].

Une technique plus récente, développée par Auroux et Blum [START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF] est le Nudging direct et rétrograde (Back and Forth Nudging, BFN). Le Nudging direct consiste simplement à ajouter un terme de rétroaction dans l'équation gouvernante, comme le fait l'observateur de Luenberger, mais aussi à supposer des observations complètes dans le (1) (c-à-d y(t) = Cx obs (t), avec C inversible) et aucune entrée u(t). En utilisant les mêmes idées, le Nudging rétrograde considère un intervalle de temps borné [0, T ] pour approcher plutôt la condition initiale x 0 en résolvant ẋ(t) = Ax(t) -K [x obs (t) -xk (t)] dans [0, T ] x(T ) = xT où xT est une observation du système à l'instant T et K est la matrice de nudging, choisie symétrique et positive définie pour assurer la convergence asymptotique (∀t ∈ (0, T ]) lim

min λ∈σ(K) λ→+∞ x(t) = x obs (t). (2)

Méthodes variationnelles

D'autre part, Sasaki [START_REF] Sasaki | An objective analysis based on the variational method[END_REF] a proposé d'appliquer une approche différente aux problèmes de météorologie : ici, l'équation gouvernante contraint une fonctionnelle de coût J, représentant l'écart entre l'état réel x(t) et les données disponibles. L'objectif de cette formulation est de minimiser J = J(u), où u(t) agit comme une entrée de l'équation gouvernante (par exemple, la condition initiale ou limite). Autrement dit, le fait de trouver la variable de contrôle u donne l'état réel x = x(u). Les problèmes de ce type relèvent de la théorie du contrôle optimal appliquée aux EDP, dont le fondement théorique a été initialement développé par J.-L. Lions [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF].

À titre indicatif, nous considérons le problème de minimisation continue

min u∈U J(u) = 1 2 Ω (u -x b 0 ) B -1 (u -x b 0 )dx + 1 2 T 0 Ω (H(x) -y) R -1 (H(x) -y) dxdt s.t. ẋ(t) = M(x(t), t) dans Ω × [0, T ] x(0) = u dans Ω (3)
La condition initiale u appartient à un espace fonctionnel U qui résume les propriétés souhaitables du contrôle, tandis que la définition de J(u) fait intervenir plusieurs éléments : les matrices de covariance des erreurs de prévision et d'observation B et R, qui mesurent l'incertitude autour de l'estimation préalable de la condition initiale x b 0 et des observations y(t), respectivement ; et un opérateur différentiel H qui décrit les observations prévues du système de contrôle. Ce dernier est également modélisé par un opérateur différentiel M.

Nous avons besoin de calculer ∇J(u), soit pour dériver les conditions d'optimalité pour (3), soit pour résoudre numériquement ce problème. Dans ce qui suit, nous décrivons différentes approches pour le faire [START_REF] Allaire | Numerical analysis and optimization: An introduction to mathematical modelling and numerical simulation[END_REF][START_REF] Nodet | Inverse problems for the environment: tools, methods and applications[END_REF][START_REF] Ide | Unified notation for data assimilation : Operational, sequential and variational[END_REF].

Un calcul direct. Puisque x dépend implicitement de u, nous devons d'abord déterminer son comportement lorsque la condition initiale est légèrement perturbée dans une direction v. Cette variation, désignée par X , satisfait

     Ẋ (t) = ∂M ∂x (x(t), t)X dans Ω × [0, T ] X (0) = v dans Ω (4)
État de l'art et ensuite, un calcul direct donne

∇J(u), v = Ω B -1 (u -x b 0 ) + T 0 X ∂H ∂x (x(t), t) * R -1 (H(x) -y) dt v dx.
où * désigne l'opérateur adjoint. L'inconvénient de cette méthode est qu'elle nécessite de résoudre (4) pour chaque direction v.

La méthode adjointe. Cette procédure réduit (3) à un problème d'optimisation sans contrainte avec des variables supplémentaires, en définissant le Lagrangien

L(u, x, λ, µ) = J(u) + T 0 Ω λ [ ẋ(t) -M(x(t), t)] dxdt + Ω µ (x(0) -u) dx,
où λ et µ sont connus sous le nom des multiplicateurs de Lagrange associés à chaque équation du système principal.

Au lieu de traiter l'équation linéaire tangente (4), la fonction Lagrangienne nous permet de dériver une équation pour la variable duale de X de la façon suivante. La mise à zéro de la dérivée du Lagrangien par rapport à x dans la direction z, à (u, x, λ, µ), conduit à

∇ x L(u, x, λ, µ), z = 0 T 0 Ω z ∂H ∂x (x(t), t) * R -1 (H(x) -y) dxdt T 0 Ω λ ż(t) - ∂M ∂x (x(t), t)z dxdt + Ω µ z(0) dx = 0,
et après avoir intégré par parties la deuxième intégrale et en réorganisant les termes, nous obtenons

T 0 Ω z ∂H ∂x (x(t), t) * R -1 (H(x) -y) + λ - ∂M ∂x (x(t), t) * λ dxdt + Ω λ (T )z(T ) + Ω [µ -λ(0)] z(0) dx = 0.
Puisque z est une direction arbitraire, λ doit nécessairement satisfaire l'équation adjointe

     λ(t) - ∂M ∂x (x(t), t) * λ = - ∂H ∂x (x(t), t) * R -1 (H(x) -y) in Ω × [0, T ] λ(T ) = 0 in Ω (5)
et la condition supplémentaire µ = λ(0). Ensuite, nous pouvons calculer le gradient pour n'importe quelle direction v par

∇J(u), v = ∇ u L(u, x, λ, µ), v = Ω B -1 (u -x b 0 ) -λ (0) v dx.
Notons que nous devons résoudre (5) une seule fois, puisque la variable adjointe ne dépend pas de v.

J(u) = 1 2 (u -x b 0 ) B -1 (u -x b 0 ) + 1 2 N n=1 (H(x n ) -y n ) R -1 (H(x n ) -y n ) s.t. x n = M [t n-1 ,tn] (x n-1 ) ∀n = 1, . . . , N x 0 = u (6) où H est l'opérateur d'observation et {M [t n-1 ,tn] (•)} N n=1
est une famille d'opérateurs qui décrivent les transitions de l'état de t n-1 à t n . Puis

x n = M [t n-1 ,tn] • • • • • M [t 0 ,t 1 ] (u) := M [t 0 ,tn] (u) et le gradient peut être calculé par ∇J(u) = B -1 (u -x b 0 ) + N n=1 M n H n • R -1 (H • M [t 0 ,tn] (u) -y n ). (7)
Les matrices M n = DM [t 0 ,tn] (u) et H n = DH(x n ) sont connues comme les opérateurs linéaires tangentes associés à M [t 0 ,tn] et H, respectivement.

D'autres variantes de cette méthode sont 3D-Var, une version indépendante du temps qui peut aussi être appliquée à des problèmes évolutifs en supposant que toutes les observations ne sont disponibles qu'au début (c-à-d qu'aucun système principale n'est nécessaire) ; 3D-FGAT, une amélioration de la dernière dans laquelle nous remplaçons seulement M n par la matrice d'identité, simplifiant le calcul du gradient ; et 4D-Var Incrémental [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF], qui consiste à approximer (6) en utilisant une séquence de problèmes de minimisation quadratique pour réduire le coût opérationnel.

Méthodes parallèles en espace-temps

Pour l'une ou l'autre des méthodes d'AD susmentionnées, le calcul numérique de l'estimation de l'état est dans de nombreuses situations aussi pertinent que sa précision, c'est-à-dire que le problème est posé sous contrainte de calcul en temps réel. Ainsi, la première exige souvent d'être réalisée dans un délai raisonnable, ce qui est possible à l'aide de méthodes parallèles dans l'espace ou dans le temps, une approche naturelle pour accélérer la résolution numérique des EDP à l'aide du calcul parallèle. En suivant Gander [START_REF] Gander | Schwarz methods over the course of time[END_REF][START_REF] Gander | 50 years of time parallel time integration[END_REF], nous en décrivons brièvement quelques-unes.

Au cours du dix-neuvième siècle, l'analyse de Fourier a été le principal outil d'étude des EDP, bien qu'elle soit limitée à des géométries simples comme les cercles ou les rectangles. Dans le but d'étendre le principe de Dirichlet à des domaines arbitraires, Schwarz [START_REF] Schwarz | Über einen grenzübergang durch alternierendes verfahren[END_REF] a proposé de résoudre l'équation de Laplace en décomposant le domaine en deux sous-domaines qui se chevauchent, où un calcul explicite peut-être utilisé, par exemple via l'analyse de Fourier, puis en résolvant alternativement un problème réduit sur chacun, dans une procédure connue aujourd'hui sous le nom de Méthode de Schwarz alternée. Cette décomposition est le principe sous-jacent des méthodes de décomposition de domaines. Cent ans plus tard, Lions [START_REF] Lions | On the Schwarz alternating method[END_REF][START_REF] Lions | On the Schwarz alternating method II: Stochastic interpretation and orders properties[END_REF][START_REF] Lions | On the Schwarz alternating method. III: a variant for nonoverlapping subdomains[END_REF] a étendu la méthode précédente à un cadre parallèle, en envisageant la possibilité d'absence de chevauchement des sous-domaines et en résolvant de manière synchrone les problèmes locaux associés. La Méthode de Schwarz parallèle était née. Depuis lors, à une époque où les ordinateurs deviennent de plus en plus performants, de nombreuses méthodes ont été développées pour tirer parti de cette stratégie. Mais qu'en est-il de la dimension temporelle ? Comme la solution d'une EDP évolutive est naturellement affectée par le passé, le temps n'est généralement pas utilisé dans le calcul parallèle, même si des méthodes parallèles en temps ont été développées depuis plus de 50 ans. Ses origines remontent à Nievergelt [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF], qui a proposé l'idée principale derrière les Méthodes de tir multiples : décomposer l'intervalle de temps en sous-intervalles disjoints et résoudre simultanément une famille de problèmes de valeur initiale, en brisant la nature séquentielle intrinsèque de l'équation différentielle dépendante du temps. D'une manière plus générale, les méthodes de parallélisation espace-temps se distinguent par le caractère itératif ou direct de la procédure et la décomposition du domaine spatio-temporel considéré. Divisées en quatre classes, les méthodes de Tir multiples parallélisent sur l'intervalle de temps, les méthodes de relaxation d'ondes et décomposition de domaine utilisent plutôt la dimension spatiale et les méthodes Multigrille travaillent simultanément sur les deux, étant toutes de nature itérative. En revanche, les méthodes parallèles à temps direct tentent de récupérer une solution en utilisant un solveur direct.

Dans ce qui suit, nous rappelons les bases de l'algorithme Parareal, une des plus récentes méthodes de tir multiple ; et quelques exemples de parallélisation de problèmes d'assimilation de données.

L'algorithme Pararéel

La parallélisation en temps du problème

u(t) = f (u(t)), t ∈ [0, T ] u(0) = u 0
exige de décomposer l'intervalle de temps en N sous-intervalles, désignés par (t n-1 , t n ); et d'introduire des cibles intermédiaires qui servent de conditions initiales sur chacune. Une façon directe de déterminer ces valeurs consiste à résoudre un système d'équations non linéaire en appliquant la méthode de Newton, une procédure coûteuse par un système de grande taille puisqu'elle repose sur le calcul d'une matrice jacobienne.

Lions, Maday et Turinici [START_REF] Lions | Résolution d'EDP par un schéma en temps «pararéel[END_REF] ont proposé l'algorithme Pararéel, une approche différente qui utilise deux solveurs F et G qui calculent une approximation numérique fine et grossière de u sur les sous-intervalles et mettent à jour les conditions initiales artificielles. Gander et Vandewalle [START_REF] Gander | Analysis of the parareal time-parallel timeintegration method[END_REF] ont prouvé que cette méthode se lit comme une méthode de tir multiple dans laquelle la matrice jacobienne est approximée par une différence finie dans une grille grossière. L'algorithme Pararéel s'approche de {u(t n )} N n=1 par une séquence {U k n } N n=1 , qui est construite comme suit :

(a) puisque les conditions initiales sont inconnues sauf sur le premier sous-intervalle, nous imposons des valeurs arbitraires sur le reste en utilisant le solveur grossier G,

U 0 n = G(t n , t n-1 , U 0 n-1 ), U 0 0 = u 0 .
où G(t n , t n-1 , U 0 n-1 ) désigne la solution obtenue avec le solveur grossier à t n , en considérant U 0 n-1 comme condition initiale à t n-1 .

(b) Puis résoudre en parallèle les problèmes restreints

   u(t) = f (u(t)), t ∈ [t n-1 , t n ] u(t n-1 ) = U k n-1
en utilisant le solveur fin F, ce qui donne les approximations

{F(t n , t n-1 , U k n-1 )} N n=1 .
(c) Enfin, lissez les discontinuités précédemment introduites en définissant la suite

U k+1 n := F(t n , t n-1 , U k n-1 ) + G(t n , t n-1 , U k+1 n-1 ) -G(t n , t n-1 , U k n-1
Theorem (Convergence du Pararéel [36]). Soit F(t n , t n-1 , U k n ) la solution exacte sur le sous-domaine temporel (t n-1 , t n ) et soit G(t n , t n-1 , U k n )
une solution approximative avec une erreur de troncature locale bornée par C 3 ∆T p+1 , et satisfaisant

F(t n , t n-1 , x) -G(t n , t n-1 , x) = c p+1 (x)∆T p+1 + c p+2 (x)∆T p+2 + ... ,
pour ∆T petit, où les coefficients {c j } j≥p+1 sont continuellement différentiables, et supposons que G satisfait la condition de Lipschitz

G(t + ∆T, t, x) -G(t + ∆T, t, y) ≤ (1 + C 2 ∆T ) x -y . ( 8
)
Ensuite, à l'itération k de l'algorithme Pararéel, nous avons la borne

u(t n ) -U k n ≤ C 3 C 1 (C 1 ∆T p+1 ) k+1 k! (1 + C 2 ∆T ) n-(k+1) k j=0 (n -j). ( 9
)
En raison du terme produit dans [START_REF] Bartels | Total variation minimization with finite elements: convergence and iterative solution[END_REF], après k itérations de l'algorithme Parareal, l'approximation est exacte sur les premiers k sous-intervalles (comme le montre la figure 2), et donc elle converge au plus en N itérations.

Même si ce n'est pas explicitement indiqué, la constante C 2 doit être positive, raison pour laquelle l'hypothèse de Lipschitz (8) ne prend pas en compte le cas décroissant, c-à-d lorsque le solveur grossier est contractif. Puisque nous sommes intéressés à coupler cet algorithme avec l'observateur de Luenberger et à tirer profit de son comportement de décroissance, nous présentons dans le Chapitre 2 une variante du Théorème précédent qui couvre cette situation.

Méthodes de parallélisation en espace-temps et l'AD

Trémolet et Le Dimet [START_REF] Trémolet | Parallel algorithms for variational data assimilation and coupling models[END_REF] ont été parmi les premiers à aborder la parallélisation des problèmes d'assimilation de données variationnelles en météorologie. Dans un cadre continu, ils ont proposé une approche de décomposition des domaines combinée avec la méthode adjointe, en assignant à chaque sous-domaine une version locale de (3), avec un terme supplémentaire sur la fonctionnelle de coût locale pour renforcer la continuité de l'état entre les domaines adjacents.

Vingt ans plus tard, des nouvelles stratégies dans ce domaine incluent également la direction temporelle. Par exemple, Rao et Sandu [START_REF] Rao | A time-parallel approach to strong-constraint fourdimensional variational data assimilation[END_REF] appliquent un solveur quasi-Newton au problème 4D-Var, mais ils parallélisent en temps le calcul du gradient. Une approche plus sophistiquée est proposée par D'Amore et Cacciapuoti [START_REF] Cacciapuoti | DD-DA PinT-based model: A domain decomposition approach in space and time, based on parareal, for solving the 4D-Var data assimilation model[END_REF], qui combine l'algorithme Pararéel avec la méthode de Schwarz multiplicatif (MPS) pour résoudre 4D-Var.

La méthodologie associée à l'algorithme Parareal a également été appliquée à des problèmes d'optimisation, mais pas nécessairement liés à l'AD. Par exemple, Maday, Salomon et Turinici [START_REF] Maday | Monotonic parareal control for quantum systems[END_REF] ont proposé un algorithme spécifique pour résoudre les systèmes d'optimalité dans le cas du contrôle quantique. En définissant des états intermédiaires et en résolvant ensuite en parallèle une famille de problèmes d'optimisation locaux, ces valeurs sont mises à jour après avoir résolu séquentiellement les équations gouvernantes et adjointes, ce qui est peu coûteux en termes de calcul par rapport à la procédure d'optimisation. Comme la fonctionnelle de coût original n'est pas parallélisable, la méthode en utilise une autre qui dépend de la variable adjointe et peut également être décomposée comme la somme des fonctionnelles de coût pour chaque sous-intervalle. En définitive, sa solution optimale coïncide avec celle du problème original.

Estimation de la bathymétrie

Récemment, une littérature considérable s'est développée autour du sujet des problèmes inverses dans les écoulements de surface libre. Une revue de Sellier identifie différentes techniques appliquées à la reconstruction bathymétrique [83, Section 4.2], qui reposent principalement sur la dérivation d'une formule explicite pour la bathymétrie, la résolution numérique des équations gouvernantes, ou l'assimilation des données variationnelles [START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF].

Une approche variationnelle est également bien adaptée à la résolution de problèmes d'ingénierie côtière, car ils doivent satisfaire plusieurs contraintes mécaniques. Par exemple, parmi les nombreux aspects à considérer lors de la conception d'un port, la construction de structures de protection est essentielle pour le protéger contre l'impact des vagues. Ces structures peuvent être optimisées pour minimiser l'énergie des vagues, en étudiant son interaction avec les vagues réfléchies [START_REF] Isebe | Optimal shape design of defense structures for minimizing short wave impact[END_REF]. Bouharguane et Mohammadi [START_REF] Mohammadi | Optimal dynamics of soft shapes in shallow waters[END_REF][START_REF] Bouharguane | Minimization principles for the evolution of a soft sea bed interacting with a shallow[END_REF] envisagent une approche temporelle pour étudier l'évolution du mouvement du sable au fond de la mer, qui pourrait également permettre à ces structures de changer dans le temps. Dans ce cas, les fonctionnelles proposées sont minimisées localement en utilisant l'analyse de sensibilité, une technique largement appliquée dans les géosciences. D'un point de vue mathématique, la résolution de ce genre de problème est surtout numérique. Les questions concernant un espace de contrôle approprié pour la bathymétrie, la contrôlabilité, la régularité de la solution ou le caractère bien posé du problème ne sont généralement pas abordées. Une approche théorique appliquée à la modélisation des bassins de surf peut être trouvée dans les publications [START_REF] Dalphin | Shape optimization of a moving bottom underwater generating solitary waves ruled by a forced KdV equation[END_REF][START_REF] Nersisyan | Generation of two-dimensional water waves by moving bottom disturbances[END_REF], où le but est de maximiser l'énergie totale de la vague prescrite. La première propose de déterminer une bathymétrie, tandis que la seconde définit la forme et le déplacement d'un objet sous-marin à une profondeur constante.

Modélisation des vagues

Dans tous les problèmes précédents, une contrainte cruciale pour déterminer la bathymétrie est son interaction avec les vagues. La pierre angulaire de sa modélisation sont les équations de Navier-Stokes

             ρ ∂u ∂t + (u • ∇) u = div (σ T ) + g dans Ω t , div (u) = 0 dans Ω t , u = u 0 dans Ω 0 , où u = (u, v,
         ∂η ∂t + ∂hu ∂x = 0 ∂u ∂t + u ∂u ∂x + g ∂η ∂x = z 2 b 3 ∂ 3 u ∂x 2 ∂t
pour une profondeur constante -z b . Ce modèle et ses variantes ultérieures sont utilisés pour décrire les mouvements de flottabilité dans les fluides, comme les vagues entrant dans la zone proche du rivage [START_REF] Barthélemy | Nonlinear shallow water theories for coastal waves[END_REF] ou les vagues générées par les glissements de terrain [START_REF] Dutykh | Boussinesq modeling of surface waves due to underwater landslides[END_REF].

Méthode de l'élément de pale

La théorie BEM repose alors sur une décomposition du système fluide/turbine en une partie macroscopique via la MT et une partie planaire locale via la BET. Cette dernière considère une décomposition radiale pour les pales et la colonne de fluide (Figure 3a), en divisant la zone du rotor en anneaux concentriques d'épaisseur infinitésimale qui n'interagissent pas entre eux.

tan ϕ = 1 -a λ(1 + a) , ( 11
)
a 1 -a = 1 4 sin 2 ϕ • Bc λ 2πr (C L (ϕ -γ λ ) cos ϕ + C D (ϕ -γ λ ) sin ϕ), ( 12
)
a 1 -a = 1 4λ sin 2 ϕ • Bc λ 2πr (C L (ϕ -γ λ ) sin ϕ -C D (ϕ -γ λ ) cos ϕ). ( 13
)
où r indique la distance entre l'élément de pale pris en compte (figure 3b) et le rotor, λ = λ(r) est le rapport de vitesse au bout et B le nombre total de pales. Le corde c λ et l'angle de torsion γ λ sont donnés par la géométrie des pales.

Diverses modifications ont été apportées au modèle pour tenir compte, par exemple, de l'écoulement autour de l'extrémité d'une pale ou d'un état de sillage turbulent. Pour un compte rendu plus détaillé, voir la section 4.1.4.

Enfin, nous évaluons le rendement de la turbine en utilisant le coefficient de puissance

C p (ϕ) = 8 λ max λmax λ min λ 3 a (1 -a) 1 - C D C L (ϕ -γ λ ) tan -1 ϕ dλ.
Inversement, nous pouvons concevoir un profil de pale qui maximise cette quantité en résolvant max

ϕ C p (ϕ) s.t (11-13) α = ϕ -γ λ ,
où α est choisi de telle sorte que le rapport C D C L soit proche de zéro. Afin de simplifier ce problème d'optimisation, on suppose généralement un coefficient de traînée nul, ce qui conduit à une solution explicite. A notre connaissance, le cas général concernant la traînée non nulle et les corrections possibles du modèle n'a pas été abordé.

General Introduction

The present thesis aims to contribute to the development of a theoretical framework for three problems in the context of renewable marine energy, namely the timeparallelization of sequential data assimilation problems, bathymetry optimization and mathematical analysis of the blade element momentum theory. Since their solving rely largely on empirical knowledge, we believe that adopting a mathematical standpoint leads to a better understanding of the different situations, this being an opportunity to encourage interdisciplinary collaboration between mathematics and applied sciences as geophysics and engineering. In what follows, we briefly introduce each of these topics.

Time-parallelization of sequential data assimilation problems

The assumptions behind a mathematical model not only determine their range of applicability, but also induce an inevitable gap between predictions and reality. In order to narrow this difference, we can sacrifice the simplicity of the model or incorporate real data instead, by following a Data assimilation (DA) procedure. Several aspects have to be considered when applying these techniques to e.g. meteorology or oceanography problems, but here we recall that due to the number of state variables and the vast amount of observations required, their numerical solving is computationally expensive. Trémolet and Le Dimet [START_REF] Trémolet | Parallel algorithms for variational data assimilation and coupling models[END_REF] were among the first to address the parallelization of variational DA problems (which are based on optimal control and use the information collected in a fixed amount of time) by using a Domain decomposition (DD) approach. Since then, the coupling between these two procedures have been widely studied.

Roughly speaking, DD methods consist in decomposing the spatial dimension into subdomains, with possible overlap, and then solve syncronously a local problem on each. This divide-and-conquer strategy seems counterintuitive when handling the time direction, due to its inherently sequential nature, which is why time-parallelization is not commonly applied to DA problems.

Parallel-in-time algorithms can be quite useful when dealing with large time intervals, as is the case of sequential DA methods, where information can arrive uninterrupted. Then a natural question arises: can we combine both procedures? We begin to answer this question in Chapter 2, by studying the Luenberger observer and its coupling with the Parareal algorithm.

Bathymetry optimization

Despite the fact that the bathymetry can be inaccurately known in many situations, wave propagation models strongly depend on this parameter to capture the flow behavior, which emphasize the importance of studying inverse problems concerning its reconstruction from observed free surface data.

These kinds of problem are usually solved by simply discretizing the governing equations or with the help of sequential DA methods. Another alternative is to consider the bathymetry as control variable of a PDE-constrained optimization problem, an approach used in coastal engineering due to mechanical constraints associated with building structures and their interaction with sea waves. However, its solving rely mostly in sensitivity analysis, numerical methods or simplifications of the model that leads to explicit solutions and then, questions regarding a suitable control space, continuity of the control-to-state mapping, regularity or well-posedness of the solution are in general not adressed. Recently, Dalphin and Barros [START_REF] Dalphin | Shape optimization of a moving bottom underwater generating solitary waves ruled by a forced KdV equation[END_REF] carried out this theoretical analysis to model a wavemaker.

In Chapter 3 we try to answer the questions above when considering a reformulation of the Helmholtz equation for modeling wave propagation and a general cost functional that can be identified, for instance, with the mismatch between the predicted and observed wave solution.

Mathematical analysis of the Blade element momentum theory

Hydrotube Énergie, enterprise devoted to nautical design and on-board electronics systems, tested in 2015 a propotype of a floating water turbine on the Garonne river, in Bordeaux. But moving towards industrial scale production requires the development of a numerical software to simulate an optimized device, which is the purpose of its partnership with team ANGE (INRIA), a research group working in modeling, analysis and simulation of geophysical flows.

Among the problems that need to be addressed, we study the optimization of the turbine efficiency via the Blade element momentum (BEM) theory. Introduced by Glauert [START_REF] Glauert | The Elements of Aerofoil and Airscrew Theory[END_REF], this method is widely used for determining the efficiency and hence the design parameters of a blade, according to their geometric and mechanical characteristics and the current to which it is exposed. It follows from the combination of two independent methods that treat the fluid/turbine system from a macroscopic and local perspective, which results in a set of equations that are solved by applying an iterative algorithm. In some cases, a correction factor is required to ensure existence of solutions of the model. Although old, this method is still used due to its relative simplicity compared to the complexity of the hydrodynamic phenomenon developed in the fluid/turbine system. Nevertheless, both the existence of solutions and the algorithm convergence has been never analyzed from a mathematical point of view, being this the aim of Chapter 4.

Contributions of this thesis

This work is divided into three independent chapters regarding the aforementioned topics, in which I had the pleasure to collaborate with Julien Salomon (INRIA Paris), Felix Kwok (Hong Kong Baptist University), Pierre-Henri Cocquet (Université de Pau et des Pays de l'Adour) and Jérémy Ledoux (Hydrotube Énergie). In the following, we summarize our main contributions to each.

Time-parallelization of sequential data assimilation problems

We start in Chapter 2 by proposing a procedure to couple sequential data assimilation methods with parallel-in-time algorithms, that consists in splitting the unbounded time interval into subintervals of the same length (windows) and then apply, following a sequential order, the time-parallel solver on each. By considering the Luenberger observer as assimilation method, we provide a stopping criterion that preserves its exponential rate of convergence, which yields an a posteriori estimate of the accuracy of the solver.

In order to go further, we set the Parareal algorithm as parallel-in-time solver. This allows us to design an alternative algorithm that provides an a priori estimate of the number of iterations required on each window, which also enables us to investigate the theoretical efficiency of the entire procedure. These results are based on a new convergence estimate that we derive for Parareal when the coarse solver is contractive.

Bathymetry optimization

We move on to study in Chapter 3 the determination of a bathymetry from an optimization problem, where a reformulation of the Helmholtz equation acts as a constraint. Even though this equation is limited to describe waves of small amplitude, it is often used in engineering due to its simplicity, which leads to explicit solutions when a flat bathymetry is assumed. By lifting this hypothesis, we obtain a different formulation in which this variable plays the role of a scatterer.

Under suitable assumptions on the cost functional and the admissible set of bathymetries, we are able to prove the continuity of the control-to-state mapping and the existence of an optimal solution, in addition to the continuity and boundedness of the resulting wave. The discrete optimization problem is also addressed, studying the convergence to the discrete optimal solution as well as the convergence of a finite element approximation.

Mathematical analysis of the Blade element momentum theory

We finally discuss in Chapter 4 the existence of solutions and convergence of solving procedures for the Blade element momentum theory. The key point of our work consists in showing that the decomposition proposed by Glauert can be used to reformulate its original set of equations into a single expression containing two very distinct terms: a universal one, independent of the turbine considered and associated with macroscopic aspects of the model; and an experimental term, concerning the characteristics of the turbine and related with the local part of the model.

The advantage of our approach is that it explicitly identifies assumptions about the turbine parameters that guarantee the existence of a solution, but it also helps us to present convergence criteria for different solving algorithms. Glauert's model is also used to optimize blade geometries, in the sense that it maximizes the turbine efficiency. We recall the details of the usual design procedure and briefly discuss the case when a correction of the model is introduced.

Chapter 1

State of the art

This chapter is devoted to an overview of the different subjects discussed in this thesis. We begin with two topics studied in Chapter 2: Data assimilation methods, both in continuous and discrete setting; followed by a summary of space-time parallel methods, in particular the Parareal algorithm. We move on to discuss the key elements of Chapter 3, Bathymetry estimation and Wave modeling, with emphasis on the main assumptions needed to derive different models for wave propagation. It is followed by a brief review of the Blade Element Momentum theory, which is treated in Chapter 4.

Data assimilation (DA)

Mathematical models are widely used to describe complex systems. They rely on approximations and simplifications of a real phenomenon, which ultimately defines its range of applicability. After proper validation, we would like to apply them to describe or even predict a real situation, as monitoring earthquakes, forecast seasonal influenza or estimate the state of charge of a battery [START_REF] Hu | Estimation of state of charge of a lithium-ion battery pack for electric vehicles using an adaptive Luenberger observer[END_REF]. In order to do so, we need to incorporate real data into our framework. The different techniques that combines mathematical models with available observations to improve the knowledge of a system are known as Data assimilation (DA). Among these approaches, here we recall the Sequential and Variational methods.

Sequential methods

Let us suppose that for some reason, e.g. physical or budgetary constraints, we only have access to partial information of the system. Even in some cases the initial condition is not precisely known, as it happens in climate science [START_REF] Werndl | Initial-condition dependence and initial-condition uncertainty in climate science[END_REF]. We can treat this lack of information and uncertainty by constructing a new system which uses the available observations to approximate the true state. In a deterministic context, this device is called an observer.

We assume that the system is governed by

       ẋ(t) = Ax(t) + Bu(t) t ∈ [0, +∞) x(0) = x 0 , y(t) = Cx(t), (1.1)
where x ∈ R m is the true state vector, y ∈ R q represents the observations (with q < m, m, q ∈ N * ), u is an input and x(0) = x 0 an unknown initial condition. The matrices A, B and C are known and their dimensions are consistent.

The Luenberger observer [START_REF] Luenberger | Introduction to Dynamic Systems: Theory, Models, and Applications[END_REF] imitates the previous model, but includes an extra term in the dynamic that measures the misfit between the observations and its own predictions. It produces a state estimate x satisfying

       ẋ(t) = Ax(t) + Bu(t) + L [y(t) -ŷ(t)] t ∈ [0, +∞) x(0) = x0 , ŷ(t) = C x(t),
with x0 an arbitrary initial condition. As long as the original model (1.1) is observable (i.e. the initial state x(0) can be uniquely determined from the observations in [0, T ], for any T ), the resulting error x(t) -x(t) can be driven to zero at exponential rate by properly choosing the observer gain matrix L, and then the true state is recovered asymptotically.

Another alternative is the Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], which takes into account measurement errors and model innacuracies represented by Gaussian white noises (both stationary and mutually uncorrelated), in order to construct a state estimate that minimizes the mean square error. Note that extensions to the nonlinear case have been developed, e.g. the nonlinear Luenberger observer [START_REF] Afri | State and parameter estimation: A nonlinear Luenberger observer approach[END_REF] and Extended Kalman filter [START_REF] Ide | Unified notation for data assimilation : Operational, sequential and variational[END_REF].

A more recent technique, developed by Auroux and Blum [START_REF] Auroux | Back and forth nudging algorithm for data assimilation problems[END_REF] is Back and forth nudging (BFN). Nudging simply consists in adding a feedback term in the governing system, as the Luenberger observer does, but also assuming full observations in (1.1) (i.e. y(t) = Cx obs (t), with C invertible) and no input u(t). Using the same ideas, Backward nudging considers a bounded time interval [0, T ] to approximate instead the initial condition x 0 by solving

ẋ(t) = Ax(t) -K [x obs (t) -xk (t)] in [0, T ] x(T ) = xT
where xT is an observation of the system at time T and K is the nudging matrix, chosen symmetric and positive definite to ensure the asymptotic convergence

(∀t ∈ (0, T ]) lim min λ∈σ(K) λ→+∞ x(t) = x obs (t). (1.2)
The BFN algorithm combines these procedures by defining the iterative method

ẋk (t) = Ax k (t) + K [x obs (t) -xk (t)] in [0, T ] xk (0) = xk-1 (0) ẋk (t) = Ax k (t) -K [x obs (t) -xk (t)] in [0, T ] xk (T ) = xk (T )
with x0 (0) = x 0 . Both sequences {x k (t)} k≥1 and {x k (t)} k≥1 converge to x∞ (t) and x∞ (t), respectively. Moreover, these limit functions also exhibit the asymptotic behavior described in (1.2), even for t = 0.

Variational methods

On the other hand, Sasaki [START_REF] Sasaki | An objective analysis based on the variational method[END_REF] proposed to apply a different approach to meteorology problems: here the governing system constraints a cost functional J, representing the mismatch between the true state x(t) and available data. The goal of this formulation is to minimize J = J(u), where u(t) acts as an input of the governing system (for instance, the initial or boundary condition). In other words, finding the control variable u yields the true state x = x(u). Problems of these kinds fall under the Optimal control theory applied to PDEs, whose theoretical background was originally developed by J.-L. Lions [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF].

For illustrative purposes, we consider the continuous minimization problem

min u∈U J(u) = 1 2 Ω (u -x b 0 ) B -1 (u -x b 0 )dx + 1 2 T 0 Ω (H(x) -y) R -1 (H(x) -y) dxdt s.t. ẋ(t) = M(x(t), t) in Ω × [0, T ] x(0) = u in Ω (1.3)
The initial condition u belongs to a functional space U that summarizes desirable properties of the control, whereas the definition of J(u) involves several elements: the background and observation error-covariances matrices B and R, which measures the uncertainty around the prior estimate of the initial condition x b 0 and observations y(t), respectively; and a differential operator H that describes the predicted observations of the governing system. The latter is also modeled by a differential operator M.

We require to compute ∇J(u), either to derive optimality conditions for (1.3) or to solve numerically this problem. In what follows, we outline different approaches to do so [START_REF] Allaire | Numerical analysis and optimization: An introduction to mathematical modelling and numerical simulation[END_REF][START_REF] Nodet | Inverse problems for the environment: tools, methods and applications[END_REF][START_REF] Ide | Unified notation for data assimilation : Operational, sequential and variational[END_REF].

A direct computation

Since x implicitly depends on u, we must first determine its behavior when the initial condition is slightly perturbed in a direction v. This variation, denoted by X , it satisfies

     Ẋ (t) = ∂M ∂x (x(t), t)X in Ω × [0, T ] X (0) = v in Ω (1.4)
and then, a direct calculation yields

∇J(u), v = Ω B -1 (u -x b 0 ) + T 0 X ∂H ∂x (x(t), t) * R -1 (H(x) -y) dt v dx.
where * denotes the adjoint operator. The drawback of this method is that requires to solve (1.4) for each direction v.

The Adjoint method

This procedure reduces (1.3) to an unconstrained optimization problem with additional variables, by defining the Lagrangian

L(u, x, λ, µ) = J(u) + T 0 Ω λ [ ẋ(t) -M(x(t), t)] dxdt + Ω µ (x(0) -u) dx,
where λ and µ are known as the Lagrange multipliers associated with each equation of the governing system.

Instead of dealing with the Tangent-linear equation (1.4), the Lagrangian function allow us to derive an equation for the dual variable of X in the following way. Setting to zero the derivative of the Lagrangian w.r.t. x in the direction z, at (u, x, λ, µ), leads to

∇ x L(u, x, λ, µ), z = 0 T 0 Ω z ∂H ∂x (x(t), t) * R -1 (H(x) -y) dxdt T 0 Ω λ ż(t) - ∂M ∂x (x(t), t)z dxdt + Ω µ z(0) dx = 0,
and after integrating by parts the second integral and rearranging terms, we obtain

T 0 Ω z ∂H ∂x (x(t), t) * R -1 (H(x) -y) + λ - ∂M ∂x (x(t), t) * λ dxdt + Ω λ (T )z(T ) + Ω [µ -λ(0)] z(0) dx = 0.
Since z is an arbitrary direction, λ must necessarily satisfy the Adjoint equation

     λ(t) - ∂M ∂x (x(t), t) * λ = - ∂H ∂x (x(t), t) * R -1 (H(x) -y) in Ω × [0, T ] λ(T ) = 0 in Ω (1.5)
and the supplementary condition µ = λ(0). Then we can compute the gradient for any direction v by

∇J(u), v = ∇ u L(u, x, λ, µ), v = Ω B -1 (u -x b 0 ) -λ (0) v dx.
Note that we need to solve (1.5) only once, since the adjoint variable does not depend on v.

Discrete variational methods

In contrast to the continuous setting, which requires to solve the tangent-linear equation (1.4) or the adjoint equation to obtain the gradient, discretizing (1.3) leads to its direct computation. Then, discrete variational methods differ in the way they describe the governing system and the treatment of possible nonlinearities.

Given a discretization {t n } N n=0 of [0, T ], we denote by x i and y i the state and observation vector at time t i , respectively. Using the same notation as before on the variables J, u, x b 0 , B and R, even if they could depend on the time-discretization considered, the 4D-Var algorithm [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] reads

min u∈U J(u) = 1 2 (u -x b 0 ) B -1 (u -x b 0 ) + 1 2 N n=1 (H(x n ) -y n ) R -1 (H(x n ) -y n ) s.t. x n = M [t n-1 ,tn] (x n-1 ) ∀n = 1, . . . , N x 0 = u (1.6)
where H is the observation operator and

{M [t n-1 ,tn] (•)} N n=1
is a family of operators that describe the transitions of the state from t n-1 to t n . Then

x n = M [t n-1 ,tn] • • • • • M [t 0 ,t 1 ] (u) := M [t 0 ,tn] (u)
and the gradient can be computed by

∇J(u) = B -1 (u -x b 0 ) + N n=1 M n H n • R -1 (H • M [t 0 ,tn] (u) -y n ). (1.7)
The matrices M n = DM [t 0 ,tn] (u) and H n = DH(x n ) are known as the tangent linear operators associated with M [t 0 ,tn] and H, respectively.

Space-time parallel methods

Other variants of this method are 3D-Var, a time independent version that can be also applied to evolutionary problems by assuming that all the observations are only available at the beginning (i.e. no governing system is needed); the 3D-FGAT, an improvement of the last in which we replace only M n by the identity matrix, simplifying the computation of the gradient; and the Incremental 4D-Var [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF], which consists in approximate (1.6) using a sequence of quadratic minimization problems to reduce the operative cost.

Space-time parallel methods

For any of the aforementioned DA methods, numerical computation of the state estimate is as relevant as its accuracy. However, the former also requires to be carry out in a reasonable amount of time, which is possible with the help of space or time parallel methods, a natural approach to speed-up the numerical resolution of PDEs using parallel computing. Following Gander [START_REF] Gander | Schwarz methods over the course of time[END_REF][START_REF] Gander | 50 years of time parallel time integration[END_REF], we briefly describe some of them.

During the nineteenth century, Fourier analysis was the main tool for studying PDEs, though it is restricted to simple geometries as circles or rectangles. With the purpose of extending Dirichlet's principle to arbitrary domains, Schwarz [START_REF] Schwarz | Über einen grenzübergang durch alternierendes verfahren[END_REF] proposed to solve the Laplace equation by decomposing the domain into two overlapping subdomains where Fourier analysis apply, and then solving alternately a reduced problem on each, in a procedure known nowadays as the Alternating Schwarz method. This decomposition is the underlying principle of Domain decomposition methods. A hundred years later, Lions [START_REF] Lions | On the Schwarz alternating method[END_REF][START_REF] Lions | On the Schwarz alternating method II: Stochastic interpretation and orders properties[END_REF][START_REF] Lions | On the Schwarz alternating method. III: a variant for nonoverlapping subdomains[END_REF] extended the previous method to a parallel framework, by considering possibly overlapping subdomains and solving syncronously the associated local problems. The Parallel Schwarz method was born. Since then, in a time when computers become more and more efficient, numerous methods have been developed to take advantage of this strategy. But what about the time dimension? Since the solution of an evolutionary PDE is naturally affected by the past, time is not usually used in parallel computing, even if parallel-in-time methods have been developed for more than 50 years. Its origins can be traced back to Nievergelt [START_REF] Nievergelt | Parallel methods for integrating ordinary differential equations[END_REF], who proposed the main idea behind Multiple shooting methods: decompose the time interval into disjoint subintervals and solve simultaneously a family of initial-value problems, breaking the intrinsic sequential nature of the time dependent differential equation.

In a more general way, space-time methods differ in the iterative of direct nature of the procedure and the decomposition of the space-time domain considered. Divided into four classes, Multiple shooting methods parallelize along the time interval, Waveform relaxation and Domain decomposition methods use the space dimension instead and Space-time multigrid methods work simultaneously on both, being all of them of iterative nature. In contrast, Direct time parallel methods attempt to retrieve a solution by using a direct solver.

In what follows we recall the basics of the Parareal algorithm, one of the most recent Multiple shooting methods; and some examples of parallelization of data assimilation problems.

The Parareal algorithm

The time parallelization of the problem

u(t) = f (u(t)), t ∈ [0, T ] u(0) = u 0
requires to decompose the time interval on N subintervals, denoted by (t n-1 , t n ); and introduce intermediate targets that act as initial conditions on each. A direct way to determine these values is solving a nonlinear system of equations by applying Newton's method, an expensive procedure for large systems since it relies on the computation of a Jacobian matrix.

Lions, Maday and Turinici [START_REF] Lions | Résolution d'EDP par un schéma en temps «pararéel[END_REF] proposed the Parareal algorithm, a different approach which uses two solvers F and G that compute a fine and coarse numerical approximation of u on the subintervals and update the artificial initial conditions. Gander and Vandewalle [START_REF] Gander | Analysis of the parareal time-parallel timeintegration method[END_REF] proved that this method reads as a Multiple shooting method in which the Jacobian matrix is approximated by a finite difference in a coarse grid.

The Parareal algorithm approximates {u(t n )} N n=1 by a sequence {U k n } N n=1 , which is contructed as follows:

(a) since the initial conditions are unknown except on the first subinterval, impose arbitrary values on the rest by using the coarse solver G,

U 0 n = G(t n , t n-1 , U 0 n-1 ), U 0 0 = u 0 . where G(t n , t n-1 , U 0 n-1
) denotes the solution obtained with the coarse at solver t n , considering U 0 n-1 as initial condition at t n-1 . (b) Then solve in parallel the restricted problems

   u(t) = f (u(t)), t ∈ [t n-1 , t n ] u(t n-1 ) = U k n-1
using the fine solver F, which yields the approximations {F(t n , t n-1 , U k n-1 )} N n=1 . (c) Finally, smooth the discontinuities previously introduced by defining the sequence

U k+1 n := F(t n , t n-1 , U k n-1 ) + G(t n , t n-1 , U k+1 n-1 ) -G(t n , t n-1 , U k n-1 )
, where the superscript k denotes the current number of iterations. Note that on the right-hand side, the first and third term were already computed, whereas the second shows that the update must be done sequentially. Hopefully, this is not computationally expensive since it only depends on the coarse solver G.

The advantage of this algorithm is its superlinear convergence rate. Indeed, due to Gander and Hairer, we have the a priori estimate: Theorem 1.2.1 (Convergence of Parareal [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF]). Let F(t n , t n-1 , U k n ) be the exact solution on the time subdomain (t n-1 , t n ) and let G(t n , t n-1 , U k n ) be an approximate solution with local truncation error bounded by C 3 ∆T p+1 , and satisfying

F(t n , t n-1 , x) -G(t n , t n-1 , x) = c p+1 (x)∆T p+1 + c p+2 (x)∆T p+2 + ... ,
for ∆T small, where the coefficients {c j } j≥p+1 are continuously differentiable, and assume that G satisfies the Lipschitz condition

G(t + ∆T, t, x) -G(t + ∆T, t, y) ≤ (1 + C 2 ∆T ) x -y . (1.8)
Then, at iteration k of the Parareal algorithm, we have the bound

u(t n ) -U k n ≤ C 3 C 1 (C 1 ∆T p+1 ) k+1 k! (1 + C 2 ∆T ) n-(k+1) k j=0 (n -j). (1.9)
As a result of the product term in (1.9), after k iterations of the Parareal algorithm the approximation is exact on the first k subintervals (as shown in Figure 1.2), and hence it converges in at most N iterations.

(t) = -i 2 u in [0, 20]
Even though it is not explicitly stated, the constant C 2 must be positive, which is why the Lipschitz assumption (1.8) does not take into account the decaying case, i.e. when the coarse solver is contracting. Since we are interested in coupling this algorithm with the Luenberger observer and take advantage of its decaying behavior, in Chapter 2 we present a variant of Theorem 1.2.1 that covers this situation.

Space-time parallel methods and DA

Trémolet and Le Dimet [START_REF] Trémolet | Parallel algorithms for variational data assimilation and coupling models[END_REF] were among the first to address the parallelization of Variational data assimilation problems in meteorology. In a continuous setting, they proposed a Domain decomposition approach combined with the Adjoint method, by assigning to each subdomain a local version of (1.3), with an extra term on the local cost functional to enforce the continuity of the state between adjacent domains. Twenty years later, new strategies in this topic also include the time direction. For instance, Rao and Sandu [START_REF] Rao | A time-parallel approach to strong-constraint fourdimensional variational data assimilation[END_REF] apply a quasi-Newton solver to the 4D-Var problem, but they time-parallelize first the computation of the gradient. A more sophisticated approach is proposed by D'Amore and Cacciapuoti [START_REF] Cacciapuoti | DD-DA PinT-based model: A domain decomposition approach in space and time, based on parareal, for solving the 4D-Var data assimilation model[END_REF], who combines the Parareal algorithm with the Multiplicative Parallel Schwarz method (MPS) to solve 4D-Var.

The methodology associated with the Parareal algorithm has also been applied to optimization problems, but not necessarily related to DA. For instance, Maday, Salomon and Turinici [START_REF] Maday | Monotonic parareal control for quantum systems[END_REF] proposed a specific algorithm to solve optimality systems in the case of quantum control. By defining intermediate states and then solving a family of local optimization problems in parallel, these values are updated after solving sequentially the forward and adjoint equations, which is computationally cheap compared with the optimization procedure. Since the original cost functional is not parallelizable, the method uses a different one which depends on the adjoint variable and can be also decomposed as the sum of cost functionals for each subinterval. Ultimately, its optimal solution coincides with that of the original problem.

Bathymetry estimation

Recently, considerable literature has grown up around the subject of inverse problems in free surface flows. A review from Sellier identifies different techniques applied for bathymetry reconstruction [83, Section 4.2], which rely mostly on the derivation of an explicit formula for the bathymetry, numerical resolution of a governing system, or variational data assimilation [START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF].

A variational approach is also well suited for solving coastal engineering problems, since they must satisfy several mechanical constraints. For instance, among the several aspects to consider when designing a harbor, building defense structures is essential to protect it against wave impact. These structures can be optimized to minimize the wave energy, by studying its interaction with the reflected waves [START_REF] Isebe | Optimal shape design of defense structures for minimizing short wave impact[END_REF]. Bouharguane and Mohammadi [START_REF] Mohammadi | Optimal dynamics of soft shapes in shallow waters[END_REF][START_REF] Bouharguane | Minimization principles for the evolution of a soft sea bed interacting with a shallow[END_REF] consider a time-dependent approach to study the evolution of sand motion at the seabed, which could also allow these structures to change in time. In this case, the proposed functionals are locally minimized using sensitivity analysis, a technique broadly applied in geosciences.

From a mathematical point of view, the solving of these kinds of problem is mostly numerical. Questions about a suitable control space for the bathymetry, controllability, regularity or well-posedness of the solution are not usually adressed. A theoretical approach applied to the modeling of surfing pools can be found in [START_REF] Dalphin | Shape optimization of a moving bottom underwater generating solitary waves ruled by a forced KdV equation[END_REF][START_REF] Nersisyan | Generation of two-dimensional water waves by moving bottom disturbances[END_REF], where the goal is to maximize the total energy of the prescribed wave. The former proposes to determine a bathymetry, whereas the latter sets the shape and displacement of an underwater object along a constant depth.

Wave modeling

In all the preceding problems, a crucial constraint for determining the bathymetry is its interaction with sea waves. The cornerstone of its modeling is the Navier-Stokes system

             ρ ∂u ∂t + (u • ∇) u = div (σ T ) + g in Ω t , div (u) = 0 in Ω t , u = u 0 in Ω 0 ,
where u = (u, v, w) denotes the velocity of the fluid in the bounded and timedependent region Ω t ; ρ, g and σ T represent the density coefficient, gravity and total stress tensor, respectively. Here we omit the boundary conditions on the water level η(x, t) and at the bottom -z b (x), possibly depending on the pressure, viscosity and friction effects, see Section 3.1.1 for more details. In what follows, for simplicity we neglect the last two.

Certain simplifications of the Navier-Stokes system gives rise to different models for wave propagation [START_REF] Bristeau | Derivation of a non-hydrostatic shallow water model; comparison with Saint-Venant and Boussinesq systems[END_REF][START_REF] Méhauté | An Introduction to Hydrodynamics and Water Waves[END_REF], all of them depending on the relationship between three characteristic parameters known as the relative depth H, the horizontal lenght L and the maximum vertical amplitude A. To recall some of them, we introduce the ratios

ε = H L , δ = A H .
Vertically averaged fluids as rivers, coastal domains and oceans can be modeled by using the shallow water assumption ε 1, whereas small amplitude waves are described by models where δ 1.

Assuming a shallow water regime, as well as different values of δ and approximations of the pressure, leads to the following models.

Saint-Venant system

The hydrostatic pressure approximation (i.e. the contribution of the vertical acceleration of the fluid in the pressure is negligible) combined with δ = O(1) yields

                     ∂η ∂t + ∂hu ∂x + ∂hv ∂y = 0 ∂hu ∂t + ∂hu 2 ∂x + ∂hu v ∂y + gh ∂η ∂x = 0 ∂hv ∂t + ∂hu v ∂x + ∂hv 2 ∂x + gh ∂η ∂y = 0 (1.10)
where h = η + z b is the water height and u, v are depth-averaged velocities. This model is commonly used for modeling dam breaks [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF], hydraulic jumps and tidal flows in estuaries.

Wave equation

Keeping the hydrostatic pressure but assuming instead δ 1 is equivalent to neglect the convective acceleration terms in (1.10) and linearize the new system around the sea level η = 0. Combining the resulting expressions brings

∂ 2 η ∂t 2 -∇ (gz b ∇η) = 0.
In particular, its solution has the form η(x, t) = Re{ψ tot (x)e -iωt }, where ω is a given frequency and the amplitude ψ tot satisfies the Helmholtz equation

ψ tot + ω 2 gz b ψ tot = 0,
when constant depth -z b is assumed. Under suitable boundary conditions, this equation is used for studying water-wave scattering problems, such as the one presented in Chapter 3, where a variable bathymetry acts as a scatterer.

Blade element momentum (BEM) theory

Classical Boussinesq system

On the other hand, a non-hydrostatic pressure introduces dispersive terms in the governing equations. A third ratio, called the Ursell number

U r = δ ε 2 ,
measures their strength with respect to the nonlinear terms [START_REF] Ursell | The long-wave paradox in the theory of gravity waves[END_REF]. In two-dimensional motion, a specific choice of this parameter (U r = O(1) and δ 1) leads to

         ∂η ∂t + ∂hu ∂x = 0 ∂u ∂t + u ∂u ∂x + g ∂η ∂x = z 2 b 3 ∂ 3 u ∂x 2 ∂t
for a constant depth -z b . This model and its subsequent variants are used to describe buoyant motion in fluids, as waves entering the near-shore zone [START_REF] Barthélemy | Nonlinear shallow water theories for coastal waves[END_REF] or landslide generated waves [START_REF] Dutykh | Boussinesq modeling of surface waves due to underwater landslides[END_REF].

Blade element momentum (BEM) theory

Blade element momentum (BEM) theory is a mechanical model widely used to evaluate turbine performance, according to the mechanical/geometric characteristics of their blades and the current to which they are exposed. Developed by Glauert [START_REF] Glauert | The Elements of Aerofoil and Airscrew Theory[END_REF], it follows from the combination of two different models: Blade element theory (BET) and Momentum theory (MT).

By cutting the blade into sections that are treated according to a planar model, Blade element theory studies the turbine behavior from a local point of view [START_REF] Froude | On the elementary relation between pitch, slip and propulsive efficiency[END_REF]. The fundamental quantities of this model are the coefficients C L and C D called drag and lift, which are introduced to account for the drag and lift forces expressed in the cut plane. The results are then integrated along the blade to obtain the quantities of global interest. In constrast, Momentum theory [START_REF] Rankine | On the mechanical principles of the action of propellers[END_REF] (also known as disk actuator theory or Axial momentum theory) is a global theory that macroscopically studies the behavior of a fluid column passing through a turbine.

BEM theory then relies on a decomposition of the fluid/turbine system into a macroscopic part via the MT and a local planar part via the BET. The latter considers a radial decomposition for the blades and fluid column (Figure 1.3a), by splitting the rotor area into concentric rings of infinitesimal thickness that do not interact with each other. Glauert's model ultimately links three variables associated with the ring under consideration: the axial induction factor a, the tangential induction factor a and the relative angle deviation ϕ. Given a blade profile and assuming a fixed rotation speed, the resulting equations are

tan ϕ = 1 -a λ(1 + a) , (1.11) a 1 -a = 1 4 sin 2 ϕ • Bc λ 2πr (C L (ϕ -γ λ ) cos ϕ + C D (ϕ -γ λ ) sin ϕ), (1.12) a 1 -a = 1 4λ sin 2 ϕ • Bc λ 2πr (C L (ϕ -γ λ ) sin ϕ -C D (ϕ -γ λ ) cos ϕ). (1.13)
where r denotes the distance between the blade element taken into account (Figure 1.3b) and the rotor, λ = λ(r) is the tip-speed ratio and B the total number of blades.

The chord c λ and twist angle γ λ are given by the blade geometry.

Various modifications have been introduced to the model to consider e.g. the flow around the tip of a blade or a turbulent wake state. For a more detailed account, see Section 4.1.4.

Finally, we evaluate the turbine efficiency using the power coefficient

C p (ϕ) = 8 λ max λmax λ min λ 3 a (1 -a) 1 - C D C L (ϕ -γ λ ) tan -1 ϕ dλ.
Conversely, we can design a blade profile that maximizes this quantity by solving

max ϕ C p (ϕ) s.t (1.11-1.13) α = ϕ -γ λ ,
where α is chosen such that the ratio C D C L is close to zero. In order to simplify this optimization problem, it is commonly assumed a zero drag coefficient, which leads to an explicit solution. To the best of our knowledge, the general case concerning non-zero drag and possible corrections of the model has not been addressed.

Chapter 2

Time-parallelization of sequential data assimilation problems

We present in this chapter a first algorithm to parallelize in time a Luenberger observer. By dividing the time interval into a sequence of data processing windows, we then perform the time-parallelization on each, following a stopping criterion that preserves the exponential rate of convergence of the observer. In the particular case of the Parareal algorithm, we obtain an a priori estimate of the number of parareal iterations required on each window. The efficiency of the procedure is also studied.

We present experiments to confirm the results obtained theoretically.

This is a joint work with Felix Kwok (Hong Kong Baptist University) and Julien Salomon (ANGE, INRIA Paris).

The Luenberger observer

Control theory usually requires total knowledge of the state vector. However, due to certain limitations related to a problem, for instance the number of available measurements, one can often have access only to partial information. A dynamic which fits into this setting is given by

ẋ(t) = Ax(t) + Bu(t), x(0) = x 0 y(t) = Cx(t) (2.1)
where A ∈ M m×m (R), B ∈ M m×p (R) and C ∈ M q×m (R) are assumed to be known (with m, p, q ∈ N * , p, q < m); x ∈ R m is the state vector, y ∈ R q is the measured output, u ∈ R p and t ∈ (0, +∞). The initial condition x(0) = x 0 is unknown.

We are interested in computing an estimate x(t) of x(t), knowing only the imput u(t) and output y(t). To tackle this problem, Luenberger [65, pp.300-307]

proposed the model ẋ(t) = Ax(t) + Bu(t) + L [y(t) -ŷ(t)] , x(0) = x0 ŷ(t) = C x(t) (2.2)
with L ∈ M m×q (R) the observer gain. System (2.2) is known as the Luenberger observer or the Identity observer.

The Luenberger observer

The matrix L needs to be specified, but let us already note that it plays an important role in the estimation error x(t) -x(t). Indeed, substracting (2.1) and (2.2), and then solving the resultant ODE, we get

x(t) -x(t) = e (A-LC)t (x(0) -x(0)) (2.3)
If the eigenvalues of A -LC lie in the open left half-plane {z ∈ C : Re{z} < 0}, then the error will decay to zero. Due to the following result, known as the Identity observer Theorem [65, p.303], we can construct a matrix L such that this property holds:

Theorem 2.1.1. Given a completely observable system (2.1), an identity observer of the form (2.2) can be constructed, and the coefficients of the characteristic polynomial of the observer can be selected arbitrarily.

To be more precise, let us recall that System (2.1) is observable if the rank of the matrix Note that for a single-input single-output system, i.e. p = q = 1, we could determine a unique L ∈ R m by equating the m coefficients of both polynomials in (2.4). However, this approach leads to highly nonlinear equations that are in practice not tractable. Another way to proceed is the Bass-Gura method [START_REF] Bass | High-order system design via state-space considerations[END_REF], which requires the first companion form of A and the coefficients of φ(s). An even more direct method consists in using the Ackermann formula [START_REF] Ackermann | Der entwurf linearer regelungssysteme im zustandsraum[END_REF] for an observable system, to define L by

C :=        C CA . . . CA m-1        is m. Then, given a set (µ i ) i=1,...,m , Theorem 2.1.1 guarantees the existence of L that satisfies det (sI -(A -LC)) = φ(s) (2.
L = φ(A)C -1 (0 • • • 0 1)
which is a consequence of the Cayley-Hamilton Theorem. For its multi-input multioutput extension, see [START_REF] Ackermann | On the synthesis of linear control systems with specified characteristics[END_REF].

Due to Theorem 2.1.1, we obtain Proposition 2.1.2. We assume System (2.1) is observable and the eigenvalues of A -LC are negative and simple. Then, we have e (A-LC)t ≤ γe -µt with µ := min ν∈σ(A-LC) |ν| and γ := cond(V ) = V -1 V , where V is the matrix whose rows are the eigenvectors of A-LC and • represents the induced 2-norm of a matrix.

In particular, combining the latter with Equation (2.3) yields

x(t) -x(t) ≤ γ x(0) -x(0) e -µt .
(2.5)

In practice, the term x(0) -x(0) is unknown, whereas µ is constructed (and known) by the procedure that designs L. Consequently, this equation only provides a rate of convergence for the Luenberger observer.

Time-parallelization setting

In the last decades many techniques have been designed to combine space-time parallelization procedures with data assimilation or optimal control algorithms on bounded time intervals, for instance [START_REF] Trémolet | Parallel algorithms for variational data assimilation and coupling models[END_REF][START_REF] Rao | A time-parallel approach to strong-constraint fourdimensional variational data assimilation[END_REF][START_REF] Cacciapuoti | DD-DA PinT-based model: A domain decomposition approach in space and time, based on parareal, for solving the 4D-Var data assimilation model[END_REF][START_REF] Maday | Monotonic parareal control for quantum systems[END_REF]. In what follows, we propose to extend this notion to the unbounded case and at the same time exploit the exponential rate of convergence of the problem, by an approach that we call the Diamond strategy.

Let us describe briefly our approach. We proceed by partitioning (0, +∞) on intervals of the same length that we call windows. Following a sequential order, we apply a parallel-in-time solver on each, up to some level of accuracy related to a specific stopping criterion. Since our analysis allows us to decompose the estimation error into two terms, corresponding respectively to the Luenberger observer and the parallelization error, we propose a suitable bound on the latter, so that our stopping criterion preserves Luenberger's rate of convergence.

Framework

In order to accelerate the assimilation and take advantage of a time-parallelization procedure, we propose to divide the time interval into windows of a given length T > 0, denoted by

W := (T , T +1 ), ∈ N
where T = •T . Then, we solve Equation (2.2) on each window, in a sequential order, using a time-parallel algorithm.

Let us describe how the latter is applied. Given ∈ N and a fixed window W , we decompose it into N subintervals of length ∆T

W = N -1 n=0 (t n , t n+1 )
with t n = T + n∆T and N ∆T = T . Since time moves forward, parallelizing in this direction requires the introduction of initial conditions Xh

,n on each subinterval, which are obtained from the time-parallelization procedure under consideration. In this setting, the parameter h is related with the accuracy of the method.

Time-parallelization setting

When n = 0, we consider as initial condition Xh 0,0 = x0 in the first window, and Xh ,0 = x (T -) for > 0. Finally, we construct a parallel version of Equation (2.2) in each subinterval (t n , t n+1 ) through the equation

   ẋ (t) = Ax (t) + Bu(t) + L y(t) -C x (t) x (t n + ) = Xh ,n (2.6)
where x (t) denotes the approximation of x(t) obtained by the parallel-in-time solver under consideration.

The Diamond strategy

Imposing initial conditions induces discontinuities at the interfaces t n of the subintervals. In what follows, we call these differences jumps, defined by J h ,n := Xh ,nx (t n -). The success of the parallel method relies on their decay to zero. Note that jumps can be computed without knowing the true solution and hence be used to design an a posteriori estimate of the method.

In this setting, we can express the relation between the solution of (2.1) and its parallel version (2.6) by Lemma 2.2.1. Under the assumptions of Proposition 2.1.2, we have

ε (t n -) ≤ γ e -µt n x(0) -x(0) + e -µ∆T J h ,n-1
where ε (t) = x(t) -x (t) is the error of approximation (2.6) and t n = (N + n)∆T .

Proof. Since this result relies only on the subintervals and does not depend on the windows, we simplify the notation using the sequence (t i ) i∈N , with i = N + n.

Let i ∈ N * . Substracting Equation (2.6) and (2.1), the latter restricted to the sub-interval (t i , t i+1 ), we get

   ε (t) = (A -LC)ε(t) ε (t + i ) = Xh ,i -x(t i )
Its solution is given by ε (t) = e (A-LC)(t-t i ) ε (t + i ). Setting t = t - i+1 in the equation above, we get

ε (t - i+1 ) = e (A-LC)∆T (ε (t + i ) -J h i ) + e (A-LC)∆T J h i (2.7)
On the other hand, from the definition of jump, we have

ε (t + i ) -J h i = x (t - i ) -x(t i ) = e (A-LC)∆T e -(A-LC)∆T (x (t - i ) -x(t i )) = e (A-LC)∆T (x (t - i-1 ) -x(t i-1 )) = e (A-LC)∆T (ε (t + i-1 ) -J h i-1 ).
Applying the previous formula recursively in (2.18) leads to

ε (t - i+1 ) = e (A-LC)(i+1)∆T ε (0 + ) + e (A-LC)∆T J h i ,
since J h 0 = 0. The result follows from taking the norm and using Proposition 2.1.2. We recall now that our approach aims to preserve Luenberger's rate of convergence. Thanks to the previous lemma, we can achieve this goal by defining on each window W the stopping criterion max

1≤n≤N J h ,n ≤ γe -µ T ,
|ν| ∆t := T N = 0 repeat T = T (∀n ∈ {1, . . . , N }) t n = T + n∆T if = 0 then Xh ,0 = x0 else Xh ,0 = x (T -) end /* Using a generic time-parallelization procedure (GTP): */ Construct the initial conditions { Xh ,n } n=0,...,N -1 = GTP( Xh ,0 ) J h ,n := Xh ,n -x (t n -) Assign ← + 1 until max 1≤n≤N J h ,n ≤ γe -µ T Proposition 2.

Let us assume that h is obtained from the stopping criterion in

W max 1≤n≤N J h ,n ≤ γe -µ T (2.8)
where γ is an arbitrary parameter. Then, the rate of convergence of x (t) to x(t) is bounded by µ, i.e.

ε (t n -) ≤ γe -µ∆T ( x(0) -x(0) + γ) e -µ T .

Parallelization

Independent to the choice of the parallel-in-time solver, Algorithm 2.1 is well-defined since the jumps are computed a posteriori. However, by specifying it, we can study in more detail the stopping criterion presented in (2.8) and the complexity of the overall procedure. In order to obtain an a priori estimate for the required accuracy h on each window and the efficiency of the Diamond strategy, throughout this section we consider the Parareal algorithm as the time-parallel method.

The Parareal algorithm

Introduced by Lions, Maday and Turinici [START_REF] Lions | Résolution d'EDP par un schéma en temps «pararéel[END_REF], the goal of the Parareal algorithm is to solve evolution problems by partitioning a bounded time interval into subintervals, and after assigning to each of them a processor, updates iteratively the initial conditions and solve a series of independent and smaller problems in parallel, reducing the time-computation of the solution.

To be more precise, given the problem

u(t) = f (u(t)), t ∈ [0, T ] u(0) = u 0 (2.9)
we divide [0, T ] in M subintervals denoted by (t n-1 , t n ). Then we consider two solvers F and G, that compute a fine and a coarse numerical approximation of u. The former is computationally expensive and consequently restricted to solve initial-value problems with high accuracy in each subinterval (t n-1 , t n ), whereas the latter is faster and can be used for solving on large intervals as [0, T ]. For an arbitrary initial condition ũ given in t = t n-1 , we denote these local approximations of u(t n ) by F(t n , t n-1 , ũ) and G(t n , t n-1 , ũ), respectively.

In this framework, (u(t n )) n=1,...,M is approximated by (U k n ) n=1,...,M , which is build as follows: Algorithm 2.2: Parareal algorithm Input: u 0 , T, M, Tol Output:

(t n ) n=1,...,M , (U k * n ) n=1,...,M ∆t := T N , t 0 = 0 /* Initialization of the initial conditions */ U 0 0 = u 0 for 1 ≤ n ≤ M do t n = n∆T U 0 n = G(t n , t n-1 , U 0 n-1 ) end k = 0 repeat U k 0 = u 0 for 1 ≤ n ≤ M do U k+1 n = F(t n , t n-1 , U k n-1 ) + G(t n , t n-1 , U k+1 n-1 ) -G(t n , t n-1 , U k n-1 ) (2.10) end J k n := U k n -u(t - n ) Assign k ← k + 1 until max 1≤n≤M J k n ≤ Tol k * = k -1
Notice that the superscript k in Algorithm 2.2 plays the role of the parameter h, introduced in the previous section.

As a remark, Gander and Vandewalle [START_REF] Gander | Analysis of the parareal time-parallel timeintegration method[END_REF] showed that the parareal algorithm can be presented as a multi-shooting algorithm applied to (2.9). Solving the multiple shooting equations with the Newton's method yields

U k+1 n = u n-1 (t n , U k n-1 ) + ∂u n-1 ∂U n-1 (t n , U k n-1 )(U k+1 n-1 -U k n-1 ),
where u n-1 (t n , U k n-1 ) denotes the exact solution of (2.9) at t n , with initial condition U k n-1 at t n-1 . If we approximate the exact solution using the fine solver and the Jacobian term by a difference on a coarse grid, we obtain Equation (2.10).

Concerning its convergence, they proved that it requires at most M iterations. Moreover, after k iterations the algorithm is exact on the first k subintervals. An improvement of their result, due to Gander and Hairer [START_REF] Gander | Nonlinear convergence analysis for the parareal algorithm[END_REF], assumes that the coarse solver must satisfy

G(t n , t n-1 , y) -G(t n , t n-1 , z) ≤ (1 + C 2 ∆T ) y -z ,
for a positive constant C 2 . However, this result does not cover the decaying case, i.e. when the Lipschitz constant is smaller than one. Since we are interested in coupling this algorithm with the Luenberger observer and take advantage of its decaying behavior, we provide a result adapted to this case, which follows from [START_REF] Gander | Multigrid interpretations of the parareal algorithm leading to an overlapping variant and MGRIT[END_REF].

Theorem 2.3.1 (Convergence of Parareal for decaying problems). Let F(t n , t n-1 , U k n-1 ) be the exact solution on the time subdomain (t n-1 , t n ) and let G(t n , t n-1 , U k n-1 ) be a coarse integrator, such that the local truncation error τ (t n , z)

:= F(t n , t n-1 , z) - G(t n , t n-1 , z) satisfies for all z τ (t n , z) ≤ α, τ (t n , y) -τ (t n , z) ≤ β y -z ,
where α, β > 0 are constants. We also assume that F and G are Lipschitz with respect to the initial conditions:

max { F(t n , t n-1 , y) -F(t n , t n-1 , z) , G(t n , t n-1 , y) -G(t n , t n-1 , z) } ≤ ε y -z ,
for a constant ε ∈ (0, 1).

Then, after k iterations of the Parareal algorithm, the error U k n -u(t n ) is bounded above by E k n , defined as

E k n =      0 n ≤ k αβ k n-k-1 i=0 k+i k ε i n > k.
(2.11)

Proof. From the definition of the Parareal algorithm in (2.10) and since F is the exact solution in (t n-1 , t n ), we obtain, after adding and substracting G(t n , t n-1 , u(t n-1 ))

U k n -u(t n ) = F(t n , t n-1 , U k-1 n-1 ) + G(t n , t n-1 , U k n-1 ) -G(t n , t n-1 , U k-1 n-1 ) -F(t n , t n-1 , u(t n-1 )) = τ (t n , U k-1 n-1 ) -τ (t n , u(t n-1 )) + G(t n , t n-1 , U k n-1 ) -G(t n , t n-1 , u(t n-1 )),
and then, taking norms and using the first assumption on τ , combined with the Lipschitz condition on G, leads to

U k n -u(t n ) ≤ β U k-1 n-1 -u(t n-1 ) + ε U k n-1 -u(t n-1 ) .
Moreover, the error in the initial guess can be estimated in the same way, although using the coarse integration U 0 n = G(t n , t n-1 , U 0 n-1 ), which gives

U 0 n -u(t n ) ≤ α + ε U 0 n-1 -u(t n-1 ) .
Therefore, an upper bound E k n for U k n -u(t n ) satisfies the recurrence relation

E k n = βE k-1 n-1 + εE k n-1 , (2.12) E 0 n = α + εE 0 n-1 , (2.13)
with E k 0 = 0 for all k. We solve this recurrence using a generating function, by defining the formal power series

ρ k (ζ) = n≥1 E k n ζ n .
Multiplying (2.12) and (2.13) by ζ n and summing over n ≥ 1 gives

ρ k (ζ) = βζρ k-1 (ζ) + εζρ k (ζ), ρ 0 (ζ) = αζ 1 -ζ + εζρ 0 (ζ),
and then, solving by induction yields the explicit formula

ρ k (ζ) = αβ k ζ k+1 (1 -ζ)(1 -εζ) k+1 .
Expanding ρ k (ζ) in a power series leads to

ρ k (ζ) = αβ k ζ k+1 i≥0 ζ i j≥0 k+j k (εζ) j = αβ k ζ k+1 n≥0 n i=0 k+i k ε i ζ n = n≥0 αβ k n i=0 k+i k ε i ζ n+k+1 .
Then, for n ≤ k we have E k 0 = . . . = E k k = 0; whereas for n > k, we obtain

E k n = αβ k n-k-1 i=0 k+i k ε i .

Combination with Luenberger observer

For a fixed window W , we approximate x (t n ) using the sequence ( Xk ,n ) k=1,...,n , constructed with the recursive formula

   Xk ,n = F(t n , t n-1 , Xk-1 ,n-1 ) + G(t n , t n-1 , Xk ,n-1 ) -G(t n , t n-1 , Xk-1 ,n-1 ) X0 ,n = G(t n , t n-1 , X0 ,n-1 ), X0 ,0 = x (T -).
(2.14)

In the framework of the Diamond strategy, jumps are computed explicitly and by Proposition 2.2.2, we can determine the number of iterations k obs := k (required to guarantee a Luenberger's rate of convergence of Algorithm 2.1) a posteriori. Instead, we propose to combine the stopping criterion with estimate (2.11) to derive an a priori upper bound of k obs , denoted by k , which is given by our next theorem.

Theorem 2.3.2. We keep the assumptions of Proposition 2.1.2 and Theorem 2.3.1. For a window W and γ > 0, we define

k =    min S S = ∅ k -1 S = ∅ (2.15)
where

S = k ∈ N * , k ≤ N -1 : β k N -k-1 i=0 k+i k ε i ≤ γ K e -µ T , K = α 2 1 -ε N 1 -ε . ( 2
max 1≤n≤N J k ,n ≤ αβ k n-k -1 i=0 k +i k ε i • max 1≤n≤N X0 ,n -x (t n ) ≤ α γ K e -µ T • max 1≤n≤N X0 ,n -x (t n ) (2.17)
The term X0

,n -x (t n ) can be bounded as well thanks to Theorem 2.3.1, and then we have

X0 ,n -x (t n ) ≤ α n-1 i=0 ε i = α 1 -ε n 1 -ε ≤ α 1 -ε N 1 -ε , since ε < 1.
Thus, from (2.17) and the definition of K in (2.16), we get

max 1≤n≤N J k ,n ≤ γe -µ T .
Since we introduce an approximation Xk ,n of the parallel solution x (t n ), we provide a bound for the estimation error between the solution x(t n ) of (2.1) and the parareal sequence.

Complexity analysis

Given a tolerance parameter Tol, we define the efficiency of the algorithm by

E = τ s N τ p (2.19)
where τ s is the CPU time required to reach the tolerance Tol using a sequential solver. For a parallel solver, we denote this quantity by τ p , and N represents the number of available processors (and hence, subintervals).

In what follows, we make a distinction between theoretical and observed efficiency. The former comes from the analysis of our algorithm developped in the previous section, whereas the latter is a result of empirical measurements obtained when using the a posteriori estimator (2.19).

The efficiency of the Diamond strategy, combined with the Parareal algorithm, is estimated by our next theorem. Theorem 2.3.4. Let Tol be the tolerance parameter on the error in Corollary 2.3.3. We define * as the number of windows required to reach this tolerance, i.e. * = min ∈ N : γ(1 + e -µ∆T ) + γe -µ∆T x(0) -x(0) e -µ T ≤ Tol . (2.20) Then, the estimated efficiency of Algorithm 2.1 is given by

E th = * τ F ∆T τ F ∆T + N τ G ∆T * -1 =0 k -1 , (2.21)
where τ F ∆T , τ G ∆T represents the amount of time spent in solving (2.2) over an interval of size ∆T with F and G, respectively.

Proof. The decay rate of the Parareal algorithm is provided by Corollary 2.3.3. For a fixed length window T , we define * by (2.20). Then, the CPU time for the parallel solver is given by

τ p = * -1 =0 k N τ F ∆T T ∆T + k τ G ∆T T ∆T = τ F ∆T + N τ G ∆T * -1 =0 k
On the other hand, we consider a sequential solver on the interval [0, * T ], which requires * T ∆T iterations. Assuming F as the solver, its CPU time is given by

τ s = τ F ∆T • * T ∆T = τ F ∆T • * N.
Finally, from Definition (2.19) of efficiency we get (2.21).

Numerical experiments

The present section is devoted to some numerical experiments for the Luenberger observer. For this purpose, we use

A = 0 1 -1 -2 , B = 0 1 , C = 0 1 , u(t) = 3 + 0.5 sin(0.75t).
We remark that the initial condition on System (2.1) is unknown, but we perform the experiments with x(0) = (0 0) . We then construct the observer x(t) by setting as initial condition x(0) = (2 1) and the eigenvalues of A -LC. For the latter, we consider {-0.8, -1} and {-0.2, -0.25} as possible choices.

To introduce the parareal procedure, given N available processors, we set

T = 1, δT = ∆T = T N , Tol = 10 -8 ,
where δT denotes the time step associated with G, chosen as a one step solver for the sake of simplicity. We use the Backward Euler method to define both propagators F and G.

Diagonalized system

We recall that the essential part of Theorem 2.3.1 is the contracting factor ε. For the Luenberger observer (2.2), we have

ε = max [I -δt(A -LC)] -∆T /δt , [I -∆T (A -LC)] -1 .
where δt is the time step associated with F, assumed to be constant. Even if we choose the eigenvalues of A -LC to guarantee a decaying rate of convergence, ε is not necessarily smaller than one. For this reason, we consider instead a diagonalized observer ż(t) = Dẑ(t) + V -1 (Bu(t) + Ly(t))

ẑ(0) = V -1 x0 (2.22)
where ẑ = V -1 x and

D = V -1 (A -LC)V .
Due to the change of variables, γ = 1. We determine the constants α, β and ε by Proposition 2.4.1. Let F and G be defined by the Backward Euler scheme, with time steps δt and δT , respectively. We assume that ∆T K ≤ 1 and Equation (2.22) satisfies

M := sup (ẑ,t) Dẑ + V -1 (Bu(t) + Ly(t)) < ∞ K := max D , sup t>0 V -1 (B u(t) + L ẏ(t)) < ∞.
Then, the constants associated with both propagators in Theorem 2.3.2 are given by

α = ∆T 2 K(M + 1) 2(1 -∆T K) , β = [I -δtD] -∆T /δt -[I -∆T D] -1 , (2.23) ε = max [I -δtD] -∆T /δt , [I -∆T D] -1 . (2.24)
The proof is standard, but we present it for the sake of completeness.

Proof. Let {t n } N n=0 be a regular partition of the interval [0, T ], with ∆T = T /N . Given ẑn-1 an approximation of ẑ(t n-1 ), we recall that the Backward Euler method applied to (2.22) is given by ẑnẑn-1 ∆T = f (ẑ n , t n ) where f (s, t) = Ds + V -1 (Bv(t) + Ly(t)).

Since δt is assumed to be constant, we then define F by F(t n , t n-1 , ẑn-1 ) = (I -δtD) -∆T /δt δtV -1 (Bv(t n-1 ) + Ly(t n-1 )) + ẑn-1 and then, a direct calculation yields

F(t n , t n-1 , y) -F(t n , t n-1 , z) = (I -δtD) -∆T /δt (y -z). (2.25)
On the other hand, G is defined as a one-step solver, which allows us to replace δt by ∆T in the previous expressions to obtain To bound the local truncation error, we proceed as follows. Starting at the exact solution z n-1 = ẑ(t n-1 ), we define z n = G(t n , t n-1 , z n-1 ) and then

G(t n , t n-1 , y) -G(t n , t n-1 , z) = (I -∆tD) -1 (y -z). ( 2
τ (t n , z n-1 ) = F(t n , t n-1 , z n-1 ) -G(t n , t n-1 , z n-1 ) = ẑ(t n ) -z n
since F is an exact solver. We use that z n = ẑ(t n-1 ) + ∆T f (z n , t n ) and then expand ẑ(t n-1 ) around t n to get

τ (t n , z n-1 ) = ∆T ż(ẑ(t n ), t n ) -f (z n , t n ) - (∆T ) 2 2 z(ẑ(ξ), ξ) (2.27)
where ξ ∈ (t n-1 , t n ). Since ż = f (z, t), we can get rid of the derivatives of z. In particular, the definition of f (s, t) shows that

ż(ẑ(t n ), t n ) -f (z n , t n ) = f (ẑ(t n ), t n ) -f (z n , t n ) = Dτ (t n , z n-1 ), z(ẑ(ξ), ξ) = df dt (ẑ(ξ), ξ) = ∂f ∂s (ẑ(ξ), ξ) • f (ẑ(ξ), ξ) + ∂f ∂t (ẑ(ξ), ξ) = Df (ẑ(ξ), ξ) + V -1 (B v(ξ) + L ẏ(ξ)).
Replacing these expressions in (2.27) and rearranging terms yields

τ (t n , z n-1 ) = - (∆T ) 2 2 (I -∆T D) -1 Df (x(ξ), ξ) + V -1 (B v(ξ) + L ẏ(ξ)) .
Finally, assuming that ∆T K < 1, we take norms and use the definitions of K and M to obtain α.

Evolution of k

As a first experiment, since the jumps involved in Equation (2.8) allows us to compute the sequence k obs = {k obs } , we propose to compare its behavior with its a priori estimate k th = {k } , provided by Theorem 2.3.2.

We observe in Figure 2.1 that increasing γ leads to enlarge the number of windows in which the algorithm requires only 1 iteration. This is expected, due to the term γe -µ T present in Proposition 2.2.2 and Theorem 2.3.2. We also notice an asymptotical behavior of both sequences in Figure 2.1a, which it is hidden after, because the total number of windows increases at a lower rate.

Observed efficiency

Our second experiment consists in comparing the observed efficiencies for both sequences k obs and k th , using different values of γ, N and δt. To evaluate E obs , the execution time for the parallel and sequential solvers was computed with the functions tic and toc of MATLAB (version 9.4.0.813654 (R2018a)).

As we notice previously, increasing γ improves the algorithm performance, but the difference between E obs (k obs ) and E obs (k th ) still remains, as observed in Figure 2.2a. In Figure 2.2c, the gap between these values varies slightly, showing that δt small enough does not affect the efficiency significantly. Increasing the number of processors N makes this difference smaller and also improves the efficiency of the algorithm, as shown in Figure 2.2b. Another way to narrow this gap is choosing smaller eigenvalues for A -LC. As Figure 2.2 suggests, the comparison between {-0.8, -1} and {-0.2, -0.25} shows that E obs (k th ) increases, whereas E obs (k obs ) becomes smaller.

E obs (k th ) ≤ E obs (k obs ),
which is simply because the sequence k th underperforms k obs . However, in Figure 2.2b (left) we note at some point the opposite behavior, which could be explained as follows: since β decreases dramatically as N increases, it leads to underestimate the number of iterations. This suggests a relation of the sort γ ≤ f (µ, ∆T, δt).

Finally, we recall that k th is useful to determine the theoretical efficiency (2.21).

Assuming that τ G

∆T is negligible, we denote it by

E th 0 = * * -1 =0 k -1
.

with * given by (2.20). We note that this value predicts quite well E obs (k th ) in all the tests.

Perspectives

The key element of our strategy is the design of a stopping criterion for the parallelin-time solver that preserves the Luenberger rate of convergence on each assimilation window. Nonetheless, we believe that several questions remains to be answered, starting with using other time-parallelization algorithms, possible extensions to Kalman filters and considering a variable window size.

Another known parallel-in-time algorithm is ParaExp [START_REF] Gander | ParaExp: A parallel integrator for linear initial-value problems[END_REF]. This method, proposed by Gander and Güttel, decomposes a linear initial-value problem into two subproblems that can be solved in parallel, to then superpose their respective solutions and retrieve the original one. In contrast to Parareal, which often shows a slow convergence when applied to hyperbolic systems, it is well suited for solving these kinds of problem. Its main feature is the approximation of the matrix exponential [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF] and consequently, to combine it with our strategy, we propose to determine the number of terms required for a given approximation on each window.

On the other hand, an extension to a stochastic framework should consider a continuous Kalman filter, which is known for being an optimal a posteriori estimator, in the sense that minimizes the variance of the state estimation error. However, an explicit upper bound for this quantity is not often provided. A second point to take into account is that time parallel algorithms are designed for deterministic differential equations and not stochastic ones. The parallelization of the latter was first addressed by Bal [START_REF]Parallelization in time of (stochastic) ordinary differential equations[END_REF], followed by applications to e.g. chemical kinetics [START_REF] Engblom | Parallel in time simulation of multiscale stochastic chemical kinetics[END_REF] or finance [START_REF] Pagès | The parareal algorithm for american options[END_REF].

Finally, the present work proposes to use an assimilation window of fixed lenght and then determine the number of iterations. The opposite idea, i.e. fixing the latter in a window and determining its size, could be useful when a large number of processors is available.

Chapter 3

Bathymetry optimization

The current chapter focuses on the determination of a bathymetry from an optimization perspective. We consider a PDE-constrained optimization problem in which the wave motion, dependent on the bathymetry and modeled by the weak formulation of the Helmholtz equation, acts as constraint. The cost functional is assumed to be general.

We begin by deriving the Helmholtz equation from the Navier stokes system. Afterwards, a C 0 -bound for its weak solution is obtained. Concerning the optimization problem, we study its well-posedness and continuity of the control-to-state mapping. The discrete optimization problem is also addressed, studying the convergence to the discrete optimal solution as well as the convergence of a finite element approximation. We illustrate our results by solving numerically two examples that describe wave damping and bathymetry reconstruction. This is a joint work with Pierre-Henri Cocquet (Université de Pau et des Pays de l'Adour) and Julien Salomon (ANGE, INRIA Paris).

Derivation of the wave model

We start from the Navier-Stokes equations to derive the governing PDE. However, due to its complexity, we introduce two approximations [START_REF] Méhauté | An Introduction to Hydrodynamics and Water Waves[END_REF]: a small relative depth (Long wave theory) combined with an infinitesimal wave amplitude (Small amplitude wave theory). An asymptotic analysis on the relative depth shows that the vertical component of the depth-averaged velocity is negligible, obtaining the Saint-Venant equations. After neglecting its convective inertia terms and linearizing around the sea level, it results a wave equation which depends on the bathymetry. Since a variable sea bottom can be seen as an obstacle, we reformulate the equation as a Scattering problem involving the Helmholtz equation.

From Navier-Stokes system to Saint-Venant equations

For t ≥ 0, we define the time-dependent region

Ω t = {(x, z) ∈ Ω × R | -z b (x) ≤ z ≤ η(x, t)}
where Ω is a bounded open set with Lipschitz boundary, η(x, t) represents the water level and -z b (x) is the bathymetry (also called bottom topography), a time independent and negative function. The water height is denoted by h

= η + z b . x z Free surface η(x, t) -z b (x) h Bottom
In what follows, we consider an incompressible fluid of constant density (assumed to be equal to 1), governed by the Navier-Stokes system

           ∂u ∂t + (u • ∇) u = div (σ T ) + g in Ω t , div (u) = 0 in Ω t , u = u 0 in Ω 0 , (3.1)
where u = (u, v, w) denotes the velocity of the fluid, g = (0, 0, -g) is the gravity and σ T is the total stress tensor, given by σ T = -pI + µ ∇u + ∇u with p the pressure and µ the coefficient of viscosity.

To complete (3.1), we require suitable boundary conditions. Given the outward normals

n s = 1 1 + |∇η| 2 -∇η 1 , n b = 1 1 + |∇z b | 2 ∇z b 1 ,
to the free surface and bottom, respectively, we recall that the velocity of the two must be equal to that of the fluid:

     ∂η ∂t -u • n s = 0 on (x, η(x, t), t), u • n b = 0 on (x, -z b (x), t). (3.2)
On the other hand, the stress at the free surface is continuous, whereas at the bottom we assume a no-slip condition

σ T • n s = -p a n s on (x, η(x, t), t), (σ T n b ) • t b = 0 on (x, -z b (x), t), ( 3.3)
with p a the atmospheric pressure and t b an unitary tangent vector to n b .

A long wave theory approach

For the sake of completeness and following the standard procedure described in [START_REF] Bristeau | Derivation of a non-hydrostatic shallow water model; comparison with Saint-Venant and Boussinesq systems[END_REF][START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF][START_REF] Sainte-Marie | Vertically averaged models for the free surface Euler system. derivation and kinetic interpretation[END_REF], we derive the Saint-Venant equations from the Navier-Stokes system. For simplicity of presentation, system (3.1) is restricted to two dimensions, but a more detailed derivation of the three-dimensional case can be found in [START_REF] Decoene | Asymptotic derivation of the section-averaged shallow water equations for river hydraulics[END_REF].

Since our analysis focuses on the shallow water regime, we introduce the parameter

ε = H L ,
where H denotes the relative depth and L is the characteristic dimension along the horizontal axis. The importance of the nonlinear terms is represented by the ratio

δ = A H ,
with A the maximum vertical amplitude.

We use the change of variables

x = x L , z = z H , t = C 0 L t, and u = u δC 0 , w = w δεC 0 , η = η A , z b = z b H , p = p gH ,
where C 0 = √ gH is the characteristic dimension for the horizontal velocity. Assuming the viscosity and atmospheric pressure constants, we define their respective dimensionless versions by

µ = µ C 0 L , p a = p a gH .
Dropping primes after rescaling, the dimensionless system (3.

-δµ ∂u ∂z + ε 2 ∂w ∂x ∂z b ∂x 2 + ε 2δµ ∂w ∂z -p ∂z b ∂x = 0. (3.10)
To derive the Saint-Venant equations, we use an asymptotic analysis in ε. In addition, we assume a small viscosity coefficient µ = εµ 0 .

A first simplification of the system consists in deriving an explicit expression for p, known as the hydrostatic pressure. After rearranging the terms of order ε 2 in (3.5) and integrating in the vertical direction, we get

p(x, z, t) = O(ε 2 δ) + (δη -z) + εδµ 0 ∂u ∂x + 2 ∂w ∂z - ∂u ∂x (x,
∂x + u(x, δη, t) • w(x, δη, t) -u(x, -z b , t) • w(x, -z b , t) = ∂(h δ u 2 ) ∂x + u(x, δη, t) ∂η ∂t 1 + (εδ) 2 ∂η ∂x 2 ,
and then, the vertical integration of the left-hand side of (3.4) brings

δη -z b δ ∂u ∂t + δ 2 u ∂u ∂x + w ∂u ∂z dz = δ ∂(h δ u) ∂t + δ 2 ∂(h δ u 2 ) ∂x + δ 2 u(x, δη, t) ∂η ∂t 1 + (εδ) 2 ∂η ∂x 2 -1 .
To deal with the term h δ u 2 , we start from (3.12) which shows that ∂p ∂x = O(δ). Plugging this expression into (3.4) yields

∂ 2 u ∂z 2 = O(ε).
From boundary conditions (3.8) and (3.10) we obtain

∂u ∂z (x, δη, t) = O(ε 2 ), ∂u ∂z (x, z b , t) = O(ε).
Consequently, u(x, z, t) = u(x, 0, t) + O(ε) and then u(x, z, t) -u(x, t) = O(ε). Hence, we have the approximation

h δ u 2 = h δ u 2 + δη -z b (u -u) 2 dz = h δ u 2 + O(ε 2 )
and finally

δη -z b δ ∂u ∂t + δ 2 u ∂u ∂x + w ∂u ∂z dz = δ ∂(h δ u) ∂t + δ 2 ∂(h δ u 2 ) ∂x + O(ε 2 δ 2 ) (3.13) + δ 2 u(x, δη, t) ∂η ∂t 1 + (εδ) 2 ∂η ∂x 2 -1 .
We then integrate the right-hand side of Equation (3.4)

δη -z b - ∂p ∂x + δ µ 0 ε ∂ ∂z ∂u ∂z + εδµ 0 2 ∂ ∂x ∂u ∂x + ∂ ∂z ∂w ∂x dz = -δh δ ∂η ∂x + O(εδ) + δ µ 0 ε ∂u ∂z (x, δη, t) - µ 0 ε ∂u ∂z (x, -z b , t) .
Combining this expression with (3.13), we get the vertical integration of the momentum equation.

In summary, we have the system

∂η ∂t 1 + (εδ) 2 ∂η ∂x 2 + ∂(h δ u) ∂x = 0 (3.14) ∂(h δ u) ∂t + δ ∂(h δ u 2 ) ∂x = -h δ ∂η ∂x + O(ε) (3.15) + µ 0 ε ∂u ∂z (x, δη, t) - µ 0 ε ∂u ∂z (x, -z b , t) + δu(x, δη, t) ∂η ∂t 1 + (εδ) 2 ∂η ∂x 2 -1 .
If δ = O(1) and ε → 0, then we recover the classical derivation of the one-dimensional Saint-Venant equations. The convergence of (3.15) is guaranteed by the boundary equations (3.8) and (3.10), from which we get

µ 0 ε ∂u ∂z (x, δη, t) = O(εδ), µ 0 ε ∂u ∂z (x, -z b , t) = O(ε).

Small amplitudes

With respect to the classical Saint-Venant formulation, passing to the limit δ → 0 is equivalent to neglect the convective acceleration terms and linearize the system (3.14-3.15) around the sea level η = 0. In order to do so, we rewrite the derivatives as

∂(h δ u) ∂t = h δ ∂u ∂t + δ ∂η ∂t u, ∂(h δ u) ∂x = δ ∂(ηu) ∂x + ∂(z b u) ∂x ,
and then, taking ε, δ → 0 in (3.14-3.15) yields

         ∂η ∂t + ∂(z b u) ∂x = 0, - ∂(z b u) ∂t + z b ∂η ∂x = 0.
Finally, after deriving the first equation with respect to t and replacing the second into the new expression, we obtain the wave equation for a variable bathymetry.

All the previous computations hold for the three-dimensional system (3.1). In this case, we obtain

∂ 2 η ∂t 2 -div (gz b ∇η) = 0. (3.16)

Helmholtz formulation

Equation (3.16) defines a time-harmonic field, whose solution has the form η(x, t) = Re{ψ tot (x)e -iωt }, where the amplitude ψ tot satisfies

ω 2 ψ tot + div (gz b ∇ψ tot ) = 0. (3.17)
We wish to rewrite the equation above as a scattering problem. Since a variable bottom z b (x) := z 0 + δz b (x) (with z 0 a constant describing a flat bathymetry and δz b a perturbation term) can be considered as an obstacle, we thus assume that δz b has a compact support in Ω and ψ tot satisfies the so-called Sommerfeld radiation condition. In a bounded domain as Ω, we impose the latter thanks to an impedance boundary condition (also known as first-order absorbing boundary condition), which ensures the existence and uniqueness of the solution [75, p.108]. We then reformulate (3.17) as div

((1 + q)∇ψ tot ) + k 2 0 ψ tot = 0 in Ω, ∇(ψ tot -ψ 0 ) • n -ik 0 (ψ tot -ψ 0 ) = 0 on ∂Ω, (3.18)
where q(x) := δz b (x) z 0 is compactly supported in Ω, k 0 := ω √ gz 0 , n is the unit normal to ∂Ω and ψ 0 (x) = e ik 0 x• d is an incident plane wave propagating in the direction d (such that | d| = 1).

Decomposing the total wave as ψ tot = ψ 0 + ψ sc , where ψ sc represents an unknown scattered wave, we obtain the Helmholtz formulation div ((1 + q)∇ψ sc ) + k 2 0 ψ sc = -div (q∇ψ 0 ) in Ω,

∇ψ sc • n -ik 0 ψ sc = 0 on ∂Ω. (3.19)
Its structure will be useful to prove the existence of a minimizer for a PDE-constrained functional, as discussed in the next section.

Description of the optimization problem

We are interested in studying the problem of a cost functional constrained by the weak formulation of a Helmholtz equation. The latter intends to generalize the equations considered so far, whereas the former indirectly affects the choice of the set of admissible controls. These can be discontinuous since they are included in the space of functions of bounded variations. In this framework, we treat the continuity and regularity of the associated control-to-state mapping, and the existence of an optimal solution to the optimization problem.

Weak formulation

Let Ω ⊂ R 2 be a bounded open set with Lipschitz boundary. We consider the following general Helmholtz equation

-div ((1 + q)∇ψ) -k 2 0 ψ = div (q∇ψ 0 ) in Ω, (1 + q)∇ψ • n -ik 0 ψ = g -q∇ψ 0 • n on ∂Ω, (3.20)
where g is a source term. We assume that q ∈ L ∞ (Ω) and that there exists α > 0 such that a.e. x ∈ Ω, 1 + q(x) ≥ α.

(3.21) Remark 3.2.1. Here we have generalized the models described in the previous section: if q has a fixed compact support in Ω, we have that the total wave ψ tot satisfying (3.18) is a solution to (3.20) with g = ∇ψ 0 • nik 0 ψ 0 and no volumic right-hand side; whereas the scattered wave ψ sc satisfying (3.19) is a solution to (3.20) with g = 0. All the results obtained in this broader setting still hold true for both problems.

A weak formulation for (3.20) is given by a(q; ψ, φ) = b(q; φ), ∀φ ∈ H 1 (Ω), (3.22) where a(q; ψ, φ) :=

Ω (1 + q)∇ψ • ∇φ -k 2 0 ψφ dx -ik 0 ∂Ω ψφ dσ, b(q; φ) := - Ω q∇ψ 0 • ∇φ dx + g, φ H -1/2 ,H 1/2 .
Note that, thanks to Cauchy-Schwarz inequality, the sesquilinear form a is continuous

|a(q; ψ, φ)| ≤ C(Ω, q, α)(1 + q L ∞ (Ω) ) ψ 1,k 0 φ 1,k 0 , ψ 2 1,k 0 := k 2 0 ψ 2 L 2 (Ω) + α ∇ψ 2 L 2 (Ω) ,
where C(Ω, q, α) > 0 is a generic constant. In addition, taking φ = ψ in the definition of a, it satisfies a Garding inequality

Re{a(q; ψ, ψ)} + 2k 2 0 ψ 2 L 2 (Ω) ≥ ψ 2 1,k 0 , ( 3

Continuous optimization problem

We are interested in solving the next PDE-constrained optimization problem min

(q,ψ)∈U Λ ×H 1 (Ω) J(q, ψ) s.t. (3.22). (3.24)
We now define the set of admissible q. We wish to find optimal q that can have discontinuities and we thus cannot look for q in some Sobolev spaces that are continuously embedded into C 0 (Ω), even if such regularity is useful for proving existence of a minimizer (see e.g. [8, Chapter VI], [START_REF] Bastide | Penalization model for Navier-Stokes-Darcy equation with application to porosity-oriented topology optimization[END_REF]Theorem 4.1]). To be able to find optimal q satisfying (3.21) and having possible discontinuities, we follow [START_REF] Chen | An augmented Lagrangian method for identifying discontinuous parameters in elliptic systems[END_REF] and introduce the set

U Λ = {q ∈ BV (Ω) | α -1 ≤ q(x) ≤ Λ a.e. x ∈ Ω} .
Above Λ ≥ max{α -1, 0} and BV (Ω) is the set of functions with bounded variations [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF] whose distributional gradient belong to the set M b (Ω, R N ) of bounded Radon measures. Note that piecewise constant functions over Ω belong to U Λ .

Some useful properties of BV (Ω) can be found in [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF] and are recalled below for the sake of completeness. This is a Banach space for the norm (see [5, p.120

, Proposition 3.2]) q BV (Ω) := q L 1 (Ω) + |Dq|(Ω),
where D is the distributional gradient and

|Dq|(Ω) = sup Ω q div (ϕ) dx ϕ ∈ C 1 c (Ω, R 2 ) and ϕ L ∞ (Ω) ≤ 1 ,
is the variation of q (see [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]p.119,Definition 3.4]).

The weak * convergence in BV (Ω), denoted by

q n q, weak * in BV (Ω), means that q n → q in L 1 (Ω) and Dq n Dq in M b (Ω, R N ).
Also, in a two-dimensional setting, the continuous embedding BV (Ω) ⊂ L 1 (Ω) is compact. We finally recall that the application q ∈ BV (Ω) → |Dq|(Ω) ∈ R + is lower semi-continuous with respect to the weak * topology of BV . Hence, for any sequence q n q in BV (Ω), we have

|Dq|(Ω) ≤ lim inf n→+∞ |Dq n |(Ω).
The set U Λ is a closed, weakly * closed and convex subset of BV (Ω). However, since its elements are not necessarily bounded in the BV -norm, in order to prove the existence of a minimizer to Problem (3.24) we need to add a penalizing distributional gradient term to the cost functional J(q, ψ). To avoid this situation, we introduce the set of admissible parameters

U Λ,κ = {q ∈ U Λ | |Dq|(Ω) ≤ κ}
which also possesses the aforementioned properties. The choice between these two sets also affects the convergence analysis of the discrete optimization problem, topic discussed in Section 3.4.

Remark 3.2.2.

In this chapter, we are interested in computing either the total wave satisfying (3.18) or the scattered wave solution to Equation (3.19). Since this require to work with q having a fixed compact support in Ω, we also introduce the following set of admissible parameters

U ε := {q ∈ U |q(x) = 0 for a.e x ∈ O ε } , O ε = {x ∈ Ω | dist(x, ∂Ω) ≤ ε} ,
which is a set of bounded functions with bounded variations that have a fixed support in Ω. We emphasize that it is a convex, closed and weak- * closed subset of BV (Ω). As a result, all the theorems we are going to prove also hold for this set of admissible parameters.

Continuity of the control-to-state mapping

This section is devoted to prove the continuity of the application q ∈ U → ψ(q) ∈ H 1 (Ω) where ψ(q) satisfies Problem (3.22). We assume that U ⊂ BV (Ω) is weakly * closed and ∀q ∈ U, a.e. x ∈ Ω, α -1 ≤ q(x) ≤ Λ.

Note that both U Λ , U Λ,κ and U ε (see Remark 3.2.2) also satisfy these two assumptions. The next result consider the dependance of the stability constant with respect to the optimization parameter q.

Theorem 3.2.3. Assume that q ∈ U . Then there exists a constant C s (k 0 ) > 0 that does not depend on q such that

ψ 1,k 0 ≤ C s (k 0 ) sup φ 1,k 0 =1 |a(q; ψ, φ)|, (3.25)
where the constant C s (k 0 ) > 0 only depend on the wavenumber and on Ω. In addition, the solution to (3.22) satisfies the bound The proof of (3.25) proceed by contradiction assuming the above inequality is false. Therefore, we suppose there exist sequences (q n ) n ⊂ U and (ψ n

ψ 1,k 0 ≤ C s (k 0 )C(Ω) max{k -1 0 , α -1/2 } q L ∞ (Ω) ∇ψ 0 L 2 (Ω) + g H -1/2 (∂Ω) , ( 3
) n ⊂ H 1 (Ω) such that q n BV (Ω) ≤ M , ψ n 1,k 0 = 1 and lim n→+∞ sup φ 1,k 0 =1 |a(q n ; ψ n , φ)| = 0. (3.27)
The compactness of the embeddings BV (Ω) ⊂ L 1 (Ω) and H 1 (Ω) ⊂ L 2 (Ω) yields the existence of a subsequence (still denoted (q n , ψ n )) such that

ψ n ψ ∞ in H 1 (Ω), ψ n → ψ ∞ in L 2 (Ω) and q n → q ∞ ∈ U in L 1 (Ω). (3.28)
Compactness of the trace operator gives that lim n→+∞ ψ n | ∂Ω = ψ ∞ | ∂Ω strongly in L 2 (∂Ω) and thus, from (3.28), we get

lim n→+∞ Ω k 2 0 ψ n φ dx + ik 0 ∂Ω ψ n φ dσ = Ω k 2 0 ψ ∞ φ dx + ik 0 ∂Ω ψ ∞ φ dσ, ∀ v ∈ H 1 (Ω), lim n→+∞ Ω ∇ψ n • ∇φ dx = Ω ∇ψ ∞ • ∇φ dx.
It only remains to pass to the limit in the term involving q n . We start from

(q n ∇ψ n , ∇φ) L 2 (Ω) -(q ∞ ∇ψ ∞ , ∇φ) L 2 (Ω) = ((q n -q ∞ )∇ψ n , ∇φ) L 2 (Ω) + (q ∞ ∇(ψ n -ψ ∞ ), ∇φ) L 2 (Ω) ,
and use Cauchy-Schwarz inequality to get

Ω q n ∇ψ n •∇φ dx - Ω q ∞ ∇ψ ∞ • ∇φ dx ≤ ((q n -q ∞ )∇ψ n , ∇φ) L 2 (Ω) + (q ∞ ∇(ψ n -ψ ∞ ), ∇φ) L 2 (Ω) ≤ |q n -q ∞ |∇φ L 2 (Ω) |q n -q ∞ |∇ψ n L 2 (Ω) + (q ∞ ∇(ψ n -ψ ∞ ), ∇φ) L 2 (Ω) ≤ 2 √ Λ √ α ψ n 1,k 0 |q n -q ∞ |∇φ L 2 (Ω) + (∇(ψ n -ψ ∞ ), q ∞ ∇φ) L 2 (Ω) .

Description of the optimization problem

The right term above goes to 0 owning to q ∞ ∈ L ∞ (Ω) and (3.28). For the other term, since q n → q ∞ strongly in L 1 , we can extract another subsequence (q n k ) k such that q n k → q ∞ pointwise almost everywhere in Ω. Also, |q n -q ∞ ||∇φ| 2 ≤ 2 √ Λ|∇φ| 2 ∈ L 1 (Ω) and the Lebesgue's Dominated Convergence Theorem then yields

lim k→+∞ |q n k -q ∞ |∇φ L 2 (Ω) = 0.
This gives that (see also [21,

Equation (2.4)]) lim k→+∞ (q n k ∇ψ n k , ∇φ) L 2 (Ω) = (q ∞ ∇ψ ∞ , ∇φ) L 2 (Ω) , ∀φ ∈ H 1 (Ω). ( 3
n k ; ψ n k , φ) = a(q ∞ , ψ ∞ , φ), ∀φ ∈ H 1 (Ω),
and the uniqueness result [START_REF] Graham | Stability and error analysis for the Helmholtz equation with variable coefficients[END_REF]Theorems 2.1,2.4] shows that ψ ∞ = 0 thus the whole sequence actually converges to 0. To get our contradiction, it remains to show that ∇ψ n L 2 (Ω) converges to 0 as well. From the Garding inequality (3.23), we have

ψ n 2 1,k 0 ≤ Re{a(q n ; ψ n , ψ n )} + 2k 2 0 ψ n 2 L 2 (Ω) ----→ n→+∞ 0,
where we used (3.27) and the strong L 2 convergence of ψ n towards ψ ∞ = 0. Finally, we get lim n→+∞ ψ n 1,k 0 = 0 which contradicts ψ n 1,k 0 = 1 and gives the desired result.

Applying then (3.25) to the solution to (3.22) finally yields

ψ 1,k 0 ≤ C s (k 0 ) sup φ 1,k 0 =1 |a(q; ψ, φ)| ≤ C s (k 0 ) sup φ 1,k 0 =1 |b(q; φ)| ≤ C s (k 0 ) sup φ 1,k 0 =1 q L ∞ (Ω) ∇ψ 0 L 2 (Ω) ∇φ L 2 (Ω) + g H -1/2 (∂Ω) φ H 1/2 (∂Ω) ≤ C s (k 0 )C(Ω) max{k -1 0 , α -1/2 } q L ∞ (Ω) ∇ψ 0 L 2 (Ω) + g H -1/2 (∂Ω)
, where C(Ω) > 0 comes from the trace inequality.

Remark 3.2.4. Let us consider a more general version of Problem (3.20), given by

-div ((1 + q)∇ψ) -k 2 0 ψ = F in Ω, (1 + q)∇ψ • n -ik 0 ψ = G on ∂Ω.
We emphasize that the estimation of the stability constant C s (k 0 ) with respect to the wavenumber have been obtained for (F, G) ∈ L 2 (Ω) × L 2 (∂Ω) for q = 0 in [START_REF] Hetmaniuk | Stability estimates for a class of Helmholtz problems[END_REF] and for q ∈ Lip(Ω) satisfying (3.21) in [START_REF] Barucq | Stability analysis of heterogeneous Helmholtz problems and finite element solution based on propagation media approximation[END_REF][START_REF] Graham | Stability and error analysis for the Helmholtz equation with variable coefficients[END_REF][START_REF] Graham | The Helmholtz equation in heterogeneous media: a priori bounds, well-posedness, and resonances[END_REF]. Since their proofs rely on Green, Rellich and Morawetz identities, they do not extend to the case (F, G) ∈ (H 1 (Ω)) ×H -1/2 (∂Ω) but such cases can be tackled as it is done in [START_REF] Esterhazy | On stability of discretizations of the Helmholtz equation[END_REF]p.10,Theorem 2.5]. The case of Lipschitz q has been studied in [START_REF] Brown | Multiscale Petrov-Galerkin method for high-frequency heterogeneous Helmholtz equations[END_REF]. As a result, the dependance of the stability constant with respect to q, in the case q ∈ U and (F, G) ∈ (H 1 (Ω)) × H -1/2 (∂Ω), does not seem to have been tackled so far to the best of our knowledge. Remark 3.2.5 (H 1 -bounds for the total and scattered waves). From Remark 3.2.1, we obtain that the total wave ψ tot and the scattered wave ψ sc are solutions to (3.22), with respective right-hand sides

b tot (q; φ) = ∂Ω (∇ψ 0 • n -ik 0 ψ 0 )φ dσ, b sc (q; φ) = - Ω q∇ψ 0 • ∇φ dx.
As a result of Theorem 3.2.3 and the continuity of the trace, we have

ψ tot 1,k 0 ≤ C(Ω)C s (k 0 )k 0 max{k -1 0 , α -1/2 }, ψ sc 1,k 0 ≤ C s (k 0 )α -1/2 q L ∞ (Ω) ∇ψ 0 L 2 (Ω) ≤ k 0 C s (k 0 )α -1/2 q L ∞ (Ω) |Ω|.
We can now prove some regularity for the control-to-state mapping. Theorem 3.2.6. Let (q n ) n ⊂ U be a sequence satisfying q n BV (Ω) ≤ M and whose weak * limit in BV (Ω) is denoted by q ∞ . Let (ψ(q n )) n be the sequence of weak solution to Problem (3.22). Then ψ(q n ) converges strongly in H 1 (Ω) towards ψ(q ∞ ). In other words, the mapping

q ∈ (U Λ , weak * ) → ψ(q) ∈ (H 1 (Ω), strong), is continuous.
Proof. Since q n BV (Ω) ≤ M and (q n ) n ⊂ U , there exists q ∞ such that q n q ∞ , weak * in BV (Ω). Using that U is weak * closed, we obtain that q ∞ ∈ U . Note that the sequence (ψ(q n )) n of solution to Problem (3.22) satisfies estimate (3.26) uniformly with respect to n. As a result, the convergences (3.28) hold and, from (3.29) and the unicity of the limiting problem, we get that a(q n ; ψ n , φ) → a(q ∞ ; ψ ∞ , φ). Since b(q n , φ) → b(q ∞ , φ) for all φ ∈ H 1 (Ω), this proves that a(q ∞ ; ψ ∞ , φ) = b(q; φ) for all φ ∈ H 1 (Ω), ψ(q ∞ ) is a weak solution to (3.22) and ψ(q n ) ψ(q ∞ ) in H 1 (Ω).

We now show that ψ(q n ) → ψ(q ∞ ) strongly in H 1 . To see this, we start by noting that, up to extracting a subsequence (still denoted q n ), we can use (3.29) to get that

lim n→+∞ b(q n ; ψ(q n )) = b(q ∞ ; ψ(q ∞ )).
Since ψ(q n ), ψ(q ∞ ) satisfy the variational problem (3.22), we infer

lim n→+∞ a(q n ; ψ(q n ), ψ(q n )) = a(q ∞ ; ψ ∞ , ψ(q ∞ )), (3.30)
where the whole sequence actually converges owing to the uniqueness of the limit. Using then that ψ(q n ) ψ(q ∞ ) in H 1 (Ω) together with (3.30), we get

√ 1 + q n ∇ψ(q n ) 2 L 2 (Ω) = a(q n ; ψ(q n ), ψ(q n )) + k 0 ψ(q n ) 2 L 2 (Ω) + ik 0 ψ(q n ) 2 L 2 (∂Ω) ----→ n→+∞ a(q ∞ ; ψ(q ∞ ), ψ(q ∞ )) + k 0 ψ(q ∞ ) 2 L 2 (Ω) + ik 0 ψ(q ∞ ) 2 L 2 (∂Ω) = √ 1 + q ∞ ∇ψ(q ∞ ) 2 L 2 (Ω)
.

Boundedness/Continuity of solution to Helmholtz problem

We prove in this section that, even if the parameter q is not smooth enough for the solution to (3.20) to be in H s (Ω) for some s > 1, we can still have continuous solution.

In order to prove such regularity for ψ, we are going to rely on the De Giorgi-Nash-Moser theory [

               Lv := div (A(x)∇v) = f 0 - N j=1 ∂f j ∂x j , ∇v • n = h + N j=1 f j n j , (3.31)
where

A ∈ L ∞ (Ω, R N ×N ) satisfy the standard ellipticity condition A(x)ξ • ξ ≥ γ|ξ| 2 for a.e. x ∈ Ω. Let p > N and assume that f 0 ∈ L p/2 (Ω), f j ∈ L p (Ω) for all j = 1, • • • , N and h ∈ L p-1 (∂Ω). Then the weak solution v to (3.31) satisfies v C 0 (Ω) ≤ C(N, p, Ω, γ)   v L 2 (Ω) + f 0 L p/2 (Ω) + N j=1 f j L p (Ω) + h L p-1 (∂Ω)   .

C 0 -bound for the general Helmholtz problem

Using Theorem 3.3.1, we can prove some L ∞ bound for the weak solution to Helmholtz equation with bounded coefficients.

Theorem 3.3.2. Assume that q ∈ L ∞ (Ω) and satisfies (3.21) and g ∈ L 2 (∂Ω). Then the solution to Problem (3.22) satisfies

ψ C 0 (Ω) ≤ C(Ω) C s (k 0 , α) q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g L 2 (∂Ω) , (3.32)
where andC(Ω) > 0 does not depend on k nor q.

C s (k 0 , α) = 1 + (1 + k 2 0 )k -1 0 + α -1/2 max{k -1 0 , α -1/2 }C s (k 0 ),
Proof. We cannot readily apply Theorem 3.3.1 to the weak solution of Problem (3.20) since the latter involves a complex valued operator. We therefore consider the Problem satisfied by ν = Re{u} and ζ = Im{u} which is given by

             -div ((1 + q)∇ν) -k 2 0 ν = div (q∇ Re{ψ 0 }) in Ω, -div ((1 + q)∇ζ) -k 2 0 ζ = div (q∇ Im{ψ 0 }) in Ω, (1 + q)∇ν • n = Re{g} -k 0 ζ -q∇ Re{ψ 0 } • n, on ∂Ω, (1 + q)∇ζ • n = Im{g} + k 0 ν -q∇ Im{ψ 0 } • n on ∂Ω. (3.33)
Since Problem (3.33) is equivalent to Problem (3.20), we get that the weak solution (ν, ζ) ∈ H 1 (Ω) to (3.33) satisfies the inequality (3.26). Assuming that g ∈ L 2 (∂Ω) and using the continuous Sobolev embedding H 1 (Ω) ⊂ L 6 (Ω), the (compact) embedding H 1/2 (∂Ω) ⊂ L 2 (∂Ω), that q ∈ L ∞ (Ω) satisfies (3.21) and the fact that ψ 0 is smooth we get the next regularities

f 0,1 = k 2 0 ν ∈ L 6 (Ω), f j,1 = q ∂ Re{ψ 0 } ∂x j ∈ L ∞ (Ω), h 1 = Re{g} -k 0 ζ ∈ L 2 (∂Ω), f 0,2 = k 2 0 ζ ∈ L 6 (Ω), f j,2 = q ∂ Im{ψ 0 } ∂x j ∈ L ∞ (Ω), h 2 = Im{g} + k 0 ν ∈ L 2 (∂Ω).
Applying now Theorem 3.3.1 to (3.33) twice with p = 3 and N = 2, we get C 0 bounds for ν and ζ

ν C 0 (Ω) ≤ C(2, 3, Ω, γ)   ν L 2 (Ω) + f 0,1 L 3/2 (Ω) + 2 j=1 f j,1 L 3 (Ω) + h 1 L 2 (∂Ω)   , ζ C 0 (Ω) ≤ C(2, 3, Ω, γ)   ζ L 2 (Ω) + f 0,2 L 3/2 (Ω) + 2 j=1 f j,2 L 3 (Ω) + h 2 L 2 (∂Ω)   .
Some computations with the Holder and multiplicative trace inequalities then give

( ν L 2 (Ω) + ζ L 2 (Ω) ) ≤ 2 ψ L 2 (Ω) , f 0,1 L 3/2 (Ω) + f 0,2 L 3/2 (Ω) ≤ k 2 0 ψ L 3/2 (Ω) ≤ |Ω| 1/6 k 2 0 ψ L 2 (Ω) , f j,l L 3 (Ω) ≤ q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) , j = 1, 2, h 1 L 2 (∂Ω) + h 2 L 2 (∂Ω) ≤ g L 2 (∂Ω) + k 0 ψ L 2 (∂Ω) , ≤ g L 2 (∂Ω) + k 0 C(Ω) ψ L 2 (Ω) ψ H 1 (Ω) .

Using then Young inequality yields

k 0 ψ L 2 (Ω) ψ H 1 (Ω) ≤ C ψ H 1 (Ω) + k 2 0 ψ L 2 (Ω) ≤ C ∇ψ L 2 (Ω) + (1 + k 2 0 ) ψ L 2 (Ω)
where C > 0 is a generic constant. We obtain the bound

ψ C 0 (Ω) = ν C 0 (Ω) + ζ C 0 (Ω) ≤ C(Ω) 1 + k 2 0 ψ L 2 (Ω) + ∇ψ L 2 (Ω) + q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g L 2 (∂Ω) .
Using the definition of ψ 1,k 0 on the estimate above, we have

ψ C 0 (Ω) ≤ C(Ω) (1 + k 2 0 )k -1 0 + α -1/2 ψ 1,k 0 + q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g L 2 (∂Ω) (3.34)
To apply the a priori estimate (3.26), we recall that the H -1/2 norm can be replaced by a L 2 norm (since g ∈ L 2 (∂Ω)) and then,

ψ 1,k 0 ≤ C(Ω) max{k -1 0 , α -1/2 }C s (k 0 ) q L ∞ (Ω) ∇ψ 0 L 2 (Ω) + g L 2 (∂Ω) ≤ C(Ω) max{k -1 0 , α -1/2 }C s (k 0 ) max{1, |Ω|} × q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g L 2 (∂Ω)
Finally, combining the last expression with Equation (3.34), we obtain that the weak solution to the Helmholtz equation satisfies

ψ C 0 (Ω) ≤ C(Ω) 1 + (1 + k 2 0 )k -1 0 + α -1/2 max{k -1 0 , α -1/2 }C s (k 0 ) × q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g L 2 (∂Ω) ,
where C(Ω) > 0.

Remark 3.3.3.

1. For the one-dimensional Helmholtz problem, the a priori estimate (3.26) and the continuous embedding H 1 (I) ⊂ C 0 (I) directly gives the continuity of u over

I ψ C 0 (I) ≤ C ψ 1,k 0 ≤ C(k) q L ∞ (Ω) ∇ψ 0 L ∞ (Ω) + g H -1/2 (∂Ω) .
It is worth noting that we do not need to assume that g ∈ L 2 (∂Ω).

For the two-dimensional

Helmholtz problem with q = 0, we can get the above C 0 estimate from the embedding H 2 (Ω) → C 0 (Ω) since

ψ C 0 (Ω) ≤ C ψ H 2 (Ω) ,
for a generic constant C. We can then see that the estimate (3.32) actually has the same dependance with respect to k 0 as the H 2 -estimate [49, p.677, Proposition 3.6].

C 0 -bounds for the total and scattered waves

Thanks to Remark 3.2.1 and following the proof of Theorem 3.3.2, these bounds can be roughly obtained by setting g = ∇ψ 0 • nik 0 ψ 0 and omitting the L ∞ -norms in (3.34) for the total wave ψ tot , and simply by setting g = 0 in the case the scattered wave ψ sc . Using after the H 1 -bounds from Remark 3.2.5, we actually get

ψ tot C 0 (Ω) ≤ C(Ω)k 0 (1 + k 2 0 )k -1 0 + α -1/2 max{k -1 0 , α -1/2 }C s (k 0 ) + 1 ψ sc C 0 (Ω) ≤ C(Ω)k 0 (1 + k 2 0 )k -1 0 + α -1/2 α -1/2 C s (k 0 ) + 1 q L ∞ (Ω) .
We emphasize that the previous estimates show that the scattered wave ψ sc vanishes in Ω if q → 0. This is expected since, if q = 0, there is no obstacle to scatter the incident wave, which amount to saying that ψ tot = ψ 0 .

Discrete optimization problem

This section is devoted to the finite element discretization of the optimization problem (3.24). We consider a quasi-uniform family of triangulations (see [30, p.76, Definition 1.140]) {T h } h>0 of Ω and the corresponding finite element spaces

V h = φ h ∈ C(Ω) | φ h | T ∈ P 1 (T ), ∀T ∈ T h .
Note that thanks to Theorem 3.3.2, the solution to the general Helmholtz equation (3.20) is continuous, which motivates to use continuous piecewise linear finite elements. We are going to look for discrete optimal design that belong to some finite element spaces K h and we thus introduce the following set of discrete admissible parameters

U h = U ∩ K h .
The full discretization of the optimization problem (3.24) then reads

Find q * h ∈ U h such that J(q * h ) ≤ J(q h ), ∀q h ∈ U h , (3.35)
where J(q h ) = J(q h , ψ h (q h )) is the reduced cost-functional and

ψ h := ψ h (q h ) ∈ V h satisfies the discrete Helmholtz problem a(q h ; ψ h , φ h ) = b(q h ; φ h ), ∀φ h ∈ V h . ( 3.36)
The existence of solution to Problem (3.36) follows from uniqueness since we are in a finite dimensional setting [46, Theorems 2.1, 2.4].

Before giving the definition of K h , we would like to discuss briefly the strategy for proving that the discrete optimal solution converges toward the continuous ones. To achieve this, we need to pass to the limit in inequality (3.35). Since J is only lower-semi-continuous with respect to the weak * topology of BV , we can only pass to the limit on one side of the inequality and the continuity of J is then going to be needed to pass to the limit on the other side to keep this inequality valid as h → 0. We discuss first the case U = U Λ for which Theorem 3.2.7 gives the existence of optimal q but only if β > 0. Since we have to pass to the limit in (3.35), we need that lim h→0 |Dq h |(Ω) = |Dq|(Ω). Since the total variation is only continuous with respect to the strong topology of BV , we have to approximate any q ∈ U Λ by some q h ∈ U h such that lim h→0 q -q h BV (Ω) = 0.

However, from [9, p.8, Example 4.1] there is an example of a BV -function v that cannot be approximated by piecewise constant function v h over a given mesh in such a way that lim h→0 |Dv h |(Ω) = |Dv|(Ω). Nevertheless, if one consider an adapted mesh that depends on a given function v ∈ BV (Ω)∩L ∞ (Ω), we get the existence of piecewise constant function on this specific mesh that strongly converges in BV toward v (see [START_REF] Bȇlík | Approximation by piecewise constant functions in a BV metric[END_REF]p.11,Theorem 4.2]). As a result, when considering U = U Λ , we use the following discrete set of admissible parameters

K h,1 = {q h ∈ L ∞ (Ω) | q h | T ∈ P 1 (T ), ∀T ∈ T h } .
Note that, from Theorem [18, p.10, Theorem 4.1 and Remark 4.2], the set U h = U Λ ∩ K h,1 defined above has the required density property hence motivated its introduction as a discrete set of admissible parameter.

In the case, U = U Λ,κ , we will not need the density of U h in the strong topology of BV but only for the weak * topology. The discrete set of admissible parameter is then going to be

U h = U Λ,κ ∩ K h,0 with K h,0 = {q h ∈ L ∞ (Ω) | q h | T ∈ P 0 (T ), ∀T ∈ T h } .
We show below the convergence of discrete optimal solution to the continuous one for both cases highlighted above.

Convergence of the Finite element approximation

We prove here some useful approximations results for any U h defined above. We have the following convergence result whose proof can be found in [START_REF] Esterhazy | On stability of discretizations of the Helmholtz equation[END_REF]p.22,Lemma 4.1] (see also [START_REF] Graham | Stability and error analysis for the Helmholtz equation with variable coefficients[END_REF]p.10,Theorem 4.1]). Theorem 3.4.1. Let q h ∈ U h and ψ(q h ) ∈ H 1 (Ω) be the solution to the variational problem a(q h ; ψ(q h ), φ) = b(q h , φ), ∀φ ∈ H 1 (Ω).

Let S * : (q h , f ) ∈ U h × L 2 (Ω) → S * (q h , f ) = ψ * ∈ H 1 (Ω)
be the solution operator associated with the following problem

Find ψ * ∈ H 1 (Ω) such that a(q h ; φ, ψ * ) = (φ, f ) L 2 (Ω) , ∀φ ∈ H 1 (Ω).
Denote by C a the continuity constant of the bilinear form a(q h ; •, •), which does not depend on h since q h ∈ U h , and define the adjoint approximation property by

δ(V h ) := sup f ∈L 2 (Ω) inf φ h ∈V h S * (q h , f ) -φ h 1,k 0 f L 2 (Ω) .
Assume that the spaces V h satisfies

2C a k 0 δ(V h ) ≤ 1, (3.37)
then the solution ψ h (q h ) to Problem (3.36) satisfies

ψ(q h ) -ψ h (q h ) 1,k 0 ≤ 2C a inf φ h ∈V h ψ(q h ) -φ h 1,k 0 .
In the case q ∈ C 0,1 (Ω) where Ω is a convex Lipschitz domain, the assumption (3.37) has been discussed in [START_REF] Graham | Stability and error analysis for the Helmholtz equation with variable coefficients[END_REF]p.11,Theorem 4.3] and, broadly speaking, (3.37) holds if k 2 0 h is small enough. Since the proof rely on H 2 -regularity for a Poisson problem, we cannot readily extend the argument here since we can only expect to have ψ ∈ H 1 (Ω) and that S * also depend on the meshsize. We can still show that (3.37) is satisfied for small enough h. Lemma 3.4.2. Assume that q h ∈ U h weak * converges toward q ∈ BV (Ω). Then (3.37) is satisfied for small enough h.

Proof. Note first that Theorem 3.2.6 also holds for the adjoint problem and thus

lim h→0 S * (q h , f ) -S * (q, f ) 1,k 0 = 0.
Using the density of smooth function in H 1 and the properties of the piecewise linear interpolant [30, p.66, Corollary 1.122], we have that

lim h→0   sup f ∈L 2 (Ω) inf φ h ∈V h S * (q, f ) -φ h 1,k 0 f L 2 (Ω)   = 0,
and thus a triangular inequality shows that (3.37) holds for small enough h.

We can now prove a discrete counterpart to Theorem 3.2.6.

Theorem 3.4.3. Let (q h ) h ⊂ U h be a sequence satisfying q h BV (Ω) ≤ M and whose weak * limit in BV (Ω) is denoted by q. Let (ψ h (q h )) h be the sequence of discrete solutions to Problem (3.36). Then ψ(q h ) converges, as h goes to 0, strongly in H 1 (Ω) towards ψ(q) satisfying Problem (3.22).

Proof. For h small enough, Lemma 3.4.2 ensures that (3.37) holds and a triangular inequality then yields

ψ h (q h ) -ψ(q) 1,k 0 ≤ ψ h (q h ) -ψ(q h ) 1,k 0 + ψ(q h ) -ψ(q) 1,k 0 ≤ 2C a inf φ h ∈V h ψ(q h ) -φ h 1,k 0 + ψ(q h ) -ψ(q) 1,k 0 ≤ (1 + 2C a ) ψ(q h ) -ψ(q) 1,k 0 + 2C a inf φ h ∈V h ψ(q) -φ h 1,k 0 .
Theorem 3.2.6 gives that the first term above goes to zero as h → 0. For the second one, we can use the density of smooth function in H 1 to get that it goes to zero as well.

Convergence of the discrete optimal solution

We are now in position to prove the convergence of discrete optimal design toward continuous one in the case

U = U Λ , U h = U Λ ∩ K h,1 .
Hence the set of discrete control is composed of piecewise linear function on T h . Theorem 3.4.4. Assume that (A1) -(A2) -(A3) from Theorem 3.2.7 hold and that the cost function (3.35). Then the sequence (q * h ) h ⊂ U Λ is bounded and there exists q * ∈ U Λ such that q * h q * weakly * in BV (Ω), ψ(q * h ) → ψ(q * ) strongly in H 1 (Ω) and

J 0 : (q, ψ) ∈ U Λ × H 1 (Ω) → J 0 (q, ψ) ∈ R is continuous with respect to the (weak * , strong) topology of BV (Ω) × H 1 (Ω). Let (q * h , ψ h (q * h )) ∈ U Λ,h × V h be an optimal pair of
J(q * ) ≤ J(q), ∀q ∈ U Λ .
Hence any accumulation point of (q * h , ψ h (q * h )) is an optimal pair for Problem (3.24).

Proof. Let q Λ ∈ U Λ,h be given as

q Λ (x) = Λ, ∀x ∈ Ω.
Then Dq Λ = 0. Since ψ h (q Λ ) is well-defined and converges toward ψ(q Λ ) strongly in H 1 (see Theorem 3.4.4), we have that

J(q Λ ) = J(q Λ , ψ h (q Λ )) = J 0 (q Λ , ψ h (q Λ )) --→ h→0 J 0 (q Λ , ψ(q Λ )).
As a result, using that (q * h , ψ h (q * h )) is an optimal pair to Problem (3.36), we get that

β|D(q * h )|(Ω) ≤ -J 0 (q * h , ψ h (q * h )) + J(q Λ , ψ h (q Λ )) ≤ -m + J 0 (q Λ , ψ h (q Λ )),
and thus the sequence (q * h ) h ⊂ U Λ,h ⊂ U Λ is bounded in BV (Ω) uniformly with respect to h. We denote by q * ∈ U Λ its weak * limit and Theorem 3.4.3 then shows that ψ h (q * h ) → ψ(q * ) strongly in H 1 (Ω). The lower semi-continuity of J ensures that

J(q * , ψ(q * )) = J(q * ) ≤ lim inf h→0 J(q * h ) = lim inf h→0 J(q * h , ψ h (q * h )).
Now, let q ∈ U Λ , using the density of smooth function in BV , we get that there exists a sequence q h ∈ U Λ,h such that q h -q * BV (Ω) → 0 (see also [START_REF] Bartels | Total variation minimization with finite elements: convergence and iterative solution[END_REF]p.10,Remark 4.2]). From Theorem 3.4.3, we get ψ h (q h ) → ψ(q) strongly in H 1 (Ω) and the continuity of J ensure that J(q h ) → J(q). Since J(q * h ) ≤ J(q h ) for all q h ∈ U Λ,h , we get by passing to the inf-limit that

J(q * ) ≤ lim inf h→0 J(q * h ) ≤ lim inf h→0 J(q h ) = J(q), ∀q ∈ U Λ ,
and the proof is therefore finished.

When U = U Λ,κ , U h = U Λ,κ ∩ K h,0 ,
the set of discrete control is composed of piecewise constant function on T h that satisfy

∀q h ∈ U h , q h BV (Ω) ≤ 2 max(Λ, κ, |α -1|).
We can compute explicitly the previous norm by integrating by parts the total variation (see e.g. [9, p.7, Lemma 4.1]). This reads

∀q h ∈ U h , |Dq h |(Ω) = F ∈F i |F ||[q h ]| F |,
where F i is the set of interior faces and

|[q h ]| F = |q h | T 1 -q h | T 2 |
is the jump of q h on the interior face F = ∂T 1 ∩ ∂T 2 . Note then that any q h ∈ U h can only have either a finite number of discontinuity or jumps that are not too large.

Theorem 3.4.5. Assume that β = 0 and (A2) -(A3) from Theorem 3.2.7 hold and that the cost function J : (q, ψ) (3.35). Then the sequence (q * h ) h ⊂ U Λ,κ is bounded and there exists q * ∈ U Λ,κ such that q * h q * weakly * in BV (Ω), ψ(q * h ) → ψ(q * ) strongly in H 1 (Ω) and J(q * ) ≤ J(q), ∀q ∈ U Λ .

∈ U Λ × H 1 (Ω) → J(q, ψ) ∈ R is continuous with respect to the (weak * , weak) topology of BV (Ω) × H 1 (Ω). Let (q * h , ψ h (q * h )) ∈ U h × V h be an optimal pair of
Hence any accumulation point of (q * h , ψ h (q * h )) is an optimal pair for Problem (3.24).

Proof. Since (q * h ) h belong to U h , it satisfies q h BV (Ω) ≤ 2 max(Λ, κ, |α -1|) and is thus bounded uniformly with respect to h. We denote by q * ∈ U Λ,κ its weak * limit. Theorem 3.4.4 then shows that ψ h (q * h ) converges strongly in H 1 (Ω) toward ψ(q * ). Now, let q ∈ U Λ,κ , using the density of smooth function in BV , one gets that there exists a sequence q h ∈ U h such that q h q weak * in BV (Ω) (see also [START_REF] Bartels | Total variation minimization with finite elements: convergence and iterative solution[END_REF]Introduction]). From Theorem 3.4.3, we get ψ h (q h ) → ψ(q) strongly in H 1 (Ω) and the continuity of J ensure that J(q h ) → J(q). The proof can then be done as in Theorem 3.4.4.

Numerical experiments

In this section, we tackle numerically the optimization problem (3.24), when it is constrained to the total amplitude ψ tot described by (3.18). We focus on two examples: a damping problem, where the computed bathymetry optimally reduces the magnitude of the incoming waves; and an inverse problem, in which we recover the bathymetry by minimizing the mismatch between the observed and predicted magnitude of the waves.

In what follows, we consider an incident plane wave ψ 0 (x) = e ik 0 x• d propagating in the direction d = (0 1) , with

k 0 = ω 0 √ gz 0 , ω 0 = 2π T 0 , T 0 = 20, g = 9.81, z 0 = 3.
For the space domain, we set Ω = [0, L] 2 , where L = 10π k 0 . We also impose a L ∞constraint on the variable q, namely that q ≥ -0.9.

We observe in Figure 3.1 that the optimal bathymetry we obtain is highly oscillating. In our experiments, this oscillation remained at every level of space discretization we have tested. This could be related to the fact that in all our results, q ∈ BV (Ω). Note also that the damping is more efficient over Ω q . This fact is coherent with the results of the next experiment. Chapter 4

Mathematical analysis of the Blade element momentum theory

The aim of this chapter is to determine suitable conditions for both the existence of solutions and the numerical solving of the Blade element momentum theory. We recall first the main features of Glauert's modeling, which is followed by the mathematical analysis of the questions above. The key element in our approach is the reformulation of the original set of equations into a single expression, which allows us to split the model into two parts associated with the blade geometry and fluid-turbine dynamic. Last of all, we study the optimization problem concerning turbine efficiency.

This is a joint work with Jérémy Ledoux (Hydrotube Énergie) and Julien Salomon (ANGE, INRIA Paris), via the Grant ANR HyFloEFlu (Hydroliennes flottantes et énergie, ANR-10-IEED-0006-04).

The Blade element momentum theory

We present the model proposed by Glauert to describe the interaction between a turbine and a flow. After having introduced the relevant variables, we recall the main steps of the reasoning leading to the governing equations. We then detail the simple and complex versions of the model that we will consider.

Variables

Glauert's theory aims at establishing algebraic relations that represent the interaction between a stream and several rotating blades, named turbine in what follows. In order to do this, it couples two models: a macroscopic one that describes the evolution of fluids rings crossing the turbine; and a local model that characterizes the behavior of a planar section of a blade (a blade element) for various angles of attack and, in some cases, for various Reynolds numbers.

The stream is supposed to be horizontal, constant in time and incompressible. The latter assumption implies that the stream velocities in the left and right neighborhoods of the turbine are equal. We denote them by U 0 and by U -∞ and U +∞ the upstream and downstream velocities, respectively.

As BEM models do not take into account interactions between blade elements, w.l.o.g. we restrict ourselves to a fixed blade element and a constant rotation speed Ω, i.e. a fixed value of the parameter

λ := Ωr U -∞ ,
where Ω denotes the rotation speed of the blades and r is the distance of the element under consideration to the rotation axis. It should also be noted that in practical cases, the turbine works at constant λ: Ω is indeed often controlled through the torque exerted by a generator in such a way that the ratio Ω U -∞ is kept constant for various values of U -∞ . It follows that the value of λ associated with one element only depends on r. In the sequel, we consequently only use the variable λ to describe the location of a blade element.

Macroscopic variables and unknowns

Glauert's model ultimately consists in linking three variables a, a and ϕ associated with the ring of fluid under consideration. Among these three unknowns, the axial induction factor a and the tangential induction factor a are defined by

a = U -∞ -U 0 U -∞ , ( 4.1)
a = ω 2Ω , ( 4.2)
where ω is the rotation speed of the ring of fluid at location of the turbine. The angle ϕ is the relative angle deviation [68, p.120] of the ring so that

tan ϕ = 1 -a λ(1 + a ) . ( 4.3)
For the sake of simplicity, and to emphasize their role of unknowns in Glauert's model, we omit in this paper the dependence of a, a , ϕ (and α in what follows) on λ in the notations.

Local variables

Given a blade profile, the lift and drag coefficients C L and C D are defined through the relations

dL = C L (α) ρ 2 U 2 rel c λ dr, dD = C D (α) ρ 2 U 2 rel c λ dr,
where ρ is the mass density of the fluid, dL and dD are the elementary lift and drag forces applying on a blade element of thickness dr and chord c λ = c λ(r) , and U rel is the relative fluid speed (also called apparent fluid speed) perceived from this blade element while rotating, that is

U rel = U 0 sin ϕ . (4.4)
The parameter α is called angle of attack and is defined as the angle between the chord and flow direction. It satisfies the relation The coefficients C L and C D correspond to the ratio of lift and drag forces with respect to the dynamic force, which is associated with the observed kinetic energy. They are assumed to be known, both depending on α and in some cases of the Reynolds number. In the following, we do not take into account their dependence on the latter, however, all the results established in this chapter can be extended without difficulty to that framework. On the contrary to α, the Reynolds number is known since it depends on α as well as U -∞ , r, Ω, c λ (see [87, p.374]), variables that can be controlled e.g. in a wind tunnel.

α = ϕ -γ λ , ( 4
Thought varying from one profile to another, the variations of C L and C D with respect to α can be described qualitatively in a general way. The coefficient C L usually increases linearly with respect to α up to a given critical angle α s , with 0 < α ≤ π/2, where the so-called stall phenomenon occurs: C L then decreases rapidly (see, e.g., [15, pp.93-94]), causing a sudden loss of lift, so that angles of attack larger that α s are not desirable. Since C D is associated with a drag force, it is always positive and defined for all angles in concrete cases. Its usual behavior consists in slightly increasing with α up to α = α s where it then becomes very large. As a consequence, the condition ϕ < α s + γ λ is always required in the blade design phase.

We summarize the properties of C L and C D required for our study in the following assumption. It should be noted that as soon as we set the blade profile, that is, C L and C D , the chord c λ and twist γ λ are the main design parameters. Their optimization will be discussed in Section 4.4.

Glauert's modeling

For the sake of completeness, we now shortly recall here the reasoning proposed by Glauert to model the interaction between a turbine and a stream. We refer to [START_REF] Burton | The Wind Energy Handbook[END_REF]Chapter 3] for a more extended presentation of the theory. We denote by dT the infinitesimal thrust and dQ the infinitesimal torque that apply on the blade element under consideration.

Macroscopic approach

The first part of the model is related to the Momentum Theory, dealing with the macroscopic evolution of the ring of fluid. It aims at expressing dT and dQ in terms of a, a and ϕ. Denoting by p -and p + the fluid pressure left and right neighborhoods of the turbine, we apply twice Bernouilli's relation between -∞ and 0 -and between 0 + and +∞ to get

p --p + = 1 2 ρ(U 2 -∞ -U 2 +∞
). Considering then the rate of change of momentum through the turbine, we obtain a second expression for the variation in the pressure, namely

p --p + = ρ(U -∞ -U +∞ )U 0 .
Combining the two latter equations and using (4.1), we get

U +∞ = (1 -2a)U -∞ .
Since dT = (p --p + )πrdr and dQ = ωρU 0 2πr 3 dr, we finally obtain

dT = 4a(1 -a)U 2 -∞ ρπrdr, (4.6) dQ = 4a (1 -a)λU 2 -∞ ρπr 2 dr. (4.7)
The quantity

C T = dT 1 2 U 2 -∞ ρ2πrdr , ( 4.8)
is often called local thrust coefficient [START_REF] Wilson | Aerodynamic performance of wind turbines[END_REF].

Local expressions

Another set of equations can be obtained through the Blade Element Theory, where local expressions for infinitesimal thrust and torque are considered. The approach can be summarized as follows: by (4.5), we have

dT = B 2 U 2 rel (C L (ϕ -γ λ ) cos ϕ + C D (ϕ -γ λ ) sin ϕ) ρc λ dr, dQ = B 2 U 2 rel (C L (ϕ -γ λ ) sin ϕ -C D (ϕ -γ λ ) cos ϕ) ρc λ rdr,
where B is the number of blades of the turbine under consideration. These equations can then be combined with (4.4) to give

dT = σ λ (1 -a) 2 sin 2 ϕ (C L (ϕ -γ λ ) cos ϕ + C D (ϕ -γ λ ) sin ϕ)U 2 -∞ ρπrdr, ( 4.9
)

dQ = σ λ (1 -a) 2 sin 2 ϕ (C L (ϕ -γ λ ) sin ϕ -C D (ϕ -γ λ ) cos ϕ)U 2 -∞ ρπr 2 dr, ( 4.10)
with σ λ = Bc λ 2πr .

Glauert's relations

To close the system of equations, Glauert combined the previous results. More precisely, equating (4.6) and (4.7) with (4.9) and (4.10) respectively, and dividing both resulting equations by 4(1 -a) 2 gives

a 1 -a = σ λ 4 sin 2 ϕ (C L (ϕ -γ λ ) cos ϕ + C D (ϕ -γ λ ) sin ϕ), (4.11) a 1 -a = σ λ 4λ sin 2 ϕ (C L (ϕ -γ λ ) sin ϕ -C D (ϕ -γ λ ) cos ϕ).
(4.12)

The system obtained by assembling (4.3), (4.11) and (4.12) is the basis of Glauert's Blade element momentum theory.

Simplified model

In the monographs devoted to aerodynamics of wind turbines, the contribution of C D is ocasionally set to zero. This point is discussed in [86, p.135], where it is particular stated that "Since the drag force does not contribute to the induced velocity physically, C D is usually omitted when calculating induced velocities."

In the same way, Manwell mentions in [68, p. This assumption is actually justified in many cases, since the procedures used to design profiles aim at minimizing their drag. Moreover, as explained in Section 4.4.1, the usual blade design procedure starts by selecting a twist angle γ λ that minimizes the ratio C D C L .

We also consider this case, referred as simplified model in the following. It corresponds to the three equations:

tan ϕ = 1 -a λ(1 + a ) , (4.13) a 1 -a = 1 sin 2 ϕ µ L (ϕ) cos ϕ, (4.14) a 1 -a = 1 λ sin ϕ µ L (ϕ), (4.15)
where we have introduced the dimensionless function µ L (ϕ) :

= σ λ 4 C L (ϕ -γ λ ), that is defined on I β,γ λ := [-β + γ λ , β + γ λ ] (4.16)
by virtue of Assumption 4.1.1.

Corrected model

To get closer to the experimental results, many modifications of the model defined by (4.3), (4.11) and (4.12) have been introduced, see e.g. [START_REF] Sørensen | General Momentum Theory for Horizontal Axis Wind Turbines[END_REF]Section 3]. Hereafter, we present three important corrections, namely non-zero drag coefficient C D , tip loss correction and a specific treatment of a for cases where its values become large. The first and last will modify significantly the analysis developed for the simplified Glauert's model.

Weak drag

In the corrected model that we will consider, we assume that C D is strictly positive. However, in our analysis, we will assume small values of this parameter and/or a slow increasing before the occurrence of the stall phenomenon.

Tip loss correction

The independence of the infinitesimal rings of fluid, that is assumed in the model of Glauert is valid for rotors with infinite many blades. In real situations, a modification of the flow at the tip of a blade has to be included to take into account that the circulation of the fluid around the blade must go down (exponentially) to zero. In this way, given a number of blades B and a radius R of the considered turbine, Glauert (see [43, p.268]) has used the Prandtl tip function F λ [START_REF] Prandtl | Vier abhandlungen zur hydrodynamik und aerodynamik[END_REF]:

F λ (ϕ) := 2 π cos -1 exp(- B 2 (1 -r R ) r R sin ϕ ) ,
as a new factor in Equations (4.6) and (4.7), which give rise to

dT = 4a(1 -a)F λ (ϕ)U 2 -∞ ρπrdr, (4.17) dQ = 4a (1 -a)F λ (ϕ)U -∞ ρπr 3 Ωdr. (4.18)
Further models of tip loss corrections have been introduced in between. We refer to [START_REF] Shen | Tip loss corrections for wind research turbine computations[END_REF] and [START_REF] Branlard | Wind turbine aerodynamics and vorticity-based methods: fundamentals and recent applications[END_REF]Chapter 13], for reviews of these models.

Correction for high values of a

For induction factors a larger than about 0.4, a turbulent wake state generally occurs and it is broadly consider that momentum theory does not apply [87, p.297]. This fact was already noted by Glauert [START_REF] Glauert | The analysis of experimental results in the windmill brake and vortex ring states of an airscrew[END_REF], who proposed, for large values of a, to modify the thrust expression (4.8) to fit with experimental observations. Subsequently, many other expressions have been proposed to improve this fitting, see [START_REF] Branlard | Wind turbine aerodynamics and vorticity-based methods: fundamentals and recent applications[END_REF]Section 10.2.2].

All these corrections lead to modify the infinitesimal thrust dT and consequently the term a(1 -a) in Definition (4.17 where (a -a c ) + := max(0, a -a c ) and ψ is a given function defined on R + . Various choices of corrections are presented via the function ψ in Table 4.1. Note that Glauert's empirical correction leads to a discontinuity at a = a c when F λ (ϕ) = 1 (see [15, p.195]). Buhl proposed in [START_REF] Buhl | New empirical relationship between thrust coefficient and induction factor for the turbulent windmill state[END_REF] a slight modification to fix this issue.

Corrected system

We now repeat the reasoning used to obtain Equations (4. 13-4.15), that is, we equalize (4. [START_REF] Buhl | New empirical relationship between thrust coefficient and induction factor for the turbulent windmill state[END_REF]) and (4.18) respectively with (4.9) and (4.10). This gives, using (4.20) and simplifying: where we have introduced the dimensionless functions

tan ϕ = 1 -a λ(1 + a ) , (4.21) a 1 -a = 1 sin 2 ϕ (µ c L (ϕ) cos ϕ + µ c D (ϕ) sin ϕ) - ψ ((a -a c ) + ) (1 -a) 2 , (4.22) a 1 -a = 1 λ sin 2 ϕ (µ c L (ϕ) sin ϕ -µ c D (ϕ) cos ϕ), ( 4
µ c L (ϕ) := σ λ 4F λ (ϕ) C L (ϕ -γ λ ), µ c D (ϕ) := σ λ 4F λ (ϕ) C D (ϕ -γ λ ), (4.24)
defined respectively on I β,γ λ and R. Since in usual applications a ≤ 1, the corrected model coincides with the simplified model when F λ (ϕ) = 1, a c = 1 and C D = 0.

Analysis of Glauert's model and existence of solution

In this section, we reduce each of the two previous versions of Glauert's model to a scalar equation, that explicitly separates local and macroscopic variables (up to one term for the corrected model). This lead us to formulate assumptions that only concern the characteristics of the turbine, that is, on the functions µ L , µ D , µ c L , µ c D , or equivalently, in C L and C D .

We define the angle θ λ ∈ (0, π 2 ) by

tan θ λ := 1 λ . ( 4.25)
and

I := I β,γ λ ∩ - π 2 + θ λ , π 2 + θ λ , ( 4.26)
I + := I ∩ (0, θ λ ].
This result shows that Glauert's model -here in its simplified version -essentially boils down to an only scalar equation: indeed, suppose that ϕ satisfies Equation (4.27), then a and a can be post-computed thanks to Equations (4.14-4.15). These quantities are consequently a by-product of the determination of ϕ.

An important property of Equation (4.27) is that its left-hand side corresponds to the local description of the problem, whereas the right-hand side is rather related to the macroscopic modeling introduced by Glauert's theory, where µ G appears to be a universal function in fluid-turbine dynamics. Another way to formulate this fact consists in associating the left-hand side with 2D experimental data (through C L and C D ) and design variables (as γ λ or σ λ ), and the right-hand side with a 3D model and the parameter θ λ , that does not depend on the blade geometry. This splitting is the heart of Glauert's model. The formulation given in Theorem 4.2.1 gives rise to many criteria to ensure existence of solution of (4.27): the existence indeed holds as soon as the graphs of µ G and µ L intersect. As an illustration, we only give a simple condition in the case of symmetric profiles. Proof. The set I is non-empty as soon as

θ λ -γ λ - π 2 ≤ β.
Since γ λ > 0 and β > 0, the equation above holds, so that max I is well defined. As we consider a symmetric profile, we have µ

L (γ λ ) = C L (0) = 0 whereas µ G (γ λ ) > 0.
Due to the continuity assumption on µ L and (4.28), the existence of solution of (4.27) in [γ λ , max I] follows from the Intermediate Value Theorem. The positivity of µ G on the interval [γ λ , max I] implies that the resulting lift is positive.

If max I = θ λ , then µ G (0) = 0 and µ L is positive on [γ λ , max I], the last assertion follows.

In the case where µ L is supplementary assumed to be increasing on [γ λ , β + γ λ ], then, the solution defined in Theorem 4.2.2 is unique.

Corrected model

We now consider the corrected model defined by Equations (4.21-4.23), for a given value a c ∈ (0, 1). The algebraic manipulations performed in the previous section to get Lemma 4.2.1 cannot be pushed as far as in the simplified model and lead to expressions still containing the unknown a. Hence, before stating a reformulation of this model and an existence result, we need to clarify the dependence of a on the variable ϕ. Again, we express our assumptions in terms of µ c L and µ c D , but the translation in terms of C L , C D , σ λ and F λ (ϕ) is straightforward.

In all this section, we suppose that 0 ∈ I, i.e. |γ λ | ≤ β, which means in particular that

I + = (0, min{θ λ , β + γ λ }], (4.29
∈ I + g(ϕ) := tan -1 ϕ tan(θ λ -ϕ) + µ c D (ϕ) sin ϕ 1 + tan -1 ϕ tan(θ λ -ϕ) . ( 4

.30)

Let ψ be one of the functions given in Table 4.1, with F λ (ϕ) = 1 in the case of Glauert empirical correction. Then, the expression

a 1 -a + 1 - cos θ λ cos ϕ cos(θ λ -ϕ) ψ ((a -a c ) + ) (1 -a) 2 = g(ϕ) (4.31)
defines a continuous decreasing mapping τ : ϕ ∈ I + → a ∈ [0, 1).

Moreover, if g is decreasing and µ c D differentiable, then τ is decreasing and differentiable in all points with a possible exception of one point ϕ c .

As a by-product of the properties of C D stated in Assumption 4.1.1, the function µ c D is always positive and defined for all angles in concrete cases so that g is well defined on I + . Note also that the only obstruction for g to be decreasing would come from the term C D which is often increasing in a neighborhood of 0. But for usual profiles, its variations are negligible when compared to the other (decreasing) terms in (4.30).

Proof. To simplify the notations, let us rewrite Equation (4.31) under the form

u(a) + v(ϕ)w(a) = g(ϕ), (4.32) with u(a) := a 1 -a , v(ϕ) := 1 - cos θ λ cos ϕ cos(θ λ -ϕ) , w(a) := ψ ((a -a c ) + ) (1 -a) 2 .
Let us first consider the left-hand side of (4.32). We see that u is positive and increasing on [0, 1) as well as w for any function ψ given in Table 4.1 (with F λ (ϕ) = 1 in the case of Glauert empirical correction). In the same way, v is positive and increasing on (0, θ λ ]. Fixing now ϕ ∈ (0, θ λ ], it is fairly easy to see that the mapping a ∈ [0, 1) → u(a) + v(ϕ)w(a) is continuous, strictly increasing (actually strictly convex), strictly positive and goes from 0 to +∞. Since g is assumed to be positive on I + , there exists an only a in [0, 1) such that (4.32) holds. Hence the existence of the mapping τ .

Suppose now that g is decreasing and µ c D differentiable. Except w in the point a = a c , all the functions involved in Equation (4.32) are differentiable. Consider ϕ ∈ I + , such that τ (ϕ) = a c . The functions u and w are differentiable in a = τ (ϕ) ∈ (0, 1) and u (a) + v(ϕ)w (a) = 0 so that we can differentiate Equation (4.32) with respect to ϕ. We get

τ (ϕ) = g (ϕ) -v (ϕ)w (τ (ϕ)) u (τ (ϕ)) + v(ϕ)w (τ (ϕ))
.

Combining the fact that g is decreasing with the above properties of v, w, u and their derivatives implies that τ (ϕ) ≤ 0. As a consequence, the mapping τ is either differentiable on the whole interval I + , or differentiable on a set of the form I + -{ϕ c } where ϕ c is the only value in I + such that τ (ϕ c ) = a c The result follows.

Remark 4.2.4. The quantity a = τ (ϕ) can generally be computed explicitly provided that the function ψ is specified analytically, as e.g. in Table 4.1. In the last case, the computation consists in solving a low order polynomial equation in a.

We can now state a result similar to Theorem 4.2.1 in the case of the corrected model. Theorem 4.2.5. Let the assumptions of Lemma 4.2.3 hold so that τ is defined on I + and ψ be one of the functions given in Table 4.1, with F λ (ϕ) = 1 in the case of Glauert empirical correction. Suppose also that (ϕ, a, a (4.21-4.23). Then ϕ satisfies

) ∈ I + -{ π 2 } × R -{1} × R satisfies Equations
µ c L (ϕ) -tan(θ λ -ϕ)µ c D (ϕ) = µ c G (ϕ), (4.33)
where

µ c G (ϕ) := µ G (ϕ) + cos θ λ sin 2 ϕ cos(θ λ -ϕ) ψ ((τ (ϕ) -a c ) + ) (1 -τ (ϕ)) 2 . (4.34)
Reciprocally, suppose that ϕ ∈ I + -{ π 2 } satisfies (4.33) and define a and a by a =τ (ϕ), (4.35) 21-4.23). As it was the case with the simplified model, the value ϕ = π 2 is excluded only for the technical reason that (4.21) is not defined for this angle.

a = 1 -τ (ϕ) λ sin 2 ϕ (µ c L sin ϕ -µ c D cos ϕ). (4.36) Then (ϕ, a, a ) ∈ I + ∈ I + -{ π 2 } × R -{1} × R satisfies Equations (4.
Proof. Let (ϕ, a, a ) ∈ I + -{ π 2 }×[0, 1)×R + satisfiying Equations (4.21-4.23). Because of (4.21), we get:

tan -1 ϕ = λ 1 + a 1 -a = λ(1 + a 1 -a ) + λ a 1 -a .
In the latter, the terms a 1-a and λ a 1-a can be eliminated thanks to (4.22) and (4.23), respectively. After some algebraic manipulations, we end up with a corrected version of (4.27): As in the simplified model, Glauert's model boils down to an only scalar equation, with ϕ as an unknown. However, on the contrary to the simplified model, the formulation (4.33) does not completely decompose the terms into a local part and macroscopic modelling part: much as the left-hand side of Equation (4.33) still relies on local features and experimental 2D data of the problem, its right-hand side now also contains an experimental term, namely µ c D through the definition of τ .

µ c L = (sin ϕ + µ c D ) tan(θ λ -ϕ) + cos θ λ sin 2 ϕ cos(θ λ -ϕ) ψ ((a -a c ) + ) (1 -a)
Before going further, let us give more details about the behavior of τ in ϕ = 0. Lemma 4.2.6. Let Assumption 4.1.1 hold so that τ is defined on I + and ψ be one of the functions given in Table 4.1, with F λ (ϕ) = 1 in the case of Glauert empirical correction. Then

τ (ϕ) = 1 - ψ(1 -a c ) µ c D (0) ϕ 3/2 + o ϕ=0 (ϕ 3/2 ).
Proof. Let us first prove that lim ϕ→0 + τ (ϕ) = 1 -. From Equation (4.30), we see that

lim ϕ→0 + g(ϕ) = +∞. Given ϕ ∈ (0, θ λ ], we have a = τ (ϕ) ∈ [0, 1) and 1 - cos θ λ cos ϕ cos(θ λ -ϕ) = sin θ λ sin ϕ cos(θ λ -ϕ) ≥ 0,
so that all the terms of the left-hand side of (4.31) are positive. As a consequence, the only possibility for the sum of these terms to go to +∞ is that lim ϕ→0 + τ (ϕ) = 1 -.

Proof. Because of Lemma 4.2.6 and Definition (4.34) of µ c G , we get

µ c G (ϕ) ≈ ϕ→0 + µ c D (0) ϕ , so that lim ϕ→0 + µ c G (ϕ) = +∞. ( 4.39)
This implies that there exists a small enough ϕ (4.38), the existence of a solution of Equation (4.33) follows from the Intermediate Value Theorem. The positivity of µ c G on I + implies that the resulting lift is positive.

0 > 0 such that µ c G (ϕ 0 ) ≥ µ c L (ϕ 0 ) - tan(θ λ -ϕ 0 )µ c D (ϕ 0 ). Because of
Suppose now that g is decreasing, max I + = θ λ and ϕ c < θ λ . Because of the last assertion of Lemma 4.2.6, the mapping τ is a decreasing function of ϕ, the correction associated with ψ is not any more active on [ϕ c , θ λ ). We then have

µ c G (max I + ) = µ G (max I + ) = µ G (θ λ ) ≤ 0 ≤ µ c L (max I + ) = µ c L (max I + ) -tan(θ λ -max I + )µ c D (max I + ).
As a consequence, Equation (4.38) holds.

Unlike the simplified model, no condition on γ λ or µ c L (γ λ ) is assumed to get the previous (and following) result, but the alternative (non restrictive) Assumption (4.29) is required. This makes the corrected model much better posed than its simplified version.

Remark 4.2.9. In the case µ c D (0) = 0, a similar reasoning gives:

µ c G (ϕ) ≈ ϕ→0 + (1 + tan θ λ )ϕ.
As a consequence, µ c G (0) = 0, so that as in the simplified model, we need an assumption about, e.g. µ c L (γ λ ) to get an existence result similar to Corollary 4.2.2.

Multiple solutions

The results presented on the previous sections can be completed with some additional remarks about multiple solutions. We identify at least three independent situations.

Simplified model

Since lim ϕ→θ λ ±π/2 µ G (ϕ) = -∞, there shall be two intersections between the graphs of µ G and µ L , e.g. in the case where µ L is affine, C L (0) = 0 and γ λ ∈ (0, θ λ ]. In this framework, one of the two roots gives rise to a negative lift.

for some ε ∈ (0, 1) and the initial value

ϕ 0 = θ λ (4.42)
converges to ϕ , the largest solution in (0, θ λ ] .

Proof. We first prove that f is increasing. Thanks to the choice of ρ and the concavity of ϕ → µ G (ϕ) on [0, θ λ ], we have:

f (ϕ) = 1 + ρ µ G (ϕ) -µ c L (ϕ) -(1 + tan 2 (θ λ -ϕ))µ c D (ϕ) + tan(θ λ -ϕ)µ c D (ϕ) ≥ 1 -ρ max ϕ∈I + µ c L + (1 + tan 2 (θ λ -ϕ)) max ϕ∈I + µ c D + sin θ λ = 1 -ε ≥ 0.
Let us then show that [ϕ , θ λ ] is stable by f . Since f is increasing and f (ϕ ) = ϕ , it remains to show that f (θ λ ) ≤ θ λ . The latter statement follows from:

f (θ λ ) = θ λ + ρ (µ G (θ λ ) -µ c L (θ λ )) = θ λ -ρµ c L (θ λ ) ≤ θ λ .
Since ϕ 0 = θ λ , (ϕ k ) k∈N is bounded and monotonically decreasing, the result follows from the definition of ϕ .

The efficiency of this algorithm depends on the choice of ε and more generally on the value assigned to ρ. In some cases, we can estimate the rate of convergence of (ϕ k ) k∈N . Theorem 4.3.2. In addition to the assumptions of Theorem 4.3.1, suppose that

1 λ = tan θ λ < min ϕ∈[γ λ ,θ λ ] µ c L (ϕ) + µ c D (ϕ). ( 4

.43)

Then the sequence (ϕ k ) k∈N defined by (4.40-4.42) satisfies

|ϕ k -ϕ | ≤ 1 -ρ min ϕ∈[γ λ ,θ λ ] µ c L (ϕ) + µ c D (ϕ) -tan θ λ k |θ λ -ρµ c L (θ λ )|.
Proof. Since we have already proved that f (ϕ) ≥ 0 on [ϕ , θ λ ], it remains to determine an upper bound for f (ϕ). To do this, we use (4.43) and µ G (ϕ) ≤ µ G (0) = tan θ λ to get:

f (ϕ) = 1 + ρ µ G (ϕ) -µ c L (ϕ) -(1 + tan 2 (θ λ -ϕ))µ c D (ϕ) + tan(θ λ -ϕ)µ c D (ϕ) ≤ 1 -ρ min ϕ∈[γ λ ,θ λ ] µ c L (ϕ) + µ c D (ϕ) -tan θ λ .
The result follows by induction.

We can actually obtain a quadratic convergence, i.e. |ϕ k -ϕ | ≤ ν|ϕ 0 -ϕ | 2 k for some ν > 0 by making use of a Newton procedure, that is by replacing ρ given in (4.41) by the sequence (ρ k ) k∈N defined by

ρ k := - 1 µ G (ϕ k ) -µ c L (ϕ k ) -(1 + tan 2 (θ λ -ϕ k ))µ c D (ϕ k ) + tan(θ λ -ϕ k )µ c D (ϕ k )
.

Though the term µ G (ϕ k ) can be computed exactly as well as most of the terms in the denominator, the latter expression remains difficult to evaluate since the functions µ c L and µ c D are only known experimentally, i.e. pointwise.

The general case

If the correction for high values of a is considered and applies, then the framework of Theorem 4.2.8 implies that there exists a solution of the corrected model in I + . As a consequence, a bisection algorithm applied on Equation (4.33) converges to such a solution.

We finally show that the solution found in the case ψ = 0 can be used to bracket the solution in the general case of the corrected model. Proof. Since ϕ satisfies (4.33) with ψ = 0, we have:

µ c L (ϕ ) -tan(θ λ -ϕ )µ c D (ϕ ) -µ c G (ϕ ) = - cos θ λ sin 2 ϕ cos(θ λ -ϕ ) ψ ((τ (ϕ ) -a c ) + ) (1 -τ (ϕ )) 2 ≤ 0.
The result follows by combining the latter with (4.38) and by applying Corollary 4.2.8.

Optimization

The BEM model does not only aim at evaluating the efficiency of a given geometry, but also provides a method to design rotors, that is, to select high-performance parameters γ λ and c λ . In this way, monographs generally consider a particular functional, often called power coefficient ([68, p.126], [44, p.328] and [START_REF] Sørensen | General Momentum Theory for Horizontal Axis Wind Turbines[END_REF]), which corresponds to the ratio between the energy received and the energy captured.

In this framework, the drag coefficient C D is taken into account (though partly neglected in the reasoning, as explained hereafter) as well as the Tip loss correction. To the best of our knowledge, no correction related to high values of a, as the one presented in Section 4.1.4, is considered in optimization procedures.

Simplified model and usual design procedure

The design procedure mainly consists in optimizing the power coefficient, defined by the relation

C p (γ λ , c λ , ϕ) = 8 λ 2 max λmax λ min λ 3 a (1 -a) 1 - C D C L tan -1 ϕ dλ, ( 4.44)
under the constraints (4.21-4.23), in the sense that only a Tip loss correction is occasionally considered, which actually does not modify the optimization procedure described in [68, p.131-137]. For the sake of completeness, we recall it here keeping the mathematical point of view adopted so far.

The first step consists in defining an angle α that minimizes the ratio C D C L . Then the coefficient C D is simply neglected: not only the factor 1 -C D C L tan -1 (ϕ) is set to 1 in (4.44), but it also vanishes from the constraints, which now become the simplified model (4. 13-4.15). Finally, Lemma 4.2.1 allows us to replace µ L by µ G , so that the simplified model (and hence C p ) depends exclusively on ϕ. Then, the optimization problem is reduced to

Asymptotical analysis of the corrected model

The optimization problem associated with the corrected model, in the case of Glauert empirical correction (F λ (ϕ) = 1), is given by max It is clear that its solving becomes more complicated than in the simplified case.

However, what if we can identify regimes where the correction is not active at the optimum? We provide a asymptotic result in a neighborhood of ϕ = 0 where we discuss this situation.

Theorem 4.4.1. Suppose that Lemma 4.2.3 holds and the functions µ D and ψ ((a -a c ) + ) are differentiable. Then, in a neighborhood [0, δ) of ϕ = 0, the cost functional satisfies

J(ϕ) = λ µ D (0) ϕ 2 + o ϕ=0 (ϕ 2 ).
Moreover, if g is decreasing and ϕ c ∈ [0, δ), then the optimal solution belongs to the region where the correction ψ((a -a c ) + ) is not active.

Before stating the proof, an improved version of Lemma 4.2.6 is required. Proof. Consider the following expansions that hold in a neighborhood of ϕ = 0:

sin θ λ sin ϕ cos(θ λ -ϕ) ψ ((a -a c ) + ) = ψ(1 -a c ) λ ϕ - ψ(1 -a c ) λ 2 ϕ 2 + o ϕ=0 (ϕ 2 ), g(ϕ) = µ D (0) λ 1 ϕ 2 + 1 λ + µ D (0) λ - µ D (0) λ 2 1 ϕ + o ϕ=0 ( 1 ϕ ),
where 1 λ = tan θ λ . Note that the first can be explained as follows: in a neighborhood of a = 1 (or equivalently around ϕ = 0), due to Lemma 4. Defining ν(ϕ) = 1 -τ (ϕ), we multiply Equation (4.31) by ν 2 (ϕ) and then use the previous expansions. After rearranging terms, we have

µ D (0) λ ν 2 (ϕ) ϕ 2 = ν(ϕ) -ν 2 (ϕ) + ψ(1 -a c ) λ ϕ - ψ(1 -a c ) λ 2 ϕ 2 + o ϕ=0 (ϕ 2 ) - 1 λ + µ D (0) λ - µ D (0) λ 2 + o ϕ=0 (1) ν 2 (ϕ) ϕ .
We define a function ξ(ϕ) such that ν 2 (ϕ) = ϕ 3 (1 + ξ(ϕ)), with = ψ(1-ac) µ D (0) . Plugging this expression on the left-hand side, and using on the right-hand side the expansions ν(ϕ) = √ ϕ 3/2 + o ϕ=0 (ϕ 3/2 ) and ν 2 (ϕ) = ϕ 3 + o ϕ=0 (ϕ 3 ), lead to

ψ(1 -a c ) λ ϕξ(ϕ) = √ ϕ 3/2 + o ϕ=0 (ϕ 3/2 ) - ψ(1 -a c ) λ 2 ϕ 2 - 1 λ + µ D (0) λ - µ D (0) λ 2 + o ϕ=0 (1) ϕ 2 .
Since ϕ 2 = o ϕ=0 (ϕ 3/2 ), we can merge all the quadratic terms with o ϕ=0 (ϕ 3/2 ). Dividing by ϕ, we finally get

ξ(ϕ) = λ √ ψ(1 -a c ) ϕ 1/2 + o ϕ=0 (ϕ 1/2 ).
Proof of Theorem 4.4.1. We begin by deriving an explicit expression for the cost functional. After replacing the constraints associated with a and a into J(ϕ), we obtain In order to determine a regime where the optimum of J(ϕ) is attained in [0, a c ), i.e. when the correction does not apply, we study the asympotical behavior of J(ϕ) around ϕ = 0 (or equivalently, a = 1). With the help of Lemma 4.2.6, we can expand the first terms of the right-hand side of (4.48), to get From the previous calculations we also obtain Finally, we denote by [0, δ) the neighborhood where all these computations hold. Since the cost functional is increasing in this interval, necessarily the optimum ϕ * is attained outside. Then, if exists ϕ c ∈ [0, δ) such that τ (ϕ c ) = a c , we have τ (ϕ * ) ≤ τ (δ) ≤ a c , due to the fact that τ is a decreasing function. It follows that a correction for high values of a is not required.

Numerical experiments

In this section, we investigate on a practical case the performance of the algorithms presented in Section 4.3 and study numerically the design optimization problem considered in Section 4.4.

A practical example

We consider a turbine consisting of three blades, designed with a NACA 4415 profile and three different blade elements of rotation speeds λ 1 = 0.5, λ 2 = 1.25, and λ 3 = 3, respectively. In all these cases, we set γ λ = γ * λ and c λ = c * λ , i.e. we use the optimal values of the simplified model given by (4.45) (where ᾱ = 0.105 radians):

Solution algorithms

We first compare the various solution algorithms, in terms of iterations. After having computed accurately a numerical solution φ ∞ , we measure the number of iterations required to reach an error satisfying

err k := |φ k -φ ∞ | ≤ Tol = 10 -5 .
We consider the three algorithms presented in Section 4.3, that is, the usual aglorithm Algorithm 4.1, the bisection algorithm detailed in Section 4.3.2 and the new fixedpoint algorithm corresponding to the iteration (4.40). For the latter, we use for ρ the value given by Equation (4.41) (though we include in our test a correction on α), with ε = 1. We see in the table above that though we do not use an optimal value for ρ, the convergence of our new fixed-point algorithm converges faster than the others. To evaluate the effect of a change in ρ, we plot the values of err k with respect to the iterations in the case where λ = λ 2 . The results are presented in The value ρ = 1.4484 is the one given by Equation (4.41) for ε = 1. The test confirms that this value is not optimal, and can be easily improved. However, using greater or smaller values can deteriorate the convergence, as is the case for ρ = 1 or ρ = 3.

Optimization

In this test, we work with three values λ 1 = 1.05, λ 2 = 2.96 and λ 3 = 4.88, corresponding to three elements of the actual turbine used by Hydrotube Energy.

We use a gradient method to compute a solution of (4.46) and compare it to the (explicit) optimal solution associated with the simplified model given by Equations (4.45). For the sake of completeness, we detail now the optimization method we use. In particular, we recall here how the introduction of an adjoint variable enables to compute the gradient of J. In this last section, we denote by C L and C D the derivatives of C L and C D respectively, and omit the dependence of µ L and µ D and their derivatives on the variable ϕ for the sake of simplicity. where The optimization algorithm is then:

We define the Lagrangian of

∂L ∂γ λ = a (1 -a) C D (ϕ -γ λ )C L (ϕ -γ λ ) -C L (ϕ -γ λ )C D (ϕ -γ λ ) C 2 L (ϕ -γ λ )

end

We apply this algorithm using as initial guess the solution (4.45) associated with the simplified model. The results are presented in Table 4.3. We see that the gradient procedure enables to improve significantly the values of J, namely by 7.4918%, 14.267% and 19.907% for the three elements under consideration. Notice that in the three cases, the correction is activated, i.e. a > a c , though the value of a is slightly lower than with the parameters values given by (4.45).

J

Perspectives

This work focuses mainly on providing conditions that ensure the existence of solutions for the BEM model. Several questions regarding the convergence of solving algorithms, as well as the associated optimization problem, remain to be answered.

Concerning the former, extending the proof to a more general framework is desirable, in particular for the simplified case, since its analysis rely on a fixed point argument. The optimization problem also offers several possibilities. An asymptotic analysis gives a general idea of the optimal solution behavior, however the required assumptions seem very restrictive. It could be useful to prove a similar result by using optimization tools, to confirm when the correction for high values of a is needed. If it is not the case, a next step concerns the negligibility of C D , focusing on the difference between the optimal solution and the obtained in the simplified case. Finally, the question of multiple optima in the corrected model remains open. A systematic use of numerical experiments on various cases could provide first answers.

4.A. Convergence in the simplified case

For ϕ ∈ (γ λ , θ λ ), µ L (ϕ)h(ϕ) ≥ 0 and µ L (ϕ)h (ϕ) ≤ 0, so that -1 1 + (λ + µ L (ϕ)h(ϕ)) 2 µ L (ϕ)h(ϕ) ≤ f (ϕ) ≤ -1 1 + (λ + µ L (ϕ)h(ϕ)) 2 µ L (ϕ)h (ϕ).

The result then follows from the fact that µ L (ϕ)h(ϕ) ≥ 0, h and h is decreasing on (γ λ , θ λ ).

We are now in the position to obtain a conditional convergence result, the condition being stated in the next assumption.

L

  'algorithme BFN combine ces procédures en définissant la méthode itérative ẋk(t) = Ax k (t) + K [x obs (t) -xk (t)] dans [0, T ] xk (0) = xk-1 (0) ẋk (t) = Ax k (t) -K [x obs (t) -xk (t)] dans [0, T ] xk (T ) = xk (T )avec x0 (0) = x 0 . Sous certaines hypothèses, les deux suites {x k (t)} k≥1 et {x k (t)} k≥1 convergent vers x∞ (t) et x∞ (t), respectivement. De plus, ces fonctions limites présentent également le comportement asymptotique décrit dans (2), même pour t = 0.

Figure 1 :

 1 Figure 1: Méthode de Schwarz (discrétisée) appliqué àd 2 u dx 2 + ηu = f in [0, 1] [37, p.8].

Figure 2 :

 2 Figure 2: L'algorithme Pararéel appliqué à l'equation de Dahlquist u(t) = -i 2 u in [0, 20]

  Décomposition de la pale en éléments annulaires

Figure 3 :

 3 Figure 3: Décompositions impliquées dans la théorie BEM [54, p.8].

Figure 1 . 1 :

 11 Figure 1.1: (Discretized) Schwarz methods applied tod 2 u dx 2 + ηu = f in [0, 1] [37, p.8].

Figure 1 . 2 :

 12 Figure 1.2: Parareal algorithm applied to the Dahlquist equation u(t) = -i 2 u in [0, 20]

Figure 1 . 3 :

 13 Figure 1.3: Decompositions involved in BEM Theory [54, p.8].

  4) where φ(s) = (s -µ 1 ) • • • (s -µ m ), i.e. (µ i ) i=1,...,m are the eigenvalues of A -LC.

Figure 2 . 1 : 5 .

 215 Figure 2.1: Comparison between k th and k obs , for N = 16 and δt = ∆T 2 5 . The eigenvalues of A -LC are {-0.8, -1} (left) and {-0.2, -0.25} (right).

Figure 2 .

 2 Figure 2.2 also shows that the observed efficiencies satisfy

(

  b) E(N ), for δt = ∆T 2 5 and γ = 210 . E(δt), for N = 16 and γ = 2 10 .

Figure 2 . 2 :

 22 Figure 2.2: Comparison between E obs (k obs ), E obs (k th ) and E th 0 . The eigenvalues of A -LC are {-0.8, -1} (left) and {-0.2, -0.25} (right).

. 23 )

 23 and the well-posedness of Problem(3.22) follows from the Fredholm Alternative. Uniqueness holds for any q ∈ L ∞ (Ω) satisfying(3.21) owning to [46, Theorems 2.1, 2.4].

. 26 )

 26 where C(Ω) > 0 only depends on the domain.Proof. We recall that the existence and uniqueness of solutions to Problem (3.22) was already stated in Subsection 3.2.1.

  (a) Norm of the numerical solution. (b) Real part of the incident wave. (c) Real part of the numerical solution.

Figure 3 . 2 :

 32 Figure 3.2: Numerical solution of a wave damping problem. The yellow part represents Ω 0 and the red part corresponds to the nodal points associated with q.

Figure 3 . 3 :

 33 Figure 3.3: Detection of a bathymetry from a wavefield. The yellow part represents Ω 0 and the red part corresponds to the nodal points associated with q.

. 5 )

 5 where -π/2 ≤ γ λ ≤ π/2 is the twist (also called local pitch) angle of the blade. The notations associated with a blade element are summarized in Figure4.1.

Figure 4 . 1 :

 41 Figure 4.1: Blade element profile and associated angles, velocities and forces

Assumption 4 . 1 . 1 .

 411 In what follows, we assume that C L is well-defined and continuous on an intervalI β := [-β, β]for some β ∈ [0, α s ) and positive on I β ∩ R + . The coefficient C D is well-defined and positive on R.

  125]: "In the calculation of induction factors,[...] accepted practice is to set C D equal to zero [...]. For airfoils with low drag coefficients, this simplification introduces negligible errors.".

  ), that becomesdT = 4χ(a, a c )F λ (ϕ)U 2 -∞ ρπrdr. (4.19)In the literature, the function χ(a, a c ) is in most cases of the form χ(a, a c ) = a(1 -a) + ψ ((a -a c ) + ) ,(4.20)

Corollary 4 . 2 . 2 .

 422 In addition to Assumption 4.1.1, suppose that the profile under consideration is symmetric with γ λ > 0, and that µ G (max I) ≤ µ L (max I), (4.28) where max I = min{θ λ , β + γ λ }. Then (4.27) admits a solution in [γ λ , max I] corresponding to a positive lift. Moreover, if max I = θ λ , i.e. θ λ ≤ β + γ λ , then (4.28) is automatically satisfied.

Lemma 4 . 3 . 3 .

 433 Keep the assumptions of Corollary 4.2.8, and denote by ϕ a solution of (4.33), when ψ = 0. Then (4.33) admits a solution in (ϕ , min{θ λ , β + γ λ }] corresponding to a positive lift.

2 ϕ

 2 Recalling that µ G (ϕ) = sin(ϕ) tan(θ λ -ϕ), we can express a and a as functions of ϕ, namely1 -a = sin ϕ cos(θ λ -ϕ) sin θ λ , a = sin ϕ sin(θ λ -ϕ) cos θ λ , so that C p reads sin (2(θ λ -ϕ)) dλ.It then remains to optimize ϕ → sin 2 ϕ sin(2(θ λ -ϕ)), for ϕ ∈ [0, θ λ ]. An easy computation shows that the maximum is attained at ϕ * = 2 3 θ λ . Finally, the design parameters γ * λ = γ λ (ϕ * ) and c * λ = c λ (ϕ * ) can then be computed from Equations (4.5) and (4.27), to get γ * λ = ϕ * -α, c * λ = 8πrµ G (ϕ * ) BC L (α) . (4.45)

ϕJµ

  (ϕ) = a (1 -a) 1 -µ D (ϕ) µ L L (ϕ) = µ c G (ϕ) + tan(θ λ -ϕ)µ D (ϕ) a = τ (ϕ) a = 1 -τ (ϕ) λ sin 2 ϕ (µ L (ϕ) sin ϕ -µ D (ϕ) cos ϕ) .

Lemma 4 . 4 . 2 . 3 / 2 D

 44232 Under the assumptions of Lemma 4.2.6, we have(1 -τ (ϕ)) 2 = ψ(1 -a c ) µ D (0) ϕ 3 + λ ψ(1 -a c ) µ (0) ϕ 3+1/2 + o ϕ=0 (ϕ 3+1/2 ).

  2.6 we haveψ ((a -a c ) + ) = ψ(1 -a c ) + dψ((a -a c ) + ) da a=1 (a -1) + o a=1 (a -1) = ψ(1 -a c ) +   dψ((a -a c ) + ) da a=1 ψ(1 -a c ) µ D (0)   ϕ 3+1/2 + o ϕ=0 (ϕ 3+1/2 ).

J(ϕ) = [( 1 - 2 λ( 1 -

 121 τ (ϕ)) 2 (µ L (ϕ) sin ϕ -µ D (ϕ) cos ϕ)] τ (ϕ)) 2 µ L (ϕ) sin 3 ϕ . (4.47)Note that Equation (4.34) allows us to write explicitely the numerator. Indeed, we have(1 -τ (ϕ)) 2 (µ L (ϕ) sin ϕ -µ D (ϕ) cos ϕ) = (1 -τ (ϕ)) 2 sin 2 ϕ tan(θ λ -ϕ) + cos θ λ cos(θ λ -ϕ) sin 3 ϕ ψ((τ (ϕ) -a c ) + ) -(1 -τ (ϕ)) 2 µ D (ϕ) .

( 1 - 3 / 2 D

 132 τ (ϕ)) 2 sin 2 ϕ tan(θ λ -ϕ) = ψ(1 -a c ) λµ D (0) ϕ 5 + o ϕ=0 (ϕ 5 ), cos θ λ cos(θ λ -ϕ) = 1 + o ϕ=0 (1), sin 3 ϕ ψ((a -a c ) + ) = ψ(1 -a c ) • ϕ 3 + o ϕ=0 (ϕ 4 ).However, applying Lemma 4.2.6 to the last term in (4.48) leads to a right-hand side of order o ϕ=0 (ϕ 3 ), and then J(ϕ) = o ϕ=0 (ϕ). Instead, we use Lemma 4.4.2, which brings(1 -τ (ϕ)) 2 µ D (ϕ) = ψ(1 -a c ) µ D (0) ϕ 3 + λ ψ(1 -a c ) µ (0) ϕ 3+1/2 + o ϕ=0 (ϕ 3+1/2 ) × (µ D (0) + µ D (0) • ϕ + o ϕ=0 (ϕ)) = ψ(1 -a c ) • ϕ 3 + λ ψ(1 -a c ) µ D (0) ϕ 3+1/2 + o ϕ=0 (ϕ 3+1/2 )and then, replacing all these expressions into (4.48) and using that ϕ 4 = o(ϕ 3+1/2 ), yields(1 -τ (ϕ)) 2 (µ L (ϕ) sin ϕ -µ D (ϕ) cos ϕ) = λ ψ(1 -a c )µ D (0) ϕ 3+1/2 + o ϕ=0 (ϕ 3+1/2 ). (4.49)

( 1 -ϕ 2 λ

 12 τ (ϕ)) 2 µ L (ϕ) sin ϕ = ψ(1 -a c ) • ϕ 3 + o ϕ=0 (ϕ 3 ). (4.50)Now we can expand the cost functional by replacing (4.49) and (4.50) in (43+1/2 + o ϕ=0 (ϕ 3+1/2 ) (ψ(1 -a c ) • ϕ 3 + o ϕ=0 (ϕ 3 )) (ϕ + o ϕ=0 (ϕ 2 )) 2 = λ µ D (0) ϕ 2 + o ϕ=0 (ϕ 2 ).

  use the correction from Wilson et al and Spera, see Table 4.1. In this example, I + = (0, θ λ ] for all elements. Graphs of the functions µ c LD : ϕ → µ c L (ϕ) -tan(θ λϕ)µ c D (ϕ), µ c G and µ G are presented below.

(a) λ 1

 1 = 0.5 (b) λ 2 = 1.25 (c) λ 3 = 3

Figure 4 . 2 :

 42 Figure 4.2: Graphs of of the functions µ c LD , µ c G and µ G for different values of λ. Note that these figures are similar to the scheme given in [68, Figure 3.27, p.126].

λ 1 =Table 4 . 2 :

 142 Algorithm 4.1 Bisection Fixed-point Algorithm λ Number of iterations required to solve Equations (4.21-4.23)

  Figure 4.3.

Figure 4 . 3 :

 43 Figure 4.3: Convergence of the new fixed-point algorithm for various values of ρ

2 -p 3 • a 1 -a - 1 λa ( 1 .

 2111 Problem 4.46 by L(ϕ, a, a , p 1 , p 2 , p 3 , c λ , γ λ ) := J(ϕ, a, a , c λ , γ λ ) -p 1 • (µ L (ϕ) -tan(θ λ -ϕ)µ D (ϕ) -µ c G (ϕ)) cos ϕ + µ c D (ϕ) sin ϕ) + ψ ((a -a c ) + ) (1 -a) sin 2 ϕ (µ c L (ϕ) sin ϕ -µ c D (ϕ) cos ϕ ,where p 1 , p 2 and p 3 are the Lagrange multipliers associated with the constraints (4.21-4.23).The optimality system is obtained by cancelling all the partial derivatives of L. Differentiating L with respect to p 1 , p 2 and p 3 and equating the resulting terms to zero gives the corrected model Equations (4.21-4.23), that can be solved using the algorithms presented in Section 4.3. Setting the derivatives of L with respect to (ϕ, a, a ) to zero gives rise to the linear systemM • p = b (4.51)where p := (p 1 p 2 p 3 ) is the Lagrange multiplier vector, and -a)C L (ϕ-γ λ )C D (ϕ-γ λ )-C D (ϕ-γ λ )C L (ϕ-γ λ ) C L (ϕ-γ λ ) 2 tan ϕ + C D (ϕ-γ λ ) C L (ϕ-γ λ ) sin 2 ϕ -a 1-C D (ϕ-γ λ ) C L (ϕ-γ λ ) tan ϕ (1 -a) 1-C D (ϕ-γ λ ) C L (ϕ-γ) tan ϕWe are in a position to detail how the gradient can be computed. Fix the values of the pair (γ λ , c λ ). If ϕ, a, a are the corresponding solutions of Equations (4.21-4.23) and p is the associated solution of Equation (4.51), then the gradient ∇J(γ λ , c λ ) of J(γ λ , c λ ) is given by ∇J(γ λ , c λ ) = ∂L ∂γ λ

Algorithm 4 . 2 : 2 . 3 . 5 .

 42235 Numerical optimizationInput: Tol > 0, κ > 0, α → C L (α), α → C D (α), λ, x → ψ(x).Initial guess: γ λ , c λ . Output: γ λ , c λ . Set err := Tol + 1. Define the functions µ c L and µ c D by (4.24) using the input data. while err > Tol do 1. Set ϕ, a, a as the solutions of Equations (4.21-4.23). Set p as the solution of Equation (4.51). Compute the gradient ∇J(γ λ , c λ ) given by Equation (4.52). κ∇J(γ λ , c λ ), Set err := ∇J(γ λ , c λ ) .

Theorem 4 .A. 3 .

 43 In addition to the assumptions of Lemma 4.A.

  Méthodes variationnelles discrètes. Contrairement au cadre continu, qui nécessite de résoudre l'équation linéaire tangente (4) ou l'équation adjointe pour obtenir le gradient, la discrétisation (3) conduit à son calcul direct. Ensuite, les méthodes variationnelles discrètes diffèrent dans la façon dont elles décrivent l'équation gouvernante et le traitement des non-linéarités possibles.Étant donné une discrétisation {t n } N n=0 de [0, T ], nous désignons par x i et y i l'état et le vecteur d'observation au temps t i , respectivement. En utilisant la même notation qu'auparavant sur les variables J, u, x b 0 , B et R, même si elles peuvent dépendre de la discrétisation temporelle considérée, l'algorithme 4D-Var se lit ainsi min u∈U

Équations de Saint-Venant. La pression

	de l'eau η(x, t) et au fond -z b (x), éventuellement en fonction des effets de la pression, quand la profondeur -z b est supposée constante. Dans des conditions aux limites
	de la viscosité et du frottement, voir la section 3.1.1 pour plus de détails. Dans ce qui appropriées, cette équation est utilisée pour étudier les problèmes de diffusion des
	suit, par souci de simplicité, nous négligeons les deux derniers. vagues, comme celle présentée au Chapitre 3, où une bathymétrie variable agit comme
	un diffuseur.
		Certaines simplifications des équations de Navier-Stokes donnent lieu à différents
	modèles de propagation des vagues [16, 59], tous dépendant de la relation entre trois Équations de Boussinesq classique. D'autre part, une pression non hydrosta-
	paramètres caractéristiques appelés profondeur relative H, longueur horizontale L et tique introduit des termes dispersifs dans les équations gouvernantes. Un troisième
	amplitude verticale maximale A. Pour rappeler certains d'entre eux, nous introduisons rapport, appelé le nombre d'Ursell
	les rapports
					ε = U r = H L , δ = δ ε 2 ,	A H	.
	Les rivières, les domaines côtiers et les océans peuvent être modélisés via une intégra-tion verticale en utilisant l'hypothèse d'eau peu profonde ε mesure leur force par rapport aux termes non linéaires [89]. En mouvement bidimen-1, alors que les vagues de petite amplitude sont décrites par des modèles où δ 1. sionnel, un choix spécifique de ce paramètre (U r = O(1) et δ 1) conduit à
		En supposant un régime d'eau peu profonde, ainsi que différentes valeurs de δ et
	des approximations de la pression, on obtient les modèles suivants.
					hydrostatique, c-à-d que la contribution
	de l'accélération verticale du fluide dans la pression est négligeable, combinée avec
	δ = O(1) donne
		                    	∂hv ∂t ∂hu ∂t	+ +	∂η ∂t ∂hu v ∂x + ∂hu 2 ∂x + ∂hu v + ∂hu ∂x ∂x + gh + ∂hv ∂y ∂η ∂y ∂hv 2 ∂y + gh ∂η ∂x	= 0 = 0 = 0	(10)
	où h = η + z b est la hauteur d'eau et u, v sont les vitesses moyennes sur la profondeur.
	Ce modèle est couramment utilisé pour modéliser les ruptures de barrage [40], les
	sauts hydrauliques et les flux de marée dans les estuaires.
	Équation des ondes. Maintenir la pression hydrostatique mais en supposant que
	δ	1 est équivalent à négliger les termes d'accélération convective dans (10) et
	linéariser le nouveau système autour du niveau de la mer η = 0. En combinant les
	expressions résultantes, on obtient
				∂ 2 η ∂t 2 -∇ (gz b ∇η) = 0.
		En particulier, sa solution a pour forme η(x, t) = Re{ψ tot (x)e -iωt }, où ω est une
	w) représente la vitesse du fluide dans la région bornée et dépendante fréquence donnée et l'amplitude ψ tot satisfait l'équation de Helmholtz
	du temps Ω t ; ρ, g et σ T représentent respectivement le coefficient de densité, la grav-ité et le tenseur de contrainte. Nous omettons ici les conditions aux limites au niveau ψ tot + ω 2 gz b ψ tot = 0,

  Combining Theorem 2.3.1 with the definition of k in (2.15), we obtain

.

[START_REF] Bristeau | Derivation of a non-hydrostatic shallow water model; comparison with Saint-Venant and Boussinesq systems[END_REF]

) Suppose that we apply Algorithm 2.1 using k iterations of the Parareal algorithm, i.e. h = k . Then, the stopping criterion (2.8) is satisfied.

Proof.

Table 4 .

 4 1: Various corrections proposed in the literature. The order corresponds to the degree of a (as a polynomial) in χ(a, a c )

	.23)

  2 , so that (4.33) is satisfied. Using now the last expression to eliminate µ c L in (4.22) gives (4.31). Consequently, Lemma 4.2.3 implies that a and ϕ satisfy (4.35). Finally, (4.36) is a direct consequence of (4.35) and (4.23). Suppose now that (ϕ, a, a ) ∈ I + -{ π 2 }×[0, 1)×R + satisfies Equations (4.33-4.36). Replacing τ (ϕ) by a in (4.36) gives immediately (4.23). Combining (4.35) with the definition of µ c G give (4.22). Finally, (4.21) is obtained by introducing a and a thus defined into (4.33).

Table 4 .

 4 3: Optimal values obtained with Algorithm 4.2

	c λ	γ λ

  1, suppose that µ L is differentiable and satisfies Then, if ϕ 0 belongs to [γ λ , θ λ ], the sequence (ϕ k ) k∈N defined by Equation (4.53) converges to the unique solution of Equation (4.27). Proof. As a consequence of Lemma 4.A.1, the function f maps [γ λ , θ λ ] onto itself. From Equations (4.55), (4.56) and Lemma 4.A.2 we deduce f is contracting. The result follows from the Banach Fixed Point Theorem.

	max I µ L h(γ λ ) 1 + λ 2	≤ 1	(4.55)
	max I µ L |h (γ λ )| 1 + λ 2	≤ 1.	(4.56)

La méthode de l'élément de pale (Blade element momentum theory, BEM) est un modèle mécanique largement utilisé pour évaluer la performance des turbines, en fonction des caractéristiques mécaniques et géométriques de leurs pales et du courant auquel elles sont exposées. Développé par Glauert[START_REF] Glauert | The Elements of Aerofoil and Airscrew Theory[END_REF], il résulte de la combinaison de deux modèles différents : Blade element theory (BET) et la Théorie de Froude (Momentum theory, MT).En découpant la pale en sections qui sont traitées selon un modèle planaire, BET étudie le comportement de la turbine d'un point de vue local[START_REF] Froude | On the elementary relation between pitch, slip and propulsive efficiency[END_REF]. Les quantités fondamentales de ce modèle sont les coefficients C L et C D appelés portance et traînée, qui sont introduits pour tenir compte ces forces exprimées dans le plan de coupe. Les résultats sont ensuite intégrés le long de la pale pour obtenir les quantités d'intérêt global. En contraste, la Théorie de Froude[START_REF] Rankine | On the mechanical principles of the action of propellers[END_REF] est une théorie globale qui étudie de façon macroscopique le comportement d'une colonne de fluide traversant une turbine.

Acknowledgments

To show that lim n→+∞ ∇ψ(q n ) 2 L 2 (Ω) = ∇ψ(q ∞ ) 2 L 2 (Ω) , note that

Using the same arguments as those to prove (3.29), we have a subsequence (same notation used) such that q n → q ∞ pointwise a.e. in Ω and thus √ 1 + q n -1 → √ 1 + q ∞ -1 pointwise a.e. in Ω. Due to Lebesgue's Dominated Convergence Theorem and √ 1 + q n ∇ψ(q n ) → √ 1 + q ∞ ∇ψ(q ∞ ) strongly in L 2 (Ω), we have

The latter, together with Theorem 3.2.6, shows that ψ(q n ) → ψ(q ∞ ) strongly in H 1 .

Existence of optimal solution

We are now in position to prove the existence of a minimizer to Problem (3.24) for U = U Λ . Theorem 3.2.7. Assume that the cost function (q, ψ) ∈ U Λ → J(q, ψ) ∈ R satisfies: (A1) There exists β > 0 such that J(q, ψ) = J 0 (q, ψ) + β|Dq|(Ω).

(A2) ∀(q, ψ) ∈ U Λ × H 1 (Ω), J 0 (q, ψ) ≥ m > -∞.

(A3) (q, ψ) → J 0 (q, ψ) is lower-semi-continuous with respect to the (weak * ,weak) topology of BV (Ω) × H 1 (Ω).

Then the optimization problem (3.24) has at least one optimal solution in U Λ ×H 1 (Ω).

Proof. The existence of a minimizer to Problem (3.24) can be obtained with standard techniques by combining Theorem 3.2.6 with weak-compactness arguments as done in [21, Lemma 2.1], [START_REF] Bastide | Penalization model for Navier-Stokes-Darcy equation with application to porosity-oriented topology optimization[END_REF]Theorem 4.1] or [START_REF] Haslinger | On a topology optimization problem governed by two-dimensional Helmholtz equation[END_REF]Theorem 1]. We still give the proof for the sake of completeness.

We introduce the following set A = (q, ψ) ∈ U Λ × H 1 (Ω) a(q; ψ, φ) = b(q; φ) ∀φ ∈ H 1 (Ω) .

The existence and uniqueness of solution to Problem (3.22) ensure that A is nonempty. In addition, combining assumptions (A1) and (A2) we have that J(q, ψ) is bounded from below on A. We thus have a minimizing sequence (q n , ψ n ) ∈ A such that lim n→+∞ J(q n , ψ n ) = inf (q,ψ)∈A J(q, ψ).

Theorem 3.2.3 and (A1) then gives that the sequence (q n , ψ n ) ∈ BV (Ω) × H 1 (Ω) is uniformly bounded with respect to n and thus admits a subsequence that converges towards (q * , ψ * ) in the (weak * ,weak) topology of BV (Ω)×H 1 (Ω). Using now Theorem 3.2.6 and the weak * lower semi-continuity of q → |Dq|(Ω), we end up with (q * , ψ * ) ∈ A and

It is worth noting that β has been introduced only to obtain a uniform bound in the BV -norm for the minimizing sequence.

When U = U Λ,κ , we note that any q ∈ U Λ,κ is actually bounded in BV since

With this property at hand, we can get a similar result to Theorem 3.2.7 without adding a penalization term in the cost function, hence β = 0. Theorem 3.2.8. Assume that the cost function (q, ψ) ∈ U Λ,κ → J(q, u) ∈ R satisfies (A2) -(A3) given in Theorem 3.2.7 and that β = 0. Then the optimization problem (3.24) with U = U Λ,κ has at least one optimal solution.

Proof. We introduce the following non-empty set

From (A2), J(q, ψ) is bounded from below on A. We thus have a minimizing sequence

) and thus admits a convergent subsequence toward some q ∈ U Λ,κ . Theorem 3.2.6 then gives that ψ(q n ) → ψ(q) strongly in H 1 (Ω) and the proof can be finished as the proof of Theorem 3.2.7.

Numerical methods

We discretize the space domain by using a structured triangular mesh of 8192 elements, that is a space step of ∆x = ∆y = 8.476472.

For the discretization of ψ sc , we use a P 1 -finite element method. The optimized parameter q is discretized through a P 0 -finite element method. Hence, on each triangle, the approximation of ψ sc is determined by three nodal values, located at the edges of the triangle, and the approximation of q is determined by one nodal value, placed at the center of gravity of the triangle.

On the other hand, we perform the optimization through a subspace trust-region method, based on the interior-reflective Newton method described in [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF] and [START_REF] Coleman | On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds[END_REF]. Each iteration involves the solving of a linear system using the method of preconditioned conjugate gradients, for which we supply the Hessian multiply function. The computations are achieved with MATLAB (version 9.4.0.813654 (R2018a)).

Example 1: a wave damping problem

We first consider the minimization of the cost functional

where Ω 0 = [ L 6 , 5L 6 ] 2 is the domain where the waves are to be damped. The bathymetry is only optimized on a subset Ω q = [ L 4 , 3L 4 ] 2 ⊂ Ω 0 .

Example 2: an inverse problem

Given the bathymetry

, where τ = 10 -3 , we use the previous methodology to reconstruct it on the domain

, by minimizing the cost functional

where ψ ref is the amplitude associated with

Note that in this case, Ω q is not contained in Ω 0 .

In Figure 3.3, we observe that the part of the bathymetry that does not belong to the observed domain Ω 0 is not recovered by the procedure. On the contrary, the bathymetry is well reconstructed in the part of the domain corresponding to Ω 0 .

Perspectives

The numerical examples highlight two situations in which our approach could be applied. Concerning batrymetry reconstruction, we believe that our standpoint would help in obtaining more precise estimates from observed free surface data. On the other hand, a damped flow is desirable when designing structures as harbors, since minimizes the damage generated by the waves and improves navigation conditions. Then, such a problem a problem could certainly be studied by our analysis.

Finally, the method presented here uses a monochromatic incident wave, focusing entirely on its amplitude. A natural extension is then consider a polychromatic wave, that can be splitted into a sum of monochromatic waves via Fourier transform. Ideally, this should lead us to a family of optimization problems associated with each frequency (that can be solved in parallel), but several questions have to be addressed in between, regarding first a possible decomposition of the cost functional and then the convergence of the whole procedure. Since we work with the Helmholtz equation, this idea cannot be extended to nonlinear wave propagation models as Saint-Venant or Boussinesq. In this case, we suggest to consider the time-dependent problem, associated with time-periodic boundary conditions.

Simplified model

In this simple setting, a reformulation of Equations (4.13-4.15) can be obtained after a short algebraic manipulation. Theorem 4.2.1. Suppose that Assumption 4.1.1 holds and that (ϕ, a, a 13-4.15). Then ϕ satisfies

where

Reciprocally, suppose that ϕ ∈ I -{0, π 2 } satisfies Eq. (4.27) and define a and a as the corresponding solutions of Eqs. (4.14) and (4.15), respectively. Then (ϕ, a, a 13-4.15). Equation (4.27) appears -up to a factor -in [68, p.128, Equation (3.85a)]. We see that the intervals I β,γ λ and (-π 2 + θ λ , π 2 + θ λ ) play similar roles in the local and macroscopic descriptions, respectively, as they both correspond to domains of definition. In the same way, I β,γ λ ∩ R + and (0, θ λ ], whose intersection is I + , correspond to angles associated with positive lift in the two descriptions.

We have excluded the angles ϕ = 0 and ϕ = π 2 for the sole reason that Equations (4. 13-4.15) are not defined correctly for these values. Actually, ϕ = 0 is naturally associated with the case a = 1. On the contrary, the value ϕ = π 2 (that belongs to

2 ) is neither a solution of Equations (4.13-4.15) nor of (4.27): it leads indeed to a = -1, a = 0 and -λ = µ L ( π 2 ) which corresponds to a negative lift, contradicting Assumption 4.1.1. Consequently, both values can be included in the interval of admissible angles without introducing a spurious solution. In addition, these angles do not pose any particular problems concerning the definition of µ G , so that the formulation (4.27) will be considered on I in the rest of this chapter. Note also that ϕ = θ λ ± π 2 are neither solutions of Equations (4.13-4.15) nor of (4.27), since they lead to an absurd.

Proof. Suppose that (ϕ, a, a ) ∈ I ×R-{1}×R-{-1} satisfies Equations (4.21-4.23). We have to prove that ϕ satisfies (4.27). Eliminating a and a in (4.13), using (4.14) and (4.15), we get

so that (4.27) follows from Definition (4.25) of θ λ . Repeating these steps backward ends the proof of the equivalence.

Defining ν(ϕ) = 1 -τ (ϕ), we expand (4.31) in a neighborhood of ϕ = 0 + to get

from which we obtain, after some computations,

We consider a sequence (ϕ n ) n∈N satisfying lim

we obtain that the sequence

Back to (4.37), we find a contradiction since the left-hand side goes to +∞ whereas the right-hand side is constant. It follows that, up to a subsequence, we can assume that lim Remark 4.2.7. In the case C D (0) = 0, a similar reasoning gives:

The quantity µ c D (0) has no specific physical meaning in the applications. We have introduced it as a constant (that can be written explicitly), for simplicity of presentation. The angle ϕ = 0 has rather a meaning from the macroscopic point of view, as is it appears when considering µ G that cancels in 0 and µ c G , see (4.39) hereafter.

We are now in the position to give an existence result about the corrected model.

Corollary 4.2.8 (of Theorem 4.2.5). Suppose that assumptions of Lemma 4.2.6 hold and that

.38)

Then Equation (4.33) admits a solution in I + corresponding to a positive lift. Moreover if g is decreasing, max I + = θ λ and ϕ c < θ λ , then the assumption (4.38) is automatically satisfied.

Stall

As mentioned in Section 4.1.1, the stall phenomenon is generally associated with a sudden decrease in C L . It follows that if the stall angle α s satisfies α s + γ λ ∈ I ,the graph of µ L shall cross the graph of µ G at an angle in ϕ ≥ α s + γ λ . This fact is quoted in [68, p.139]: "In the stall region [...] there may be multiple solutions for C L . Each of these solutions is possible. The correct solution should be that which maintains the continuity of the angle of attack along the blade span."

Corrected model

Though lim ϕ→0 + µ c G (ϕ) = +∞, the graph of µ c G may no longer be concave when a correction is used for values of a close to 1. Hence possible multiple solutions, e.g. in the case µ L is affine.

Solution algorithms

In this section, we focus on the numerical solving of Glauert's model. In the literature, one algorithm is particularly highlighted: it consists in a fixed point iteration applied on the three equations of the model, i.e. either Equations (4. 13-4.15) or Equations (4.21-4.23), that we describe in Section 4.3.1. Many articles note that this algorithm is unstable in some cases, preventing convergence and ultimately the solving of the problem. In Section 4.3.2, we propose a new algorithm, for which we prove convergence in a less restrictive framework.

Usual algorithm

The following procedure is broadly used to solve Equations (4. 21-4.23). An early presentation of this algorithm is given in [START_REF] Wilson | Aerodynamic performance of wind turbines[END_REF], but can also be found in several monographs [START_REF] Branlard | Wind turbine aerodynamics and vorticity-based methods: fundamentals and recent applications[END_REF][START_REF] Manwell | Wind Energy Explained: Theory, Design and Application[END_REF][START_REF] Sørensen | General Momentum Theory for Horizontal Axis Wind Turbines[END_REF][START_REF] Hansen | Aerodynamics of Wind Turbines[END_REF]. In particular, the version given in [START_REF] Song | BEM simulation and performance analysis of a small wind turbine rotor[END_REF] includes a correction for high values of a and corresponds to the following algorithm. Note that the stopping criterium in Step 4 is arbitrary and usually not mentioned in monographs. The convergence of this algorithm is problematic. Instabilities are often observed in practice, as quoted in [START_REF] Song | BEM simulation and performance analysis of a small wind turbine rotor[END_REF]: "Note that this set of equations must be solved simultaneously, and in practice, numerical instability can occur" and "When local angle of attack is around the stall point, or becomes negative, getting the BEM code to converge can become difficult."

We also refer to [START_REF] Maniaci | AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, chapter An Investigation of WT_Perf Convergence Issues[END_REF] for a specific study on the convergence issues, as well as Appendix 4.A for additional assumptions that guarantee its convergence.

Algorithm 4.1: BEM, usual procedure

Input:

Initial guess: a, a .

Output: a, a , ϕ.

Set err := Tol + 1. Define the functions µ c L and µ c D by (4.24) using the input data.

while err > Tol do

2. Set a = τ (ϕ), i.e. the solution of Equation (4.31).

Set

4. Set err := tan ϕ -1-a λ(1+a ) . end

Alternative algorithms

The corrected model involves many cases that make the description and analysis of a general algorithm difficult. We can however propose some strategies and systematic approaches to compute a solution of different versions of the corrected model.

Fixed-point algorithm associated with (4.33)

If the correction for high values of a is not considered, µ c G = µ G and an alternative fixed-point iteration can be proposed: assuming that µ c L is differentiable and given an initial value ϕ 0 , we define the sequence (ϕ k ) k∈N by

with

, where ρ > 0.

We then have the following result of convergence in the case ϕ = 0.

Appendix 4.A Convergence in the simplified case

In the simplified case, Algorithm 4.1 can be summarized as the calculation of the terms of a sequence (ϕ k ) k∈N associated with the following recursion:

with f (x) := π 2 -atan (λ + µ L (x)h(x)) and h(x) := λ tan -1 x+1 sin x

.

The stability can actually be obtained in the simplified case with additional assumptions. If the initial value ϕ 0 belongs to [γ λ , θ λ ], then the sequence defined by (4.53) satisfies:

Proof. The result can be obtained by induction. Assume that for some k ∈ N,

Because of (4.53), we have tan -1 ϕ k+1 := λ + µ L (ϕ k )h(ϕ k ), so that, since µ L ≥ 0 and by Definition (4.25) of θ λ , then ϕ k+1 ≤ θ λ . On the other hand, thanks to the assumption (4.54), we obtain

Because of (4.54), the left-hand side of the last inequality is bounded by tan -1 γ λ . The result follows.

To get a sufficient condition for convergence, the last result must be completed by a contraction property. This is the object of the following result. Lemma 4.A.2. Suppose that µ L is differentiable and denote by µ L its derivative. The derivative of f satisfies

Proof. Differentiating f , we find that

MOTS CLÉS

Algorithme pararéel, Assimilation parallèle de données, Optimisation d'une bathymétrie, Équation de Helmholtz, Méthode de l'élément de pale, Design des pales.

RÉSUMÉ

La présente thèse vise à contribuer à l'élaboration d'un cadre théorique pour trois problèmes dans le contexte des énergies marines renouvelables. Dans sa première partie, nous proposons une procédure pour coupler des méthodes d'assimilation de données temporelles non limitées avec des algorithmes parallèles en temps. La combinaison entre l'observateur de Luenberger et l'algorithme Pararéel est étudiée, ce qui permet d'estimer le nombre d'itérations pararéelles nécessaires pour préserver le taux de convergence de l'observateur et d'obtenir une estimation de l'efficacité théorique de l'ensemble de la procédure.

Nous discutons ensuite la détermination d'une bathymétrie dans une perspective d'optimisation. En imposant que la propagation des vagues optimise un certain critère associé à une fonctionnelle de coût, nous considérons un problème d'optimisation sous contrainte d'EDP où la bathymétrie joue le rôle de contrôle et la propagation des vagues est décrite par une équation de type Helmholtz. Nous sommes en mesure de prouver, sur la base d'hypothèses appropriées, la continuité de la fonction contrôle-état et l'existence d'une solution optimale, incluant aussi quelques résultats sur les solutions au problème de Helmholtz et la convergence dans un cadre discret. Ce travail est complété par des expériences numériques.

La dernière partie de ce travail est consacrée à l'analyse de la méthode de l'élément de pale (BEM), une méthode classique utilisée pour déterminer les performances d'une hélice ainsi que des paramètres de design. Nous proposons une reformulation de la méthode qui permet d'obtenir des conditions d'existence des solutions et d'établir la convergence de certains algorithmes de résolution. Nous étudions également le problème d'optimisation associé dans certains contextes.

ABSTRACT

The present thesis aims to contribute to the development of a theoretical framework for three problems in the context of renewable marine energy. In the first part, we propose a procedure to couple unbounded in time data assimilation methods with time-parallel algorithms. The combination between the Luenberger observer and Parareal algorithm is studied, providing a way to estimate the number of parareal iterations required to preserve the observer rate of convergence, as well as an estimation of the theoretical efficiency of the entire procedure.

We then discuss the determination of a bathymetry from an optimization perspective. Imposing that wave propagation must fulfill a certain criterion associated with a cost functional, we consider a PDE-constrained optimization problem where the bathymetry plays the role of control and wave propagation is described by the Helmholtz equation. We are able to prove, under suitable assumptions, the continuity of the control-to-state mapping and the existence of an optimal solution, including also some results about solutions to Helmholtz problem and convergence in a discrete framework. This work is complemented by numerical experiments.

The last part of this work is devoted to analyze the convergence of the Blade element momentum (BEM) theory, a classical method used to determine the propeller efficiency as well as its design parameters. We propose a reformulation of the method that allows to obtain conditions for existence of solutions and establish the convergence of some solving algorithms. We also study the associated optimization problem in certain contexts.

KEYWORDS

Parareal algorithm, Parallel data assimilation, Bathymetry optimization, Helmholtz equation, Blade element method, Blade design.