2 Copyright (C) 2008-2018 EDF R&D
4 This file is part of SALOME ADAO module.
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
11 This library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 .. _section_ref_operator_requirements:
26 Exigences pour les fonctions décrivant un opérateur
27 ---------------------------------------------------
29 Les opérateurs d'observation et d'évolution sont nécessaires pour mettre en
30 oeuvre les procédures d'assimilation de données ou d'optimisation. Ils
31 comprennent la simulation physique par des calculs numériques, mais aussi le
32 filtrage et de restriction pour comparer la simulation à l'observation.
33 L'opérateur d'évolution est ici considéré dans sa forme incrémentale, qui
34 représente la transition entre deux états successifs, et il est alors similaire
35 à l'opérateur d'observation.
37 Schématiquement, un opérateur doit donner une solution étant donné les
38 paramètres d'entrée. Une partie des paramètres d'entrée peut être modifiée au
39 cours de la procédure d'optimisation. Ainsi, la représentation mathématique d'un
40 tel processus est une fonction. Il a été brièvement décrit dans la section
41 :ref:`section_theory` et il est généralisée ici par la relation:
43 .. math:: \mathbf{y} = O( \mathbf{x} )
45 entre les pseudo-observations :math:`\mathbf{y}` et les paramètres
46 :math:`\mathbf{x}` en utilisant l'opérateur d'observation ou d'évolution
47 :math:`O`. La même représentation fonctionnelle peut être utilisée pour le
48 modèle linéaire tangent :math:`\mathbf{O}` de :math:`O` et son adjoint
49 :math:`\mathbf{O}^*`, qui sont aussi requis par certains algorithmes
50 d'assimilation de données ou d'optimisation.
52 En entrée et en sortie de ces opérateurs, les variables :math:`\mathbf{x}` et
53 :math:`\mathbf{y}` ou leurs incréments sont mathématiquement des vecteurs, et
54 ils sont donc passés comme des vecteurs non-orientés (de type liste ou vecteur
55 Numpy) ou orientés (de type matrice Numpy).
57 Ensuite, **pour décrire complètement un opérateur, l'utilisateur n'a qu'à
58 fournir une fonction qui réalise uniquement l'opération fonctionnelle de manière
61 Cette fonction est généralement donnée comme un script qui peut être exécuté
62 dans un noeud YACS. Ce script peut aussi, sans différences, lancer des codes
63 externes ou utiliser des appels et des méthodes internes SALOME. Si l'algorithme
64 nécessite les 3 aspects de l'opérateur (forme directe, forme tangente et forme
65 adjointe), l'utilisateur doit donner les 3 fonctions ou les approximer.
67 Il existe 3 méthodes effectives pour l'utilisateur de fournir une représentation
68 fonctionnelle de l'opérateur. Ces méthodes sont choisies dans le champ "*FROM*"
69 de chaque opérateur ayant une valeur "*Function*" comme "*INPUT_TYPE*", comme le
70 montre la figure suivante:
72 .. eficas_operator_function:
73 .. image:: images/eficas_operator_function.png
77 **Choisir une représentation fonctionnelle de l'opérateur**
79 Première forme fonctionnelle : utiliser "*ScriptWithOneFunction*"
80 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
82 .. index:: single: ScriptWithOneFunction
83 .. index:: single: DirectOperator
84 .. index:: single: DifferentialIncrement
85 .. index:: single: CenteredFiniteDifference
87 La première consiste à ne fournir qu'une seule fonction potentiellement non
88 linéaire, et d'approximer les opérateurs tangent et adjoint. Ceci est fait en
89 utilisant le mot-clé "*ScriptWithOneFunction*" pour la description de
90 l'opérateur choisi dans l'interface graphique ADAO. L'utilisateur doit fournir
91 la fonction dans un script, avec un nom obligatoire "*DirectOperator*". Par
92 exemple, le script peut suivre le modèle suivant::
94 def DirectOperator( X ):
95 """ Opérateur direct de simulation non-linéaire """
101 Dans ce cas, l'utilisateur doit aussi fournir une valeur pour l'incrément
102 différentiel (ou conserver la valeur par défaut), en utilisant dans l'interface
103 graphique (GUI) le mot-clé "*DifferentialIncrement*", qui a une valeur par
104 défaut de 1%. Ce coefficient est utilisé dans l'approximation différences finies
105 pour construire les opérateurs tangent et adjoint. L'ordre de l'approximation
106 différences finies peut aussi être choisi à travers l'interface, en utilisant le
107 mot-clé "*CenteredFiniteDifference*", avec 0 pour un schéma non centré du
108 premier ordre (qui est la valeur par défaut), et avec 1 pour un schéma centré du
109 second ordre (qui coûte numériquement deux fois plus cher que le premier ordre).
110 Si nécessaire et si possible, on peut :ref:`subsection_ref_parallel_df`. Dans
111 tous les cas, un mécanisme de cache interne permet de limiter le nombre
112 d'évaluations de l'opérateur au minimum possible du point de vue de l'exécution
113 séquentielle ou parallèle des approximations numériques des opérateurs tangent
114 et adjoint, pour éviter des calculs redondants.
116 Cette première forme de définition de l'opérateur permet aisément de tester la
117 forme fonctionnelle avant son usage dans un cas ADAO, réduisant notablement la
118 complexité de l'implémentation de l'opérateur. On peut ainsi utiliser
119 l'algorithme ADAO de vérification "*FunctionTest*" (voir la section sur
120 l':ref:`section_ref_algorithm_FunctionTest`) pour ce test.
122 **Avertissement important :** le nom "*DirectOperator*" est obligatoire, et le
123 type de l'argument ``X`` peut être une liste, un vecteur ou une matrice Numpy.
124 La fonction utilisateur doit accepter et traiter tous ces cas.
126 Seconde forme fonctionnelle : utiliser "*ScriptWithFunctions*"
127 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
129 .. index:: single: ScriptWithFunctions
130 .. index:: single: DirectOperator
131 .. index:: single: TangentOperator
132 .. index:: single: AdjointOperator
134 **En général, il est recommandé d'utiliser la première forme fonctionnelle
135 plutôt que la seconde. Un petit accroissement de performances n'est pas une
136 bonne raison pour utiliser l'implémentation détaillée de cette seconde forme
139 La seconde consiste à fournir directement les trois opérateurs liés :math:`O`,
140 :math:`\mathbf{O}` et :math:`\mathbf{O}^*`. C'est effectué en utilisant le
141 mot-clé "*ScriptWithFunctions*" pour la description de l'opérateur choisi dans
142 l'interface graphique (GUI) d'ADAO. L'utilisateur doit fournir trois fonctions
143 dans un script, avec trois noms obligatoires "*DirectOperator*",
144 "*TangentOperator*" et "*AdjointOperator*". Par exemple, le script peut suivre
145 le squelette suivant::
147 def DirectOperator( X ):
148 """ Opérateur direct de simulation non-linéaire """
152 return quelque chose comme Y
154 def TangentOperator( (X, dX) ):
155 """ Opérateur linéaire tangent, autour de X, appliqué à dX """
159 return quelque chose comme Y
161 def AdjointOperator( (X, Y) ):
162 """ Opérateur adjoint, autour de X, appliqué à Y """
166 return quelque chose comme X
168 Un nouvelle fois, cette seconde définition d'opérateur permet aisément de tester
169 les formes fonctionnelles avant de les utiliser dans le cas ADAO, réduisant la
170 complexité de l'implémentation de l'opérateur.
172 Pour certains algorithmes, il faut que les fonctions tangente et adjointe puisse
173 renvoyer les matrices équivalentes à l'opérateur linéaire. Dans ce cas, lorsque,
174 respectivement, les arguments ``dX`` ou ``Y`` valent ``None``, l'utilisateur
175 doit renvoyer la matrice associée.
177 **Avertissement important :** les noms "*DirectOperator*", "*TangentOperator*"
178 et "*AdjointOperator*" sont obligatoires, et le type des arguments ``X``,
179 ``Y``, ``dX`` peut être une liste, un vecteur ou une matrice Numpy.
180 La fonction utilisateur doit accepter et traiter tous ces cas.
182 Troisième forme fonctionnelle : utiliser "*ScriptWithSwitch*"
183 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185 .. index:: single: ScriptWithSwitch
186 .. index:: single: DirectOperator
187 .. index:: single: TangentOperator
188 .. index:: single: AdjointOperator
190 **Il est recommandé de ne pas utiliser cette troisième forme fonctionnelle sans
191 une solide raison numérique ou physique. Un accroissement de performances n'est
192 pas une bonne raison pour utiliser la complexité de cette troisième forme
193 fonctionnelle. Seule une impossibilité à utiliser les première ou seconde formes
194 justifie l'usage de la troisième.**
196 La troisième forme donne de plus grandes possibilités de contrôle de l'exécution
197 des trois fonctions représentant l'opérateur, permettant un usage et un contrôle
198 avancés sur chaque exécution du code de simulation. C'est réalisable en
199 utilisant le mot-clé "*ScriptWithSwitch*" pour la description de l'opérateur à
200 travers l'interface graphique (GUI) d'ADAO. L'utilisateur doit fournir un script
201 unique aiguillant, selon un contrôle, l'exécution des formes directe, tangente
202 et adjointe du code de simulation. L'utilisateur peut alors, par exemple,
203 utiliser des approximations pour les codes tangent et adjoint, ou introduire une
204 plus grande complexité du traitement des arguments des fonctions. Mais cette
205 démarche sera plus difficile à implémenter et à déboguer.
207 Toutefois, si vous souhaitez utiliser cette troisième forme, on recommande de se
208 baser sur le modèle suivant pour le script d'aiguillage. Il nécessite un fichier
209 script ou un code externe nommé ici "*Physical_simulation_functions.py*",
210 contenant trois fonctions nommées "*DirectOperator*", "*TangentOperator*" et
211 "*AdjointOperator*" comme précédemment. Voici le squelette d'aiguillage::
213 import Physical_simulation_functions
214 import numpy, logging, codecs, pickle
216 return pickle.loads(codecs.decode(data.encode(), "base64"))
219 for param in computation["specificParameters"]:
220 if param["name"] == "method":
221 method = loads(param["value"])
222 if method not in ["Direct", "Tangent", "Adjoint"]:
223 raise ValueError("No valid computation method is given")
224 logging.info("Found method is \'%s\'"%method)
226 logging.info("Loading operator functions")
227 Function = Physical_simulation_functions.DirectOperator
228 Tangent = Physical_simulation_functions.TangentOperator
229 Adjoint = Physical_simulation_functions.AdjointOperator
231 logging.info("Executing the possible computations")
233 if method == "Direct":
234 logging.info("Direct computation")
235 Xcurrent = computation["inputValues"][0][0][0]
236 data = Function(numpy.matrix( Xcurrent ).T)
237 if method == "Tangent":
238 logging.info("Tangent computation")
239 Xcurrent = computation["inputValues"][0][0][0]
240 dXcurrent = computation["inputValues"][0][0][1]
241 data = Tangent(numpy.matrix(Xcurrent).T, numpy.matrix(dXcurrent).T)
242 if method == "Adjoint":
243 logging.info("Adjoint computation")
244 Xcurrent = computation["inputValues"][0][0][0]
245 Ycurrent = computation["inputValues"][0][0][1]
246 data = Adjoint((numpy.matrix(Xcurrent).T, numpy.matrix(Ycurrent).T))
248 logging.info("Formatting the output")
249 it = numpy.ravel(data)
250 outputValues = [[[[]]]]
252 outputValues[0][0][0].append(val)
255 result["outputValues"] = outputValues
256 result["specificOutputInfos"] = []
257 result["returnCode"] = 0
258 result["errorMessage"] = ""
260 Toutes les modifications envisageables peuvent être faites à partir de cette
261 hypothèse de squelette.
263 .. _section_ref_operator_control:
265 Cas spécial d'un opérateur d'évolution avec contrôle
266 ++++++++++++++++++++++++++++++++++++++++++++++++++++
268 Dans certains cas, l'opérateur d'évolution ou d'observation doit être contrôlé
269 par un contrôle d'entrée externe, qui est donné *a priori*. Dans ce cas, la
270 forme générique du modèle incrémental est légèrement modifié comme suit:
272 .. math:: \mathbf{y} = O( \mathbf{x}, \mathbf{u})
274 où :math:`\mathbf{u}` est le contrôle sur l'incrément d'état. En effet,
275 l'opérateur direct doit être appliqué à une paire de variables :math:`(X,U)`.
276 Schématiquement, l'opérateur doit être construit comme suit::
278 def DirectOperator( (X, U) ):
279 """ Opérateur direct de simulation non-linéaire """
283 return quelque chose comme X(n+1) (évolution) ou Y(n+1) (observation)
285 Les opérateurs tangent et adjoint ont la même signature que précédemment, en
286 notant que les dérivées doivent être faites seulement partiellement par rapport
287 à :math:`\mathbf{x}`. Dans un tel cas de contrôle explicite, seule la deuxième
288 forme fonctionnelle (en utilisant "*ScriptWithFunctions*") et la troisième forme
289 fonctionnelle (en utilisant "*ScriptWithSwitch*") peuvent être utilisées.
291 Remarques complémentaires sur l'adimensionnement des opérateurs
292 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
294 .. index:: single: Adimensionnement
295 .. index:: single: Sans dimension
297 Il est fréquent que les grandeurs physiques, en entrée ou en sortie des
298 opérateurs, présentent des différences notables d'ordre de grandeur ou de taux
299 de variation. Une manière d'éviter des difficultés numériques est d'utiliser, ou
300 d'établir, un adimensionnement des calculs menés dans les opérateurs
301 [WikipediaND]_. Par principe, dans la mesure où la simulation de la physique
302 devrait être la plus adimensionnée possible, il est en premier lieu recommandé
303 d'utiliser les capacités existantes d'adimensionnement du code de calcul.
305 Néanmoins, dans le cas courant où l'on ne peut en disposer, il est souvent utile
306 d'environner le calcul pour l'adimensionner en entrée ou en sortie. Une manière
307 simple de faire cela en entrée consiste à transformer les paramètres
308 :math:`\mathbf{x}` en argument d'une fonction comme "*DirectOperator*". On
309 utilise le plus souvent les valeurs par défaut :math:`\mathbf{x}^b` (ébauche, ou
310 valeur nominale). Pourvu que chaque composante de :math:`\mathbf{x}^b` soit non
311 nulle, on peut en effet poser:
313 .. math:: \mathbf{x} = \mathbf{\alpha}\mathbf{x}^b
315 et optimiser ensuite le paramètre multiplicatif :math:`\mathbf{\alpha}`. Ce
316 paramètre a pour valeur par défaut (ou pour ébauche) un vecteur de 1. Attention,
317 l'application d'une démarche d'adimensionnement nécessite aussi la modification
318 des covariances d'erreurs associées dans la formulation ADAO du problème
321 Une telle démarche suffit rarement à éviter tous les problèmes numériques, mais
322 permet souvent d'améliorer beaucoup le conditionnement numérique de
