%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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.


