%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Le cas des algorithmes ne respectant pas la symétrie\label{sec:inversion}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lorsque la discrétisation ne préserve pas la symétrie du transport, on a le choix pour les applications entre l'algorithme adjoint ou l'algorithme de rétro-transport. Les possibles implications sont illustrées et discutées dans le cadre de calculs d'advection unidimensionnelle avec le schéma de Van Leer. \subsection{Calculs de sensibilité} \def\NT{N} \def\dt{\delta_t} \def\bq{{\breve{c}}} % q au bord de la maille \def\q{c} Pour l'advection unidimensionnelle par un champ de vent non divergeant $u>0$, sur une grille régulière, l'évolution temporelle de la concentration du traceur dans la maille $i$ est donnée par \begin{equation} \label{eq:up1d} \q_i^{n+1}-\q_i^n=\alpha \dep{\bq_{i-\dem}^n-\bq_{i+\dem}^n} \mbox{\ \ \ \ avec\ } \alpha=U/m \end{equation} où $\bq_{i-\dem}^n$ est la concentration moyenne dans la partie de l'air transférée au cours du pas de temps entre les mailles $i-1$ et $i$. Pour le schéma I de Van Leer, on utilise une approximation linéaire de la distribution sous-maille avec une pente $\dq_i$ correspondant à des concentrations $\q_i^\pm=\q_i\pm\dq_i/2$ aux bords de la maille (se reporter à la \fig{schema2b} pour le principe du schéma et les notations). $\bq_{i+\dem}$ vaut alors \begin{eqnarray}\label{eq:bqim} \bq_{i+1/2} & = & \q_i + \frac{1}{2}\dep{1-\alpha} \dq_i \mbox{\ si\ } \alpha >0\\\label{eq:bqip} & = & \q_{i+1} - \frac{1}{2}\dep{1+\alpha} \dq_{i+1} \mbox{ , sinon.} \end{eqnarray} Dans le schéma I de Van Leer, la pente est d'abord calculée par différences finies avant d'appliquer un limiteur de pente de sorte que le calcul de la pente peut s'écrire sous la forme de la séquence d'opérations et de tests suivante~: \begin{eqnarray} \mbox{si } & \depc{\ (c_{i+1}-c_{i})\times(c_{i}-c_{i-1}) >0\ } \mbox{\ alors\ } \label{eq:nomax} \\ & \dq_i=(\q_{i+1}-\q_{i-1})/2 &\label{eq:slope}\\ & \mbox{si } \depc{\ |\dq_i| > 2 | \q_{i}-\q_{i-1}|\ } \mbox{alors\ } &\dq_i=2 ( \q_{i}-\q_{i-1} )\label{eq:gauche}\\ & \mbox{si } \depc{\ |\dq_i| > 2 | \q_{i+1}-\q_i|\ } \mbox{alors\ } &\dq_i=2 ( \q_{i+1}-\q_i)\label{eq:droite}\\ \mbox{sinon}& \dq_i=0 \end{eqnarray} \def\EPS1D{/home/hourdin/TEX/ARTICLES/RETRO/1DTESTS/INFLUENCEB} \def\EPS1D{\local/FIGURES} \def\tmpsz{6.cm} \begin{figure} \centerline{ \begin{tabular}{|c|c|} \multicolumn{2}{c}{{\bf A) TRACEUR DIRECT }}\\\hline Distribution initiale sinusoïdale & Distribution initiale carrée \\\hline \includegraphics[width=\tmpsz]{\EPS1D/dirvl2E1mu2.eps} & \includegraphics[width=\tmpsz]{\EPS1D/dirvl2E2mu2.eps} \\\hline \multicolumn{2}{c}{}\\ \multicolumn{2}{c}{}\\ \multicolumn{2}{c}{{\bf B) FONCTION D'INFLUENCE AU PAS DE TEMPS INITIAL}}\\\hline Mesure au pas de temps 30 & Mesure au pas de temps 150 \\\hline &\\ \multicolumn{2}{|c|}{Cas 1: Etat initial plat, schéma Van Leer I}\\ \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt30.E0.eps} & \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt150.E0.eps} \\\hline &\\ \multicolumn{2}{|c|}{Cas 2: Etat initial sinusoïdal, Schéma Van Leer I}\\ \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt30.E1.eps} & \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt150.E1.eps} \\\hline &\\ \multicolumn{2}{|c|}{Cas 3: Etat initial carré, Schéma Van Leer I}\\ \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt30.E2.eps} & \includegraphics[width=\tmpsz]{\EPS1D/vl2.mu2.nt150.E2.eps} \\\hline &\\ \multicolumn{2}{|c|}{Cas 4: Etat initial carré, Schéma Van Leer II}\\ \includegraphics[width=\tmpsz]{\EPS1D/vl1.mu2.nt30.E2.eps} & \includegraphics[width=\tmpsz]{\EPS1D/vl1.mu2.nt150.E2.eps} \\\hline \end{tabular} } \caption{ Calculs de sensibilités ou fonctions d'influence. \label{fg:1d}} Partie A. Distribution initiale de traceur (courbe continue) et distribution obtenue après 150 pas d'advection (courbe discontinue) pour des distributions initiales sinusoïdale (à gauche) et carrée (à droite). Partie B. Fonction d'influence à l'instant initial pour une mesure effectuée au pas de temps 30 (colonne de gauche) ou 150 (colonne de droite). Pour chaque graphique, la courbe épaisse continue montre la distribution de mesure. Les autres courbes montrent la sensibilité calculée avec l'algorithme de rétro-transport (courbe épaisse grise), avec l'adjoint du code direct (croix) et par perturbations de simulations directes (carrés). Les différents cas (de haut en bas) diffèrent par la distribution initiale de traceur ou par le schéma utilisé pour l'advection. \end{figure} Nous montrons ci-dessous les résultats obtenus avec trois méthodes différentes pour calculer la sensibilité d'une mesure par rapport à la distribution initiale de traceur. Les deux premières méthodes consistent en l'intégration à rebours de l'algorithme de rétro-transport ou de l'algorithme adjoint. Le troisième calcul est obtenu par des perturbations successives de l'état initial suivies de simulations directes. Cette méthode requière autant d'intégrations que de points de grille. La même méthode de perturbation a aussi été testée avec le modèle linéaire tangent. Comme il se doit, les résultats obtenus ainsi sont identiques à ceux de l'intégration adjointe et ces résultats ne sont donc pas montrés. Pour un modèle linéaire et symétrique en temps, les trois estimations de la sensibilité doivent être égales. Ce point a été vérifié avec le schéma de Godunov (résultats non montrés). Pour un algorithme non symétrique mais linéaire, l'estimation par perturbations et l'estimation adjointe sont identiques mais diffèrent du calcul de rétro-transport. A cause des limiteurs de pentes (\ref{eq:nomax}, \ref{eq:gauche} et \ref{eq:droite}), le schéma I de Van Leer n'est pas linéaire. Dans ce cas, les trois estimations sont différentes. L'estimation par rétro-transport ne dépend que de la distribution de mesure alors que l'intégration adjointe dépend du calcul de base au voisinage duquel on calcule les sensibilités. Quant au calcul direct par perturbation, il dépend à la fois de la solution de base et de l'amplitude de la perturbation initiale. A cause des non-linéarités, l'intégration de l'adjoint du modèle nécessite le stockage de la solution de base qui est ensuite relue et utilisée pour activer les propositions conditionnelles dans les calculs (dans \ref{eq:nomax}, \ref{eq:gauche} et \ref{eq:droite}) au cours du calcul de sensibilité. Les résultats sont montrés sur la \fig{1d}. Le domaine est périodique avec 60 points et on a choisi un nombre de Courant $\alpha=$0,2. Trois distributions initiales de traceurs sont testées~: pas de traceur (état plat) ainsi qu'une onde sinusoïdale et une distribution carrée, toutes deux d'amplitude 1. Les deux graphiques du haut de la figure (partie A) montrent la distribution initiale et les résultats de l'intégration directe de base au pas de temps 150 pour les distributions sinusoïdale (à gauche) et carrée (à droite). La solution exacte de l'advection au pas 150 est une translation de $\alpha\NT$=30 points de grille. La partie B de la même figure montre des calculs de sensibilité ou fonctions d'influence. La mesure consiste en un prélèvement uniforme sur un intervalle de 6 points de grille. La distribution de mesure correspondante (\ie\ la fonction $\mu$ de l'\eq{cout}) correspond à la courbe noire épaisse sur tous les graphiques. La sensibilité est calculée par rapport à la concentration initiale de traceur pour une mesure effectuée au pas de temps 30 (pour la colonne de gauche) et 150 (colonne de droite). La sensibilité exacte serait obtenue en translatant la distribution de mesure de -6 et -30 points respectivement. Pour l'estimation directe, une perturbation d'amplitude $g=10^{-8}$ est ajoutée successivement à chaque point de grille. Pour le cas 1 (graphiques du haut de la partie B), qui correspond à un état initial sans traceur, la condition (\ref{eq:nomax}) n'est jamais satisfaite. La pente est donc nulle dans le calcul adjoint, qui équivaut alors à un calcul avec le schéma de Godunov (Eqs \ref{eq:up1d}, \ref{eq:bqim} et \ref{eq:bqip} avec $\dq=0$). La sensibilité correspondante est beaucoup plus diffuse que celle obtenue avec l'algorithme rétro ou le calcul par perturbations, ces deux derniers étant à la fois plus proches l'un de l'autre et plus proche de la solution exacte. Ce cas classique, où l'on essaie de reconstituer une source en partant d'une concentration supposée nulle a priori, est en fait ici un cas limite singulier dans lequel la sensibilité dépend de choix arbitraires sur l'écriture du schéma. En changeant le signe $>$ en $\ge$ dans l'\eq{nomax}, on obtiendrait par exemple des sensibilités différentes. Le cas 2, qui correspond à une onde sinusoïdale, est moins pathologique. Pour un pas de temps donné, la condition (\ref{eq:nomax}) est satisfaite presque partout et le limiteur de pentes est actif seulement pour un petit nombre de points de grille. Les trois estimations de la sensibilité sont affectées par la diffusion numérique mais sont relativement proches les unes des autres (celle calculée par rétro-transport, qui ne dépend que de la distribution de mesure, est la même que dans le cas 1). A noter que les estimations directes par perturbations et adjointes produisent de faibles valeurs négatives (on le voit particulièrement pour la sensibilité d'une mesure effectuée au pas de temps 150). Le cas 3 correspond à une distribution carrée. La sensibilité calculée par rétro-transport est à nouveau la même que précédemment. Pour une mesure au pas 30, la solution adjointe commence à montrer des oscillations numériques alors que les deux autres restent proches l'une de l'autre. Pour une mesure au pas 150, les intégrations adjointes et par perturbations produisent des oscillations bien plus grandes que le signal lui-même avec des sensibilités négatives très importantes. L'instabilité rencontrée ici peut être expliquée de la façon suivante. Après un certain nombre de pas de temps pour la simulation directe, la distribution n'est plus vraiment carrée mais décroît rapidement de part et d'autre du pic de traceur, typiquement par un ordre de grandeur d'un point de grille au suivant. De ce fait, en amont du pic, dans la région où le pic de la sensibilité est advecté (dans notre cas particulier), la condition \ref{eq:droite} est atteinte partout. Du coup, le schéma adjoint correspond en fait à un schéma direct associé à une pente $\dq_i=2 ( \q_{i+1}-\q_{i} )$. Ce schéma est pathologique comme schéma de transport à la fois parce que l'estimation de la pente est systématiquement décentrée et parce que cette pente est double d'une estimation classique par différences finies. Afin d'illustrer ce cas plus en détail, on montre, dans la partie de gauche de la \fig{1dbis}, les sensibilités directes calculées pour une mesure au pas 150 (comme pour la partie de droite de la \fig{1d}) mais en utilisant différentes valeurs pour $g$ (rappelons que les résultats montrés dans la \fig{1d} correspondent tous à $g=10^{-8}$). Pour des valeurs petites de $g$, comme on peut s'y attendre, les résultats sont proches du calcul adjoint. Quand $g$ croît, les calculs directs par perturbations tendent à se lisser et deviennent plus proches de la solution obtenue par rétro-transport (ainsi que de la solution physiquement correcte). Olivier Talagrand a proposé l'explication suivante pour ce phénomène intéressant. A cause de la présence des opérations conditionnelles associées aux limiteurs de pentes (\ref{eq:gauche}) et (\ref{eq:droite}), chaque intégration du modèle direct traverse une série de bifurcations. Pour des valeurs de perturbations suffisamment petites, la séquence des bifurcations n'est pas modifiée (une petite perturbation n'est pas ``vue" par le limiteur) et le calcul adjoint montre que la sensibilité ainsi obtenue est très grande. Quand on fait croître l'amplitude de cette perturbation, la séquence des bifurcations est de plus en plus modifiée, avec des fluctuations très fortes dans la sensibilité. Ces fluctuations tendent à s'annuler les unes les autres, ce qui n'est pas surprenant dans la mesure où les limiteurs de pentes sont justement introduit pour garantir le bon comportement physique du schéma de transport. Il est remarquable que le rétro-transport (qui inclue bien sûr les limiteurs de pente) atteint le même but pour un coût bien moindre. Le cas 4 correspond aussi à une distribution initiale carrée, mais l'intégration est effectuée avec le second schéma de Van Leer, déjà décrit dans la \sec{limiteurs}, dans lequel la pente (\ref{eq:slope}) est remplacée par $\dq_i = 2\dep{\q_{i+1}-\q_i}\dep{\q_i-\q_{i-1}}/\dep{\q_{i+1}-\q_{i-1}}$. Dans ce cas, les conditions \ref{eq:gauche} et \ref{eq:droite} ne sont jamais atteintes. Ce schéma est légèrement plus diffusif que le schéma I mais son comportement est moins pathologique car sa non-linéarité est moins forte. On observe encore cependant le même effet que pour le cas précédent bien que dans une moindre mesure. Pour $g\le 0,0001$ et à la précision de la figure, la sensibilité adjointe et celle obtenue par perturbation sont confondues. Pour $g=1$, la sensibilité directe est en revanche plus proche du calcul par rétro-transport. D'autres résultats (non montrés) confirment que le calcul de rétro-transport fournit en générale une bonne estimation de la sensibilité calculée par perturbation du calcul direct avec des perturbations de relativement grande amplitude. \begin{figure} \centerline{ \begin{tabular}{cc} Van Leer I, pas de temps 150 & Van Leer II, pas de temps 150 \\ \includegraphics[width=\tmpsz]{\EPS1D/zoomvl2.mu2.nt150.E2.eps} & \includegraphics[width=\tmpsz]{\EPS1D/zoomvl1.mu2.nt150.E2.eps}\\ \end{tabular}} \caption{Calcul de sensibilités par perturbations directes de la distribution initiale pour les schéma I (à gauche) et II (à droite) de Van Leer (mêmes conditions que pour les deux graphiques les plus bas de la colonne de droite de la \fig{1d}). Les sensibilités sont calculées pour différentes amplitudes de la perturbation $g$. Les courbes pour $g=10^{-8}$ ainsi que les calculs adjoint et rétro sont aussi montrés sur la \fig{1d}. \label{fg:1dbis}} \end{figure} On voit que dans le cas où la sensibilité de l'algorithme direct est définie de manière non ambiguë, le calcul adjoint en fournit, comme il se doit, une bonne estimation. Cependant, pour des algorithmes fortement non-linéaires, cette estimation, exacte d'un point de vue numérique, peut s'avérer non physique (oscillations, valeurs négatives). En comparaison, le calcul par rétro-transport, robuste et préservant la positivité (ainsi que la monotonie), fournit des sensibilités qui restent réalistes et proches de la solution exacte même pour des schémas fortement non linéaires. \subsection{Expériences de minimisation} \def\ct{y} \def\obs{o} On compare à présent les différentes approches dans le contexte de l'assimilation variationnelle, a priori plus favorable au calcul adjoint, avec pour but de reconstruire la distribution initiale de traceur à partir d'observations distribuées dans le temps. On effectue des tests académiques classiques de type ``expériences jumelles", standards pour l'évaluation des algorithmes d'assimilation. On commence par effectuer une simulation de référence, notée ${\ct}_i^n$, sur $\NT$ pas de temps. Cette expérience est considérée comme la réalité à reconstruire. Des observations synthétiques sont générées à chaque pas de temps $n$ à partir de cette simulation sous la forme $\sum_i\mu_i{\ct}_i^n$, où $\mu=(\mu_i)$ est à nouveau la distribution de mesure. Pour toute solution ${c}_i^n$ de l'équation de transport, la fonction objective (scalaire) \begin{equation}\label{eq:J} J\dep{\ct,c}=\sum_n \depb{\sum_i{\mu_i\ct}_i^n-\sum_i\mu_i{c}_i^n}^2 \end{equation} mesure l'écart entre cette solution et l'observation. On minimise alors cette fonction coût par rapport à la concentration initiale ${c}_i^0$. La minimisation est effectuée grâce au code M1QN3 développé par \cite{Gilb:89}. Il s'agit d'une procédure itérative qui nécessite, à chaque pas, au moins une approximation du gradient de la fonction objective. L'algorithme est de type quasi-Newtonien, ce qui signifie que le gradient local est utilisé pour bâtir petit à petit une approximation de l'inverse du Hessien (matrice des dérivées secondes) de la fonction objective. L'utilisation de cet inverse rend la minimisation particulièrement efficace, au moins dans le cas où le gradient varie de façon relativement douce. Deux séries d'expériences ont été effectuées ici, en utilisant, pour calculer le gradient, dans le premier cas l'adjoint exact du modèle direct et dans le second cas l'adjoint approché obtenu par rétro-transport. Dans toutes les expériences, effectuées dans le même cadre unidimensionnel que pour la section précédente, la longueur de la simulation de référence est de $\NT=300$ pas de temps, ce qui correspond à une révolution complète sur le domaine d'advection. La fonction de mesure est également comme avant un prélèvement uniforme sur 6 points de grille (courbe épaisse discontinue sur la \fig{min}). Avec ces choix, le minimum de la fonction objective est unique et égal à zéro. \def\EPSM{/home/hourdin/TEX/ARTICLES/RETRO/1DTESTS/MINIMISATION/MIN4} \def\EPSM{\local/FIGURES} \begin{figure} \centerline{\begin{tabular}{cc} {\bf a)} Fonction objective & {\bf b)} Erreur quadratique moyenne \\ \includegraphics[angle=-90,width=7cm]{\EPSM/cost.eps} & \includegraphics[angle=-90,width=7cm]{\EPSM/err.eps} \\ {\bf c)} Concentration initiale après minimisation & {\bf d)} Distribution de traceur au pas de temps 50 \\ \includegraphics[width=7.2cm]{\EPSM/qi.eps} & \includegraphics[width=7.2cm]{\EPSM/q50.eps} \end{tabular}} \caption{ Résultats d'expérience d'assimilation variationnelle avec des observations synthétiques. L'état initial ``réel" à reconstruire correspond à la courbe épaisse continue du graphique {\bf c}. La fonction objective (\eq{J}) à minimiser est calculée avec la distribution de mesure représentée par une courbe épaisse discontinue sur les graphiques {\bf c} et {\bf d}. La minimisation est commencée à partir d'un état initial sans traceurs. \label{fg:min}} {\bf a}~: Evolution de la valeur de la fonction objective au fil de la minimisation.\\ {\bf b}~: Evolution de l'erreur d'estimation des états initial et final de l'intervalle d'assimilation.\\ {\bf c}~: Concentrations initiales ``réelle" et reconstituée à la fin de la minimisation.\\ {\bf d}~: Les mêmes états que pour le graphique {\bf c} mais au pas de temps $n=50$ de l'intervalle d'assimilation. \end{figure} Les résultats sont présentés sur la \fig{min}. L'état initial de la simulation de référence utilisé comme réalité correspond à la distribution carrée déjà utilisée précédemment (courbe épaisse sur les graphiques du bas). On utilise le schéma I de Van Leer. La minimisation est initiée à partir d'un état sans traceurs. La minimisation est poursuivie jusqu'à ce que la décroissance de la fonction objective entre deux états successifs de la minimisation soit plus petite qu'un seuil prescrit. Le graphique {\bf a)} montre que, comme on peut s'y attendre, la minimisation aboutit à une valeur plus faible de la fonction objective quand on calcule le gradient avec l'adjoint exact plutôt qu'avec l'algorithme de rétro-transport. Cependant, la minimisation est plus rapide au début avec le rétro-transport. La reconstruction de la concentration initiale (graphique {\bf c)} est également très légèrement meilleure avec l'algorithme de rétro-transport (ce qui se voit aussi sur l'erreur quadratique moyenne montrée sur le graphique {\bf b)}. En fait, si les deux reconstructions obtenues en fin de minimisation diffèrent à l'instant initial, elles sont en revanche quasiment identiques au pas d'advection $n=50$ (graphique {\bf d)}. C'est ici la diffusion numérique qui rend difficile la reconstruction de l'état initial, et qui explique que les oscillations de l'état initial (graphique {\bf c)} ne soit pas ``vues" par les mesures. Le même effet se produirait dans la réalité avec la diffusion turbulente. D'autres résultats (non montrés) confirment les résultats obtenus ici à savoir que la minimisation est généralement plus poussée avec l'adjoint mais plus rapide au début avec l'algorithme de rétro-transport. Les différences concernant la reconstruction de l'état initial sont souvent non significatives, mais à nouveau, cette reconstruction est souvent légèrement meilleure avec l'algorithme de rétro-transport. Ces résultats sont cohérents avec ceux présentés sur la \fig{1dbis}. Le gradient approché mais lisse obtenu avec l'algorithme de rétro-transport permet une localisation imparfaite mais rapide du minimum de la fonction objective. L'utilisation du gradient exact permet une localisation plus précise de ce minimum. Mais ce gain ne semble pas utile en pratique. De plus, dans certains cas (non montrés), ce gradient exact peut osciller tellement que la minimisation échoue. On peut en tirer la conclusion que, au moins sur les exemples académiques présentés ici, l'utilisation d'un gradient certes approché, mais lisse et présentant un bon comportement physique, est préférable à l'utilisation d'un gradient exact mais présentant de fortes oscillations.