diff --git a/CMakeLists.txt b/CMakeLists.txt index 6c994cd9e0f18c1f7e06666ec24a88d744283d45..dfa4f8de7a565ad4cab40c41d0b2cc015319f6a0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,5 @@ target_include_directories(test_write_eddy PRIVATE add_subdirectory(Inst_eddies) add_subdirectory(Overlap) - add_subdirectory(Common) include(TAGS) diff --git a/Inst_eddies/Documentation_texfol/documentation.tex b/Inst_eddies/Documentation_texfol/documentation.tex index c255ed8a2ec73b44cd0c89e903d1773451d979f1..78929e46ebc8b06ccd759151297bb1c9548ad9f0 100644 --- a/Inst_eddies/Documentation_texfol/documentation.tex +++ b/Inst_eddies/Documentation_texfol/documentation.tex @@ -59,15 +59,39 @@ Chaque ligne de niveau de la hauteur correspond à une ligne de courant de la vitesse. -But du programme : produire une base de données de tourbillons. Le -programme ne traite qu'une date. Il écrit un ensemble de tourbillons -instantanés dans un slice, numérotés à partir de 1, sans saut. Ce slice -va peut-être recevoir d'autres tourbillons lors d'autres -exécutions. Il est plus robuste que le programme marque les +But du programme : produire une base de données de tourbillons +instantanés. Le programme ne traite qu'une date. Il écrit un ensemble +de tourbillons instantanés dans un slice, numérotés à partir de 1, +sans saut. Ce slice va peut-être recevoir d'autres tourbillons lors +d'autres exécutions. Il est plus robuste que le programme marque les tourbillons avec un indice de date. Le programme ne peut pas en général définir cet indice à partir du fichier d'entrée. L'indice de date est donc une entrée supplémentaire du programme. +Cf. algorithme \ref{alg:principal}. +\begin{algorithm}[htbp] + \begin{algorithmic} + \STATE entrer(corner, step, ssh, u, v) + \STATE appel de set\_all\_outerm(s, step, ssh, corner) + \FOR{e dans s\%list} + \IF{e\%valid} + \STATE llc, urc = bbox de e\%out\_cont, dans l'espace des indices + \STATE outside\_points = coordonnées des extremums dans llc, urc à éviter + + \STATE appel de set\_max\_speed(e, outside\_points, llc, urc, ssh, + u, v, step) + + \ELSE + \STATE e\%speed\_cont = null + \STATE e\%max\_speed = manquant + \ENDIF + \ENDFOR + \STATE écrire(s) + \end{algorithmic} + \caption{Algorithme principal} + \label{alg:principal} +\end{algorithm} + \section{Identification of instantaneous eddies} \label{sec:identification} @@ -682,6 +706,7 @@ au contour extérieur, avec un contour vide dans \verb+max_speed_contour_$m+. \section{Structures de données} +\label{sec:structures} On suppose que le maillage est régulier. @@ -832,7 +857,6 @@ Les procédures géométriques utilisées : \item polygon\_area\_2d (bibliothèque Geometry) \item polygon\_contains\_point (bibliothèque Geometry) \item find\_contours\_reg\_grid (bibliothèque Contour\_531) -\item write\_overlap \item mean\_speed \end{itemize} travaillent dans le plan et non sur la sphère. Ces procédures sont @@ -969,15 +993,6 @@ plus indépendants les uns des autres, pour ne pas chercher inutilement des contours de maximum de vitesse. Dans le nouvel algorithme, il n'y a donc plus de procédure get\_eddy. -Si le domaine considéré n'est pas global, comment traiter un -tourbillon au bord du domaine ? La taille du contour le plus extérieur -peut être artificiellement diminuée. (Ce n'est pas le même cas qu'un -tourbillon près d'une côte.) C'est-à -dire qu'on pourrait trouver un -contour plus extérieur avec un domaine élargi. Il pourrait être -intéressant d'éliminer un tel tourbillon. Comment le détecter ? Il -faudrait tester après l'appel de \verb+get_1_outerm+ si le cadre le -plus extérieur touche le bord du domaine. - Notons ici : \begin{align*} & u = \mathtt{urc(1)} \\ @@ -1047,6 +1062,41 @@ Ou encore : k = \left\lceil \frac{w - \lambda}{2 \pi} \right\rceil \end{equation*} +Cf. algorithme \ref{alg:set_all_outerm}. +\begin{algorithm}[htbp] + Déclarations : + \begin{algorithmic} + \STATE résultat s + \STATE données step, ssh, corner + \end{algorithmic} + Instructions: + \begin{algorithmic} + \STATE appel de local\_extrema(s\%extr\_map, innermost\_level, + s\%list\%cyclone, ssh) + + \FOR{e dans s\%list} + + \STATE * calculer e\%coord\_extr en utilisant s\%extr\_map, corner + et step * + + \ENDFOR + \FOR{e dans s\%list par ordre d'orientation et de ssh} + + \STATE llc, urc = bbox dans l'espace des indices : max\_radius + autour de l'extremum + + \STATE outside\_points = coordonnées des extremums dans llc, urc à éviter + + \STATE e\%out\_cont = get\_1\_outerm(e\%cyclone, e\%coord\_extr, + innermost\_level, outside\_points, ssh, llc, urc, corner, step) + + \STATE e\%valid = e\%out\_cont\%area >= min\_area + \ENDFOR + \end{algorithmic} + \caption{subroutine set\_all\_outerm(s, step, ssh, corner)} + \label{alg:set_all_outerm} +\end{algorithm} + \subsection{set\_max\_speed} Notons $d$ la distance dans l'espace des indices entre l'extremum @@ -1064,10 +1114,8 @@ contour venant de \verb+max_speed_contour_ssh+. Il est écrit dans le shapefile extremum comme simple diagnostic. Dans le cas où radius vaut 1, je calcule \verb+max_speed+ et cette -valeur n'est pas utilisée dans \verb+set_max_speed+. Ce calcul ne -serait utile que dans weight si le tourbillon a une intersection avec -un tourbillon à une autre date. Donc, pour certains tourbillons, ce -calcul est complètement inutile. +valeur n'est pas utilisée dans \verb+set_max_speed+. Donc ce calcul +est complètement inutile. \subsection{spher\_polyline\_area} @@ -1400,9 +1448,20 @@ pourquoi. \section{Améliorations, prolongements} +Dans \verb+set_all_outerm+, si le domaine considéré n'est pas global, +comment traiter un tourbillon au bord du domaine ? La taille du +contour le plus extérieur peut être artificiellement diminuée. (Ce +n'est pas le même cas qu'un tourbillon près d'une côte.) C'est-à -dire +qu'on pourrait trouver un contour plus extérieur avec un domaine +élargi. Il pourrait être intéressant d'éliminer un tel +tourbillon. Comment le détecter ? Il faudrait tester après l'appel de +\verb+get_1_outerm+ si le cadre le plus extérieur touche le bord du +domaine. + Idée de chercher les contours extérieurs dans un rayon de 10 points -d'abord et d'agmenter éventuellement le rayon jusqu'à max\_radius si -le contour extérieur trouvé est à un point du bord du domaine. +d'abord et d'augmenter éventuellement le rayon jusqu'à +\verb+max_radius+ si le contour extérieur trouvé est à un point du +bord du domaine. Créer une carte de présence. \`A chaque point du maillage, enregistrer le numéro du tourbillon présent ou 0. Méthode : pour chaque contour @@ -1423,6 +1482,18 @@ maximale non nulle. \subsection{Grille non longitude-latitude} Comment faire si la grille n'est pas une grille longitude-latitude ? + +Dans les données Aviso, les champs SSH, u et v sont sur une même +grille rectangulaire longitude, latitude. Dans les sorties du modèle +utilisé par Jonathan Gula, le champ $\zeta$ et les composantes de la +vitesse dans le plan de projection sont chacun sur une grille +propre. Les trois grilles sont rectangulaires mais non +longitude-latitude. On peut se ramener à deux grilles : une pour +$\zeta$ et une pour $(u, v)$. (Il faut de toutes façons se ramener à +la même grille $\psi$ pour passer des composantes de la vitesse dans +l'espace projeté aux composantes zonales et méridiennes dans l'espace +physique.) + La position en longitude et latitude intervient essentiellement à travers les variables corner, step et \verb+min_area+. Idée de faire tout le traitement dans l'espace des indices mais la latitude et la @@ -1438,29 +1509,36 @@ contour. Comme les points d'un contour sont par construction sur les arêtes du maillage, l'interpolation est en fait linéaire à une seule dimension. -Le passage $(i, j) \in \mathbb{N}^2 \to (\lambda, \phi)$ est ok : +Le passage $(i, j) \in \mathbb{N}^2 \mapsto (\lambda, \phi)$ est ok : c'est juste un accès direct aux tableaux venant du fichier NetCDF. Le -passage $(i, j) \in \mathbb{R}^2 \to (\lambda, \phi)$ est ok aussi : +passage $(i, j) \in \mathbb{R}^2 \mapsto (\lambda, \phi)$ est ok aussi : c'est une interpolation bilinéaire. Le problème est le passage -$(\lambda, \phi) \to (i, j) \in \mathbb{R}^2$. +$(\lambda, \phi) \mapsto (i, j) \in \mathbb{R}^2$. Idée : faire passer \verb+snapshot%ind_extr+ dans eddy et l'écrire -dans le shapefile extr dbf (pour \verb+eddy_graph+). Ajouter dans eddy -une composante polyline \verb+out_cont_no_coord+. Dans -\verb+set_max_speed+, calculer \verb+mean_speed+ sur (radius4 - 1) -$\times$ 4 contours. Mettre \verb+outside_points+ vide dans -\verb+good_contour+. Faire la recherche avec \verb+good_contour+ dans -l'espace des indices. +dans le shapefile extr dbf (pour \verb+eddy_graph+). Mais inconvénient +: cf. § \ref{sec:structures}. Donc garder \verb+snapshot%ind_extr+ et +l'écrire dans un fichier indépendant (binaire Fortran à accès direct ? +HDF5 ?) ou écrire un fichier NetCDF contenant \verb+extr_map+ ? +Ajouter dans eddy une composante polyline +\verb+out_cont_no_coord+. Dans \verb+set_max_speed+, calculer +\verb+mean_speed+ sur (radius4 - 1) $\times$ 4 contours. Mettre +\verb+outside_points+ vide dans \verb+good_contour+. Faire la +recherche avec \verb+good_contour+ dans l'espace des indices. \begin{algorithmic} - \STATE lire ($\lambda$, $\phi$, u, v, ssh, angle) - \STATE définition de extr map, ind extr, innermost level, cyclone + \STATE lire ($\lambda$, $\phi$, ssh, $\lambda_{uv}$, $\phi_{uv}$, u, v) + + \COMMENT{set\_all\_outerm} + \STATE définition de extr\_map, ind\_extr, + innermost\_level, cyclone + \FORALL{extremum} - \STATE calculer coord extr (pas d'interpolation) + \STATE calculer coord\_extr (pas d'interpolation) \ENDFOR \FORALL{extremum} - \STATE définir outside points dans l'espace des indices - \STATE get 1 outerm + \STATE définir outside\_points dans l'espace des indices + \STATE get\_1\_outerm \ENDFOR \end{algorithmic} diff --git a/Inst_eddies/local_extrema.f90 b/Inst_eddies/local_extrema.f90 index 78e98a55e0055e127efcc0c9c70cd57e09226627..2ae835f132e4ed1746d32dd6f2633a7ec4c1ace4 100644 --- a/Inst_eddies/local_extrema.f90 +++ b/Inst_eddies/local_extrema.f90 @@ -9,12 +9,12 @@ contains subroutine local_extrema(extr_map, ind_extr, innermost_level, local_min, & field, periodic, my_lbound) - ! Finds local extrema in a field and closest value around each - ! extremum. The extrema are output both in a map and as a - ! list. The input field may be periodic in the first - ! dimension. Then it is assumed that it has one duplicate column - ! at each bound. If periodic then some extrema may be duplicated - ! in extr_map but there is no duplication in ind_extr, + ! This procedure finds the local extrema in a field and the + ! closest value around each extremum. The extrema are output both + ! in a map and as a list. The input field may be periodic in the + ! first dimension. Then it is assumed that it has one duplicate + ! column at each bound. If periodic then some extrema may be + ! duplicated in extr_map but there is no duplication in ind_extr, ! innermost_level and local_min. ! Note that no coordinate grid is used here so there is no diff --git a/Inst_eddies/nearby_extr.f90 b/Inst_eddies/nearby_extr.f90 index 3558e8083ae547246e5bbc409006dfd95a434110..80ba49aca882900983093f401436d8fe57c33162 100644 --- a/Inst_eddies/nearby_extr.f90 +++ b/Inst_eddies/nearby_extr.f90 @@ -6,8 +6,8 @@ contains pure function nearby_extr(extr_map, list, i) - ! Returns a list of extrema that cannot be engulfed in a good - ! contour around the target extremum. + ! This procedure returns a list of extrema that cannot be engulfed + ! in a good contour around the target extremum. use derived_types, only: eddy diff --git a/Inst_eddies/set_all_outerm.f90 b/Inst_eddies/set_all_outerm.f90 index 40b95e8de1caba35409be417ba7f0cac78e8cbd9..85f99ab92f8326a69d1c416eb7ad9c4aa10edfe9 100644 --- a/Inst_eddies/set_all_outerm.f90 +++ b/Inst_eddies/set_all_outerm.f90 @@ -122,6 +122,8 @@ contains selection = selection(indexx(s%list(selection)%ssh_extr)) sorted_extr(n_cycl + 1:) = selection + ! Done sorting + do l = 1, s%number_extr i = sorted_extr(l) diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90 index 00ac4703bc257c7b8942479e17a2025d06d3db71..c4d3259be44ce9f3c5268806171eceabd4df2d20 100644 --- a/Inst_eddies/set_max_speed.f90 +++ b/Inst_eddies/set_max_speed.f90 @@ -7,9 +7,9 @@ contains subroutine set_max_speed(e, ind_targ_extr, outside_points, ssh, u, v, & corner, step) - ! Defines the components speed_cont, max_speed and radius4 of - ! argument e. On return, e%speed_cont may be a null ssh contour - ! but e%max_speed is never set to missing_speed. However, + ! This procedure defines the components speed_cont, max_speed and + ! radius4 of argument e. On return, e%speed_cont may be a null ssh + ! contour but e%max_speed is never set to missing_speed. However, ! e%max_speed might be NaN if speed_outerm is NaN (with any value ! of e%radius4). The value of e%radius4 set here is always >= ! 1. On return, if e%speed_cont is not a null ssh contour then diff --git a/Overlap/Documentation_texfol/Graphiques/input_output.odg b/Overlap/Documentation_texfol/Graphiques/input_output.odg index 3e15a961bdfd46c7ec09cb9451dbdb148e085ff7..e1124f93a1a751c481837cf3b4cc34e729fa9b27 100644 Binary files a/Overlap/Documentation_texfol/Graphiques/input_output.odg and b/Overlap/Documentation_texfol/Graphiques/input_output.odg differ diff --git a/Overlap/Documentation_texfol/documentation.tex b/Overlap/Documentation_texfol/documentation.tex index 0911711e99a99e212ef28ebff288606152c445a2..5f1f9cebbb5410c17d107ca258aa865a05cfcfa8 100644 --- a/Overlap/Documentation_texfol/documentation.tex +++ b/Overlap/Documentation_texfol/documentation.tex @@ -70,7 +70,7 @@ numérotation par dates successives. Le graphe non orienté associé n'est pas une forêt. Cf. figure \ref{fig:cycle}. \begin{figure}[htbp] \centering - \includegraphics{cycle} + \includegraphics[scale=0.5]{cycle} \caption[Exemple possible de cycle dans le graphe non orienté associé]{Exemple possible de cycle dans le graphe non orienté associé. (2, 3, 4, 6, 2) est un cycle donc le graphe n'est pas un @@ -152,7 +152,7 @@ prédécesseurs visibles et s'il a des successeurs visibles, mais aussi $m$ est le numéro du processus MPI, compris entre 0 et $n_p - 1$. On pourrait écrire l'algorithme principal en explicitant les -\verb+get_snpashot+, en explicitant le cas $m < n_p - 1$ et en +\verb+get_snapshot+, en explicitant le cas $m < n_p - 1$ et en supprimant la variable \verb+k_end_main_loop+. Ce serait mieux du point de vue de la performance, mais je pense que ce serait finalement moins clair. @@ -235,7 +235,7 @@ Le programme lit une liste de SHPC. Nous supposons que la première date présente dans un SHPC est strictement supérieure à la dernière date du SHPC précédent. -\section{Déclarations pour l'algorithme principal} +\section{Algorithme principal, séquentiel} On suppose que les extremums sont sur les n\oe{}uds d'un maillage rectangulaire uniforme. On suppose que ce maillage est le même à @@ -766,13 +766,14 @@ croissant des $\delta$. Plusieurs dates dans un triplet de shapefiles. Mais toutes les dates ne tiennent pas a priori dans un seul triplet de shapefiles (à cause -de la limite sur la taille des shapefiles). Ce serait plutôt à -\verb+get_snapshot+ de gérer le choix du bon triplet. Mais il ne faudrait -pas ouvrir et fermer plusieurs fois le même triplet. Pour un processus -donné, \verb+get_snapshot+ est appelé avec des dates -croissantes. \verb+get_snapshot+ n'a donc besoin d'avoir à un instant donné -qu'un seul triplet ouvert. Je fais le pari que plusieurs processus MPI -pourront accéder en même temps en lecture à un même fichier. +de la limite sur la taille des shapefiles). \verb+read_snapshot+ gère +le choix du bon triplet. Il ne faut pas ouvrir et fermer plusieurs +fois le même triplet. Pour un processus donné, \verb+get_snapshot+ est +appelé avec des dates croissantes. \verb+get_snapshot+ n'a donc besoin +d'avoir à un instant donné qu'un seul triplet ouvert. Mais chaque +processus commence par ouvrir toutes les tranches. Je fais le pari que +plusieurs processus MPI pourront accéder en même temps en lecture à un +même fichier. \subsection{overlap} diff --git a/Trajectories/cost_function.py b/Trajectories/cost_function.py index 52bc7f3f94b6939abd2adc0715c6c73fab7c5e07..a8d1dbbc5570c37c292428d5de5cd51728c7db2a 100755 --- a/Trajectories/cost_function.py +++ b/Trajectories/cost_function.py @@ -187,8 +187,8 @@ for edge in g.edges(): # (latitude needed for conversion of degrees to kilometers) lat_for_conv = math.radians(lat_for_conv) # need to convert to radians - # because of the wrapping issue (360° wrapping incorrectly to 0°), - # we check for that here + # Because of the wrapping issue (360° wrapping incorrectly to 0°), + # we check for that here: lon_diff = abs(g.vp.pos_last[source_node][0] - g.vp.pos_first[target_node][0]) if lon_diff > 300: lon_diff = 360 - lon_diff