I.9 : Petit r�sum� :

        L'�tude qui fait l'objet de ce m�moire a plusieurs caract�ristiques et contraintes :

        Ces crit�res expliquent pourquoi ce projet ne ressemble pas � des projets classiques : l'unique param�tre int�ressant �tant la vitesse, de nombreux probl�mes �tant pr�cedemment r�solus de mani�re th�oriques ou conceptuels, le d�fi � relever est la gestion de la complexit� croissante du projet.

II.1 : Introduction :

        Cette partie ne peut pas et n'a pas la pr�tention d'�tre exhaustive : le sujet est trop vaste pour pouvoir envisager d'en faire le tour complet. Il n'est pas non plus possible ici d'expliquer (encore et encore) ce que sont les LGA : d'autres articles ou livres le font tr�s bien. Le lecteur novice et press� peut consulter l'annexe B, les explications les plus claires en fran�ais sont donn�es en [18], [14] et [27]. Ce m�moire discute de la programmation d'un mod�le physique, nous allons donc nous concentrer sur les notions essentielles qui permettent de comprendre les choix de conception.

II.2 : Nomenclature :

        Le terme de "gaz sur r�seau", ou "lattice gas automata" en anglais (LGA) n'a pas fait l'objet d'une �tude par l'Acad�mie Fran�aise et la premi�re question est de savoir s'il faut mettre "r�seau" au pluriel. La r�ponse la plus satisfaisante consiste � consid�rer le r�seau comme un m�dium, un "ether" virtuel et digital sur lequel prennent vie des ph�nom�nes, gr�ce au support physique de l'ordinateur. Un "gaz sur r�seau" d�signe donc un gaz particulier (par exemple "FHP-3") alors que l'expression "les gaz sur r�seaux" d�signe un ensemble de gaz aux propri�t�s diff�rentes (par exemple l'ensemble des LGA en deux dimensions).

        On peut aussi trouver dans la litt�rature diverses d�nominations pour d�signer un ou des gaz sur r�seaux : "LGA", "LG" pour "Lattice Gas", "LGCA" pour "Lattice Gas Cellular Automata", GR ou GSR pour "gaz sur r�seau(x)"... Les nombreuses variantes des mod�les apportent encore d'autres noms ou acronymes.

II.3 : G�n�se :

        L'existence des LGA est li�e � plusieurs courants d'id�es ou techniques. D'abord, les LGA sont une simplification extr�me des lois et des programmes de "Dynamique Mol�culaire" o� chaque particule (atomes ou mol�cules) d'une mati�re est simul�e avec sa vitesse, sa masse, sa direction, avec autant de valeurs en virgule flottante pour chaque grandeur et dans chaque dimension. Ainsi, contrairement aux m�thodes d'�l�ments finis, les particules sont simul�es afin de laisser ressortir leur comportement individuel, au lieu de les enfermer dans des moyennes et des interpolations. Les Gaz sur R�seaux se situent entre le monde de la dynamique mol�culaire ("MD") d'une part et les techniques volumiques classiques (Euler ou Navier-Stokes par exemple) d'autre part, car ils apportent la capacit� de simuler de plus grands ph�nom�nes qu'avec de la MD classique. Ils peuvent faire �merger des comportements de plus grande �chelle, accessibles seulement par les techniques d'�l�ments finis mais avec un niveau de d�tail sup�rieur.

        L'autre facteur d�terminant pour l'existence des LGA est la d�mocratisation des ordinateurs. Il suffit en fait de peu de ressources pour commencer � exp�rimenter sur ce type de "m�dium". C'est cette d�pendance technologique qui explique en partie l'historique de la m�thode : tout comme pour les Automates Cellulaires classiques, les premi�res exp�riences ont �t� r�alis�es � la fin des ann�es 60 (Kadanoff & Swift, 1968) et le premier mod�le intens�ment �tudi� (HPP) date officiellement de 1973. Le domaine a ensuite explos� au d�but des ann�es 80 gr�ce � des ressources informatiques r�pandues et l'int�r�t du public. Les LGA ont des liens de parent� tr�s forts avec d'autres mod�les comme les Automates Cellulaires classiques ou le mod�le d'Ising, utilis�s en thermodynamique ou pour �tudier les s�parations de phases sans hydrodynamique. Le mod�le HPP repr�sente une convergence naturelle de toutes ces influences et a cr�� un nouveau domaine hybride et fascinant.

II.4 : Un premier mod�le : HPP

        Les LGA sont issus de recherches de m�canique statistique, un des objectifs �tant de simplifier les lourds calculs de dynamique mol�culaire pour en extraire les composantes fondamentales. Des mod�les � v�locit� finie, � temps fini puis des simplifications de plus en plus radicales (binarisation et discr�tisation totale) ont donn� le jour en 1973 au mod�le dit "HPP" [37] (initiales de Jean Hardy, Olivier de Pazzy et Yves Pomeau). Ce mod�le n'est pas utilisable en pratique mais sa compr�hension est importante pour des raisons historiques et techniques car il n'est pas possible de faire plus simple (en deux dimensions). Il cristalise donc tous les d�fauts et toutes les caract�ristiques que l'on retrouve dans les mod�les plus �volu�s et son �tude permet de g�n�raliser des techniques aux autres mod�les. Par exemple, les �quations de collision ou la complexit� algorithmique (comme dans [31]) sont �tudi�es d'abord sur HPP avant d'�tre adapt�es aux autres mod�les. HPP ressemble � un Automate Cellulaire � voisinage de Von Neumann (carr�) simple, synchrone et homog�ne. Il diff�re d'un Automate Cellulaire classique car chaque cellule ne poss�de pas d'�tat interne : comme not� pr�c�demment, le gaz sur r�seau est un "m�dium" sur lequel transitent des particules bool�ennes qui s'entrechoquent aux intersections du grillage. Les cellules sont ici appel�es "noeuds" et sont le si�ge des collision. Comme le tapis d'un billard sur lequel roulent les boules, les particules se d�placent sans friction sur cet "ether" informatique.

        Les frictions que l'on veut faire apparaitre sont du type hydrodynamique (par exemple la viscosit�) et concernent les interactions entre particules. Les interactions avec le m�dium servent � d�terminer les "conditions aux limites" (les parois ou les objets qui agissent sur le fluide). Les collisions entre les particules d�terminent le comportement du fluide et avec HPP il n'y a pas beaucoup de choix : les r�gles de conservation sont strictes et il n'y a pas beaucoup de degr�s de libert�. Ainsi :

• Toute collision doit conserver le nombre de particules : il y a autant de particules qui entrent dans le noeud (aff�rentes) que de particules sortantes. Le contraire ne serait pas logique.
• Toute collision doit conserver l'impulsion g�n�rale des particules : en "clair", la somme vectorielle des vecteurs mouvements de toutes les particules entrantes doit �tre identique � la somme vectorielle des vecteurs mouvements de toutes les particules sortantes. Nous aurons l'occasion de revenir plus en d�tail sur ce sujet bient�t.
• Il n'y a que 2^4 (16) configurations possibles pour chaque site. Les seuls quatre voisins limitent le nombre de combinaisons int�ressantes.

        Parmi toutes les combinaisons possibles, le choc frontal est la seule qui permette de r�organiser les particules aff�rentes tout en respectant les lois de conservation �nonc�es pr�c�demment.

Description d'une collision frontale.

Description des vecteurs mouvements d'une collision frontale.

        L'impulsion est la somme des vecteurs mouvements des particules. Le cas de la collision frontale avec HPP se r�sume � additionner deux vecteurs (unitaires) de direction oppos�e, ce qui donne un vecteur nul. L'impulsion d'autres configuration donnera d'autres vecteurs mais le choc frontal est la seule configuration o� toutes les grandeurs �nonc�es sont conserv�es et o� l'on puisse avoir une autre configuration en sortie. Ce type d'�change de configuration est la base de l'�mergence de ph�nom�nes d'hydrodynamique � l'�chelle macroscopique et sera raffin� plus tard.

II.5 : Caract�ristiques du mod�le HPP :

        Le choc frontal repr�sente 12,5% du champ de collision (2 combinaisons sur 16) et se produit � une densit� de 0.5 (2 particules sur 4). Par rapport aux mod�les r�cents, le premier rapport montre une faible efficacit� (faible rapport de nombre de Reynolds par site) mais HPP est inutilis� surtout parcequ'il a quatre "invariants" qui sont nuisibles dans la plupart des cas. Les "invariants" sont des quantit�s conserv�es dans le mod�le et qui ne correspondent pas au comportement d'un fluide r�el.

• temp�rature : HPP est un mod�le o� toutes les particules ont la m�me vitesse. Or la vitesse des particules (dans l'air par exemple) est fonction de leur temp�rature (et vice versa). Avec HPP la masse de toutes les particules est identiques, donc leur �nergie cin�tique ne change pas. Il n'est pas possible de simuler des transferts de chaleur ou tout ph�nom�ne o� la temp�rature change : tout le fluide simul� est � une temp�rature uniforme, homog�ne. Cela simplifie l'�tude des �quations caract�ristiques mais interdit d'explorer les pans les plus passionnants de la m�canique des fluides et de la thermodynamique.
• invariance de parit� lin�aire  (d�sol� pour le n�ologisme) : un ph�nom�ne curieux de non-r�partition homog�ne est implicitement port� par HPP. Il n'est pas possible de diffuser uniform�ment une masse de particules dans un fluide, bien que la diffusion ne fasse pas apparaitre d'onde carr�e prononc�e. Cela est d� au fait que la collision frontale, seule permise par le mod�le, ne traite qu'un nombre pair de particules par lignes. Pour poser le probl�me simplement, on peut consid�rer que la parit� du nombre de particules sur une ligne ne change pas. L'effet ind�sirable n'est sensible qu'avec des g�om�tries tr�s petites mais reste g�nant car il est conserv� dans les �quations � grande �chelle.
• anisotropisme : tout ph�nom�ne ne se produira pas de la m�me mani�re selon l'orientation par rapport � la grille du gaz. Plus simplement, le fluide a un "axe de pr�f�rence", l'orientation du m�dium agit indirectement sur le fluide. C'est un d�faut majeur qui emp�che d'�tudier correctement des tourbillons par exemple.
• invariance galil�enne : l'aspect monocin�tique du fluide emp�che les tourbillons de se d�placer � la m�me vitesse que le fluide. Si un tourbillon apparait dans le fluide, il sera advect� (emport�) plus vite que celui-ci (selon les cas).

II.6 : Le mod�le FHP :

        Les Automates Cellulaires peuvent-ils r�soudre des �quations diff�rentielles partielles, comme celles de Navier-Stokes ?
Telle �tait la question cruciale pos�e par de nombreux scientifiques depuis l'apparition de ce domaine d'�tude. Rappelons que les premiers calculateurs automatiques ont �t� con�us dans cette optique, comme la machine � diff�rences finies de Babbage ou le calculateur �lectronique de John Vincent Atanasoff (1937-1942). Stanislaw Ulam et Konrad Zuse dans les ann�es 50, Stepen Wolfram et Richard Feynman dans les ann�es 80, ont milit� pour r�soudre cet �pineux probl�me. Pourtant, ce n'est que dix ans apr�s l'introduction du mod�le HPP que la connexion entre les deux sujets a �t� comprise. Uriel Frisch, Brosl Hasslacher et Yves Pomeau (d'o� FHP) ont propos� en 1986 une l�g�re modification qui permet de retrouver les �quations diff�rentielles de Navier-Stokes [19].

        La premi�re modification du mod�le HPP porte sur le "medium" anisotropique. Au lieu d'un maillage carr�, un maillage hexagonal est n�cessaire et suffisant pour r�soudre ce probl�me. Cette l�g�re modification permet en fait d'augmenter le nombre de degr�s de libert� lors des collisions et d'avoir plus de possibilit�s de sortie diff�rentes. Le fluide a donc plus d'opportunit�s de diffuser ses particules dans chaque direction et de diminuer sa viscosit�. Les simulations deviennent soudain plus int�ressantes ...

        La deuxi�me am�lioration du mod�le HPP r�soud "l'invariance de parit� lin�aire" en exploitant ces nouveaux degr�s de libert�s. Le nombre de collisions �quivalentes augmente et il suffit d'ajouter une nouvelle loi de collision : la "collision triangulaire".

        L'invariance galil�enne sera pour l'instant corrig�e par un adimensionnement d'une grandeur. Si l'advection d'un vortex a lieu � quatre fois la vitesse du fluide, on divisera par quatre les vitesses mesur�es pour retrouver l'invariance galil�enne (en monophasique). Cette m�thode n'est plus valable avec diff�rentes phases mais cela sort du sujet du m�moire.

        La temp�rature n'est pas un probl�me en soi, elle n'intervient pas dans les exp�riences qui int�ressent les utilisateurs de LGA. Si elle est n�cessaire, la solution est simple : avoir plusieurs vitesses de particules. Deux moyens existent : les particules elles-m�mes sautent plusieurs sites � chaque pas de temps, ou/et les liens ont plusieurs longueurs. Les r�gles de collision doivent refl�ter ces changements tout en respectant les r�gles de conservation initiales. En deux dimensions, un tel r�seau est g�n�ralement le voisinage de Moore � 8 voisins : les diagonales sont 1.41 fois plus longues que les liens horizontaux ou verticaux. Le r�seau hexagonal reste toutefois plus int�ressant pour son plus grand nombre d'isom�tries et donc par un plus grand nombre d'opportunit�s de collisions �quivalentes (D2Q12 ?).

II.7 : Les r�gles et les nouvelles propri�t�s de FHP-1 :

        Dans leur fameux article [19], Frisch, Hasslacher et Pomeau vont d�crire une premi�re version du mod�le FHP qui r�pond � la contrainte d'�tre n�cessaire et suffisante pour retrouver les �quations de Navier-Stokes � grande �chelle. Seule la partie th�orique est r�solue dans l'�tude, l'efficacit� n'est pas au rendez-vous : les am�liorations viendront plus tard, la d�monstration �tant d�j� un grand pas.

Les r�gles de collision du mod�le FHP-1.

        Le deuxi�me effet du changement de g�om�trie, apr�s avoir bris� l'anisotropie, fut d'apporter d'autres axes de sym�trie : 3 au lieu de 2. Il y en a maintenant 64 possibilit�s de configuration en entr�e, au lieu de 16. Les auteurs vont introduire 2 types de collisions au lieu d'une seule, soit 5 collisions au lieu de 2 pour HPP, si on tient compte des sym�tries et des rotations. On passe de 12,5 % � 12,8 % d'exploitation du champ des collisions mais l'am�lioration se situe autre part.

        D'abord, la collision triangulaire am�liore sensiblement la qualit� de l'�coulement et permet � elle seule de faire apparaitre des comportements hydrodynamiques r�alistes. Elle permet de faire "communiquer" entre elles des lignes et de mieux r�partir les particules sur les diff�rents axes du r�seaux.

        Ensuite, et c'est tout aussi remarquable : les collisions frontales n'ont plus seulement une, mais deux "voies de sortie". Il existe trois configurations d'entr�e �quivalentes, �quiprobables, au lieu de deux pour HPP. Il faut donc choisir, lorsqu'un choc frontal a lieu, la configuration de sortie.

• La premi�re solution (na�ve) est de fixer statiquement un ordre : par exemple effectuer une rotation de 60 degr�s pour chaque configuration.
   
  => Cette solution introduit un nouvel effet parasite : la chiralit�, ou pr�f�rence d'un sens de rotation, qui est n�faste lors de l'�mergence de ph�nom�nes tr�s turbulents. Il faut �viter de favoriser une direction autant que possible afin de garder le fluide "pur". La chiralit� peut �tre un moyen de simuler l'effet de Coriolis mais ce n'est pas coh�rent avec notre �chelle macroscopique.
   
  
• Deuxi�me solution, pour les programmeurs fatigu�s ou press�s (comme Haarlan Stockman) : changer de chiralit� une ligne sur deux. La parit� du num�ro de la ligne indique si le sens de rotation est de +60 ou -60 degr�s.
   
  => Des �tudes [d'Humi�res ?] ont montr� que les effets de la chiralit� ne sont pas totalement effac�s, principalement au niveau des parois. Ce mauvais compromis n'est pas recommand� en pratique, m�me si les effets ne sont pas directement visibles.
   
  
• La solution recommand�e est d'effectuer la rotation une fois sur deux, au hasard. Le g�n�rateur de nombres al�atoires n'a pas besoin d'�tre de qualit� indiscutable, mais suffisante pour briser la chiralit�. Il doit donc simplement �tre �quiprobable et avoir une longue p�riode de r�p�tition. Dans notre programme d'exp�rimentation, le nombre al�atoire est tout simplement mis � jour en fonction des sites pr�c�dents (un simple ADD avec les sites de direction A) ce qui fournit des donn�es de nature suffisamment ind�pendantes pour briser la chiralit�, apr�s une p�riode d'amor�age n�gligeable.

        Pour illustrer cette propri�t� incroyable, on proc�de parfois � l'exp�rience suivante : soit un programme FHP avec une configuration des particules initiales choisies arbitrairement. Le g�n�rateur de nombres al�atoires est d�terministe et r�versible.

1) Le programme est lanc� et calcule T pas de temps, autant que n�cessaires pour obtenir une relaxation suffisante du fluide (une soupe brownienne bien homog�ne).

2) Les direction des particules sont toutes invers�es et le programme est relanc� avec le g�n�rateur de nombres al�atoires tournant dans le sens inverse.

3) Apr�s un m�me nombre T de pas de temps, le fluide retrouve sa configuration initiale, telle qu'au d�but du 1).

4) Ensuite, on restaure l'�tat du 2) et on change un bit, une direction de particule ou la "graine" du g�n�rateur de nombres al�atoires. Quel que soit le nombre de pas de temps de calcul, la configuration initiale ne sera pas retrouv�e et elle restera certainement � l'�tat de "soupe brownienne".

        Le fluide FHP a des propri�t�s fortement non lin�aires et dissipatives, ainsi que de nombreuses autres propri�t�s curieuses, malgr� sa simplicit� �tonnante. Il ob�it � la loi de Mariotte et des gaz parfaits. Il est un milieu compressible et peut donc propager une onde "sonore" de mani�re circulaire comme dans la nature. Son calcul est par construction inconditionnellement stable et exact. Mais son efficacit� est mauvaise en pratique : le chemin parcouru par une particule entre deux collisions (mean free path en anglais) est tr�s long en moyenne, la viscosit� du fluide est forte et il faut donc des millions de sites pour simuler des ph�nom�nes int�ressants.

II.8 : Les am�liorations de FHP-2 :

        La premi�re am�lioration du mod�le FHP original, ensuite renomm� FHP-1, fut d'ajouter une "particule immobile". Cet ajout permet de porter le nombre de configurations � 128 et de rajouter de nombreuses opportunit�s inesp�r�es pour r�arranger les particules en sorties. La viscosit� chute et le mod�le devient plus efficace. On commence aussi � s'int�resser aux propri�t�s de dualit� du mod�le : dans le cas du choc frontal, remplacer une particule par un "trou" et vice versa est tout � fait valable.

Les r�gles de collision additionnelles du mod�le FHP-2.

        La viscosit� diminue et le nombre de collisions augmente : 20 entr�es mais l'occupation du champ de collision reste faible : 15 %. Bien que les possibilit�s augmentent et que la d�finition du mod�le FHP-2 ne soit pas tr�s pr�cise, il reste encore des efforts � faire.

II.9 : Caract�ristiques particuli�res du mod�le FHP-3 :

        Le mod�le sur lequel nous allons nous attarder est le mod�le FHP-3. Il est une extension du mod�le FHP-2 avec une jeu de collisions satur� : 76 combinaisons sur 128 donnent lieu � un r�arrangement des particules � la sortie. 59 % : c'est le maximum d'occupation du champ de collisions que l'on puisse obtenir avec 7 bits. La table des collisions sera analys�e en d�tail dans la partie IV alors attardons-nous ici sur les raisons d'�tudier ce mod�le particulier.

        Tout d'abord, ce mod�le offre un compromis acceptable de simplicit� et d'efficacit�. Il n'est pas plus avantageux de faire plus simple : les deux mod�les FHP d�crits pr�c�demment sont conceptuellement plus simples mais la surcharge en calcul d�e aux collisions plus complexes de FHP-3 est largement compens�e par la r�duction du nombre de donn�es � traiter, donc le temps de calcul est r�duit. Nous verrons aussi que les mod�les FHP sont memory bound sur la plateforme qui nous concerne, une acc�l�ration des calculs n'a donc pas d'influence radicale sur le temps total de calcul car la m�moire est trop lente. La programmation de FHP-3 est plus avantageuse que FHP-1 ou FHP-2 malgr� une plus grande complexit�.

        Il n'est pas non plus envisageable � notre niveau de programmer un mod�le plus efficace : apr�s FHP3, aucun mod�le ne fait autorit� et n'est suffisamment connu pour �tre utilis� dans le cadre de notre �tude de cas. Apr�s FHP-3, le foisonnement d'am�liorations a dispers� les efforts et aucun nouveau mod�le discret n'est assez maitris�e par le public novice. De plus, la plus grande partie des nouveaux mod�les implique une r�organisation compl�te des donn�es et du programme : nouveaux r�seaux, nouvelles g�om�tries, nouvelles lois et nouveaux artefacts � maitriser. Enfin, notre travail doit rester une �tude de cas simple et notre but n'est pas d'�tudier un nouveau mod�le en profondeur : nous voulons simplement en faire fonctionner un le plus vite possible.

        FHP-3 est donc l'un des derniers mod�les "stables" et connus avant la diversification des mod�les. On peut ainsi esp�rer que des techniques d�velopp�es pour FHP-3 sont r�utilisables facilement dans les autres mod�les. Ces techniques (notamment le strip mining) et leurs implications seront ainsi facilement comprises par les utilisateurs qui les r�utiliseront et les am�lioreront selon les cas.

II.10 : Propri�t�s physiques des mod�les FHP :

        Pour illustrer les propos du chapitre pr�c�dent, nous �tudierons le tabelau suivant :

mod�le	FHP-1	FHP-2	FHP-3
Cs (vitesse du son en site par cycle)
g
f* (densit� id�ale, en occupation du site)	0,187	0,179	0,285
Re* (Nombre maximal de Reynolds par site � la densit� id�ale)	0,387	1,08	2,22

        Pour fonctionner correctement, � un nombre de Reynolds maximum, les Gaz sur R�seaux doivent �tre utilis�s autour d'une densit� particuli�re qui d�pend du mod�le. Pour FHP-3, cette densit� est heureusement simple : deux particules par site (2/7 = 0,285). Si l'invariance galil�enne est n�cessaire, la densit� peut �tre diff�rente. La vitesse des �coulements doit aussi �tre limit�e, en pratique � Mach 0.3 soit au maximum 0,6*0,3 = 0,2 sites par pas de temps.

        Puisque nous parlons ici du nombre de Reynolds, attardons-nous dans ce paragraphe sur sa d�finition. Pour simplifier un peu arbitrairement, il met en relation la taille d'un objet ou d'un ph�nom�ne avec la viscosit� du fluide. Pour un gaz sur r�seau, la viscosit� correspond � une sorte de r�sistance � la diffusion d'une force ou d'une perturbation dans toutes les directions. Ainsi, plus le nombre de Reynolds est �lev�, plus les ph�nom�nes se diffusent et les g�om�tries sont complexes. Il est ainsi possible de comparer des ph�nom�nes qui ont lieu � des vitesses diff�rentes, � des �chelles diff�rentes, � des pressions diff�rentes dans des mati�res diff�rentes, la complexit� du ph�nom�ne est ainsi r�duite � un seul nombre sans unit�. Il correspond aux �quivalences suivantes, de la plus th�orique � la plus pratique pour nous :

avec

        Le tableau montre que le gain en efficacit� en nombre de Reynolds entre FHP-1 et FHP-3 est environ d'un facteur 6. Il faut donc en th�orie 36 fois moins de sites et 6^3=216 fois moins de temps de calcul avec FHP-3 (puisque le temps de calcul cro�t environ � la puissance troisi�me de la longueur caract�ristique). Les param�tres id�aux d'utilisation changent aussi, tout comme les caract�ristiques (par exemple la vitesse du son). Le respect de ces param�tres est crucial pour utiliser FHP de mani�re optimale et sans surprise : il suffit d'un �cart pour observer des ph�nom�nes non physiques (artifacts).

        De plus, on ne reviendra jamais assez sur l'importance de l'adimensionnement li� � l'invariance galil�enne ! A densit� id�ale, la vitesse du fluide doit �tre multiplit�e par 3/10 (voir la formule du tableau avec une densit� de 2/7).

        Pour les Gaz sur R�seaux du type FHP, d'autres ph�nom�nes curieux, inattendus et inqui�tants apparaissent. Pour des �coulements � pression constante, les r�gions � plus grande vitesse ont une plus grande densit� :

avec

ce qui ne correspond pas � la formule normale :

        L'�quation de Navier-Stokes pour le cas g�n�ral est celle-ci :

        Mais pour un Gaz sur R�seau (hexagonal) � basse pression et � basse vitesse, elle devient :

        Le terme g(rho) est ajout� pour corriger l'invariance galil�enne. Il est donn� dans le tableau pr�c�dent, o� il d�pend de la densit�. Selon les cas, il peut �tre int�gr� dans l'adimensionnement d'une grandeur (vitesse, pression, temps...). Il faut alors choisir la densit� en fonction du crit�re d'exp�rimentation, soit un grand nombre de Reynolds ou l'invariance galil�enne, restaur�e d'une mani�re ou d'une autre. La pression, la densit�, la vitesse, la viscosit� et le temps sont ainsi enchev�tr�s de mani�re complexe et leur �tude sort du cadre de ce m�moire puisqu'elle est disponible dans la litt�rature sp�cialis�e. L'utilisateur est mis en garde du fait que la validit� des r�sultats d�pend de la qualit� de l'analyse du probl�me et de la compr�hension intime des m�canismes mis en jeu. Le programme seul ne pourra donner de r�sultat correct que dans des conditions contr�l�es par l'utilisateur.

        Les gaz sur r�seaux FHP ne sont pas parfaits, nous le savons maintenant. Mais il y a surtout un d�tail � rendre insomniaque car aucune solution simple n'est satisfaisante. Le probl�me peut �tre �nonc� de cette mani�re : si deux particules voyagent sur le m�me lien en sens inverse, elles ont moins d'une chance sur deux pour �tre d�vi�es dans un choc frontal. En cons�qunce, la viscosit� de FHP-3 est inf�rieure � ce que le mod�le pourrait fournir dans un cas id�al o� tout �v�nement d'une particule croisant une autre donnerait lieu � une r�organisation des configurations. Le probl�me touche la nature m�me du r�seau car un lien entre deux noeuds a deux canaux unidirectionnels et ind�pendants. Il n'est pas possible de diviser ce canal en deux tron�ons et de tester les collisions � cet endroit car il faudrait disposer d'une autre dimension ou direction pour r�ordonner les configurations � chaque temps t+1/2. Le probl�me a �t� abandonn�, les caract�ristiques des r�seaux existants �tant estim�s suffisamment satisfaisantes. Mais surtout, il faut respecter le principe d'exclusion de Fermi qui stipule que deux particules ne peuvent pas avoir la m�me position et la m�me vitesse au m�me instant.

        Le dernier d�faut reproch� aux LGA bool�ens est la quantisation et le fort niveau de bruit. Il faut effectuer des moyennes, ce qui emp�che quasiment toute mesure ponctuelle. Un physicien pr�f�re souvent 10% d'erreur � 10% de bruit lorsqu'il veut effectuer un calcul. Notons toutefois pour terminer que la pr�cision d'un gaz sur r�seau approche 1% dans des conditions id�ales, soit environ la pr�cision de la mesure. Les m�thodes traditionnelles sont qualitativement et quantitativement plus impr�cises en pratique et n�cessitent de nombreuses mises � l'�chelle ainsi que des mesures r�elles.

II.11 : Extensions diverses :

        Pour clore cette partie, nous allons voir quelques variations sur le th�me des Gaz sur R�seaux. La premi�re famille conserve l'aspect binaire, ou dicret, de HPP/FHP : des mod�les thermiques, comme �voqu�s pr�c�demment, ont �t� �tudi�s (g�om�trie D2Q9) [13]. Diverses �tudes ont permis de mieux contr�ler l'invariance galil�enne du mod�le : par exemple le mod�le FHP-4 avec plusieurs particules immobiles [20][28][29]. Des mod�les � plusieurs phases (mati�res) miscibles ou immiscibles permettent d'injecter des "traceurs" dans des �coulements, ou de simuler la s�paration de deux fluides (comme de l'huile se s�parant de l'eau) [9][15][16]. En 1988, Jean-Pierre Rivet [11] programme un Gaz sur R�seau en 3D � grande �chelle sur un r�seau FCHC (HyperCube � Faces Centr�es). En 1990, C. Appert et S. Zaleski ajoutent des forces non locales entre particules pour calculer la s�paration de phases.

        Les LGA discrets ont vite �t� supplant�s par la famille des LGA en virgule flottante car leur rapport efficacit�/temps de calcul est plus int�ressant. La disponibilit� de plateformes suffisamment puissantes a permis � ce domaine de cro�tre au point de canaliser la plupart des efforts de recherche. Les utilisateurs de CAM-8 sont probablement les derniers � utiliser des mod�les discrets. La premi�re raison d'utiliser les nouveaux mod�les est simple : il n'y a pas besoin d'int�grer de nombreux points pour effectuer une mesure. Les mod�les discrets sont intrins�quement bruit�s et les physiciens d�testent le bruit. Ce bruit est pourtant int�ressant pour faire �merger des ph�nom�nes spontan�ment alors qu'ils sont forc�s dans la pratique (par exemple : dissym�trie dans un tunnel pour forcer des all�es de Von Karman � apparaitre rapidement, voir la partie III.9).

        Pour r�duire la viscosit�, les nombres en virgule flottante sont utilis�s car ils ont une plage dynamique bien sup�rieure � un seul bit. Toutefois on ne peut plus utiliser les techniques de collision classiques ! On utilise alors les m�thodes de Bolzman, l'op�rateur BGK est le plus souvent utilis� actuellement. Le nombre de Reynolds est bien plus grand et les mesures sont plus faciles avec un m�me nombre de sites, m�me si ces sites prennent plus de place en m�moire qu'un site FHP. La m�thode de Bolzman est utilis�e dans presque tous les domaines actuellement, avec autant ou plus de diversit� que le mod�le FHP, m�me si le d�bogage est encore plus difficile qu'avec FHP : lorsque la densit� n'est pas exactement conserv�e, on ne peut pas savoir exactement si c'est une erreur d'arrondi, de frappe ou de formule...

       Enfin, puisque le bruit est n�cessaire dans certains cas, Boghossian et al. [35] ont introduit les ILG : Integer Lattice Gases, sens�s avoir � la fois les avantages des mod�les discrets et continus. Les premi�res exp�riences ont calm� l'optimisme initial, mais tout espoir n'est pas perdu.

II.12 : Conclusion :

        Le mod�le FHP-3 allie une simplicit� relative � des caract�ristiques suffisamment int�ressantes pour justifier les efforts de programmation qui sont effectu�s dans ce domaine. Nous avons toutefois constat� qu'il ne peut �tre utilis� que dans un nombre restreint de cas et dans des conditions s�v�rement contr�l�es mais notre �tude porte uniquement sur les aspects architecturaux et algorithmiques de ce type de mod�les. Cette deuxi�me partie est la seule qui porte sur les aspects purement th�oriques et nous pouvons maintenant nous int�resser aux algorithmes dans tout le reste de ce m�moire, tout en �tant conscient des possibilit�s et des limites du mod�le.

III.1 : Introduction :

        Apr�s la partie pr�c�dente, qui rappelle les points cl�s de m�canique qui nous concernent, cette partie pr�sente les techniques classiques et basiques pour programmer les Gaz sur R�seaux. Elle est la base et le point de d�part du travail r�alis� pour ce m�moire. La premi�re impl�mentation de r�f�rence est donn�e en annexe B : elle a fait l'objet d'un article dans le journal Pascalissime et permet de comprendre les probl�mes de base pos�s par le mod�le FHP3. Nous �tudierons la structure et la conception du programme puis les probl�mes rencontr�s pour arriver � la description d'une deuxi�me impl�mentation de r�f�rence plus sophistiqu�e. Enfin nous analyserons plus en profondeur les caract�ristiques de la programmation des LGA pour pr�parer la version de ce m�moire.

III.2 : Id�es de base :

        Le premier programme, �crit � l'origine en Turbo Pascal, a �t� inspir� par l'article de Pierre Lallemand, paru Revue du Palais de la D�couverte [18] et d�crivant approximativement le mod�le FHP-2. Le seul paragraphe dans lequel la programmation est abord�e est celui-ci :

"
Nous donnons d'abord quelques br�ves indications sur ces exp�riences, pour susciter d'�ventuelles vocations aupr�s des lecteurs disposant d'un ordinateur personnel. Chaque noeud du r�seau est repr�sent� par un nombre form� par la juxtaposition de sept nombres �gaux � 0 ou 1 (bits) : un octet de la m�moire d'ordinateur suffira donc pour d�crire chaque site. L'�volution temporelle se fait en deux �tapes :

        Mais ces indications m'ont longtemps laiss� perplexe. J'ai ensuite cherch� de l'aide dans les th�ses de Val�rie Pot [15] et d'Umberto d'Ortona [14] � Jussieu mais aucun code et aucune indication ne sont fournis. Pourtant, les questions suivantes sont simples :

• Que contient la table, et comment la fabrique-t-on ?
• Que repr�sentent les bits ? Un lien ou un noeud ?

        La repr�sentation des donn�es est un probl�me tout aussi compliqu� lorqu'aucune r�f�rence n'est disponible. En effet, m�me s'il est clair qu'un bit peut repr�senter une particule et qu'un octet permet de repr�senter toutes les directions, un pas de temps comporte de nombreux pas interm�diaires et le r�seau repr�sente l'�tat fig� des particules � un instant qui n'est pas pr�cis� par le mod�le. En effet, un pas de temps commence-t-il par une collision ou un d�placement ? Repr�sentons-nous les particules qui entrent ou qui sortent ? Un bit repr�sente-t-il un lien ou un noeud ? Comment organiser les donn�es pour qu'elles soient facilement manipul�es ?

        Une fois que la transposition d'un r�seau carr� � un r�seau hexagonal est comprise, nous pouvons partir d'un "ether" simple et b�tir l'algorithme de calcul � partir du code de d�placement. Tout d'abord, le calcul est effectu� cycle apr�s cycle, et l'�tat du r�seau repr�sente les particules � chaque cycle, sans se pr�occuper de la sous-�tape (avant ou apr�s collision ou d�placement). Ce qui importe est le changement entre chaque cycle : une particule se propage cycle apr�s cycle en sautant d'un noeud � un autre. Pour effectuer ce changement, la fonction de "calcul" effectue les deux op�rations (propagation et collision), une cellule apr�s l'autre. Effectuer les deux op�rations l'une apr�s l'autre sur tout le r�seau impliquerait la programmation de deux boucles ind�pendantes et augmenterait le nombre d'acc�s � la m�moire. Dans tous les algorithmes qui suivront, les deux op�rations sont effectu�es � l'int�rieur de la m�me boucle pour b�n�ficier de la localit� des registres.

        La technique d'�laboration du programme est simple :
    - D'abord mettre en place la boucle externe avec la pr�paration des pointeurs vers les donn�es.
    - Ensuite, programmer la propagation des particules : charger un octet de m�moire et distribuer un par un tous les bits vers les noeuds voisins, avec la gestion des variables temporaires (pour �viter le recouvrement dans certaines directions, comme expliqu� dans le chapitre suivant).
    - Enfin, lorsque le d�placement est fonctionnel (tous les bits se d�placent comme pr�vu sur le r�seau), il ne reste plus qu'� inclure la consultation de la table, juste apr�s l'endroit o� le noeud courant est lu.

        Il devient clair aussi que la table, puisqu'elle est consult�e entre le chargement et la distribution, doit contenir une repr�sentation aussi simple et compl�te que possible des op�rations � effectuer. La table peut effectuer des op�rations complexes et r�duire la taille du programme, elle concentre donc toute l'intelligence du code de collision. L'annexe B montre deux mani�res de la programmer : par code explicite ou par constante � la compilation. La table exploite ainsi le huiti�me bit pour effectuer la collision avec des "murs" virtuels sans ajouter une seule ligne de code. Des effets de pesanteur ou d'attraction peuvent m�me �tre ajout�s en modifiant l�g�rement certaines entr�es de la table. C'en est presque trop facile ...

III.3 : Les plans temporaires :

        Un aspect m�connu des programmes de ce type concerne les "plans temporaires" : le probl�me ne peut �tre appr�hend� dans toute sa complexit� que lors de la programmation, lorsqu'il est d�j� trop tard. Ce probl�me deviendra encore plus pr�pond�rant avec les codes de strip mining et il est important pour la suite du travail de maitriser parfaitement l'algorithme et les donn�es associ�es.

        Les "plans temporaires" deviennent tr�s important pour les g�om�tries tr�s larges car m�me avec un ordinateur disposant de toute la m�moire vive du monde, il est important de l'utiliser correctement. Les plans temporaires sont d�cris succintement dans l'annexe B et dans la partie V.6 mais �tudions ici leur th�orie g�n�rale en partant d'un cas de dimension 1 (par exemple un automate cellulaire lin�aire) :

        Dans un ordinateur "s�quenciel", les noeuds sont trait�s un par un. Supposons que le sens de balayage soit le m�me que le sens de x et regardons ce que ferait un algorithme simple mais mauvais :

1) d'abord il lit la valeur du noeud courant en x,
2) ensuite il calcule la valeur suivante,
     par consultation de la table par exemple,
3) enfin il �crit une partie du r�sultat dans chaque
     partie correspondante : x-1, x et x+1
4) il incr�mente x et retourne en 1) si la ligne n'est pas termin�e.

        La solution la plus simple est d'utiliser deux plans de travail : un plan "source" et un plan "destination", permettant une complexit� arbitraire dans le voisinage avec l'argument choc que le "transfert" s'effectue simplement en �changeant les deux pointeurs vers les tableaux :
 

        Toutefois, lorsque la taille totale des tableaux atteint les limites de la m�moire de l'ordinateur, il est clair que seule une moiti� est vraiment utile car les informations sont redondantes ou inutilis�es. Cette consid�ration devient incontournable pour les simulations en 3D car elles utilisent une quantit� ph�nom�nale de m�moire (de 16MO � 4GO dans certains cas). L'argument de simplicit� de la technique de l'�change de pointeurs perd toute cr�dibilit� devant ce probl�me.

        La solution privil�gi�e, bien que plus complexe, utilise un "plan temporaire" pour m�moriser le r�sultat de chaque cycle et �viter le "court circuit temporel" d�crit plus haut. Les caract�ristiques de cette m�moire d�pendent de la quantit� d'informations qui doit voyager dans le m�me sens que celui du balayage. Pour l'exemple de l'automate cellulaire lin�aire d�crit plus haut, il faut m�moriser un mot ou un bit car le d�placement se fait vers le voisin imm�diat. Au niveau de l'algorithme, cela se traduit par la n�cessit� de retarder l'�criture du r�sultat. Le programme est dot� d'une valeur temporaire (plan ponctuel) et devient ainsi :

1) lire la valeur du noeud courant en x,
2) �crire en x la valeur temporaire,
3) calculer la valeur suivante de x,
4) �crire chaque partie du r�sultat dans chaque partie correspondante :
   x-1, x mais x+1 va dans la valeur temporaire.
5) incr�menter x et retourner en 1) si la ligne n'est pas termin�e.

        Une difficult� suppl�mentaire est d'amorcer le programme en fournissant une valeur initiale correspondant � l'extr�mit� de la ligne : cela peut �tre astucieusement utilis� pour injecter des particules dans le domaine d'�tude et cr�er un vent artificiel. De m�me, la valeur finale de la variable temporaire (� la fin de la ligne) peut �tre ignor�e pour faire disparaitre les particules et cr�er une sorte de bouche de sortie pour le fluide. Les particules peuvent ainsi �tre cr��es puis effac�es � des extr�mit�s oppos�es du tunnel, ce qui g�n�re naturellement un vent dans le tunnel. Des techniques plus sophistiqu�es sont cependant recommand�es car cette m�thode est aussi simple que limit�e dans la pratique.

        Toutefois nous traitons des tableaux en deux dimensions et les choses ont tendance � s'emm�ler et rendre la programmation tr�s d�licate. Il y a pourtant une r�gle simple � retenir : il faut un buffer temporaire par dimension (que la dimension soit temporelle ou spatiale) et l'information � m�moriser correspond au voyage qu'effectuent les particules dans le sens du balayage.

        Pour appliquer l'algorithme en 2 dimensions, il faut un plan temporaire ind�pendant pour chaque dimension. Les colonnes sont trait�es exactement de la m�me mani�re qu'une ligne mais apr�s projection sur la dimension perpendiculaire : ce n'est plus un noeud de m�moire qui est n�cessaire pour le plan temporaire mais toute une ligne. Pour un tableau de x * y noeuds, il faudra en tout M = 1 + x + xy noeuds en m�moire. En r�gle g�n�rale, on compte (x+1)*y noeuds pour l'allocation de la m�moire d'un tableau 2D.

        La formule se g�n�ralise facilement � toute dimension N>1 et elle se r�duit approximativement � la un polyn�me si toutes les dimensions sont similaires. On peut ainsi prouver que la quantit� totale de m�moire ne s'approche pas du double de la taille du tunnel comme dans l'alternative pr�c�dente. Par exemple, pour un r�seau en 3D de dimensions (x,y,z) avec un voyage d'un noeud par pas de temps, il faudra en tout M = 1 + x + xy + xyz noeuds de m�moire. Si x, y et z sont des valeur rapproch�es, la formule devient le polyn�me suivant : M = x^0 + x^1 + x^2 + x^3 = 1 + x + x^2 + x^3. Elle tend ainsi vers M = (x+1)^N  et permet, pour un cas donn�, d'utiliser moins de m�moire ou d'avoir un tableau plus grand, par rapport � la technique de l'�change de pointeurs (M = 2x^N).

        Naturellement, la complexit� du programme augmente mais comme pour le reste cela d�pend de l'exp�rience, des ressources, de la patience et de la compr�hension de la technique. Un cas interm�diaire (compromis complexit�/espace m�moire) serait de diviser le tunnel en de nombreuses parties et de disposer d'un plan temporaire : chaque bloc ayant la m�me taille, il est facile d'utiliser la technique d'�change de pointeurs sans pour autant doubler l'occupation de la m�moire. Cette technique ne permettant pas de r�duire directement le temps de calcul, elle n'est pas �tudi�e ici. De plus, les fortes contraintes en m�moire sur les PC ainsi que les m�moires caches (principalement les modes write through et write back qui ne peuvent pas toujours �tre contr�l�s) favorisent la technique de plan temporaire minimal, d�crite plus haut.

III.4 : Premier code de r�f�rence :

        Nous allons ici �tudier succintement un premier morceau de code qui servira de r�f�rence pour estimer les performances, les contraintes et les limites des algorithmes r�alis�s. Il est extrait de l'annexe B et c'est la version �crite en assembleur pour i286 en 1995 :

BEGIN
(*...*)
 asm
(*...*)

(* boucle principale *)

@BOUCLE_EXTERIEURE:
  mov seg_ligne,$A000
  mov y,99
@BOUCLE_Y:
   mov es,seg_ligne

(* PARTIE IMPAIRE: *)

(* force la particule D *)
    xor cl,cl
    rol rand,1
    jnc @pas_retenue1
    mov cl,1
@pas_retenue1:

    mov bp,1            (* BP est le registre de contrôle de la boucle *)
    mov di,xmax+2               (* DI pointe sur le noeud courant *)
    mov si,offset temp_impair+1   (* SI pointe sur les bits E et F *)
@BOUCLE_IMPAIRE:

(* Al collectera les bits du noeud courant *)
     lodsb                         [1]
     mov byte[si-1],0              [2]

(* Ah désignera les bits à envoyer *)
     mov ah,byte ptr es:[di]       [3]
     or ah,ah                      [4]
     jz @vide1                     [5]

(* consulte le tableau: *)
     mov bl,ah                     [6]
     xor bh,bh                     [7]
     rol rand,1                    [8]
     jnc @pas_rol1                 [9]
     inc bh                        [10]
@pas_rol1:
     mov ah,byte[bx+offset p]      [11]

(* distribue les bits: *)
     shr ah,1                      [12]
     jnc @pas_A1                   [13]
     or byte ptr es:[di-1],A       [14]
@pas_A1:
     shr ah,1                      [15]
     jnc @pas_B1                   [16]
     or byte ptr es:[di-xmax-1],B  [17]
@pas_B1:
     shr ah,1                      [18]
     jnc @pas_C1                   [19]
     or byte ptr es:[di-xmax],C    [20]
@pas_C1:
     shr ah,1                      [21]
     jnc @pas_D1                   [22]
     or cl,2                       [23]
@pas_D1:
     shr ah,1                      [24]
     jnc @pas_E1                   [25]
     or byte ptr ds:[bp+offset temp_pair],E [26]
@pas_E1:
     shr ah,1                      [27]
     jnc @pas_F1                   [28]
     or byte ptr ds:[bp+offset temp_pair-1],F [29]
@pas_F1:
     shl ah,6                      [30]
     or al,ah                      [31]
@vide1:
     shr cl,1                      [32]
     jnc @pas_retenue_D1           [33]
     or al,8                       [34]
@pas_retenue_D1:
     stosb                         [35]
     inc bp                        [36]
    cmp bp,xmax                    [37]
    jbe @BOUCLE_IMPAIRE            [38]


(* PARTIE PAIRE: *)

(* coup�e : elle est similaire � la partie impaire *)

   add seg_ligne,40
   dec y
  jnz @BOUCLE_Y

  mov ah,1
  int 016h
  jz @BOUCLE_EXTERIEURE
 end;
(* boucle principale *)

END.

        or byte ptr ds:[bp+offset temp_pair],E : nous atteignons ici des sommets de programmation CISC. Rappel : le nombre r�duit de registres nous oblige � utiliser comme compteur/pointeur le pointeur de trame BP (frame Base Pointer), alors qu'il est sens� servir pour le passage de param�tres sur la pile pour les langages de haut niveau. Le probl�me est que tout adressage utilisant BP utilise par d�faut le segment de pile SS alors que les donn�es sont adress�es par ES (le tableau dans la m�moire vid�o) et DS (les donn�es statiques et le plan temporaire). Nous devons donc ajouter un pr�fixe d'adressage ("DS:") pour restaurer le segment, ce qui consomme inutilement des cycles. Rappelons aussi que sur le i286, bien qu'il existe de nombreux modes d'adressage, peu de registres peuvent �tre effectivement utilis�s comme pointeur. Il est donc facile de comprendre que c'est l'architecture bancale de la machine qui fait dire � certains que "les compilateurs peuvent g�n�rer du code aussi efficace qu'un code �crit en assembleur", puisque la flexibilit� r�duite diminue artificiellement les �carts possibles de performance.

        Nous voyons aussi, dans le reste du code, de nombreuses boucles, des manipulations de segments, des calculs de pointeurs et la gestion du plan temporaire horizontal : de nombreuses parties du programme seraient susceptibles d'�tre �limin�es en retravaillant la structure globale. Bien que le programme en assembleur soit environ deux fois plus rapide que le programme en Pascal, il semble effectuer des op�rations trop complexes pour ne bouger que quelques bits. Il faut environ 38 instructions pour traiter une seule cellule et le traitement est ralenti pour de nombreuses raisons, dont :

• l'acc�s aux instructions qui ne sont pas "cach�es" et utilisent une partie de la bande passante de la m�moire (le probl�me n'existe plus depuis le i486)
• les acc�s aux donn�es en m�moire vid�o qui entrent en concurrence avec le balayage de la carte
• les nombreux sauts en avant de quelques octets qui r�duisent l'efficacit� de la queue de pr�d�codage interne

        Les mesures manuelles sur un PC i286 � 12MHz donnent une vitesse d'environ 100000 noeuds par seconde, la g�om�trie �tant limit�e architecturalement � 320*200 par la m�moire segment�e et la carte VGA (granularit� de 64 Ko). Il y a environ une image et demie affich�e par seconde ce qui est honnorable pour un ordinateur de cette classe (comparer � l'iPSC en D.6). La recherche des goulots d'�tranglements est difficile car le processeur i286 ex�cute les instructions avec des dur�es tr�s variables et le temps d'acc�s � la m�moire est tr�s fluctuant. Toutefois nous pouvons effectuer quelques estimations : le calcul d'un noeud n�cessite environ 120 cycles d'horloge et 38 instructions, soit 3,2 cycles d'horloge par instruction en moyenne.

III.5 : Influence de l'architecture :

        Le passage du i286 � 12MHz au Pentium � 100MHz en 1996 fut un grand bouleversement. Les r�gles de codage ont �t� profond�ment modifi�es : les instructions recommand�es ne sont pas les mêmes, la m�moire cache sur deux niveaux peut enfin contenir tout le code et toutes les donn�es (L1 : 8Ko et L2 : 256Ko), de nombreuses innovations sont incluses... Pourtant le mode r�el ne b�n�ficie pas directement de toutes les am�liorations, comme par exemple les modes d'adressage �tendus et orthogonaux du i386 qui permettent d'utiliser toutes les combinaisons de registres pour acc�der � la m�moire.

        La m�moire virtuelle en mode prot�g� reste un probl�me car la manipulation des registres de segments d�clenche de longues s�quences de microcode pour v�rifier la validit� des adresses, ce qui en plus peut d�clencher des exceptions pour g�rer les pages en m�moire et sur le disque dur s'il faut les swapper. On s'attend donc naturellement � ce que le programme fonctionne moins vite lorsqu'il est lanc� � partir de Windows. Pourtant, il s'y ex�cute visiblement plus vite ! Plus pr�cis�ment, en mode fen�tr�, il s'ex�cute beaucoup plus vite qu'en mode normal/"plein �cran".

        Le changement d'architecture, du 286 au Pentium, a chang� compl�tement les rapports entre les flux de donn�es et d'instructions. L'introduction de deux niveaux de m�moire cache balance la plus grande latence relative de la m�moire vid�o. De plus, l'acc�s � la m�moire vid�o subit de nombreuses restrictions : par exemple, elle est acc�d�e au travers du bus PCI et elle n'est pas cachable. L'exp�rience avec Windows a montr� la situation paradoxale o� un programme "optimal" pour le 286 est inefficace sur Pentium. Pour expliquer cela, il faut pr�ciser que Windows (en mode fen�tr�) �mule la m�moire vid�o : il redirige les acc�s vid�o vers la DRAM centrale gr�ce � la translation d'adresse du mode prot�g�/pagin�. Ainsi, nous n'acc�dons pas nous-m�mes � la m�moire vid�o, c'est l'�mulateur qui copie la DRAM p�riodiquement vers la vraie carte. La m�moire centrale �tant cachable, les acc�s sont beaucoup plus rapides et le transfert par bloc vers l'�cran est plus rapide qu'octet par octet. Les instructions de type read-modify-write sont encore plus lourdes puisqu'elles n�cessitent une lecture (blocante, pour �viter toute modification intempestive au milieu de l'instruction) puis une �criture et il faut une transaction compl�te sur le bus PCI (attente que le bus soit libre - envoyer l'adresse - envoyer/recevoir la donn�e) pour transf�rer un simple octet ou un bloc entier. Le calcul en m�moire cachable puis le transfert par bloc (ou en rafale, burst pour le PCI) est forc�ment plus rapide que le premier programme.

        Le diagramme ci-contre est extrait de [33]. Il repr�sente une carte m�re aux caract�ristiques suivantes :

* Processeur Pentium 100MHz avec 2 caches int�gr�es de 8Ko chacune * Bus local : donn�es 64 bits, 66 MHz, 528Mo/s th�oriques * Cache L2 de 256Ko avec des puces SRAM � 15 ns de temps d'acc�s (sans compter les latences de gestion de la LRU et de la SRAM de "tag" � 15 ns aussi) * M�moire DRAM � 70ns de temps d'acc�s (30ns en mode page), EDO possible. * Bus PCI multiplex� � 33MHz et 32 bits de large (132Mo/s th�oriques)

        Nous pouvons y voir que le nombre de circuits et de connexions � traverser est proportionnel � la latence et � la vitesse de transfert. Le circuit imprim� montre en plus que cela est proportionnel � la distance parcourue dans les fils : la m�moire cache est tr�s proche du processeur alors que les slots PCI et DRAM sont �loign�s. Des mesures r�elles sur cette plateforme, en conditions "id�ales", donnent une vitesse de transfert unidirectionnel de 28Mo/s vers la carte vid�o. En pratique, le processeur peut donc afficher un mot de 32 bits tous les 4 cycles, soit une instruction d'affichage entrelac�e avec sept instructions de calcul. Il n'est pas raisonnable dans ces conditions d'utiliser le premier code de r�f�rence qui utilise de nombreuses op�rations complexes vers la m�moire vid�o.

        La le�on est simple : en mode "normal" (MS-DOS) il ne faut plus utiliser la m�moire vid�o pour stocker le tunnel. Il faut au contraire utiliser une m�thode qui paradoxalement prendrait plus de temps si l'on ne regardait que le nombre d'instructions � ex�cuter. Le calcul doit donc s'effectuer dans la m�moire centrale et, comme l'�mulateur MSDOS de Windows le faisait, afficher r�guli�rement le r�sultat � l'�cran, tous les N pas de temps (N pouvant varier � la demande de l'utilisateur). Nous �vitons ainsi de perdre du temps dans les transactions "blocantes" du bus PCI, dont le temps d'acc�s est probablement sup�rieur � 100ns (carte vid�o comprise). Pour qu'un programme tourne vite, il est crucial de placer les donn�es au meilleur endroit.

III.6 : Deuxi�me code de r�f�rence :

        Le second code de r�f�rence date de 1997 et est inclus ci-dessous. L'id�e directrice est de r�duire le code de distribution des bits en groupant les noeuds par quatre. Le r�seau est alors tourn� � 90 degr�s et il n'y a plus de code pour les lignes paires et impaires, comme pour le circuit d�crit en V.6. Ce code est destin� aux plateformes � partir du i386 et "optimis�" pour le Pentium classique (P53C) o� nous pouvons utiliser les registres sur 32 bits au lieu de 16. La boucle externe affiche les calculs puis traite N pas de temps de cette mani�re :

;
;  BALAYAGE:
;
   mov ebx,es:[640]
   mov XB,ebx

   mov di,640
boucle_ext:
     mov bx,320
boucle_int:
       mov word ptr X,bx        ;U (1)

;
; d�placement des variables:
;
       mov esi,XC               ;*V                 [3]
       mov ebx,es:[di-316]      ;*prefixe + U (3)   [6]
       mov XD,esi               ;*V                 [8]
       mov XC,ebx               ;*U (4)             [10]

       mov edx,X1               ;*V                 [12]
       mov eax,XB               ;*U (5)             [14]
       mov ecx,es:[di+4]        ;*pr�fixe+U (7)     [17]
       and edx,010000000h   ;XE ;*V                 [19]
       mov X1,eax               ;*U (8)             [21]
       mov XB,ecx               ;*V                 [23]

       mov ebp,XA               ;*U (9)             [25]
       mov ebx,es:[di+320]      ;*pr�fixe+U (11)    [28]
       and ebp,020000000h   ;XF ;*V                 [30]
       mov XA,ebx               ;*U (12)            [32]
       or ebp,edx               ;*V                 [34]

;verticaux:
       and eax,0C0C0C0C0h       ;*U (13)            [36]
       and esi,008080808h       ;*V                 [38]
       and ebx,001010101h       ;*U (14)            [40]
       or eax,esi               ;*V                 [42]
;XL
       mov ecx,XD               ;*U (15)            [44]
       or eax,ebx ;report des verticaux *V          [46]
       mov ebx,XA               ;*U (16)            [48]
       and ecx,000100010h       ;*V                 [50]
       and ebx,000002000h       ;*U (17)            [52]
       or ebp,ecx               ;*V                 [54]
       mov ecx,X1               ;*U (18)            [56]
       or ebp,ebx               ;*V                 [58]
       and ecx,000201020h       ;*U (19)            [60]
       or ebp,ecx               ;*U (20) (d�pendance sur ECX) [62]
;XR:
       mov ebx,XB               ;*V                 [64]
       rol ebp,8  ; report XL    *U (21)            [66]
       mov ecx,XC               ;*V                 [68]
       or eax,ebp ; report XL    *U (22)            [70]
       mov ebp,XD               ;*V                 [72]
       and bx,2                 ;U (23)             [72']
       and ebp,000040000h       ;*V                 [74]
       and cx,4                 ;U (24)             [74']
       mov edx,XA               ;*V                 [76]
       or bx,cx                 ;U (25)             [76']
       mov ecx,X1               ;*V                 [78]
       or bp,bx                 ;U (26)             [78']
       and edx,002000200h       ;*V                 [80]
       or ebp,edx               ;*U (27)            [82]
       and ecx,004020400h       ;*V                 [84]
       or ebp,ecx               ;*U (28)            [86]
       ror ebp,8                ;*V                 [88]
       or eax,ebp               ;*U (29)            [90]

; d�termine quelle banque LUT on utilise:
       rol word ptr seed1,1     ;U (32) 3 cycles + non pairable  [94]
       setc bh                  ;U (35) microcode                [97]

       mov bl,ah                ;U (36) d�pendance sur EBX       [98]
       mov ah,[offset p+bx]
  ; acc�de au tableau (#1)  U (38) AGI sur EBX, pr�fixe possible [101]
       mov bl,al                ;U (39) d�pendance sur EAX       [102]
       mov al,[offset p+bx]
  ; acc�de au tableau (#2)  U (41) AGI sur EBX, pr�fixe possible [105]
       ror eax,16               ;U (42) d�pendance sur EAX       [106]

       mov bl,ah                ;U (43) d�pendance sur EAX       [107]
       mov ah,[offset p+bx]
  ; acc�de au tableau (#3)  U (45) AGI sur EBX, pr�fixe possible [110]
       mov bl,al               ;U (46) d�pendance sur EAX        [111]
       mov al,[offset p+bx]
  ; acc�de au tableau (#4)  U (48) AGI sur EBX, pr�fixe possible [114]
       ror eax,16               ;U (49) d�pendance sur EAX       [115]

;writeback:
       mov bx,word ptr X        ;V                  [115']
       mov ecx,[offset temp+bx] ;**U (50)           [117]
       mov [offset temp+bx],eax ;*V                 [119]
       mov es:[di-320],ecx      ;pr�fixe+*U (51)    [122]

       add di,4                 ;V                  [123]
       sub bx,4                 ;U (52)             [123']
       jnz boucle_int           ;V                  [124]

     cmp di,63360
     jbe boucle_ext

        Mais en r�alit�, le code a �t� con�u pour �tre ex�cut� en mode r�el (16 bits sous MS-DOS) et l'utilisation des registres en mode 32 bits ajoute un pr�fixe de taille, invisible dans la syntaxe, � chaque instruction : cela brise compl�tement le pairage des instructions ! La plupart des d�pendances de donn�es entre les registres ont �t� aplanies dans le code de mouvement afin d'ex�cuter deux instructions par cycle sur le Pentium mais les pr�fixes de taille �taient oubli�s dans l'analyse. Les instructions marqu�es d'un ast�risque prennent ainsi deux cycles au lieu d'un demi-cycle ! Cela montre bien que le nombre d'instruction, le temps d'ex�cution et le nombre de sites trait�s par seconde sont des valeurs qui ne sont plus naturellement ou simplement proportionnelles : l'architecture complexe du Pentium a des "cycles cach�s" qui n'apparaissent pas en lisant le code source, m�me en assembleur ! Alors que ce code est tr�s bon pour le i386, le Pentium a un pipeline diff�rent qui favorise certains types de codes tout en pouvant ex�cuter les autres programmes mais � une vitesse nettement inf�rieure et sans pr�venir.

        L'analyse sous ce nouvel angle (voir les nombres entre crochets) donne environ 124 cycles (comme pour la version i286, mais 31 cycles par noeud). Cette baisse incroyable des performances, malgr� le soin apport� au code, a motiv� la conception d'un DOS-extender pour la suite du projet. Le code de mouvement a �t� assez bien r�duit mais reste lourd et consomme tous les registres. L'analyse du code r�v�le aussi que le processeur effectue des op�rations inutiles ou perd du temps � manipuler des donn�es entre les registres. Ainsi, la consultation des tables, en raison du faible nombre de registres disponibles, du jeu d'instructions, des pr�fixes et des d�pendances crois�es, prend 20 cycles au moins, bien qu'en th�orie cela n�cessite 4 cycles pour un cas "simple et id�al".

La densit� globale dans le domaine d'�tude peut �tre d�termin�e en accumulant le nombre de particules � chaque pixel. Nous voyons ici une accumulation sur 8 bits avec deux d�passements (les m�mes couleurs sont utilis�es deux fois pour des densit�s diff�rentes) et un �coulement de Poiseuille caract�ristique. L'image peut �tre calcul�e interactivement (les param�tres peuvent �tre chang�s) en une minute environ.

        L'accumulation doit commencer apr�s la stabilisation du fluide et la disparition des ondes de choc. Il faut ici pouvoir g�rer trois tableaux de 64 Ko, ce qui est � la limite des possibilit�s du mode r�el du x86. Avec un code de ce type, il a fallu deux heures � un Pentium � 100 MHz pour calculer l'image suivante :

A chaque pas de temps, la densit� de chaque noeud est accumul�e sur 16 bits dont nous voyons ici les 8 bits de poids fort. Il faut donc un tableau de 128Ko en m�moire et acc�der en tout � plus de 200Ko, ce qui a n�cessit� l'emploi du mode flat (ou unreal). Ce programme est donc difficilement transportable et n�cessite une configuration tr�s pr�cise pour fonctionner.

        Le bruit au niveau de l'octet de poids fort est r�duit mais il a fallu beaucoup de temps de calcul pour que toute la dynamique de l'octet soit utilis�e et donne une image utilisant toutes les couleurs de la palette (plus de 64K*3=200000 pas de temps). La vitesse de rafraichissement pour ce type de code sur cette machine atteint environ 10Hz et la taille est fix�e � 320x200 noeuds, c'est � dire la r�solution de l'�cran.

III.7 : Conditions aux limites et effets de bords :

        Les LGA permettent de calculer les �quations diff�rentielles de Navier-Stokes dans des conditions id�ales de tr�s faible vitesse et de taille infinie du r�seau [19]. Toute autre condition ne permet d'effectuer qu'une approximation, par exemple lorsque la taille est finie (limit�e par la taille de la m�moire de l'ordinateur), ou alors la simulation ne correspond pas aux r�alit�s physiques. Entre autres exemples, la vitesse du fluide ne doit pas d�passer Mach 0,3 environ mais en plus des limites th�oriques, les limites pratiques sont aussi difficiles � connaitre et � comprendre.

        Nous allons �tudier en particulier le cas d'une �prouvette dans une soufflerie ainsi que les effets de bords li�s aux conditions aux limites. Dans les images pr�c�dentes, nous pouvons observer certains ph�nom�nes qui influencent la simulation et ses r�sultats. Tout d'abord, avant chaque pas de temps, le vent est g�n�r� par une boucle qui cr�e de nouvelles particules afin de g�n�rer un "vent" dans le tunnel. Les propri�t�s de ce flux ont des cons�quences directes sur le temps de disparition des ph�nom�nes transitoires et ind�sir�s.

        La mani�re la plus simple, comme le permet l'algorithme de plan temporaire d�crit pr�c�demment, est de cr�er une particule directement dans le sens du vent. Cela correspond bien � la th�orie mais dans la r�alit�, un vent n'est pas constitu� de particules allant toutes dans la m�me directions : ce serait un vent supersonique ! De plus, les ph�nom�nes d�passant Mach 0,4 ne sont pas fid�les � la r�alit�. Il faut donc cr�er un flux d'air, non une masse homog�ne ou homocin�tique irr�elle. La vitesse et la direction des particules doivent �tre bien choisies. La vitesse des particules ne peut pas �tre modifi�e car elle est unitaire : la vitesse du flux sera r�gl�e par la "vitesse" d'introduction et de destruction des particules de chaque c�t� du tunnel. Une autre m�thode est de "forcer" certaines particules au hasard dans le tunnel en changeant leur direction mais c'est plus compliqu� car le tunnel devra �tre boucl� sur lui-m�me comme un tore, ce qui implique qu'en se propageant, les turbulences se retrouveront en amont et se perturberont elles-m�mes. Cette m�thode est adapt�e � d'autres cas (�tude d'un cisaillement pour d�terminer la viscosit� par exemple) mais pas � celui des simulations qui nous int�ressent (m�me si ce n'est pas le domaine le plus adapt� pour les LGA).

        Pour cr�er un vent, nous introduisons des particules � un certain ryhtme, qui est contr�l� au clavier par l'utilisateur dans les exp�riences r�alis�es. Leur direction est aussi importante : elle doit �tre d�corr�l�e, aussi al�atoire que possible, pour �viter que les caract�ristiques du vent n'influencent ou ne perturbent les ph�nom�nes en aval. Pour d�corr�ler un signal, il suffit de le moduler par un signal connu pour �tre al�atoire : nous disposons d�j� de ph�nom�nes explicitement al�atoires au niveau de certaines configurations de collision. La technique est alors de cr�er des particules allant en sens inverse du vent pour, si elles en ont le temps, frapper le mur adjacent. Cela cr�e une zone de haute densit� propice � de nombreuses collisions, favorisant une r�partition al�atoire et naturelle des particules, m�me si elles sont cr��es de mani�re d�terministe. Leur surnombre local et la configuration de la paroi forcent ainsi un flux de particules al�atoires dans le sens d�sir� (en moyenne). L'entropie du mod�le leur assure une r�partition naturellement �quilibr�e dans le temps et dans l'espace.

        Le deuxi�me probl�me important concerne les autres conditions aux limites. D'abord, les parois horizontales (sup�rieure et inf�rieure) sont "rugueuses" : le mod�le le plus simple sp�cifie que les particules repartent par le lien o� elles �taient venues. Ensuite, la paroi de droite est totalement "absorbante" et fait disparaitre toutes les particules. Dans la r�alit�, cela correspondrait � un tuyau d�bouchant dans le vide sid�ral, ce qui n'est pas notre intention : une veine de soufflerie r�elle maintient une densit� assez homog�ne autour du domaine d'�tude et le vent n'est pas "aval�" d'une mani�re ou d'une autre. Les conditions sont r�unies pour faire apparaitre un �coulement parabolique de Poiseuille qui interf�re avec le domaine d'�tude comme dans l'image de densit� sur 8 bits. La densit� n'est pas homog�ne ni lin�aire autour de l'�prouvette et le ph�nom�ne de portance que l'on veut mettre en �vidence est perturb� par l'influence des parois et de la disparition des particules, ce qui est similaire � une anisotropie de la pression dans tout le tunnel.

        La solution adopt�e pour l'exp�rience suivante (l'accumulation sur 16 bits) est de r�tablir la pression dans le tunnel gr�ce � une sorte de "peigne" qui pi�ge une partie des particules qui allaient dispara�tre. De proche en proche, elles cr�ent une pression qui s'oppose en partie � leur disparition totale : le fluide devient plus homog�ne et les directions des particules sont plus diverses. Nous voyons sur l'image de densit� 16 bits que le profil n'est plus parabolique, bien que les parois rugueuses laissent subsister de l�g�res perturbations locales. La pression dans le tunnel est mieux r�partie et �quilibr�e.

        Une deuxi�me solution pour r�duire encore plus, ou anihiler, l'�coulement parabolique est de rendre les parois glissantes pour que leur vitesse par rapport au vent ne soit pas chang�e. Cette solution a par exemple �t� adopt�e pour l'all�e de von Karman dans l'exemple de David Hanon en III.9. Cela n'a pas �t� essay� dans les premiers codes car le programme ne le permettait pas encore, seules les parois rugueuses pouvant �tre cod�es dans un octet. En compl�ment de la premi�re solution, nous disposons alors d'un domaine d'�tude dont la pression et la vitesse sont homog�nes et favorables pour �tudier des ph�nom�nes turbulents comme les all�es de von Karman. Dans le cas contraire, par exemple si l'�coulement de Poiseuille subsistait, il faudrait agrandir le tunnel de mani�re � ce que les paraboles n'interviennent pas substanciellement dans les mesures, ce qui augmenterait quadratiquement les besoins de m�moire et de temps de calcul.

        Pour que ces conditions aux limites soient remplies, il faut plus de flexibilit� dans le programme. La nature du mod�le FHP permet d'ajouter facilement les fonctions d�sir�es mais leur impl�mentation est souvent une aventure plus complexe que l'on pourrait s'y attendre. De plus, le bon sens des �quations classiques contredit la r�alit� microscopique des particules : les equations d'Euler ne traitent pas explicitement du mouvement chaotique brownien. Un bon programme de calcul FHP a donc besoin de parois glissantes comme rugueuses et d'un contr�le s�r de la cr�ation et de la destruction des particules. Les algorithmes vus jusqu'� pr�sent ne permettent pas de tel contr�le sur la rugosit� des parois.

III.8 : Remplissage et redimensionnement dynamique du tunnel :

        Une fois le fluide en r�gime "stationnaire", c'est � dire lorsque tous les ph�nom�nes transitoires (comme les ondes de choc ou les ph�nom�nes anisotropiques hexagonaux) ont disparu, on effectue les mesures et on r�fl�chit aux calculs suivants. En d'autres termes, une grande partie des calculs a pour seule fonction d'�liminer les transitoires. Une formule empirique : t=4*(h+l) donne le nombre de pas de temps � calculer pour que disparaisse la plupart des ondes de choc, soit deux � trois allers et retours de l'onde de paroi � paroi. Or les simulations les plus int�ressantes sont les plus grosses et le temps de calcul augmente tr�s vite ; en utilisant la formule pr�c�dente, il faut environ : N=4*(h+l)*h*l pas de calcul en tout, ou environ N=8*l^3 si h et l sont proches. Dans ce chapitre, nous allons explorer des moyens de r�duire ce temps de calcul.

        Pour commencer, nous pouvons constater que la veine de simulation FHP est souvent vide au d�but de l'exp�rience. Il faut l pas de temps et h*l*l = hl� calculs de site pour remplir le tunnel en particules. Ce sont autant de cycles facilement gagn�s si la veine contenait d�j� des particules !

        Le probl�me est maintenant : comment remplir le tunnel ? Il est facile de le remplir uniform�ment de particules, ou avec des particules au hasard, mais d'une part c'est trop simple et d'autre part il faut g�rer le cas d�licat des �prouvettes qui seraient remplies de particules. En effet, la forme que l'on place dans le tunnel n'est pas obligatoirement "pleine" et il faut que l'algorithme de remplissage tienne bien compte de ce cas de figure. Comme dit plus haut, un algorithme simple ne convient pas et il faut utiliser un algorithme de remplissage "intelligent" tel qu'on le rencontre dans les logiciels de dessin bitmap (exploration r�cursive ou floodfill).

        La premi�re simplification oubliait aussi que ce qui nous int�resse est d'obtenir un �tat le plus proche possible du r�gime stationnaire. Or le cas le plus simple consomme beaucoup de temps dans la mise en mouvement le fluide. Une possibilit� serait d'utiliser une r�solution statique et approximative des �quations d'Euler ou de Navier-Stokes pour "guider" le remplissage pr�liminaire du tunnel. La limitation ne se situe pas dans la charge de calcul mais dans sa complexit� qui augmente avec celle de la g�om�trie des parois. De plus, ce n'est pas un sujet qui nous concerne dans cette �tude.

        Une autre m�thode serait d'utiliser les mod�les plus simples comme FHP-2 ou FHP-1 pour effectuer une premi�re passe dans le fluide. Ces mod�les sont plus a priori rapides � calculer et peuvent effectuer une premi�re approximation, ce qui est int�ressant puisque nous abandonnons la technique de consultation de table. Nous verrons pourtant que le calcul a une importance aussi grande que le d�placement des donn�es ("FHPIII est memory bound") et cette m�thode ne r�duit pas le nombre de mouvements de particules. Toutefois, nous avons vu dans la partie pr�c�dente que la principale diff�rence entre les versions du mod�le FHP concernent la viscosit�, ou le nombre de Reynolds par noeud. Si nous pouvons calculer � un nombre de Reynolds plus faible, nous avons donc int�r�t � r�duire le nombre de noeuds plut�t que d'augmenter la viscosit� : le programme de calcul reste identique mais nous pouvons r�duire le nombre de d�placements en r�duisant la taille du tunnel lors de l'amor�age du fluide.

        Pour r�duire le temps total de calcul, nous avons donc int�r�t � commencer avec une version dont la taille est une fraction de la taille r�elle. En pratique, le programme de ce m�moire limite la largeur minimale � 256 noeuds. Nous pouvons donc amorcer le tunnel avec cette largeur, durant le temps n�cessaire � la disparition des transitoires, soit environ 2000 pas de temps. Ensuite, le tunnel est agrandi : le cas le plus simple est un doublement de toutes les dimensions et un quadruplement ais� du nombre de cellules et de leur contenu (le contenu est recopi� dans 3 cellules voisines). Le calcul continue et fait disparaitre les transitoires li�es � l'agrandissement soudain du tunnel. La proc�dure de calcul/agrandissement est r�it�r�e jusqu'� obtenir la taille d�sir�e pour l'exp�rience. Il est alors possible d'acc�l�rer consid�rablement le temps de calcul n�cessaire � l'�tablissement du r�gime stationnaire.

        Etudions l'exemple d'une simulation dans un tunnel de dimension l=3h avec l'intention d'effectuer une ou des mesures br�ves lorsque le r�gime transitoire a disparu. Dans notre cas, la largeur minimum est de 256 noeuds, la hauteur est donc de 256/3=86 noeuds. Le temps de stabilisation du fluide est environ de 4*(256+86)=1400 pas de temps soit 1400*256*86=30 millions de calculs de sites. Ensuite, le tunnel peut �tre agrandi et contenir 512*172 noeuds. Le calcul red�marre alors pendant un temps n�cessaire pour que le fluide s'adapte au nouveau nombre de Reynolds, mais toutefois moins que le temps n�cessaire pour l'amor�age du fluide : environ 2*(512+172)=1400 pas de temps et 1400*512*172=123 millions de sites. Nous pouvons continuer ainsi � calculer et agrandir jusqu'� ce que la m�moire soit satur�e ou le nombre de Reynolds d�sir� soit atteint. Dans notre cas, avec une limitation � 64MO, nous obtenons le temps de calcul total suivant :

taille              pas de calcul           sites calcul�s
256*86=22016        4*(256+86)=1400          30M
512*172=88064       2*(512+172)=1400         123M
1024*344=352256     2*(1024+344)=2800        986M
2048*688=1409024    2*(2048+688)=5600        7890M
4096*1376=5636096   2*(4096+1376)=11200      63124M
8192*2752=22544384  2*(8192+2752)=22400      504994M

total : 577 Milliards de sites calcul�s

        L'agrandissement progressif permet de diviser par deux le temps d'amor�age du fluide : en commen�ant � la taille pr�vue au d�part (8192*2752) il faudrait calculer 10^12 sites au moins ! A la fin du calcul, les 22400 pas de temps sont largement suffisant pour effectuer les mesures d�sir�es, gr�ce � des moyennes dans le temps et dans l'espace.

        Il faut maintenant mentionner au moins trois limitations pratiques importantes, li�es � la complexit� des ph�nom�nes. Nous atteignons un nombre de Reynolds id�al d'environ 5000 et le type de programme pr�sent� jusqu'� pr�sent ne convient plus. La premi�re limite concerne le temps de calcul : si l'ordinateur calcule un million de sites par secondes (comme dans le deuxi�me code de r�f�rence) il faudra une semaine pour calculer les 577 Milliards de sites de l'exp�rience. Nous sommes encore loin du "temps r�el" et de l'interactivit� d�sir�s. Ensuite, les programmes actuels utilisent le PC sous MS-DOS en mode r�el, ce qui limite les exp�riences � 320*200 points, il n'est pas possible d'atteindre ainsi des nombres de Reynolds int�ressants et justifiant les efforts de programmation. Enfin, bien que le redimensionnement soit une �tape relativement simple � programmer, il devient plus compliqu� de d�finir les parois par de simples bitmaps ou par programmation au cas par cas. Il faut donc utiliser un mode de repr�sentation vectoriel qui sera r�interpr�t� � chaque fois que le tunnel change de taille. Les vecteurs sont simples � g�rer mais leur repr�sentation interm�diaire dans le code du projet actuel est tr�s complexe.

III.9 : Etude d'un cas r�el : benchmark de von Karman :

        Afin de pouvoir comparer les r�sultats et l'efficacit� de la m�thode employ�e, nous utiliserons les donn�es fournies par une �tude allemande [36], d�crivant les conditions de simulation d'all�es de von Karman. Puisque nos LGA simulent les fluides de mani�re explicitement temporelle en 2D, nous ne nous int�resserons qu'� la partie correspondante du benchmark qui comporte des exp�riences en 2D et en 3D, en statique ou en dynamique.

G�om�trie du tunnel et conditions aux limites : (extrait de [36])

        Avec un adimensionnement correct (si l'on tient compte de l'invariance galil�enne entre autre) nous avons besoin approximativement d'un site par millim�tre carr�, soit 2200*410=902000 sites ou environ 1 m�gaoctet. Les PC actuels disposent d'au moins 64 Mo actuellement et le calcul � 1Mc/s permet de soutenir un affichage � 1Hz dans ce cas. Dans des conditions id�ales (densit�, vitesse...), nous pouvons donc tr�s facilement atteindre le nombre de Reynolds requis par le benchmark. Il reste ensuite assez de marge pour corriger les artifacts du r�seau.

        Les all�es de von Karman ont �t� simul�es par plusieurs laboratoires, dont les laboratoires Fuji (Japon), l'Universit� Libre de Bruxelles, l'Universit� de Munich ou l'Observatoire de Nice pour citer quelques exemples. Dans la comparaison suivante, nous �tudions la diff�rence entre une r�solution temporelle d'�quations classiques et un calcul FHP-3 :

El�ments finis : (extrait de [36])	FHP-3: (ULB, David Hanon, r�seau de 800*200, Mach 0,45 et densit�=0,28)
	t=0   t=100   t=7000

        Dans des conditions appropri�es et avec l'adimensionnement correct, les r�sultats sont similaires. Toutefois, les gaz sur r�seaux permettent de faire appara�tre les fluctuations d�es au mouvement brownien. La simulation est donc bruit�e et n�cessite une int�gration temporelle et spatiale pour effectuer une mesure. Cependant, l'apparition de ph�nom�nes complexes est naturelle et n'a pas besoin d'�tre forc�e (ou si peu) : le cercle du benchmark est d�cal� d'un centim�tre dans le tunnel (1/41=2,5%) alors que dans l'exp�rience FHP il n'est d�cal� que d'un pixel (1/200=0,5%, mais cela est probablement involontaire et li� � l'algorithme de dessin du cercle). D'habitude, le d�centrage est n�cessaire pour forcer l'apparition rapide des tourbillons contrarotatifs car la m�thode de r�solution classique n'est pas bruit�e par nature. La recherche d'une solution exacte implique qu'une dissym�trie n'apparaitra qu'avec une erreur d'arrondi. En revanche, les gaz sur r�seaux comme FHP permettent au bruit brownien d'organiser les turbulences microscopiques (� l'�chelle de quelques sites) en turbulences macroscopiques (qui peuvent �tre �tudi�es par int�gration sur de nombreux sites) et de faire appara�tre les all�es caract�ristiques sans for�age ou dissim�trie artificiels. De plus, la complexit� de la g�om�trie du tunnel est arbitraire et n'influence pas le calcul.

        Le programme de David Hanon m�lange par codage les parois glissantes et rugueuses. En raison de l'imp�ratif du changement dynamique de la g�om�trie du tunnel, cette approche devient complexe et il faut pouvoir g�rer des parois de toute nature dans notre code final.

III.10 : Conclusion :

        Les deux codes de r�f�rence, ainsi que les nombreuses versions interm�diaires et les exp�riences qu'ils ont permis de conduire, montrent que le niveau de performance esp�r� ne peut �tre atteint qu'avec des techniques de codage et des algorithmes plus sophistiqu�s. L'approche par lookup table a atteint ses limites car il n'est pas possible de diminuer le nombre d'instructions par noeuds � cause des limitations et des contraintes architecturales (pas assez de registres, mode r�el ou unreal trop limit�s, pr�fixes divers et m�moire lente). L'apparition vers 1997 du jeu d'instructions MMX ainsi que la lecture plus attentive du livre de Michael Abrash [7] renforcent la conviction qu'il faut revoir le programme depuis le d�but.

IV.1 Introduction :

        Les premiers programmes de LGA sont relativement petits (le kernel comprend moins de 100 instructions) et commencent d�j� � n�cessiter une attention soutenue pour fonctionner correctement. Le temps n�cessaire � mettre un programme au point se mesure en semaines puis en mois pour la version en assembleur. Le projet d�crit ici est encore plus ambitieux et n�cessite une longue pr�paration pour que tous les �l�ments puissent fonctionner correctement, d'abord seuls puis en coop�ration avec les autres. Il faut concevoir chaque �l�ment ind�pendemment tout en pr�voyant leur assemblage final, il est donc n�cessaire d'avoir une vue g�n�rale du projet aussi claire que possible.

        Ce projet a �t� lanc� au d�part car il a �t� pr�par� pendant plusieurs ann�es : la plateforme est mieux maitris�e (en particulier le d�veloppement en assembleur MMX en mode prot�g�) et l'algorithme de strip mining est imagin�. Certains autres aspects algorithmiques et structurels sont imagin�s comme l'organisation des donn�es et la repr�sentation des parois. En th�orie, le projet ne devrait pas �tre compliqu� mais "le d�mon se cache toujours dans les d�tails"...

IV.2 : Intel : la plateforme id�ale malgr� elle

        Pour les nombreuses raisons expos�es pr�c�demment, la plateforme reste le PC sous MS-DOS. De plus, le d�veloppement dans un autre environnement et � ce niveau n�cessiterait l'apprentissage d'autres techniques et MS-DOS est le seul environnement permettant un contr�le total de la machine, condition n�cessaire pour trouver les goulots d'�tranglement des algorithmes test�s, sans interf�rences du syst�me d'exploitation ni d'�v�nements ext�rieurs incontr�lables. Enfin, bien que l'avenir commercial de ce syst�me soit incertain, des alternatives libres comme FreeDOS permettront de distribuer le programme avec le syst�me d'exploitation pr�configur�, afin d'en simplifier l'utilisation et la diffusion gratuite. Il n'est pas question d'utiliser d'ALPHA, de SPARC ou de PowerPC malgr� leur puissance sup�rieure. Un CRAY (T3E ou SV) conviendrait tr�s bien car les outils de profiling sont livr�s mais outre le prix et la place d�raisonnables, le but du projet n'est pas d'atteindre la plus haute performance possible mais bien de trouver des techniques afin de mieux exploiter une machine existante et facilement disponible. Le d�fi algorithmique est d�j� suffisamment complexe.

        Avec le temps, de nouveaux processeurs Intel et clones sont apparus : r�cemment l'Athlon d'AMD a surpass� le PIII d'Intel (voir "la Jihad des CPU de 7�me g�n�ration" par Paul Hsieh) et les compagnies se font une concurrence fr�n�tique, au niveau des prix, de la performance et des fonctionalit�s, pour la ma�trise du march�. Pour des raisons pratiques, nous devons nous concentrer sur une partie seulement des architectures disponibles, en esp�rant que les autres architectures ne soient pas compl�tement diff�rentes. La compatibilit� binaire entre les versions garantit cependant que le programme fonctionne partout o� les fonctions n�cessaires (souris s�rie, VESA2, MMX, MS-DOS et HIMEM.SYS en mode r�el) sont pr�sentes.

        Nous allons nous concentrer sur deux processeurs et leur architecture : le Pentium MMX � 200MHz, disponible � la maison, et le Pentium II dont plusieurs sont disponibles � l'universit�, dont deux en version biprocesseur. Ce sont deux architectures r�pandues et connues, mais suffisamment diff�rentes pour n�cessiter une �tude particuli�re. Le premier (P200MMX) est un processeur superscalaire � deux voies, proche d'un RISC classique, alors que le deuxi�me (PII) a un coeur OOO (Out Of Order execution) pouvant traiter 40 instructions � diff�rents stades de leur ex�cution, disposant de 80 registres renomm�s. L'interface avec la m�moire centrale et la "r�sistance � la charge de calcul" sont compl�tement diff�rentes et n�cessiteraient deux versions diff�rentes du programme : le P200MMX est un processeur "statique", qui reste assez pr�visible gr�ce � sa cache d'instructions, alors que le PII est enti�rement dynamique et non d�terministe ! Pour ne rien am�liorer, les r�gles de codage sont radicalement diff�rentes. Les sch�mas ci-dessous illustrent les diff�rences architecturales entre les deux types de processeurs :

Pentium MMX :	Pentium II :
* jusqu'� 2 instructions d�cod�es par cycle * 2 caches internes (donn�es et instructions) de 16 Ko * bus m�moire externe : Socket 7 � 66MHz (comme les Pentium � 100 MHz)	* jusqu'� 3 instructions d�cod�es par cycle et traduites en 6 micro-instructions (�ops). * 2 caches internes (donn�es et instructions) de 16 Ko et cache transactionnelle de 256Ko dans le module * bus m�moire externe : Slot 1 � 66MHz puis 100 et 133MHz, transactionnel

courtesy of Intel

        De nombreuses ressources sont disponibles, notamment dans la bibliographie et sur Internet. Michael Abrash [7] d�crit bien les aspects importants de la programmation du Pentium classique et l'adjonction du jeu d'instruction MMX est relativement simple :

        Le DOS-Extender (ou "loader 32 bits") est un morceau de code situ� au d�but d'un programme MS-DOS autonome et qui place le processeur en mode 32 bits : il permet ainsi d'utiliser les registres sur 32 bits sans utiliser de pr�fixe de taille, donc de faire fonctionner le programme plus vite. Son d�veloppement a commenc� vers 1997 et s'est stabilis� rapidement car il est court et con�u pour �tre tr�s compact : le code binaire occupe moins d'un kilo-octet. Le mode DPMI (DOS Protected Mode Interface) permettant d'utiliser le mode prot�g� 32 bits sous Windows a �t� essay� mais de nombreuses difficult�s ont emp�ch� de poursuivre l'effort dans cette direction : il a �t� plus simple de tout recoder � la main car cela �vite d'�tre d�pendant de conventions et d'interfaces complexes, contraignantes et inadapt�es. Le loader est minimal et ne remplit que les fonctions n�cessaires. Il a d'abord �t� programm� avec TASM mais le m�lange des codes 16 bits et 32 bits a rendu le codage tr�s difficile. Il a �t� port� tr�s facilement sous NASM et ce sera l'assembleur final : il est distribu� sous GPL, il est de plus en plus utilis�, il est stable et surtout sa syntaxe est tr�s simple et tr�s pratique, d�barrass�e de toute convention inutile. L'utilisation des directives USE16 et USE32 a fait dispara�tre de nombreuses difficult�s qui ont frein� le d�veloppement avec TASM, en particulier l'adjonction manuelle des pr�fixes de taille ("db 66h" et "db 67h"). Pour compl�ter, j'ai �crit un ensemble de fonctions permettant de se passer d'un linker et qui permet de g�n�rer l'ent�te d'un fichier EXE. Ces fonctions font actuellement partie de la biblioth�que de macros de NASM. Enfin, deux programmes (une sorte de pr�processeur et un patcheur d'ent�te) �crits en Turbo Pascal fournissent les derniers utilitaires pour programmer sans se pr�occuper de certains d�tails de bas niveau.

        Une fois le programme lanc�, nous disposons de l'environnement id�al : nous pouvons acc�der lin�airement � 64Mo (limite du standard XMS) en nous souciant beaucoup moins des segments qu'avant. Seules trois interruptions sont utilis�es : le clavier, la souris et la proc�dure d'erreur fatale (qui est d�clench�e par le processeur lorsqu'un bug se manifeste violemment). Ainsi, il n'y a pas un seul cycle que nous ne pouvons contr�ler et nous disposons de l'int�gralit� de la puissance de la machine, sans autre entrave que son architecture. Cet environnement est augment� de fonctions au fur et � mesure que leur int�gration devient possible. Par exemple, une fois que le contr�le de la carte vid�o en mode 1024*768*256/LFB est assur�, la souris est int�gr�e pour entrer des ordres. Presque toutes les fonctions sont test�es hors de l'environnement avant d'�tre int�gr�es, souvent en langage Turbo Pascal pour faciliter le d�bogage et la compr�hension des algorithmes. Par exemple, voici le programme qui est � la base du driver de souris :
 

Uses Crt, Dos;
{$F+}

const COM1INTR = $0C;
      COM1PORT = $3F8;

var bytenum : word;
    combytes : array[0..2] of byte;
    x, y : longint;
    button1, button2,changed : boolean;
    MouseHandler : procedure;

procedure MyMouseHandler; Interrupt;
var dx, dy : integer;
var inbyte : byte;
begin
  inbyte := Port[COM1PORT];  { Get the port byte }

  { Make sure we are properly "synched" }
  if (inbyte and 64) = 64 then bytenum := 0;

  { Store the byte and adjust bytenum }
  combytes[bytenum] := inbyte;
  inc(bytenum);

  { Have we received all 3 bytes? }
  if bytenum = 3 then
    begin

      { Yes, so process them }
      dx := (combytes[0] and 3) shl 6 + combytes[1];
      dy := (combytes[0] and 12) shl 4 + combytes[2];
      if dx >= 128 then dx := dx - 256;
      if dy >= 128 then dy := dy - 256;
      x := x + dx;
      y := y + dy;
      button1 := (combytes[0] And 32) <> 0;
      button2 := (combytes[0] And 16) <> 0;

      { And start on first byte again }
      bytenum := 0;
      changed := true;
    end;

  { Acknowledge the interrupt }
  Port[$20] := $20;
end;

var error: boolean;
begin
  ClrScr;

  error:=false;

  { Initialize the normal mouse handler }
  asm
    mov ax, 0
    int $33
    inc ax
    jne @the_end

    mov ax,24h
    int 33h     { get mouse parameters }
    inc ax
    je @the_end  { no mouse or other error }
    cmp ch,2
    jne @the_end { we want a serial mouse }
    cmp cl,2   { COM port }
    jae @not_the_end
@the_end:
    mov word ptr [error], -1
@not_the_end:

  end;

  if error then halt(1);

  { Initialize some of the variables we'll be using }
  bytenum := 0;
  x := 0;
  y := 0;
  button1 := false;
  button2 := false;

  { Save the current mouse handler and set up our own }
  GetIntVec(COM1INTR, @MouseHandler);
  SetIntVec(COM1INTR, Addr(MyMouseHandler));


  while not keypressed do
    if changed then
     begin
       WriteLn(x : 5, y : 5, button1 : 7, button2 : 7);
       changed:=false;
     end;

  SetIntVec(COM1INTR, @MouseHandler);
end.

        Une fois l'algorithme compris, le code est traduit en langage assembleur et retest� avec Turbo Pascal en inline. Ensuite, la partie d�velopp�e et valid�e peut �tre recopi�e dans l'environnement en assembleur, apr�s une l�g�re traduction de la syntaxe (Turbo Pascal -> NASM) et le renommage des variables. C'est � cet endoit que surviennent la plupart des probl�mes : faute de frappe, erreur de copier/coller, syntaxe erronn�e mais non d�tect�e par NASM ... Les probl�mes sont facilement �limin�s avec beaucoup de pratique, de m�thode et de patience et les cas simples permettent de se pr�parer aux cas tr�s difficiles. Par exemple, certains bugs ont n�cessit� deux semaines pour �tre r�solus alors qu'un module peut �tre int�gr� en quelques heures : la programmation est un exercice non lin�aire ! Le type de code pr�c�dent n'a pas �t� optimis� tr�s fortement : son ex�cution n'est pas absolument critique et l'acc�s aux entr�es-sorties prend la plupart des cycles. Il a donc �t� recod� pour �tre le plus compact possible et effectuer seulement les op�rations n�cessaires, gr�ce � de nombreux param�tres pr�-cod�s ou d�finis une fois pour toute (comme le curseur de la souris). A ce niveau, cela suffit pour atteindre l'objectif de performance d�sir�.

        Lorsque le code le permettait, un stub a �t� int�gr� pour supporter les plateformes bi-processeur. Cela a demand� trois jours non-stop de travail pour lire les documents d'Intel mais le r�sultat est simple, compact et bien cern�. La communication s'effectue simplement avec des s�maphores en m�moire. Le syst�me de coh�rence de caches MESI a �t� mis � l'�preuve et des mesures ont montr� que cette plateforme n'est int�ressante que pour des applications CPU bound (les deux processeurs � 266 MHz doivent se partager un unique bus similaire � celui du Pentium 100). Toutefois, apr�s quelques mois d'utilisation, le stub a �t� d�valid� afin de pouvoir se concentrer sur la partie la plus importante : le calcul. En effet, il faut programmer une version mono- et bi-processeur pour chaque version du kernel de calcul et cela ralentit le d�veloppement. Le code bi-processeur reste toujours disponible et la structure du programme est pr�serv�e pour permettre son retour lorsque le code sera abouti et parall�lisable.

IV.3 : Description de l'algorithme de strip mining

        Comme nous l'avons constat� avec les programmes de la partie III, la bande passante avec la m�moire est une "ressource" tr�s importante car elle peut faire ralentir le programme d'une mani�re dramatique. Le deuxi�me code de r�f�rence ne r�soud pas tout � fait le probl�me car la vitesse de calcul est inversement proportionnelle � la vitesse de l'affichage et il faut donc choisir entre interactivit� et temps de calcul. Le choix de l'algorithme est mieux adapt� au Pentium que pour le premier exemple destin� au i286 mais l'affichage hors de la boucle de calcul est une erreur strat�gique car elle g�n�re autant de cache misses que le calcul. L'affichage entrelac� avec le calcul r�duirait le nombre de cache misses, am�liorerait la vitesse de rendu � l'�cran et b�n�ficierait naturellement de la localit� spatiale et temporelle.

Cependant, la situation va encore se d�t�riorer dans l'avenir car malgr� les lourds investissements r�alis�s dans ce domaine par les industriels, la bande passante vers la m�moire centrale s'accro�t plus lentement que la vitesse de traitement des processeurs. Si notre algorithme est bas� sur un balayage lin�aire de toute la m�moire, l'acc�l�ration r�alis�e en utilisant une plateforme "plus rapide" ne sera pas proportionnelle � l'augmentation de la vitesse du processeur. L'entrelacement de l'affichage avec le calcul, n�cessaire pour r�duire le nombre de transactions sur le bus m�moire, n'est pas suffisant pour permettre � l'algorithme d'exploiter les plateformes modernes et celles qui leur succ�deront : nous serons toujours limit�s par la vitesse d'acc�s � la m�moire centrale, acc�d�e en lecture et en �criture.

extrait de la page d'accueil du benchmark STREAM

        Pour traiter un tableau de N sites, il faudra donc acc�der � N*2 sites ainsi qu'� la m�moire vid�o ce qui sature tragiquement le bus m�moire. Si le probl�me est li� au balayage lin�aire, il faut alors y renoncer. Il faut balayer d'une mani�re qui permette � la m�moire cache d'�tre utilis�e de mani�re optimale, c'est � dire : extraire la localit� spatiale et temporelle du mod�le. Nous devons pour cela analyser la d�pendance des donn�es � chaque pas de temps entre une cellule et ses voisins :

        Le probl�me est simplifi� si l'on s'inspire d'une technique qui traite des lignes au lieu de cellules individuelles : la d�pendance des donn�es se confond alors dans une croissance lin�aire, non hexagonale, du domaine � traiter. La technique de strip mining est connue et principalement utilis�e sur des plateformes sophistiqu�es, multi-processeurs, avec plusieurs niveaux de m�moire, c'est � dire o� l'acc�s � la m�moire n'est pas uniforme, comme dans notre cas. Le principe de cette technique est de diviser le domaine calcul� en bandes (strips) et de les balayer de mani�re non lin�aire. Par exemple, en utilisant des propri�t�s math�matiques, il est possible de r�organiser les donn�es d'une matrice afin de calculer son inverse : le programme de strip mining ne balaiera plus les donn�es de mani�re triviale et les donn�es resteront plus longtemps en m�moire cache, ce qui acc�l�re le calcul. Dans notre cas, ce ne sont pas des propri�t�s math�matiques mais g�om�triques qu'il faut �tudier.

        Nos strips correspondront aux lignes du tableau : elles seront calcul�es par une boucle simple, en distingant toutefois les lignes paires et impaires (comme pour le premier code de r�f�rence). Les lignes sont ensuite balay�es � un niveau sup�rieur avec une "fen�tre glissante" : c'est un balayage similaire � une double boucle imbriqu�e, portant sur l'axe y. L'algorithme de balayage consiste donc en trois boucles les unes dans les autres :

for y:=0 to y_max do (* balayage g�n�ral : g�n�re les cache misses *) begin for i:=0 to strip do (* boucle interne de strip mining, sans inclure l'amro�age *) for x:=0 to x_max do (* balayage d'un strip *) calcule (x,y+i); (* le d�bordement en y n'est pas trait� dans ce pseudo-code *) affiche_ligne(y); end;

        Nous pouvons voir que cette adaptation convient aux n�cessit�s de l'affichage, sans que la cache ne soit vid�e : le transfert des donn�es vers l'ext�rieur s'effectue avec la garantie que les donn�es soient d�j� en cache si le param�tre strip est bien choisi. La bande passante du bus externe du processeur est donc divis�e en 3 : lecture, �criture des �l�ments modifi�s, affichage. Cependant les �l�ments mis en jeu sont tr�s complexes : ce balayage compos� implique que l'on peut trouver sur le tableau, lors du calcul, des sites dont les valeurs correspondent � des pas de temps �loign�s, il faut un m�canisme adapt� pour g�rer les buffers n�cessaires. Le sch�ma suivant montre les d�pendances sur l'axe y � tous les pas interm�diaires du calcul d'un tableau de 7 lignes et avec un strip mining de 3 lignes :

Le sens de balayage global est sur l'axe y, le sens de balayage de la boucle interne est symbolis� par les larges fl�ches transparentes. Les cases repr�sentent diff�rentes valeurs pour chaque ligne � diff�rentes �tapes, elles sont num�rot�es dans l'ordre d'appel de la fonction de calcul. Les petites fl�ches montrent les d�pendances de donn�es entre les cellules.

        Nous pouvons apercevoir que le cas de buffer temporaire expliqu� en partie III.3 est un cas particulier de ce nouvel algorithme lorsque strip est �gal � 1. Cela signifie que le plan temporaire dans ce cas particulier ne change pas. Par contre, lorsque strip augmente, la gestion du plan temporaire devient beaucoup plus complexe : il faut strip lignes de buffer et autant de pas de temps diff�rents sont pr�sents sur le tableau. Le d�placement des donn�es dans le tableau demande alors encore plus de soin lors du d�veloppement et n�cessite une pr�paration importante. La r�alisation du programme a �t� plus difficile que pr�vu, principalement car l'algorithme de balayage initialement utilis� comportait des d�pendances spatio-temporelles cach�es qui �taient impossibles � r�soudre sans une r�vision compl�te du programme.

        Dans le programme final, tous les indices sont transform�s en pointeurs et lin�aris�s afin de r�duire � la fois le nombre n�cessaire de registres et le nombres d'instructions de calcul. La structure de boucle pr�sent�e initialement ne convient pas pour g�rer l'amor�age du buffer de strip mining : nous avons vu que le balayage doit commencer et se terminer en augmentant puis en r�duisant progressivement le nombre de strips. La fen�tre de calcul est alors cern�e par deux pointeurs, modifi�s par de simples op�rations de comparaison/assignement. La boucle est ainsi contr�l�e par des conditions complexes mais le nombre de param�tres est r�duit � deux variables : "d�but" et "fin". "fin" est "poursuivi" par "d�but", la chasse se termine lorque "d�but" d�passe la "fin" (c'est � dire : la fin du tunnel, mais pour r�duire les d�pendances entre les registres, ce param�tre est confondu avec une valeur imm�diate connue). D'autres pointeurs sont aussi n�cessaires pour g�rer les tableaux temporaires.

const strip=2;
var buf:array[0..strip-1,0..159,0..1] of byte; (* buffer temporaire *)
var j,debut,fin,debut_buf:word;
label calcul;
begin

 while not keypressed do
 begin
  debut:=1; (* initialisation des pointeurs *)
  fin:=1;
  debut_buf:=0;

calcul:
  for j:=fin downto debut do
  begin
   if (j and 1)=1 then (* parit� de la ligne *)
(* calcule une ligne impaire *)
   else
(* calcule une ligne paire *)
  end;

  if fin=199 then
   inc(debut_buf)
  else
   inc (fin);

  if fin<strip+1 then
   goto calcul;
(* affichage ligne ici *)
  inc(debut);
  if debut<=199 then
   goto calcul;
 end;

        Les m�saventures d�es au strip mining sont nombreuses mais la plus importante met en jeu la structure compl�te du programme. Nous avons vu dans les pseudo-codes pr�c�dents que nous pouvions profiter de l'algorithme pour y inclure la fonction d'affichage. Le programme r�alis� va plus loin en essayant d'utiliser la structure du buffer temporaire pour m�moriser des informations sur la fen�tre courante, en particulier afin d'effectuer une accumulation des particules et afficher la densit� dans le tunnel. L'avantage �vident est que cela �conomise beaucoup de m�moire : le buffer d'affichage n�cessite alors strip lignes au lieu de y et cela nous ram�ne au cas d�crit pour le buffer temporaire simple (au total : y+1 lignes au lieu de 2y). Pourtant la r�alisation fait appara�tre que les donn�es ne sont pas acc�d�es dans le m�me ordre que pour le buffer temporaire de d�placement : l'affichage n�cessite un pointeur suppl�mentaire ainsi qu'un algortihme nouveau, m�me si la place occup�e n'est pas un probl�me. Tous les registres utilisables sont d�j� allou�s et l'accumulation de donn�es sur la fen�tre n'est pas possible. Seules les donn�es disponibles � un intervale de strip pas de temps peuvent �tre affich�es. Aussi, l'affichage sophistiqu� pr�vu au d�part n'a pas �t� inclus dans le programme m�me s'il est d�j� con�u et test�. A posteriori, la technique simple (consommant beaucoup de m�moire) aurait �t� pr�f�rable mais la structure du code existant est trop avanc�e et fig�e pour pouvoir �tre modifi�e. Un recodage complet est n�cessaire et serait plus rapide.

        Le speedup permis par l'algorithme de strip mining d�pend du nombre de lignes contenues dans la fen�tre et sa taille d�pend de nombreux param�tres. Afin de maximiser la r�utilisation des donn�es, il faut donc que la fen�tre corresponde le mieux possible aux caract�ristiques de la m�moire de la machine, ce qui est absolument ind�terministe et impr�visible. La premi�re mesure � prendre est d'orienter le tableau dans la hauteur afin que la fen�tre soit la plus �troite possible. Par exemple, dans le cas des all�es de von Karman, il faut tourner le tunnel � 90 degr�s (hauteur plus grande que la largeur) afin que l'algorithme soit plus efficace. En effet, le speedup est proportionnel � strip, ce que confirment les mesures. Pour le cas d'une version multiprocesseur, cela veut dire que le d�coupage du domaine doit se faire dans la largeur afin de r�duire la bande passante du bus externe, m�me si cela augmente le nombre de s�maphores de synchronisation.

        Ensuite, il faut que le nombre de lignes dans la fen�tre ne soit ni trop grand ni trop petit mais il n'y a pas de formule g�n�rale permettant de d�terminer ce param�tre. La solution la plus efficace et la plus simple est de mesurer en conditions r�elles toutes les possibilit�s : c'est l'�tape de calibration qui calcule plusieurs pas de temps et mesure le temps pour chaque largeur de la fen�tre. A la fin des mesures, il d�termine le meilleur score (le temps divis� par la taille de la fen�tre) : le r�sultat, en nombre de lignes, servira pour les calculs futurs. La calibration doit �tre effectu�e � chaque fois que les param�tres de simulation changent (nombre de sites ou g�om�trie du tunnel) afin que le temps de calcul soit toujours optimal. Pour des raisons historiques et pratiques, la recherche exhaustive est limit�e � 32 valeurs de strip mais il n'est pas n�cessaire d'augmenter ce nombre et cela suffit pour de donner une bonne id�e sur l'architecture de la m�moire de la machine. Par des mesures, calculs et d�ductions, il est possible de d�terminer la bande passante r�elle de la m�moire centrale ou la taille de la m�moire cache. Les processeurs les plus r�cents (Pentium II par exemple) ont une excellente r�ponse au strip mining, utilisant souvent un nombre maximal de lignes dans la fen�tre, alors que les processeurs de cinqui�me g�n�ration avec cache L2 externe ont un score faible correspondant � la taille de la m�moire cache L1.

        A premi�re vue, le strip mining est une technique complexe dont l'efficacit� n'est pas garantie : le gain en performance est d�pendant de l'architecture de la machine et sa complexit� pratique est peu �vidente � appr�hender. Toutefois c'est un algorithme adaptatif qui ne surcharge pas la machine lorsque c'est inutile et son efficacit� est d�montr�e en pratique pour les plateformes actuelles. Nous allons voir � la fin de cette partie une application r�elle o� le speedup approche 2,5.

        Dans le cas qui nous concerne, le strip mining utilis� ici est une version tr�s sp�ciale adapt�e aux tableaux en deux dimensions. Son extension � d'autres dimensions est envisageable mais sa complexit� intrins�que doit �tre compl�tement maitris�e : les m�saventures avec la version actuelle du code (en particulier les croisements de pointeurs pour l'affichage) montrent que ce domaine particulier devrait �tre �tudi� plus en profondeur de mani�re "acad�mique" pour �tre connu et utilis� largement. Le manque de connaissances pr�liminaires a rendu la programmation difficile, malgr� la longue pr�paration. Esp�rons que dans le futur, cet algorithme sera mieux �tudi� et se banalisera car c'est une solution simple et non calculatoire (ind�pendante du "Grand O") pour acc�l�rer les calculs lourds sur les machines actuelles et futures.

IV.4 : Nouvelle structure des donn�es et m�thode de calcul

        Nous avons vu que m�me en groupant les octets par 4, il faut toujours 47 instructions de mouvement des donn�es contre 12 pour la consultation de la table : le premier programme n�cessitait 30 instructions de mouvement pour 6 de consultation, nous sommes pass�s de 5:1 � 4:1 avec une augmentation de la taille du code de 1,6. Il est clair que pour acc�l�rer le programme, il faut r�duire le temps pass� � bouger les bits. En ignorant les autres param�tres, il faudrait diviser par quatre le nombre d'instructions de mouvement de donn�es pour doubler la vitesse de calcul. Malheureusement nous ne pouvons pas agir sur le jeu d'instructions, nous devons jouer sur l'organisation des donn�es pour r�duire cet �norme probl�me.

        Si le programme passe les 4/5�mes de son temps � bouger les bits, c'est parcequ'il doit extraire puis ins�rer des bits dans d'autres mots (le processeur ne dispose pas d'instruction sp�ciale pour cette t�che), ceci sept fois par site, pour tous les sites. Nous devons donc utiliser une repr�sentation des donn�es qui r�duit au maximum le nombre d'instructions n�cessaires au mouvement. La structure la plus simple est �videmment le mot long, permettant en MMX de coder 64 particules allant dans la m�me direction.

        Le calcul est un nouveau probl�me : il faut l'effectuer par op�rations logiques explicites au lieu de consulter des tables. Nous ne disposons pas d'instruction ad�quate permettant de transformer 7 mots de 64 bits en 64 octets. Le probl�me est toutefois plus compliqu� car l'enjeu est de pouvoir suivre la croissance de la puissance des ordinateur et cette croissance passe par l'augmentation de la largeur des mots : les programmes consultant des tables sont donc condamn�s � stagner au niveau de la performance car le nombre de sites calcul�s restera proportionnel au nombre de consultations. Au contraire, en utilisant une technique bas�e sur la repr�sentation parall�le des sites, une augmentation de la largeur des mots augmentera la vitesse de calcul. Une estimation situe le "point de rupture" � des mots de 32 bits : cette technique n'est pas int�ressante pour un processeur 16 bits mais cela d�pend aussi beaucoup du nombre de registres, du jeu d'instructions et du nombre d'instructions d�cod�es par cycle. L'introduction des instructions MMX et les registres associ�s favorisent beaucoup l'utilisation du mod�le � calcul explicite, malgr� la complexit� inh�rente du programme.

        L'organisation des donn�es au niveau macroscopique joue aussi un r�le important. Traditionnellement, la technique "multi-spin" (en r�f�rence au mod�le d'Ising) regroupe les sept directions dans sept tableaux car les ordinateurs vectoriels (CRAY surtout) peuvent traiter 4096 sites en une seule instruction (8 registres vectoriels de 64 fois 64 bits). Ce type d'organisation perdure encore dans les habitudes mais se heurte comme les autres algorithmes aux architectures complexes des microprocesseurs r�cents. Par exemple, il faut 7 pointeurs pour d�signer le tableau d'origine et 7 pointeurs pour le tableau de destination (en utilisant l'organisation simple, sans buffer temporaire) et les CRAY (comme la plupart des autres processeurs courants) n'ont que 8 registres d'adresse. Ensuite, la m�moire cache, ajout�e au processeur pour pallier au manque de bande passante avec la m�moire centrale, cr�e des probl�mes d'associativit� d'ensemble : pour le Pentium par exemple (8Ko associatif � 2 voies) des acc�s cons�cutifs � des adresses identiques modulo 4096 provoquent de nombreux cache misses artificiels. Nous voyons encore ici que les algorithmes, m�me si leur programmation semble facile, survivent mal aux �volutions de la plateforme.

        L'organisation que nous allons utiliser maintenant est un compromis entre "multi-spin" et "multi-site" : la premi�re contrainte est de diminuer au maximum le nombre de registres pointeurs n�cessaires. Il faut donc "lin�ariser" les acc�s et limiter les modes d'adressage au mode registre + index imm�diat. Pour r�duire encore plus le nombre de registres inutilement utilis�s dans le kernel, l'index sera parfois modifi� � la vol�e par du code automodifiant, par exemple pour l'acc�s � des cellules d'une ligne diff�rente. Le code automodifiant est �videmment plac� hors de la boucle et n'est utilis� que lorsque la taille du tableau change, pour �viter les lourdes p�nalit�s � la premi�re ex�cution.

/* quelques d�finitions : */

#define max_x 512
#define max_y 512
#define lli long long int
#define slli (sizeof(long long int))

/* d�claration en C d'un tableau multi-site : */

char tableau[max_y][max_x];
(comme pour les programmes �tudi�s jusqu'ici)

/* d�claration en C d'un tableau multi-spin : */

lli a[max_y][max_y/slli];
lli b[max_y][max_y/slli];
lli c[max_y][max_y/slli];
lli d[max_y][max_y/slli];
lli e[max_y][max_y/slli];
lli f[max_y][max_y/slli];
lli g[max_y][max_y/slli];
(pour les ordinateurs vectoriels)

/* d�claration en C d'un tableau composite : */

struct site {
  lli a,b,c,d,e,f,g;
} tableau[max_y][max_y/slli];

        L'annexe D pr�sente diff�rentes variations des structures de donn�es qui ont �t� programm�es. Comme ces logiciels sont �crits en langage de haut niveau, les contraintes de performance sont difficilement visibles au niveau du code source. De plus, la complexit� croissante des plateformes n'est pas le souci majeur des programmeurs dont le seul but est d'obtenir un r�sultat valide. La structure composite n'est donc pas encore r�pandue dans les codes FHP et la d�nomination officielle n'existe pas encore, cette organisation peut aussi bien �tre appel�e multispin entrelac�.

IV.5 : Impl�mentation des parois

        La nouvelle organisation des donn�es permet � la fois d'�conomiser de la bande passante vers la m�moire et de cr�er des parois plus flexibles mais, comme dans le reste du projet, aux d�pends de la complexit� du code et du temps de d�veloppement.

        Tout d'abord, il faut remarquer la limitation inh�rente aux parois "rugueuses" telles que programm�es dans les programmes pr�c�dents o� un seul bit �tait utilis�. D'un c�t�, les parois sont tr�s faciles � g�rer car il suffit de mettre un bit � 1 au bon endroit pour cr�er un mur. De l'autre c�t�, moins d'un de ces bits sur cent est utilis� en pratique : nous pouvons donc gagner facilement un huiti�me de bande passante en ne transmettant pas ce bit, car maintenant il est d�coupl� du domaine des particules.

Ensuite, le mod�le � un bit est limit� par le type de parois qu'il peut impl�menter : nous avons vu que les parois rugueuses sont n�cessaires mais insuffisantes dans les cas g�n�raux. Pourtant, le mod�le FHP permet une plus grande vari�t� d'orientation des parois. M�me si pour un site, une particule peut repartir dans 5 directions diff�rentes lors d'un choc avec une paroi, le nombre r�el est de 12 directions de parois, c'est � dire une pr�cision de 30 degr�s. Nous avons donc int�r�t � r�utiliser la bande passante gagn�e pour trouver une autre repr�sentation des parois, qui permettrait par exemple de repr�senter une sph�re lisse.

        La structure et l'algorithme n�cessaires � l'impl�mentation des parois est assez simple pour la partie calcul mais cette facilit� est largement compens�e par la complexit� n�cessaire � la fabrication des donn�es repr�sentant les parois. Il suffit d'environ 20 instructions pour d�vier jusqu'� 128 particules de leur trajectoire gr�ce � un �change contr�l� par un masque et deux pointeurs log�s dans une structure de type mur. Ces pointeurs et ce masque n�cessitent cependant un m�canisme sophistiqu� pour �tre mis en place. Pour le programme de calcul, les parois sont "interpr�t�es" gr�ce � une sorte de liste chain�e en m�moire. Son adresse de d�part est fournie par un pointeur de 32 bits log� dans la structure de type site qui devient donc :

/* d�claration en C d'un �l�ment de liste de modification : */
typdef struct {
  long int cpt; /* nombre d'�l�ments � interpr�ter */
  short int offset1, offset2; /* pointeurs relatifs au site courant */
  lli masque;
} mur;

/* d�claration en C d'un tableau composite : */
struct site {
  lli a,b,c,d,e,f,g;
  mur * liste;
  long int pad; /* pour aligner la structure sur 8*64 bits */
} tableau[max_y][max_y/slli];

Dans la pratique, seuls les murs verticaux et horizontaux sont utilis�s : les parois obliques ou non rectilignes posent des probl�mes � cause de la parit� des lignes. En effet, cela d�cale les parois d'un quart d'unit� d'un c�t� ou de l'autre et l'algorithme de Bresenham ou de diff�rences finies num�rique ne peuvent �tre utilis�s directement. Toutefois ce n'est pas le probl�me le plus important : la programmation de l'algorithme de g�n�ration automatique des listes de modifications est beaucoup plus complexe que lorsqu'on g�n�re les listes � la main comme lors des premiers essais. Les parois sont m�moris�es sous forme de coordonn�es vectorielles et leur transformation en listes de modifications n'est pas �vidente. Des difficult�s suppl�mentaires concernent l'efficacit� de la transformation (format vectoriel->listes de modification, qui doit donner un r�sultat le plus compact possible) et les croisements entre des parois d'inclinaison diff�rentes. Ces deux probl�mes ont �t� trait�s mais celui de l'efficacit� n'est pas consid�r� comme r�solu, en l'absence de preuves formelles d'optimalit�. Quant au croisement de parois aux orientantions diff�rentes, le probl�me est r�solu pour l'instant en fabricant un "point chaud" : il faut cr�er un mur ponctuel aux parois non glissantes, sinon des particules peuvent disparaitre dans l'interstice situ� � l'intersection des parois.

        Au fur et � mesure des appels, la deuxi�me partie va construire un r�seau de structures en essayant de r�utiliser au maximum celles d�j� construites, tout en �vitant que des intersections ne laissent �chapper des particules. C'est un algorithme qui fonctionne bien, mais dont l'activit� est difficile � contr�ler pour v�rifier son efficacit�. De plus, l'algorithme en Pascal prend beaucoup de place et a �t� difficile � traduire en assembleur car il est difficile justement de "voir" ce que l'algorithme fait, en plus des habituelles erreurs de fatigue et de frappe. Toutefois, son fonctionnement satisfaisant prouve que le programme est possible et il a permis de d�gager les contraintes inh�rentes au mod�le. Un post-processeur, balayant les structures cr��es, devrait permettre de les comprimer encore plus en �liminant celles qui ne seraient plus utilis�es : l'algorithme actuel n'est pas parfait et souffre de memory leaks. Une �tude � un plus haut niveau et l'utilisation de techniques plus efficaces (�tude et allocation globale des ressources) devraient permettre de g�n�rer des listes optimales et ainsi de consommer moins de m�moire cache. Une autre possibilit� serait de compiler les structures de modifications au lieu d'interpr�ter des listes

        L'objectif d'agrandissement arbitraire du tableau est preque atteint car les parois peuvent �tre compl�tement redimensionn�es dynamiquement. Il reste par contre � trouver un bon algorithme pour le redimensionnement des masses de particules, mais cela ne sera n�cessaire que lorsqu'il sera possible de cr�er et de d�truire des particules afin de g�n�rer du vent dans le tableau. Un algorithme similaire � celui des listes de modifications serait probablement utilis�.

IV.6 Algorithme d'affichage

        Non seulement le calcul des collisions est plus difficile qu'avec l'organisation multisite, mais l'affichage pose des probl�mes que l'on pourrait qualifier d'int�ressants. Rappelons que l'affichage s'effectue avec un octet par site, la couleur �tant contr�l�e par une palette effectuant la traduction dans le DAC de la carte vid�o.

        Dans les programmes pr�c�dents, le site est tout simplement envoy� vers la carte vid�o et il suffit de modifier la palette pour faire apparaitre les �l�ments d�sir�s, comme la densit� ou la direction des particules. L'affichage est direct et naturel. Maintenant, il faut transformer les plans de bits en octets et le m�me type de probl�me qu'avant se pose : il faut r�duire le nombre total d'instructions ainsi que le nombre d'instructions par site. Pourtant les vieilles habitudes de codage sont tenaces : si le programme devait �tre cod� dans un langage de haut niveau, sans prendre le temps de r�fl�chir, son efficacit� serait faible. Par exemple, le code suivant en C pour effectuer la translation a une complexit� ("grand O") proportionnelle au nombre de sites et ne profite pas de la largeur croissante des registres :

struct site {
  lli a,b,c,d,e,f,g;
  mur * liste;
  graph *p;
} *s;
unsigned char c, *d=video;
int i;

/* ce pseudo-code transforme un site point� par s
 en 64 octets pour afficher la densit� : */
for (i=0; i<64; i++) {
  c= ((s->a >>i)&1) + ((s->b >>i)&1) + ((s->c >>i)&1)+ ((s->d >>i)&1)
   + ((s->e >>i)&1) + ((s->f >>i)&1) + ((s->g >>i)&1);
  *(video++)=c;
}

        Le premier gros d�faut de ce code est l'acc�s par octet � la m�moire vid�o, nous avons pourtant vu que cela ralentit consid�rablement le code qui souffrira au moins de 64 lourdes p�nalit�s. Ensuite, ce code est compl�tement lin�aire, il traite chaque bit de chaque site ind�pendamment, comme lors du d�placement des bits lors du calcul en multisite. Cet algorithme "simple" n'exploite pas le parall�lisme permis par les mots tr�s larges. De plus, nous verrons dans le prochain chapitre que le code de calcul doit �tre pr�c�d� d'un code d'accumulation qui transforme les 7 bits d'entr�e en 3 bits sous forme de repr�sentation binaire classique (voir aussi le sch�ma principal et l'annexe C, page 33, pour le code). La r�utilisation de ces donn�es temporaires r�duit le corps de la boucle interne de moiti� :

for (i=0; i<64; i++) {
  c= ((s0 >>i)&1) + (((s1 >>i)<<1)&1) + (((s2 >>i)<<2)&1);
  *(video++)=c;
}

        L'annexe C, page 40, traite de l'accumulation lors du calcul puis de l'affichage. Une mani�re plus sophistiqu�e pour effectuer ces op�rations est d'utiliser la technique d'explosion, d�couverte lors de la programmation de la routine d'affichage du curseur. En utilisant les instruction PUNPCK d'une mani�re particuli�re, il est possible de transformer 8 bits cons�cutifs en 8 octets aux valeurs choisies par des masques.

movq mm1, Sx ; S0, S1 ou S2
punpcklbw mm1,mm1
punpckldw mm1,mm1
punpckldq mm1,mm1
;ici les 8 LSB sont recopi�s 8 fois. il faut donc trier ici les bits :
pand mm1,p1 ; p1 pointe vers la donn�e 0x8040201008040201
pcmpeqb mm1,[=0, reg ou mem] ; mm1=0 ou 0FFh

; pour le cas de S0 : astuce ! x-(-1)=x+1
psubb mm0,mm1 ; mm0 est la valeur accumul�e

; pour le cas de S1 et S2: 
pand mm1,p2 ; p2 pointe vers la donn�e 0x0202020202020202
                ; ou 0x040404040x04040404, selon S1 ou S2
paddb mm0,mm1 ; et voila !

        Ce code serait donc d�s�quilibr� : il utilise trop la CPU lors du calcul et perd du temps lors de l'affichage. Pourtant l'algorithme de strip mining permet de n'afficher qu'une fraction des donn�es, en fonction du param�tre strip. Il faut donc r�duire au maximum le nombre d'op�rations n�cessaires � l'accumulation des donn�es, m�me si l'affichage est plus complexe : d'une part le strip mining compense la perte lorsque strip est �lev�, d'autre part l'entrelacement d'instructions utiles lors du transfert permet de recueillir des statistiques permettant par exemple d'optimiser le contraste de la palette.

        L'algorithme choisi pour l'accumulation est en fait tr�s simple et il utilise au maximum les propri�t�s des donn�es : il consiste � additionner 8*3 bits � 8 octets en parall�le, en masquant directement, sans "�clatement" ou r�organisation, les donn�es. Le calcul de l'accumulation prend donc beaucoup moins de temps et la r�organisation des donn�es est effectu�e lors du transfert vers la m�moire vid�o.

struct site {
  lli a,b,c,d,e,f,g;
  mur * liste;
  graph *p;
  lli v[8];
} *s;
lli mask=0x0101010101010101;

/* pseudo-code d'accumulation : */
for (i=0; i<8; i++)
  s->v[i]+=((S0>>i)&mask)
          |(((S1>>i)<<1)&mask))
          |(((S2>>i)<<2)&mask));

        Maintenant que les octets sont accumul�s, il faut les afficher quand arrive la fin de la fen�tre : le gros probl�me est que la fonction d'accumulation d�crite pr�c�demment m�lange compl�tement les octets. Nous apercevons cependant qu'ils sont ordonn�s d'une mani�re qui ressemble � celle de la sortie du butterfly d'un code de FFT :

                                      v[0]:  0  8 16 24 32 40 48 56
                                      v[1]:  1  9 17 25 33 41 49 57
                                      v[2]:  2 10 18 26 34 42 50 58
                                      v[3]:  3 11 19 27 35 43 51 59
ordre des octets apr�s accumulation:  v[4]:  4 12 20 28 36 44 52 60
                                      v[5]:  5 13 21 29 37 45 53 61
                                      v[6]:  6 14 22 30 38 46 54 62
                                      v[7]:  7 15 23 31 39 47 55 63
 

        Un code a aussi �t� con�u pour afficher la "densit� de r�organisation" (les collisions donnant lieu � une r�organisation des sorties), � d�faut de pouvoir calculer la chiralit� ponctuelle instantann�e. La s�lection entre les deux versions (densit� de particules ou densit� de collisions r�organis�es) sera effectu�e par un code automodifiant (SMC : Self-Modifying Code).

        Le sch�ma suivant r�sume les algorithmes et la structure des donn�es li�es � l'affichage des densit�s :

IV.7 : Analyse bool�enne de l'op�rateur de collision

        La r�alisation du projet de maitrise a souffert de plusieurs al�as dont le plus important est celui du code de collision. Le projet avait �t� lanc� avec la conviction que le travail serait facile gr�ce aux "outils modernes" mais en r�alit�, seule la vieille m�thode du papier et du crayon a permis d'arriver � bout du probl�me.

La table ci-contre pr�sente les collisions du mod�le FHP-3 classique. Il contient � la fois la valeur num�rique de l'entr�e et des sorties, afin de pouvoir programmer une lookup table et les vecteurs vitesses sont repr�sent�s pour comprendre le type de collision qui a lieu.         Ce mod�le utilise un seul g�n�rateur binaire de nombres al�atoires, d'une probabilit� 1/2. Selon le bit fourni � chaque site par ce g�n�rateur, on choisit la colonne de gauche ou de droite. La colonne de gauche est d�termin�e par la formule (1) et les valeurs de la colonne de droite sont donn�es par la formule (2).

        Les formules qui permettent de calculer les collisions sont extraites de [32] et reproduites de la th�se de James Buick qui a aussi indiqu� des corrections. Elles sont �nonc�es de la mani�re suivante :

Nous pouvons alors d�finir les variables temporaires suivantes :

o� , a + b, et correspondent respectivement aux op�rateurs bool�ens et, ou, ou exclusif et n�gation.

Alors , l'ensemble des variables bool�ennes repr�sentant les �tats de sortie, est donn� par

(1)

ou

(2)

selon la chiralit�.

        La premi�re id�e fut de partir du tableau fabriqu� pour les programmes pr�c�dents. Le projet a �t� bas� sur l'espoir d'utiliser des outils existants pour fabriquer l'�quation centrale :

• Donner le tableau � un compilateur VHDL
• Transformer la sortie du synth�tiseur logique en code C
• Compiler le code C et g�n�rer du code assembleur
• Copier/coller le kernel de calcul dans le programme
• Assembler le programme

• Le tableau n�cessite une traduction en VHDL : c'est possible mais ce n'est pas un langage de "programmation" courant.
• Le synth�tiseur est inadapt� � cette t�che pour de nombreuses raisons. En premier lieu, le type de cible a une architecture compl�tement diff�rente, donc les contraintes et les r�gles de synth�se ne conviennent pas. Par exemple les op�rations permises par les circuits logiques programmables ne sont pas les m�mes que celles permises par le processeur, la synth�se s'effectue par optimisation de "sommes de produits" (un grand OU de ETs logiques). Le processeur au contraire n'a que des instructions � deux op�randes et dispose du XOR mais pas du NOT. Il dispose de peu de registres alors que les circuits permettent un grand parall�lisme interne. En conclusion, la synth�se g�n�re de nombreuses op�rations tr�s parall�les (afin de r�duire le chemin critique) au lieu de r�duire le nombre total d'op�rations logiques � deux op�randes et de diminuer le nombre de variables temporaires.
• La transformation de rapport de synth�se en code C n'est pas une t�che ais�e, elle doit �tre effectu�e � la main.
• Aucun compilateur � l'heure actuelle ne permet de g�n�rer de code MMX automatiquement.
• L'allocation des registres diff�re entre un code g�n�r� par un compilateur et un code �crit en assembleur.

        Afin de programmer � la main le code de collisions, leur fonctionnement a �t� �tudi� plus en profondeur. L'analyse du code VHDL fabriqu� pr�c�demment a montr� que des collisions manquaient et que les r�gles n'�taient donc pas satur�es. La fabrication des r�gles de collisions a �t� recommenc�e et a donn� ce tableau suivant.

        Le principe de conception repose sur le fait que les r�gles de collision sont plus efficaces quand il y en a le plus possible et lorsqu'elles sollicitent au maximum la particule immobile. Une recherche exhaustive des combinaisons et des configurations int�ressantes est difficile � la main, et aucun logiciel ne peut nous aider. Une des mani�res de constituer les r�gles est de r�duire les cas � un seul vecteur de mouvement, lorsque l'"impulsion" est dirig�e dans un seul sens. Il suffit alors de chercher les r�gles pour la direction puis d'appliquer cinq rotations pour obtenir les autres. Mais ce n'est pas une r�gle tr�s formelle.

        A force de t�tonnements, je suis arriv� � une classification plus formelle avec un tableau � deux dimensions, index� par la longueur du vecteur d'impulsion et le nombre de particules mobiles.

Cette table fait appel � plusieurs astuces de simplification, dont : * La dualit� particule/trou qui permet de ne traiter que la moiti� des cas. Ainsi, lorsqu'il y a plus de 3 particules en entr�e, les particules A � F sont invers�es, ce qui ne change pas la diff�rence � appliquer en sortie. Cette d�tection du nombre de particules fait l'objet d'un code pr�liminaire avant l'application de la table. * La particule G sert pour affiner la s�lection du type de collision : E1 et E2 sont compl�tement ind�pendants.         La difficult� conceptuelle r�side dans l'utilisation correcte des valeurs al�atoires. En sortie de ce tableau, nous disposons de plusieurs sous-r�sultats dont la combinaison permet d'effectuer un XOR sur les particules du site trait�.

IV.8 M�thodes de programmation manuelle en assembleur

        Ce chapitre pr�sente la technique d'optimisation en assembleur qui a �t� utilis�e dans certaines zones du programme exp�rimental. Malheureusement, certaines parties sont tellement complexes que les bugs y sont inextricables, le programme en annexe A pr�sente la version fonctionnelle mais pas compl�tement optimis�e. Tous les d�tails importants sont consign�s en annexe C.

        La technique pr�sent�e ici est l'�volution naturelle des techniques manuelles de codage pour Pentium, ou en g�n�ral pour processeur pipelin� superscalaire. Elle n'est pas d�riv�e de techniques du "Dragon Book" [6] qui traite trop peu de l'optimisation et sous l'angle de code g�n�r� par compilateur, donc des techniques de base et �videntes pour un humain. Le Dragon Book ne traite pas non plus des processeurs superscalaires et suppose un nombre "suffisant" de registres. Nous nous trouvons dans un cas o� la technique de "coloration de registres" � la main n'est pas possible � cause de la taille des graphes � traiter et du faible nombre de registres disponibles. Nous allons voir ici les contraintes cruciales qui ont n�cessit� de concevoir cette technique manuelle.

        La plateforme cible de nos efforts est le Pentium MMX et le Pentium II, deux processeurs au nom similaire mais aux caract�ristiques fondamentalement diff�rentes. Le premier est un processeur superscalaire � deux pipelines de cinq � sept niveaux, alors que le PII est une machine complexe qui ex�cute les instructions dans le d�sordre afin de pallier les lenteurs de la m�moire. Il n'est pas possible de faire deux versions du code pour des raisons de vitesse de codage. D'ici � ce que le projet aboutisse vraiment, le Willamette sera r�pandu. Nous allons nous efforcer de retenir les caract�ristiques de codage communes pour les deux plateformes, ce qui permet au code de fonctionner "honnorablement" partout au lieu d'�tre trop cibl�. En effet, le code n'a pas �t� test� sur des processeurs concurrents, bien que Yannick Sustrac ait report� des performances int�ressantes avec son Cyrix.

        La raison principale et cruciale d'utiliser l'assembleur dans ce projet est que les compilateurs actuels ne peuvent pas g�n�rer "automatiquement" les instructions MMX. Au mieux, ils supportent des pragmas, des librairies ou l'inclusion d'opcodes, mais alors quel est l'int�r�t du compilateur s'il faut faire soi-m�me le codage en assembleur ?

        Il est tr�s peu probable que l'extension MMX du jeu d'instructions x86 ait �t� destin�e � �tre facilement utilis�e dans les compilateurs classiques. Par contre, comme nous pourrons le voir plus tard, Intel a compris l'int�r�t d'un jeu d'instructions simplifi� et plus orthogonal. M�me si les caract�ristiques sont encore �loign�es d'un jeu d'instruction RISC habituel, l'othogonalit� du codage MMX facilite beaucoup la programmation : pas de registre sp�cial ni d'opcode complexe.

Intel Architecture MMX(TM) Instruction Set (Courtesy of Intel) Packed Arithmetic Wrap Around Signed Sat Unsigned Sat Addition PADD PADDS PADDUS Subtraction PSUB PSUBS PSUBUS Multiplication PMULL/H Multiply & add PMADD Shift right Arithmetic PSRA Compare PCMPcc Conversions Regular Signed Sat Unsigned Sat Pack PACKSS PACKUS Unpack PUNPCKL/H Logical Operations Packed Full 64-bit And PAND And not PANDN Or POR Exclusive or PXOR Shift left PSLL PSLL Shift right PSRL PSRL Transfers and Memory Operations 32-bit 64-bit Register-register move MOVD MOVQ Load from memory MOVD MOVQ Store to memory MOVD MOVQ Miscellaneous Empty multimedia state EMMS

Nous voyons dans le tableau ci-contre un extrait de la pr�sentation de l'extension MMX par Intel, avec la base des mn�moniques (sans leurs nombreuses variations). Le jeu d'instruction MMX est � la fois plus organis� que les instructions normales et pourtant il n'est pas compl�tement orthogonal. Nous regr�tons par exemple que certains formats ne soient pas support�s avec certaines instructions, que certaines instructions comme la rotation ne soient pas inclues, que seules 4 op�rations logiques soient disponibles (la n�gation d'une donn�e n�cessite 2 instructions !) et d'autres d�tails qui se r�v�lent dans la pratique. A un niveau sup�rieur, les instructions sont toujours au format op reg,reg/mem ce qui limite les possibilit�s, surtout avec seulement 8 registres. Enfin, selon les architectures, les limitations sont vari�es et parfois contradictoires. Il faut pouvoir jongler entre les nombreuses contraintes et �viter tous les types de blocages, au niveau du d�codage, de la m�moire et des unit�s d'ex�cution, tout en allouant correctement le peu de registres disponibles.

        Les plus grandes limitations concernent les r�gles de d�codage et de pairage : certaines instructions ne peuvent �tre d�cod�es qu'� certaines positions du 'slots' de pairage. Pour le PMMX, les instructions vont par paire (une par pipeline, U et V), alors que sur le PII le d�codage doit suivre la r�gle 4-1-1 : le groupe doit faire moins de 16 octets, la premi�re instruction du groupe ne peut g�n�rer que 4 �ops et les deux suivantes qu'une seule (les tables de correspondances instruction->�ops sont donn�es par les manuels Intel). Les r�gles sont g�nantes dans la pratique et nous retiendrons principalement celles-ci :

• Pairage de style PMMX : deux instructions par cycle.
• Instruction d'�criture ou de lecture m�moire dans le premier 'slot' (pairage en U).
• Instructions de Shift dans le pipeline V.
• Ecriture vers la m�moire : deux cycles apr�s la modification du registre (contrainte du PMMX).

        J'ai donc d�velopp�, en programmant � la main des exemples de plus en plus complexes, une technique bas�e sur les graphes de d�pendances de donn�es. L'astuce principale r�side dans l'algorithme tr�s simple pour g�n�rer le code optimis� (en lisant le graphe s�quentiellement par la fin) mais cette facilit� est largement compens�e par le travail n�cessaire pour pr�parer le graphe. Le workflow de programmation est le suivant :

• V�rification de l'algorithme par du pseudo-code dans un langage de haut niveau, par exemple avec Turbo Pascal
• Fabrication du graphe de dataflow (voir le graphe initial)
• R�duction de la largeur du graphe si n�cessaire, d�termination des variables � stocker temporairement en m�moire
• D�termination des "faisceaux" (crit�re du plus petit, sauf avec ANDN)
• D�termination du "chemin critique" (CDP en anglais)du graphe
• R�arrangement autour du CDP
• Contrainte principale : �quilibrage des load/store
• Association des registres aux faisceaux
• D�termination des "slots" de d�codage
• Arrangement fin des instructions
• "Lecture" du graphe (par un simple balayage) et �criture du code assembleur correspondant

        Chaque �tape n�cessite � la fois un algorithme et du bon sens. Il ne s'agit pas d'appliquer aveugl�ment une m�thode mais d'�tudier au cas par cas l'influence de telle ou telle transformation. L'analyse du chemin critique du graphe peut ne pas �tre n�cessaire si le graphe est simple et petit, toutefois elle devient n�cessaire pour casser certaines d�pendances complexes. L'allocation des registres se fait en deux �tapes : les branches du graphes sont d'abord r�unies en "faisceaux" lorsqu'il faut d�terminer quel registre sert de destination (le choix n'est pas possible pour l'instruction PANDN par exemple), puis les registres sont associ�s aux faisceaux au dernier moment, lorsque le graphe est �quilibr�. Ce dernier est toujours balay� par le fin : on remonte et on alloue les registres en commen�ant par le premier registre libre. Un registre est donc "occup�" lorsqu'il est associ� � un faisceau et il est lib�r� lorsque le faisceau se termine : le registre est donc libre de servir pour un autre faisceau. Cet algorithme est presque automatique lorsque le graphe est adapt� au nombre exact de registres disponibles, c'est la raison pour laquelle la pr�paration est si importante.

Par exemple, la table de d�tection des collisions est repr�sent�e dans le graphe ci-contre, apr�s traduction du pseudo-code et explicitation temporelle des d�pendances. A ce stade d'�tude, le but est d'explorer au maximum le parall�lisme des instructions permis par le graphe. Cela permet de d�gager les axes importants et le chemin critique, les besoins en instructions d'acc�s � la m�moire, leur r�partition...

Ensuite, le nombre de registres est r�duit afin de pr�server 3 registres tournants (pour le renommage des entr�es). La r�partition des acc�s � la m�moire devient importante. Les feuilles sont repli�es, les donn�es persistantes (celles qui resteront dans les registres) sont choisies en fonction du nombre d'utilisation, les autres sont stock�es en m�moire.         Une fois les registres allou�s et le graphe pr�t, on va balayer le graphe de bas en haut et de gauche � droite, tout en consignant les noeuds rencontr�s sous forme d'opcodes dans le fichier assembleur. Lors de cette phase terminale, il faut veiller � ce que les r�gles de pairage soient respect�es. Si une des instructions pr�c�demment consign�es entre en conflit (slot ou d�pendance de registre par exemple) avec l'instruction courante, il faut r�organiser le balayage de la ligne et cela d�borde parfois sur d'autres lignes.

        Le code de collision est d�crit sch�matiquement ci-dessous, comme s'il devait �tre c�bl� :

• La premi�re partie compte le nombre de bits afin de profiter de la dualit� du mod�le et ainsi r�duire la complexit� de la table. Si le nombre de bits � 1 dans les particules A � F est sup�rieur � 3, toutes les particules sont invers�es. Cela revient � effectuer un XOR de ces particules avec le bit 2 du r�sultat.
• La deuxi�me partie est la d�tection elle-m�me : elle prend 4 entr�es et g�n�re 3 r�sultats temporaires en fonction du tableau et du code d�crits pr�c�demment. Les registres sont allou�s afin de pouvoir garder au moins trois valeurs d'entr�e et les faire "tourner" sans n�cessiter d'acc�s m�moire.
• La troisi�me partie combine les r�sultats temporaires. Une partie de ces op�rations (combinaison des Pn) est effectu�e dans la deuxi�me moiti� de la d�tection pour b�n�ficier de la localit� des donn�es. L'autre partie de la combinaison est entrelac�e avec le code de d�placement de donn�es qui utilise aussi une allocation "tournante" des registres.

        Le d�veloppement d'un code sain aussi complexe n�cessite de bonnes habitudes de codage et beaucoup de bon sens. Par exemple, lorsqu'un morceau de code fonctionne, il n'y a plus besoin dans la plupart des cas de le modifier : il ne faut pas attendre qu'un bug apparaisse pour sauver le fichier sous un autre nom. Lorsqu'un probl�me se d�clenche, cela limite le nombre de lignes � v�rifier et permet de toujours avoir un programme fonctionnel sur lequel se rabattre lorsque le d�veloppement est dans une impasse. Ainsi le loader 32 bits est stable car il est d�velopp� depuis longtemps, il remplit son r�le et la plateforme ne change pas.

        Les bugs ne sont pas des �tres malins destin�s � nous faire perdre la raison mais c'est un �cart du programme dans la r�alisation de la fonction d�sir�e. Puisqu'un ordinateur effectue exactement ce qu'on lui demande, l'�cart ne peut provenir que des ordres erron�s que le programmeur lui fournit.

        Pour chasser les bugs, les m�mes recettes sont applicables pour toutes les situations :
    1)   Relire tr�s attentivement le code, plusieurs fois si n�cessaire. Simuler mentalement le fonctionnement du processeur et l'�tat de chaque registre, comparer avec les sympt�mes et chercher des similitudes avec ceux-ci.
    2)   R�unir tous les indices, tous les sympt�mes et noter toutes les conditions d'utilisation lorsque le bug se d�clare. Cela peut aider � d�terminer la partie fautive. Par exemple, si le bug est un probl�me d'affichage, il faut d�terminer toutes les parties du code qui acc�dent � l'�cran, directement ou non.
    3)   Isoler la partie fautive : la technique la plus simple est de mettre en commentaire (ou %ifdef zorglub) certaines parties du code r�cemment modifi�es puis restreindre progressivement les recherches par dichotomie jusqu'� ce que le bug soit circonscrit. Cela peut prendre une minute ou plusieurs heures et il ne faut pas h�siter � remettre en cause toute la structure du programme.
    4)   Identifier la d�pendance entre les donn�es et les codes concern�s : on d�couvre parfois que l'erreur vient d'une autre partie que la partie incrimin�e ou d�tect�e par la dichotomie. Le bug d�couvert peut �tre la cons�quence d'une erreur en amont mais qui ne se manifeste que dans certaines conditions que la partie incrimin�e remplit. Certaines donn�es ou certaines parties du code peuvent d�clencher un comportement ind�sir� ou inattendu. Le plus souvent, il s'agit "simplement" d'un oubli ou d'une erreur de frappe. Par exemple, mettre une relocation d'adresse vers l'�cran au lieu du tunnel. Dans l'empressement, ce type d'erreur est facile � commettre mais aussi � corriger. Il en va diff�remment du renommage des registres ...
    5)   Si trop de temps s'est �coul� lors du d�buggage, il faut souvent tout recoder � partir de z�ro. Il est parfois dix fois plus rapide de recommencer que de tout analyser. Le nouveau code est m�me souvent de meilleure qualit� car les bugs se cachent souvent dans les parties r�centes et peu �prouv�es, qui auront �t� m�ries depuis le premier jet. C'est parfois la solution de la derni�re chance pour �liminer les bugs les plus r�calcitrants.

        Il n'est pas forc�ment n�cessaire dans ce cas de disposer d'un d�bogueur interactif avec suivi du code source : il faut de la m�thode, respecter quelques r�gles simples et surtout ne jamais faire confiance aveugl�ment � un code sans l'avoir test� dans les conditions d'utilisation pr�vues ...

IV.9 : R�alisation et r�sultats :

        Le code utilis� actuellement n'est pas la version d�finitive du programme : le l�ger bug dans les collisions passe presque inaper�u mais emp�che d'avoir un r�sutltat r�ellement utilisable. La version beta est toutefois suffisamment rapide pour montrer le gain substantiel apport� par l'analyse soign�e du programme. Sa vitesse d'ex�cution a �t� compar�e sur le m�me ordinateur avec un autre code en C d'allure plus classique :

plateforme : Pentium MMX 200MHz, SDRAM 66MHz, L2 : 256Ko
g�om�trie : 1000 it�rations sur 1024*640  (655M site calcul�s)

SP65.asm: 15 secondes, strip=6, zoom=4:1
vitesse : environ 43Mc/s

Programme de Benjamin Temko (voir l'annexe D)
compilation:
gcc -I/usr/X11R6/include -lX11 -L/usr/X11R6/lib -O3 fhp_0.c
run :
nice --20 time -o perf ./a.out

   58.76 user
    0.44 system
 1:00.87 elapsed
     97% CPU
(0 avgtext+0avgdata 0maxresident) k0 inputs+0outputs
 (153major+178minor)page faults
 0 swaps

vitesse :   11.2Mc/s

        La version alpha du code de collision (lorsqu'il sera d�bogu�) devrait permettre d'am�liorer la vitesse d'environ 20 ou 30% gr�ce � un meilleur scheduling des instructions dans la phase de d�tection et de d�placement. La meilleure performance absolue actuelle a �t� mesur�e au MIT en janvier 2000 sur un Pentium III � 550MHz :

taille: 860*736, cycles par pas de temps : 2,3M, 3,6 cycles par site, soit   150 Mc/s.
 

        L'exp�rience suivante a �t� r�alis�e pour fabriquer une petite animation illustrant les capacit�s du programme existant. Le temps de calcul total est de quelques minutes, en comptant le long temps de sauvegarde des images qui ont �t� choisies pour leur esth�tique, non pour leur int�r�t scientifique. Le nombre de pas entre les images n'est pas fixe mais permet de montrer l'int�r�t d'un affichage plus �volu�, non en "noir ou blanc".

Les dimensions des images suivantes tiennent compte de la projection sur le r�seau hexagonal (sin(60�)).

Ecoulement apr�s 528 pas de temps : la deuxi�me chambre se remplit lentement et les artefacts hexagonaux sont visibles.	T = 4501 : le flux commence � se distordre.	T = 6010
T = 6316	T = 6931 : le flux est moins marqu�, noy� dans le bruit.	T = 8089 : la pression et la vitesse sont favorables pour faire disparaitre les artefacts hexagonaux mais le contraste est trop faible pour distinguer les �coulements.

        Nous voyons ici les limites de l'affichage de la version de d�veloppement : le bruit rend la distinction de ph�nom�nes fins difficile. La refonte du programme permettra d'inclure le code d'affichage � niveaux de gris qui est d�crit en au chapitre IV.6.

        L'exp�rience s'est d�roul�e sur le PII-350 du d�partement MIME. Elle consiste en la diffusion d'une masse de particule de la chambre gauche � la chambre droite du tunnel. Les chambres �tant s�par�es par une paroi munie d'un "trou" de quelques pixels, l�g�rement d�centr� pour que le flux de particules oscille. Apr�s quelques dizaines de milliers de pas de temps, la densit� dans les deux chambres est similaires et le flux disparait, laissant une "bouillie brownienne homog�ne" � l'�cran. Les dimensions et les param�tres initiaux sont les suivants :
- chambre gauche : densit� : 1, hauteur : 725, largeur : 434
- chambre droite : densit� : 0, hauteur : 725, largeur : 267
- paroi : largeur du trou : 26 noeuds
- vitesse : strip=32, 2014269 cycles CPU par pas de temps, 5,75 ms pour balayer un tableau de 704*728=512512 noeuds, soit 89 millions de noeuds par seconde.

        L'image suivante montre le r�sultat de la routine de calibration du strip mining. Nous pouvons y apercevoir les caract�ristiques de la m�moire et en particulier mesurer le speedup permis par l'algorithme, par rapport � un balayage lin�aire normal (4981/2014=2,47). En particulier, le speedup varie peu au-del� d'un strip mining de 10 lignes, la d�croissance est d�e � la r�duction de la bande passante utilis�e par l'affichage et est en 1/x.

        Les algorithmes de gaz sur r�seaux ont d�j� fonctionn� sur une large vari�t� de machines dans le monde, par exemple : Alpha, Apple, CRAY-XMP, CRAY-2, Connexion Machine 200, Connexion Machine 5, FPS 164, IBM RS6000, IBM 3090, PC (Intel et clones), SUN/SPARC (jusqu'� 256 CPU), Silicon Graphics Origin 8-CPU, VAX...
La nature du mod�le des gaz sur r�seaux s'adapte � une tr�s grande vari�t� de machines, des processeurs 16 bits aux calculateurs vectoriels parall�les en passant par les machines de bureau actuelles. Elle permet aussi aux �lectroniciens de cr�er des machines d�di�es, plus ou moins configurables : c'est le sujet de ce chapitre.

        Les avantages sont nombreux : si une programmation tr�s soign�e permet de "gagner" 3 ans en performance, elle peut en faire gagner le double ou plus dans certains cas. Par exemple, d�s ses d�buts, la CAM-8 a atteint des performances encore in�gal�es par les stations de travail actuelles. Le parall�lisme inh�rent au mod�le permet aussi de cr�er des architectures modulaires qui peuvent s'adapter aux contraintes budg�taires fluctuantes, laissant esp�rer des mont�es en puissance vertigineuses. Cela permet enfin de r�soudre le plus simplement du monde des probl�mes de plus en plus difficiles sur les plateformes actuelles : un simple fil peut remplacer de nombreuses lignes de code.

        Tous les essais effectu�s ces quinze derni�res ann�es ont aussi montr� les nouveaux probl�mes qu'une architecture d�di�e entraine. Tout d'abord, le manque de d�bouch�s sur le march� du grand public marginalise les efforts : combien de personnes veulent calculer un milliard de sites FHP par seconde ? La cons�quence directe est que ces architectures suscitent peu d'int�r�t r�el en dehors du monde restreint des gaz sur r�seaux, donc les machines existantes ne sont que des prototypes. Les autres machines sont rest�es dans l'imagination des r�veurs : si concevoir une architecture "accueillante" pour les gaz sur r�seaux semble tr�s facile, les probl�mes �lectroniques et de moyens mat�riels sont bien diff�rents dans la r�alit� pour un projet de cette �chelle.

        Enfin, il ne faut pas perdre de vue qu'une fois construit, le mat�riel n'�volue pas. Il faut souvent revoir compl�tement l'architecture pour chaque g�n�ration de machine, comme pour un ordinateur normal, mais le co�t de d�veloppement est lourd et constant, � la charge du laboratoire qui lance le projet. Il n'y a donc pas de "lign�e" de calculateurs d�di�s comme il y en a pour les voitures, les logiciels, les avions ou les ordinateurs. Les changements d'organisation dans les entrailles de la machine entre chaque nouveau prototype obligent � r��crire tout le logiciel de gestion pour chaque version. Les cycles de d�veloppement sont donc tr�s longs et on perd une partie de l'avance par rapport � des logiciels sur des stations de travail.

        Les quatre premiers chapitres d�crivent des impl�mentations r�elles, les suivants sont des �tudes de cas imaginaires (architecture fiction).

V.2 : RAP-1 :

        Mis au point vers 1986 par Dominique d'Humi�res et Andr� Clouqueur � l'Ecole Normale de la rue d'Ulm, c'est une machine typique de son temps comme le montrent ses caract�ristiques. Ses performances sont de 6,5Mc/s, ou 256 sites par 512 � 50 Hz. Les sites peuvent contenir 16 bits qui sont mis � jour lin�airement par une LUT (2 SRAM de 64Ko) et l'�tude en cours peut �tre visualis�e sur un �cran VGA ind�pendant (synchronisation directe, pas de frame buffer). L'h�te est un PC et les informations peuvent �tre post-trait�es par un VAX en local [11].

        Le RAP-1 a une structure proche de la CAM-6 mais ces deux architectures resten distinctes. La m�moire centrale du RAP-1 est compos�e de VRAM, des puces DRAM de 64 Kbits con�ues pour les cartes vid�o pouvant lire et �crire une donn�e en un cycle gr�ce � deux ports s�par�s. Des circuits sp�ciaux d'adressage permettent la r��criture vers les sites voisins pour certaines lignes, ce qui permet de g�rer les voisinages hexagonaux, de Moore ou de Von Neumann. La CAM-8 permet la r��criture du r�sultat de chaque site vers n'importe quel autre site.

        Le "RAP-1" signifie "R�seau d'Automates Programmable" et n'a plus qu'un int�r�t historique actuellement. Il est limit� architecturalement et ne permet pas de contr�le fin ou d'actions complexes comme le permettent les logiciels. La performance, qui nous int�resse ici, est largement d�pass�e par les PC actuels. Voir l'article dans [20] pour plus de pr�cisions architecturales.

Une all�e de Von Karman calcul�e sur RAP-1 et post-trait�e.

V.3 : CAM-8 :

        La CAM-8 est probablement la machine la plus int�ressante actuellement. Elle est issue des recherches de Tomaso Toffoli et Norman Margolus au MIT o� ils ont con�u plusieurs g�n�rations de "CAM" (Cellular Automaton Machine). Ils ont donc une exp�rience et une renomm�e qui leur ont permis de concevoir une architecture flexible et performante. En particulier, sa flexibilit� est bien sup�rieure par rapport au RAP-1 (dont elle s'inspire un peu), ce qui lui vaut d'�tre encore en usage. Dix � quinze exemplaires ont �t� fabriqu�s depuis 1992 : un record pour le domaine !

Une boite CAM-8 avec une carte ins�r�e sur laquelle on aper�oit certains ASIC.

Description d'un ASIC (STEPchip) (extrait du manuel de la CAM-8).

        La CAM-8 est une machine complexe, d�crite dans le manuel disponible sur le site qui lui est d�di� (http://www-im.lcs.mit.edu/cam8/ps/hard_ref.ps). Pour r�sumer, la machine est con�ue sur plusieurs �chelles : syst�me, bo�te, carte et ASIC.

Connexion du fond de panier d'une boite CAM-8 pour une configuration infinie. Seule une direction est montr�e, sans bouclage aux bords.

Connexion de 4 boites pour un tore 2D.

        Comme pour le RAP-1, un moniteur externe permet de visualiser l'activit� de la machine (c'est assez spectaculaire) mais le plus complexe reste la configuration de la machine par l'h�te. Une suite de logiciels est en d�veloppement depuis de nombreuses ann�es. De nombreuses applications ont �t� d�montr�es et la CAM-8 sert pour des recherches dans des domaines tr�s vari�s mais l'interface est encore rudimentaire.

        La CAM-8 est tr�s flexible. Il y a quelques limitations mais elles s'inscrivent dans une architecture sophistiqu�e, destin�e � simuler virtuellement des espaces jusqu'� 32 dimensions. L'espace m�moire est adress� dans chaque carte en associant un ou plusieurs bits � une dimension, ce qui permet de cr�er des g�om�tries quasi arbitraires o� chaque longueur est une puisance de 2 (256*256*256, 128*512, 4096*4096...).

        Chaque carte effectue une consultation de table sur 16 bits pour chaque cellule mais on peut effectuer des op�rations plus complexes (consultation multiple, LUT virtuelle...). Par exemple, pour le mod�le FHP, les 16 bits peuvent servir pour repr�senter 2 sites � 1 bit par directions ou 1 site � 2 bits par directions. Avec 8 cartes � 25 MHz, un boite CAM-8 soutient donc 400Mc/s en FHP3 ou 200Mc/s avec 2 bits par direction (Integer Lattice Gas � 2 bits).

        Le voisinage est configurable site par site : une carte peut acc�der � des �l�ments non contigus dans la m�moire. Cela permet aussi bien de calculer dans des espaces � N dimensions ou effectuer des interactions non locales. C'est ce dernier d�tail qui fait une grande diff�rence avec les autres machines.

        Le site Internet du laboratoire de Norman Margolus n'est pas avare en d�tails techniques. On peut y trouver les applications o� la CAM-8 excelle et m�me un simulateur logiciel. Un article [23] dans la preprint archive montre une application de la CAM-8 aux simulations en 3 dimensions dans le r�seau FCHC � 24 bits. La capacit� � effectuer des consultations successives de la LUT permet de "casser" l'op�rateur de collisions en isom�tries de 16 bits.

V.4 : EXA :

        La soci�t� EXA a �t� cr��e par un professeur du MIT, Molvig. Elle utilise les techniques issues des gaz sur r�seaux dans les milieux industriels. Son nom vient de l'ordre de grandeur du m�me nom (mega, giga, tera, peta puis exa). C'est le nombre d'op�rations n�cessaires pour des simulations r�alistes en 3D, hors de port�e des ordinateurs classiques de l'�poque (vers 1990).

        EXA vend ses services aux soci�t�s qui ont besoin d'�tudes sur des projets particuliers, dans les limites permises par les mod�les utilis�s. Sa puissance de calcul consistait au d�but en une station Silicon Graphics dot�e de cartes acc�l�ratrices dot�es d'ASIC traitant les calculs bool�ens. Le premier mod�le (Molvig/Vichniac, vendu sous le nom "DigitalPhysics") �tait une extension du mod�le pseudo 4D avec 2 vitesses, soit un mod�le thermique avec 48 bits par site. Peu d'informations subsistent sur ce mat�riel qui souffrait des limitations inh�rentes aux mod�le bool�ens.

        Rapidement, la soci�t� a �tendu ce mod�le en utilisant la technique de Bolzman (BGK ?) pour b�n�ficier de la mont�e en performance des calculateurs massivement parall�les constitu�s de stations de travail (beowulfs). EXA d�veloppe, utilise et vend un logiciel appel� Powerflow. Il existe peu de papiers d�crivant le mod�le utilis� mais l'approche choisie est propri�taire (secr�te) et utilise des astuces contestables pour augmenter artificiellement le nombre de Reynolds simulable.

        L'espace simul� est divis� en "voxels" (sous-divisions de l'espace 3D en cubes) de tailles variables dans l'espace simul�. Cela permet d'adapter la quantit� de m�moire et de calculs n�cessaires en fonction de la vorticit� locale. La "granularit�" sera beaucoup plus fine pr�s des parois et dans les zones turbulentes. Certains sp�cialistes ont object� que cela brise la continuit� du milieu et de ses propri�t�s, emp�chant les turbulences de se propager d'un voxel � un autre d'une taille (donc d'une viscosit�) diff�rente. Les techniques qui y rem�dient sont secr�tes et ne peuvent pas �tre examin�es.

        On pense que certains artifacts et d�fauts sont compens�s par des marges de s�curit� et une analyse "intelligente" du partitionnement de l'espace. Une approche similaire au maillage dynamique ou progressif (en cours du calcul) permet de r�duire l'impact du probl�me dans une majorit� des cas. Le reste est noy� dans le bruit num�rique et l'int�gration lors des mesures. Enfin, cet inconv�nient et les �ventuels autres probl�mes sont souvent n�gligeables par rapport aux avantages, lorsqu'ils sont compar�s aux techniques classiques.

        Les gaz sur r�seaux n'ont pas d'avantage majeur en temps de calcul ou en utilisation m�moire par rapport aux techniques �tablies, pour un cas identique. Par contre, la qualit� des r�sultats est souvent sup�rieure : ils sont fiables � quelques pourcents alors que les techniques classiques sont fiables � quelques dizaines de pourcents. De plus, ils sont r�solus explicitement dans le domaine temporel, ils peuvent donc capturer la dynamique d'un probl�me qui n'apparaitrait pas autrement (par exemple avec une r�solution o� le terme temporel est ignor�). Ensuite, les r�sultats sont plus pr�cis. Enfin, contrairement aux autres techniques, la simulation n'a pas besoin de donn�es issues de simulations r�elles pour ajuster les r�sultats des calculs : il suffit juste de fournir les g�om�tries et les param�tres du fluide (Re et Mach) pour obtenir le r�sultat. Cette derni�re qualit� a permis � EXA de se distinguer en r�solvant des probl�mes de turbulences complexes inaccessibles aux autres m�thodes. On comprend donc que les techniques employ�es ne soient pas publiquement connues.

        Le domaine d'utilisation reste limit� au bas subsonique (inf�rieur � Mach 0.4), l'a�ronautique n'est donc pas la cible de l'entreprise. Les �tudes portent sur les car�nages et carosseries de motos, voitures, camions, jusqu'� Re=6.10^6 avec un (tr�s) gros serveur. Le domaine a r�cemment �t� �tendu aux �coulements internes et aux �changes thermiques. Il devient ainsi possible de simuler des injections de moteurs � explosion ou des syst�mes de climatisation.

        Dans la pratique, cette qualit� des calcul a un prix tr�s cher en expertise et en temps CPU. Pour donner un ordre de grandeur, un grand fabricant de voiture a d�cid� d'acheter un cluster de 256 stations de travail dans le but unique d'utiliser Powerflow. Cette initiative a �t� suivie par la plupart des constructeurs europ�ens et am�ricains. EXA loue aussi du temps CPU sur un cluster SGI Origin � 8 CPU en mode batch (partag� entre plusieurs utilisateurs moins fortun�s) et inaugure ainsi un nouveau mode de "service commercial" sur Internet (la mode est vraiment au "e-business" bien que le travail par lot ne soit pas une invention r�cente).

        Peu de d�tails de programmation filtrent sur Internet ou sur le site web de la soci�t�. Une interview avec un ing�nieur fran�ais a cependant permis de d�gager certains points. EXA utilise un r�seau FCHC avec des particules � trois vitesses (0, 1 ou 2 sites par pas de temps) en virgule fixe. L'espace est divis� en "voxels" et en "surfels" de tailles variables (mais multiples les unes des autres pour simplifier le pavage de l'espace). Chaque type de zone ("surfel" ou "voxel") est trait� par un code sp�cial. Un voxel utilise environ 130 octets de m�moire, un surfel en utilise 1300 et peut repr�senter 12 ou 24 faces. Une simulation normale utilise facilement plusieurs millions de voxels, ce qui n�cessite l'emploi d'ordinateurs multiprocesseurs disposant de plusieurs gigaoctets de m�moire. Enfin, un pas de temps �quivalent � une seconde n�cessite 22000 pas de temps de calcul.

        La phase critique de l'expertise est l'adimensionnement : la d�termination des param�tres de calcul en fonction des param�tres r�els. Le but est d'utiliser la puissance de calcul le plus efficacement possible. Les param�tres importants sont la puissance de l'ordinateur (m�moire et vitesse), la vitesse du fluide et la taille de l'�prouvette. La vitesse du fluide dans l'exp�rience doit �tre maximale (environ 0.3) pour maximiser le nombre de Reynolds. Le nombre de Reynolds d�pend aussi de la viscosit� du fluide, qui est fixe pour chaque voxel, et de la longueur caract�ristique du probl�me. L'adimensionnement d�termine donc le nombre de voxels � utiliser en fonction de la d�finition du nombre de Reynolds et en ajustant les termes de l'equation de Navier-Stockes aux �quations caract�ristiques du mod�le avec les termes non-lin�aires qui apparaissent [26]. Le temps de calcul d�pendra du nombre de voxels, du nombre de pas de temps, de la complexit� du fluide et du nombre de processeurs. C'est un algorithme facilement parall�lisable donc l'efficacit� reste satisfaisante avec plusieurs centaines de processeurs, comme le montre une fiche technique disponible sur le site web. Les plateformes de calcul utilis�es sont couramment Sun et SGI. Les r�sultats ou les vecteurs de calculs sont trait�s sur des stations de travail graphiques d�port�es.

        EXA a eu l'amabilit� de communiquer un exemple de calcul r�alis� en 1998 pour donner un ordre de grandeur de l'efficacit� du logiciel. Le calcul portait sur les turbulences autour de la carosserie d'une voiture sportive afin d'�valuer l'efficacit� de volets de d�flexion. Cela a n�cessit� 22000 pas de temps soit 4 jours sur un serveur � 8 processeur � 167 MHz avec 11 millions de voxels et 1 million de surfels. Les entr�es d'air ou les interactions des pneus n'ont pas �t� prises en compte mais le r�sultat est saisissant.

        Cette approche logicielle a beaucoup d'avenir car elle permet concevoir des v�hicules plus silencieux et moins turbulents sans louer ou construire des souffleries r�elles (sans compter le prix de la maquette). Powerflow permet de mod�liser des turbulences li�es aux asp�rit�s d'une structure et donc de v�rifier son impact sur le bruit qu'il g�n�re et sur le Cx de l'engin. Le temps de calcul peut �tre plus court que la construction et les mesures d'une maquette dans une soufflerie silencieuse. Toutefois, de nombreux probl�mes subsistent comme la limitation de la taille des fichiers sur les syst�mes UNIX : la taille limite traditionnelle de 2Go est vite atteinte et les efforts de d�veloppement sont encore en cours.

V.5 : "Fourmi" :

        Ce projet du d�partement MIME n'est pas encore op�rationnel. L'architecture est simple dans les grandes lignes et est enti�rement d�di�e aux automates cellulaires 1 et 2D, sans �tre sp�cialement pr�vue pour un usage pr�cis. Ce projet est appel� ainsi car il part du principe que l'union fait la force, comme dans une fourmili�re : une fourmi isol�e n'est pas tr�s productive alors qu'une arm�e de fourmis peut effectuer beaucoup de travail. C'est une architecture extensible par addition de modules : le prototypage s'effectue avec un seul module puis d'autres modules seront construits en nombre (donc moins chers � l'unit�) lorsque la technique est au point.

        La "fourmi" est une machine SIMD divis�e en deux : elle dispose d'un s�quenceur (partie contr�le) et d'une ou plusieurs cartes d'ex�cution. Chaque �l�ment est configurable : le s�quenceur est c�bl� et "ex�cute" un programme charg� en SRAM locale, tout en envoyant des signaux aux parties d'ex�cutions. Chaque carte "active" est compos�e essentiellement d'un FPGA et d'un bloc de DRAM adress�e par le FPGA. Le FPGA est en technologie SRAM, il peut �tre reconfigur� � tout moment par l'ordinateur h�te (Sun SPARC). Les fonctions r�alisables sont donc potentiellement infinies mais limit�es par la taille et l'architecture du FPGA ainsi que par son langage de commande : le projet Fourmi a aussi pour objet de d�velopper un langage de description des automates cellulaires. Le tout fonctionne � quelques dizaines de m�gahertz.

        La pratique est beaucoup plus complexe, comme le confirme la dur�e du d�veloppement. Le programme de commande traduit le langage de description en code VHDL, ce qui est loin d'�tre facile. Le bus VME, autour duquel est architectur�e la machine, ne permet d'envisager qu'une dizaine de cartes au maximum. Les cartes sont reli�es par un r�seau unidimensionnel. Le d�veloppement est frein� par le manque de moyens mat�riels, par le peu de personnes travaillant au projet et par la d�pendence envers les fabricants de FPGA : de nouvelles versions de FPGA et de logiciels apparaissent plus vite que le projet et il est difficile de conserver d'anciennes versions car elles ne sont plus support�es.

        La CAM-8 a �chapp� � ces probl�mes : la topologie est beaucoup plus flexible et en 3D, laissant envisager la construction de cubes de dizaines de boites de c�t� (le probl�me est alors la dissipation thermique...). Le projet CAM-8 b�n�ficie aussi de subventions du d�partement de recherche de l'arm�e am�ricaine, ce qui permet de r�soudre certains probl�mes de mani�re plus triviale que lorsque le budget est d�pass� ou inexistant.

V.6 : Carte ISA :

        Ce chapitre pr�sente la premi�re �tude personnelle de mat�riel d�di� au calcul de gaz sur r�seau hexagonal, vers 1996, apr�s la sortie de l'article de l'annexe B. Aucun acc�l�rateur de cette structure n'a jamais �t� r�alis�. Le but �tait principalement d'effectuer les calculs et les mouvements complexes de donn�es par de simples circuits �lectroniques sur une carte ISA avec un budget tr�s restreint. La plateforme cible �tait un i286 � 12MHz dot� d'un bus ISA 16 bits. La vitesse th�orique est donc limit�e par la vitesse du bus ISA et la taille des exp�riences est limit�e par le x86 en mode r�el (64Ko).

        La structure de la carte est d�riv�e d'une �tude des mouvements de donn�es dans l'algorithme de l'annexe B. La carte ne pouvait �tre dot�e de m�moire spacieuse, elle devait donc d�pendre de la m�moire centrale de l'h�te. Les composants doivent �tre simples, �conomiques et peu nombreux, sur un circuit imprim� � 2 faces dessin� � la main. Aucune exploitation des r�sultats n'�tait pr�vue : pas de sommation car les particules �taient organis�es en multisite et l'affichage �tait contr�l� par la palette de la carte VGA.

        L'�tude des d�pendances des donn�es a permis de d�finir une strat�gie : seules quelques lignes ont besoin d'�tre m�moris�es dans la carte, ce qui permet de ne contenir que quelques dizaines de kilooctets de SRAM ou de FIFO. Dans ce cas, le plus difficile est de synchroniser le contenu des diff�rentes m�moires : il faut que le logiciel collabore �troitement avec la carte pour permettre � l'algorithme de fonctionner correctement.

        Le sch�ma suivant d�crit les d�pendances de donn�es :

Description des d�pendances de donn�es pour le mod�le multisite

        Comme les mots transmis sont larges de 16 bits et puisque la gestion des lignes paires/impaires est inutilement complexe, le r�seau a �t� tourn� � 90 degr�s et les noeuds sont trait�s par paires. Il faut donc 2 LUT et chaque cycle d'horloge traite deux sites. La figure suivante montre le nommage des directions, similaire � la strat�gie du deuxi�me code de r�f�rence :

        Deux strat�gies existent : utiliser des FIFO ou des SRAM. L'�tude a port� d'abord sur les FIFO. Par convention, N sera le nombre de sites par ligne.

        L'ordinateur et la carte sont asynchrones, l'ordinateur doit partager la bande passante entre le flux entrant dans la carte et le flux sortant. Pour des raisons de simplicit�, la gestion de la DMA n'a pas �t� envisag�e. Le programme de l'ordinateur h�te est assez simple : il envoie un bloc puis en re�oit un autre � mettre au m�me endroit. Il faut v�rifer � chaque fois que le buffer que l'on va acc�der est pr�t, gr�ce � des s�maphores pr�vus � cet effet. La boucle interne doit �tre �crite en assembleur car elle utilise des instructions sp�ciale rep insw et rep outsw qui ne sont pas disponibles dans les langages de haut niveau :

        ; init de la boucle externe :
           push ds ; ds pointe vers le tableau
           pop es  ; ds = es = d�but du tableau
           xor di,di ; di = 0
           xor si,si ; si = 0
           mov dx,inout_port ; adresse de la carte ISA
 
           mov bx,nb_lignes
        loop_bx:

          ; synchro 1:
             add dx,2
        loop_dx1:
               in al, dx
               and al,1   ; v�rifie que le buffer d'entr�e est pr�t
             jz loop_dx1
             sub dx,2

             mov cx,N/2
             rep outsw dx,si ; bloc DS:SI vers port DX

        ; synchro 2:
             add dx,2
        loop_dx2:
               in al, dx
               and al,2   ; v�rifie que le buffer de sortie est pr�t
             jz loop_dx2
             sub dx,2

             mov cx,N/2
             rep insw dx,di ; port DX vers bloc ES:DI

           dec bx
           jnz loop_bx

        ; fin

        Si le processeur envoie 2 octets par cycle, s'il faut 2 cycles par site (�criture puis lecture) et si le bus envoie 16 bits par cycle en rafale, la bande passante maximale th�orique est de 12MHz/2*2 soit 12 millions de sites par seconde � 12 MHz. Les al�as du bus ISA abaissent ce d�bit dans la pratique, en particulier pour acc�der � la m�moire DRAM lente et pour garder la compatibilit� avec les cartes ISA 8 bits. Il faut en th�orie 2 cycles pour un acc�s 16 bits et 3 cycles pour deux acc�s 8 bits en mode 16 bits. Cela donne beaucoup de temps � la carte ISA pour traiter les informations par paquets 16 bits, � 6 MHz maximum en th�orie. Le chemin critique de l'acc�l�rateur dispose donc de 160 ns au minimum pour effectuer son travail, ce qui est peu contraignant. Le processeur effectue deux fois plus de travail de mouvement que la carte elle-m�me car sa bande passante est partag�e par plusieurs flux : la carte doit donc pouvoir garder en m�moire deux lignes (2*N sites) au minimum pour �viter de perdre des donn�es.

        En utilisant des FIFO toutes faites, on gagne en simplicit� de gestion : pas d'adressage de SRAM � g�n�rer. La figure pr�c�dente dit en a) qu'il faut conserver pour le buffer d'entr�e 2N sites, soit 2 FIFO de N octets. Une ligne doit �tre constament gard�e en mat�riel afin de simplifier le logiciel et les pointeurs. La deuxi�me se vide au fur et � mesure que le calcul avance.

        Le calcul donne 4 bits pour la premi�re ligne, 8 bits � la suivante et les 4 derniers pour la troisi�me ligne. De plus, de petits buffers doivent �tre correctement plac�s pour temporiser certains r�sultats de un ou deux cycles. On modifie donc la structure de la figure a) pour exploiter des FIFO sur 8 ou 16 bits : on utilise des bits ind�pendants dans la m�me FIFO puisque tous les flux sont synchrones. Il faut donc trois �tages de FIFOs � N �tages, ou en tout 4 FIFOs de N octets. Dans la versions a), on s'aper�oit aussi que le calcul s'effectue � chaque lecture d'un mot en sortie : on peut donc �conomiser un �tage de FIFO. Les parties d�licates sont la programmation des LUT, la gestion du "m�lange" des bits et les buffers associ�s, ainsi que le contr�le g�n�ral de la carte, non repr�sent� ici. Le g�n�rateur de nombres al�atoires peut �tre impl�ment� dans une simple PAL avec un registre � d�calage � r�troaction lin�aire.

        Une autre optimisation porte sur l'utilisation des circuits : on peut r�duire le nombre de circuits buffers 8 bits (des bascules latch non transparentes) � seulement deux en traitant les bits ind�pendamment. Certains bits peuvent donc �tre utilis�s � diff�rents degr�s de d�lai et il n'y a pas de broche inutilis�e. Ainsi, quatre circuits suffisent � "brasser" les bits, le reste du circuit imprim� doit g�rer les autres fonctions : acc�s aux LUT pour leur configuration, d�codage du bus ISA, synchronisation et validation des horloges de tous les circuits. Ce sont ces parties qui sont les plus complexes et qui sont susceptibles de demander du temps pour la mise au point.

        L'architecture d�crite dans ce chapitre n'a pas �t� impl�ment�e pour des raisons financi�res, de temps, de moyens mais surtout parce que les programmes qui ont �t� d�velopp�s ensuite sont beaucoup plus rapides que les 12Mc/s th�oriques que l'acc�l�rateur permettait d'atteindre. La version � base de SRAM pour �muler les FIFOs (ch�res, voir RadioSpares) est plus compliqu�e car il faut adresser des SRAM : les composants discrets sont trop nombreux, les PAL sont trop petites et les FPGA d�passent le budget.

V.7 : Anneau de strip mining :

        L'id�e de cette architecture est venue dans une discussion, � la fin d'une pr�sentation par Amal Stri du projet Fourmi en 1998. C'est une application directe du principe de "strip mining" utilis� dans le programme : il y a deux niveaux de m�moire de vitesse et de quantit� diff�rentes afin de diminuer le co�t total du syst�me. Une m�moire de type DRAM est � la fois spacieuse et �conomique, elle est ici contr�l�e par une puce sp�ciale qui effectue le rafra�chissement et le transfert de blocs, alternant lectures et �critures en rafale. Ensuite, chaque "ligne" de strip mining est impl�ment�e par un circuit qui peut �tre r�pliqu� � volont� selon les besoins et le budget. Il peut consister en un gros FPGA ou tout simplement correspondre au circuit d�crit dans le chapitre pr�c�dent.

        Ce type d'architecture a un int�r�t particulier : plus le nombre de circuits de calcul de lignes est grand, plus le calcul est rapide, ind�pendamment de la taille du tableau � calculer. La partie calcul est s�par�e de la partie stockage et chacune peut avoir la taille d�sir�e. Il devient possible d'augmenter le nombre total de sites calculables en ajoutant une barette de DRAM disponible dans le commerce, ou d'augmenter la vitesse de calcul en fabricant d'autres cartes de calcul. Le contr�leur de DRAM peut aussi g�rer des communications plus lentes afin d'inclure l'anneau dans un anneau plus grand et impl�menter une architecture de strip mining du deuxi�me ordre (avec un anneau d'anneaux).

        L'affichage du r�sultat peut �tre effectu� en lisant le trafic sur un des liens et en envoyant le flux vers un tampon d'affichage pour le synchroniser avec le balayage de l'�cran. Il faut donc que le flux de l'anneau contienne un "jeton" indiquant un retour de balayage vertical mais il est plus simple de g�rer ce cas directement avec le contr�leur de m�moire : il doit contenir les pointeurs de d�but et de fin du balayage de la DRAM et peut donc g�n�rer le signal de synchronisation sans complexit� inutile.

        Ce type de syst�me parall�le exploite les caract�ristiques du mod�le calcul� et son extensibilit� est diff�rente compar�e aux autres syst�mes cellulaires. D'habitude la puissance est augment�e en ajoutant des modules qui contiennent � la fois de la m�moire et des circuits de calcul mais le prix de la DRAM au m�gaoctet est beaucoup plus faible que de la SRAM de petite quantit�. L'acc�s par blocs dans ce mod�le est aussi un facteur qui rend cette architecture possible : des acc�s al�atoires p�naliseraient la bande passante � cause de l'interface compliqu�e des puces de DRAM. Ce syst�me peut donc �tre �tendu par ajout de composants m�moire �conomiques et par ajout de modules de calcul (circuits imprim�s ou FPGA).

        Le probl�me de cette approche, malgr� le rapport performance/prix int�ressant, est qu'il est limit� aux mod�les � interactions courtes en 2D (voisinage de Moore ou de Von Neumann par exemple). Il ne peut �tre utilis� pour des mod�les uni- ou tridimensionnels et l'int�r�t est r�duit en dehors des gaz sur r�seaux : ce n'est pas une architecture "g�n�raliste". Si le domaine � �tudier est limit� aux LGA 2D, c'est l'architecture parall�le la plus recommend�e, m�me avec des mod�les non binaires : ILG, BGK...

V.8 : tableau de PAL :

        L'exploration suivante est un "exercice d'�cole" destin� � comprendre � quel point le mod�le physique influence l'architecture mat�rielle. Les contraintes de prix, de flexibilit� et de r�alisation sont ensuites utilis�es pour modifier les choix.

        Commen�ons par un cas de figure "id�al" car il repr�sente exactement le mod�le FHP : chaque site est associ� � un circuit logique programmable, par exemple une PAL 16L8.

   \   /   \   /   \   /   \   /   \   /   \   /   \   /   \
--- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL
   /   \   /   \   /   \   /   \   /   \   /   \   /   \   /
PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL ---
   \   /   \   /   \   /   \   /   \   /   \   /   \   /   \
--- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL
   /   \   /   \   /   \   /   \   /   \   /   \   /   \   /
PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL ---
   \   /   \   /   \   /   \   /   \   /   \   /   \   /   \
--- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL
   /   \   /   \   /   \   /   \   /   \   /   \   /   \   /
PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL ---
   \   /   \   /   \   /   \   /   \   /   \   /   \   /   \
--- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL
   /   \   /   \   /   \   /   \   /   \   /   \   /   \   /
PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL --- PAL ---

     Chaque lien correspond � deux connexions (un bit dans chaque sens). Ce syst�me pourrait fonctionner � 20MHz : le parall�lisme massif permet d'atteindre environ 1Gc/s avec la configuration de circuits ci-dessus (un carr� de 8 par 8). La performance dans des applications r�elles (500*500 environ) serait donc impressionnante mais les obstacles pratiques sont tr�s nombreux :

        Finalement, cette approche est r�serv�e � des cas tr�s particuliers o� l'investissement en mat�riel est justifi� par une utilisation intensive, hors du cadre de la recherche. Une puce VLSI avec un r�seau de 16 * 16 sites pourrait servir de base � des g�n�rateurs de nombres pseudo-al�atoires (de probabilit� 2/7) contr�l�s par cl�.

        Une alternative serait d'utiliser des m�moires SRAM rapides comme celles utilis�es dans les caches L2 des PC. Elles sont relativement �conomiques, largement r�pandues et ne sont pas limit�es par la complexit� de l'�quation. Par exemple, partons de mon stock personnel : environ 100 SRAM 62256 15ns en bo�tier DIL �troit 28 broches, toutes r�cup�r�es sur des cartes m�res d�fectueuses. Comme leur contenu peut �tre chang� � tout moment, le syst�me simul� peut �tre int�ractif. Une SRAM de 32 Koctets contient plus de donn�es que n�cessaire pour une seule table (512 octets), un bo�tier peut donc en plus contenir les informations de 64 lignes diff�rentes pour lesquelles les informations ne seraient pas homog�nes (comme pour des parois plus ou moins glissantes, de informations sur la g�om�trie...). Il faudrait alors que l'h�te optimise l'allocation de ces sous-tables mais ce n'est pas un probl�me inqui�tant. La m�moire est organis�e par octets et les murs peuvent �tre trait�s s�par�ment, il reste donc un bit libre dans les donn�es, utilisable pour impl�menter FHP4 :

        Il n'y a pas assez de puces pour constituer un tableau entier (un tableau de 10 par 10 n'est pas int�ressant pour des applications normales), il faut donc balayer la surface ligne par ligne. L� encore, une ligne de 100 sites n'est pas assez pour des cas classiques, il faudra donc plusieurs cycles pour balayer une ligne r�elle. L'organisation des circuits est toutefois 1D : tous les circuits sont plac�s en ligne. Cela pose un gros probl�me : il faut m�moriser les autres lignes, une moiti� des puces m�moire sera alors d�di�e � cette fonction et il faudra les adresser s�par�ment. Le tableau peut donc contenir 50 x 32K = 1,6 millions de sites qui peuvent �tre organis�s par groupes de 50 sites. Une autre solution serait d'acqu�rir 100 puces FIFO de plusieurs milliers d'octets chacune pour simplifier l'adressage et le c�blage mais ce sont des pi�ces ch�res. Le "brassage" des bits pour les r�partir sur plusieurs ligne cons�cutives peut s'effectuer comme avec la carte ISA, avec des buffers discrets.

        Le probl�me le plus inqui�tant est la programmation les SRAM : il faudrait intercaler des transmetteurs qui s�lectionnent la provenance des donn�es dans le chemin critique, ce qui diminuerait la vitesse, augmenterait le co�t et la consommation �lectrique. De toute fa�on, la sortie de ces puces m�moires n'est pas "latch�e" (la sortie ne peut �tre maintenue lorsque que l'entr�e change). Un syst�me synchrone rapide n'est donc pas possible car il faut en plus alterner les cycles de lecture et d'�criture des SRAM de stockage. On peut consid�rer que la vitesse de consultation des LUT est de 15 ou 20 MHz, soit environ 1 milliard de sites par seconde avec 100 puces de SRAM. La consommation �lectrique est au minimum de : 100mA sous 5V x 100 = 50 Watts pour les SRAM seules.

        La question reste en suspens, mais il est �videmment plus facile de programmer un ordinateur "classique" que de mettre au point un syst�me de ce style, m�me avec des FPGA modernes, en maintenant des performances et une flexibilit� similaires.

        Il appara�t que la cr�ation d'une machine d�di�e au calcul de FHP n'est pas aussi ais�e que le mod�le le laisse penser au d�part. L'�lectronique "c�bl�e" reste pourtant le dernier recours lorsque la vitesse demand�e justifie le manque de flexibilit� et le budget consacr� au projet. L'approche parall�le en 1D ou en 2D permet de gagner un ordre de grandeur dans la vitesse de calcul. Une machine � base de dizaines d'ASIC est une solution brute � un probl�me tr�s pr�cis mais la nature des circuits int�gr�s implique une perte de vitesse d�s qu'il faut communiquer avec un circuit int�gr� voisin : une broche suivie d'une piste de circuit imprim� constituent une charge capacitive qui devient pr�pond�rante � tr�s haute vitesse. Il serait donc illusoire de croire � une augmentation lin�aire de la vitesse de traitement lorsque l'horloge est acc�l�r�e : un syst�me synchrone bas� sur l'architecture PAL ou sur un tableau d'ASIC sera limit� par la vitesse de communication entre les puces. Malgr� ses avantages, le calcul "c�bl�" sera donc toujours limit� par les constantes fondamentales de la physique (comme la vitesse de la lumi�re ou la mobilit� des �lectrons).

V.9 : Conclusion :

        Il semble que l'avenir de ce type de machines se trouve dans les r�seaux logiques reconfigurables sur site : les FPGA. Les nouvelles g�n�rations sont tr�s, tr�s rapides et peuvent m�me contenir des microprocesseurs. Leurs architectures, comme pour les microprocesseurs, deviennent de plus en plus complexes et il est difficile de les utiliser � 90% de leur capacit� th�orique, sans parler du manque de transparence des constructeurs, du prix �lev� des logiciels et de plateformes... Le projet Fourmi a explor� ces diff�rents probl�mes.

        Cependant les FPGA offrent un niveau d'abstraction qui permet de faire b�n�ficier l'utilisateur des performances d'une nouvelle plateforme au prix d'une recompilation du source VHDL ou Verilog. Des plateformes diverses existent d�j�, allant de la simple carte PCI aux syst�mes massivement parall�les. Ces syst�mes reconfigurables se d�mocratisent (lentement mais s�rement) car ils touchent beaucoup plus d'applications : ils permettent de calculer rapidement tout ce qu'un ordinateur classique peut calculer. L'ordinateur "h�te" se charge des t�ches complexes, de l'interface avec l'ext�rieur, de la configuration, pendant que la partie de "calcul brut" est prise en charge par le syst�me reconfigurable.

        Il faut toutefois remarquer qu'aucun standard n'existe actuellement. On peut sp�culer sur l'utilisation du VHDL pour programmer des automates cellulaires, Cellang n'a qu'un usage tr�s restreint et le reste reste au stade exp�rimental. Dans l'absence d'un outil d�terministe et fiable pour programmer les FPGA, le pr�cieux mat�riel risque d'�tre sous-utilis� comme un ordinateur classique. Doit-on en venir aux ASIC ? C'est ce que pr�pare Norman Margolus [30] avec une puce m�langeant circuits de m�moire et de calcul, b�n�ficiant ainsi des toutes derni�res technologies.

        Les r�seaux d'interconnexion entre les processeurs �l�mentaires sont tr�s divers. Pour calculer des gaz sur r�seaux en 2D, la m�thode la plus pratique, �conomique et extensible est l'anneau d�crit au chapitre 6. Les autres techniques semblent soumises aux caprices des �volutions technologiques. Les calculs en 3D sont moins faciles � c�bler de mani�re flexible, il faut encore attendre que la technologie �volue. Un hypercube 4D semblerait convenir mais la seule exp�rience effectu�e a utilis� une CAM-8. Les architectures adapt�es aux r�seaux en 2D ne peuvent pas �tre �tendues directement � la 3D car c'est une projection d'un r�seau 4D. De nombreuses recherches suppl�mentaires doivent �tre effectu�es pour d�terminer une architecture adapt�e mais les techniques futures seront certainement extrapol�es � partir des techniques 2D comme l'anneau (issue de l'algorithme de strip mining).

1.1 : Approche globale du problème :

        Programmer des LGA ne se limite pas en l'expression d'un algorithme dans un langage particulier : le premier problème à surpasser avant de concevoir un programme efficace est de s'affranchir du niveau conceptuel des langages textuels. Bien que le premier charme des LGA soit d'être programmables "facilement", on vérifie aujourd'hui que les algorithmes "naïfs" ne sont plus adaptés aux ordinateurs actuels (à partir de 1995 environ) et nécessitent une connaissance approfondie du modèle physique et de l'architecture de la machine.

        Pour programmer les LGA, il faut partir du modèle physique exprimé abstraitement, examiner de nombreuses options et comparer leurs intéractions dans un cadre choisi. L'interface entre le logiciel et le matériel, qui permet de contr�ler la performance du programme, s'étend sur le système d'exploitation, sur l'algorithme, le modèle physique, l'application, la machine cible... Ce sont autant de détails, parfois contradictoires, qu'il faut considérer attentivement, longtemps avant de fixer des choix définitifs ou commencer à programmer.

        La philosophie de programmation joue aussi un r�le important : selon les objectifs, les moyens et les connaissances, des niveaux différents de performance pourront être atteints. Les objectifs doivent être réalistes et les moyens doivent être préalablement maitrisés, sinon le projet risque de ne pas aboutir. Il ne faut pas hésiter à constituer une large base de connaissances (bibliographie, discussions, rencontres, butinage sur Internet) afin de maitriser le sujet dans ses nombreux détails.

        Lorsque suffisament d'éléments sont réunis et préparés, l'intransigeance et le souci du détail permettent de consolider l'édifice, qui ne manque pas de montrer des signes de faiblesse à mesure qu'il s'agrandit. La cohérence du programme devient de plus en plus difficile à gérer lorsqu'on ajoute des éléments : s'il y a N éléments dans un programme, l'ajout d'un autre élément nécessite N vérifications qui peuvent remettre en cause d'autres éléments. Le développement logiciel de ce projet a d'ailleurs été interrompu car il y a trop de parties à réécrire et les outils actuels ne permettent pas de les développer suffisamment rapidement. Ce n'est pas un problème lié au langage assembleur car des intéractions complexes apparaissent dans tous les langages.

        La programmation d'un logiciel optimisé n'est pas seulement un travail de science fondamentale : ce mémoire a fait appel à des compétences en électronique, compilation, système d'exploitation. Il emprunte des méthodes utilisées entre autres par les codeurs de jeux vidéo, qui ont comme préoccupation commune l'utilisation la plus efficace possible de chaque ressource de l'ordinateur. Il faut pourtant faire attention aux parties théoriques car le résultat risque d'être inutile : l'utilisation du modèle FHP en dehors de ses limites raisonnables donne des résultats aberrants. L'utilisation d'un modèle trop primitif semble économique en complexité de calcul mais est sous-efficace en temps de calcul total sur la machine... Un codeur de jeu vidéo sera plus intéressé par l'effet visuel que par l'exactitude d'un calcul.

        Pour résumer, il faut entre autres : être très attentif, être compétent dans de nombreux domaines, pouvoir faire le lien abstrait entre chaque partie, travailler à la fois sur le très haut et le très bas niveau. En réduisant le nombre de contraintes (par exemple : temps du projet trop court, moyens trop réduits, ressources trop faibles) il est possible de faire des programmes de meilleure qualité mais cela est rarement le cas car une contrainte est souvent compensée par une autre.

1.2 : Organisation des données :

        La refonte de la structure des données et des algorithmes (principaux et annexes) est souvent la première voie à explorer pour améliorer fondamentalement un programme. Dans le contexte d'une approche globale, c'est à dire en tenant compte des détails d'implémentation, il faut essayer et comparer les différentes possibilités qui existent. Lors de ce projet, quatre organisations ont été examinées et on peut résumer leurs propriétés dans le tableau suivant :

type	avantage(s)	inconvénient(s)
tableau 2D d'octets (LUT, multisite) (chapitre III.4)	simplicité de progammation/compréhension, adapté pour les ordinateurs 8/16 bits (ex.: i286)	sous-efficace pour les microprocesseurs récents (mots large, coeurs OOO, occupation de la cache...)
tableau 2D de mots / 4 sites sur 32 bits (multisite à traitement parallèle) (chapitre III.6)	plus efficace pour les coeurs 32 bits (i386-486-P53C), moins d'instructions pour le mouvement de bits individuels	plus de manipulations d'octets individuels, donc plus contraignants pour les coeurs récents
multispin entrelacé dans un tableau (chapitre IV.4)	convient le mieux aux microprocesseurs modernes (registres larges, OOO, cache sur la puce)	complexe (mais ce m�moire montre que c'est possible)
multispin (equation booléenne) sur plusieurs tableaux séparés	convient naturellement aux calculateurs vectoriels	nécessite trop de pointeurs (pression sur les registres et le compilateur pour les processeurs classiques, mauvaise localité spatiale, risque de cache thrashing avec certaines granularités)

        Puisque l'étude porte sur les microprocesseurs modernes (Pentium MMX et Pentium II) il est donc recommandé d'utiliser la technique "multispin entrelacé" sur ces plateformes. Les autres techniques sont sous-efficaces, mettent une pression inutile sur l'allocation des registres (pointeurs et données) et utilisent mal la hiérarchie de la mémoire.

        De plus, il est fortement recommandé d'utiliser des buffers temporaires comme décrits dans au chapitre III.3, afin d'éviter une occupation trop importante. Même si cela peut paraître peu important sur des ordinateurs actuels disposant d'au moins 64MO de mémoire, l'espace supplémentaire nécessaire à la collecte des statisques compense vite cette économie. De plus, le léger surco�t en complexité de programmation est justifié par une meilleure utilisation de la bande passante vers la mémoire. Par exemple, la stratégie de write back du Pentium est no-write-allocate : elle court-circuite la cache si la zone mémoire n'est pas déjà cachée, ce qui pénalise certains algorithmes selon leur ordre de lecture/écriture.

1.3 : Strip mining :

        L'utilisation du strip mining dans le projet permet de gagner entre 25 et 50% du temps de calcul selon la taille des tableaux, alors même que le coeur du CPU est saturé par les calculs. Plus qu'une simple "astuce", cela montre bien que les microprocesseurs sont très dépendants de la localité spatio-temporelle des données car toute leur hiérarchie mémoire est conçue selon ce postulat optimiste.

        Les premiers microprocesseurs ne disposaient pas de mémoire cache, et les machines dédiées ont des canaux spécialisés qui garantissent une certaine bande passante, mais les microprocesseurs actuels ont une bande passante vers la mémoire centrale qui est très réduite par rapport à leur fréquence de fonctionnement interne. A mesure que les architectures deviennent plus complexes, de nombreux paramètres rendent le travail plus difficile : il faut faire attention à la cachabilité des zones accédées, aux temps de latence, aux transferts par paquets, � l'ordre et � l'alignement des accès...

        En concevant un algorithme, il faut donc absolument tenir compte de ces paramètres, avec la difficulté supplémentaire de la nature des traitements à effectuer. En effet, la fenêtre de strip mining doit respecter certaines propriétés spatio-temporelles (qui étaient mal identifiées au début du projet). D'autres types d'automates cellulaires, des traitements de signaux (images, sons) ou des opérations mathématiques complexes (inversions de matrices) ont des dépendances différentes entre les données et nécessiteront un algorithme de strip mining différent, qu'il faudra analyser sérieusement pour éviter que la mémoire centrale ne ralentisse le processeur. En règle générale, on peut considérer que tout balayage d'un tableau avec des intéractions limitées aux voisins peut et doit bénéficier du strip mining.

        Il faut espérer que les prochaines générations de microprocesseurs banalisent encore plus les instructions dédiées au contr�le de la mémoire cache. Le strip mining présenté ici est au départ un palliatif au manque d'instruction spéciale et fonctionne par "cache touching" (une partie de la ligne de cache est accédée, le reste est ensuite utilisé au maximum et on attend que le mécanisme de LRU vide la ligne automatiquement). Une amélioration simple consisterait à effectuer le travail de LRU grâce à une instruction explicite et gagner ainsi quelques pourcents sur l'espace réellement utilisé dans la cache (pour être remplacée, une ligne doit être inutilisée, elle prend donc de la place inutilement) et donc accélérer la vitesse du programme. Jusqu'à maintenant, la seule instruction disponible (WBINVD) vidait indifféremment toute la cache et ne procurait aucun bénéfice net (voir les tests en annexe C). Intel a amélioré la situation lors de l'introduction du Pentium III et des instructions SSE, on peut donc espérer que les prochaines versions du programme en bénéficieront.

        Le Pentium II avait déjà changé le "paysage" de la mémoire cache par une architecture à bus séparés très efficaces et rendait ainsi la L2 presque aussi "rapide" que la L1. La mémoire principale n'avait pourtant pas été accélérée mais le schéma original de strip mining fonctionne très bien. Les choses se compliquent dans les cas où la taille des simulations nécessite l'emploi de liens encore plus lents : réseau local Ethernet ou disque dur. Il faut alors utiliser le strip mining au deuxième degré : la complexité et le gain de cette technique éprouvée au niveau de la carte mère peuvent être retrouvés au niveau d'un système à une échelle différente. Des études sont en cours pour le projet beo-kragen.

        D'une simple réorganisation des ordres d'accès à la mémoire, le point de vue a évolué vers la gestion de tampons mémoire et de flux de données à l'échelle du système entier, prouvant que le programme est actuellement étudié à un niveau d'abstraction beaucoup plus élevé qu'à l'origine.

1.4 : Parois complexes :

        L'utilisation de listes de modifications pour gérer les parois permet de bénéficier de toute la puissance intrinsèque du modèle FHP. Bien que cela soit plus difficile à programmer, les efforts sont justifiés par une grande flexibilité des parois tout en ayant une occupation raisonnable en temps CPU et en mémoire. Le projet a permis de programmer et de résoudre les couches de bas niveau, le reste de l'effort étant laissé à la discrétion d'autres programmeurs spécialisés dans d'autres disciplines (l'algorithme de Bresenham dans un repère rhombohédrique ne faisant pas l'objet de ce mémoire). Ce projet de maitrise prouve que les listes de modifications sont possibles et montre comment les réaliser (bien qu'on puisse toujours améliorer la technique). Il existe aujourd'hui moins de raisons de s'en passser.

1.5 : Programmation en assembleur :

        La pratique de l'assembleur est souvent vue comme un exercice de haute voltige, inutile et peu pratique. Ce mémoire montre que ce n'est qu'une des nombreuses difficultés que l'on rencontre dans les exercices de programmation courante car la plus grande difficulté est souvent au niveau algorithmique (le coeur du calcul utilise un nombre très réduit de types d'instructions, ce qui invalide l'argument de complexité). Le reste de la "difficulté" est au niveau de l'interface, ce qui peut être résolu par de nombreux moyens autres que celui utilisé ici.

        L'utilisation de l'assembleur ici correspond à une philosophie de contr�le total sur les instructions que le microprocesseur va exécuter, et donc maitriser étroitement la performance du programme, quasiment au cycle près. Aucun autre moyen n'existe actuellement pour garantir la "pureté" d'un code exécutable. Les meilleurs compilateurs du monde ne permettent qu'un contr�le limité sur le code généré, malgré l'existence de centaines d'options de compilation : il manque toujours celle dont on a besoin.

        Non seulement les compilateurs actuels ne permettent aucun contr�le du code à l'instruction près, mais il ne savent pas encore gérer automatiquement les instructions MMX ou SSE qui sont nécessaires pour accélérer les calculs avec les algorithmes multispin. L'utilisation d'un compilateur pour du code "final" est donc un effort superflu, car il ajoute un niveau d'abstraction parasite dans l'analyse du problème. Les espoirs de générations automatique de codes "preque parfaits" se sont évanouis devant l'ampleur du manque de réelle intelligence des outils. Ecrire le coeur du calcul à la main, à l'aide d'un simple outil disponible gratuitement sous GPL, s'est révélé finalement beaucoup plus efficace et plus facile que de forcer les compilateurs à créer du code qu'ils ne peuvent pas générer. De plus, cela a nécessité un important effort au niveau de la théorie derrière les calculs de collision, qui a permis la création d'une nouvelle formule.

        Le code actuel a atteint un point critique qui nécessite de nombreux efforts pour le dépasser. Dans l'état actuel du projet, il est alors rentable de se concentrer sur les outils de codage "assisté", dans lesquels les techniques de codage développées pour ce projet peuvent être réinvesties. Le projet "GNL" permettra de recoder entièrement le source du projet, de le porter vers d'autres architectures, tout en permettant d'atteindre et de se maintenir à la puissance de crête de la plateforme cible. Si cet outil interactif avait existé au début du projet, le programme serait fonctionnel plus rapidement et plus facilement. Beaucoup de "travail sur papier", d'efforts et de temps aurait été économisé. Toutefois, le "travail sur papier"effectué pour ce projet a permis de mettre au point et tester en grandeur réelle des techniques à la base de GNL (allocation des registres, évaluation des accès à la mémoire etc).

1.6 : Une nouvelle formule :

        L'une des découvertes les plus inattendues de ce projet est l'équation booléenne au coeur du calcul des collisions. C'est aussi un point très important du mémoire car il découle d'une analyse plus aprofondie et sous un autre angle que la formule (classique aujourd'hui) donn�e au chapitre IV.7.

        La nouvelle formule essaie de correspondre aux exigences du Pentium : peu de registres, bande passante réduite, pairage des instructions compliqué... alors que la formule classique nécessite de nombreux termes temporaires dans le graphe de dépendences, favorisant les architectures RISC avec de très nombreux registres. En comparaison, la nouvelle formule peut être utilisée sur des calculateurs vectoriels (par exemple CRAY classiques avec 8 registres vectoriels).

        La nouvelle formule ne remplace pas l'ancienne, nous pouvons remarquer qu'elle la complémente et élargit le domaine d'application et d'utilisation du modèle FHP saturé. En l'absence d'analyse booléenne plus poussée, il est difficile et peut-être inefficace d'utiliser la formule de d'Humières si moins de 32 registres sont disponibles pour les données (c'est à dire que les pointeurs, compteurs de boucles etc. doivent �tre séparés ou très peu nombreux). Au contraire, si le rapport mémoire/calcul d'une machine particulière efface le problème des variables temporaires d'une manière ou d'une autre, la formule de d'Humières est plus efficace car le nombre total d'opérations booléennes est inférieur à la nouvelle formule.

        Il reste encore à prouver que la vielle formule est optimale, ce qui est une tache difficile en raison de la spécialisation de chaque cas : quelles opérations booléennes sont admises ? XOR, ANDN, NAND, ORN ? La formule peut changer complétement en ajoutant ou en supprimant un opérateur et la pratique montre que les compilateurs VHDL sont incapables de répondre à cette question de manière satisfaisante. L'analyse brute est inefficace, il faudrait donc trouver une nouvelle méthode pour "casser" ce type de formule lourde en instructions simples.

        Loin d'avoir résolu le problème, ce travail relance la question épineuse de l'efficacité des calculs. Il apporte un nouvel élément (une autre formule) ainsi qu'une nouvelle manière d'appréhender la question des collision saturées. Le long travail qui leur a permis de voir le jour n'aurait pas été nécessaire si les microprocesseurs courants avaient plus de bande passante vers la mémoire, plus de registres, et des outils plus efficaces pour "casser" le graphe de dépendances de données de manière satisfaisante.

        Cependant, la formule de d'Humières ressemble un peu à l'algorithme de Bresenham dans le sens où avant son existence, le problème était considéré comme "difficile" (peu de personnes avaient envie de transformer le modèle en code informatique, en raison de sa complexité qui intéressait peu les physiciens). Lorsque le code a été publié, il a été réutilisé abondament et le problème n'a pas ressurgi car le code était estimé "satisfaisant" (quand il était correct...). Aujourd'hui, comme l'algorithme de Bresenham, la formule "canonique" ne satisfait plus les contraintes complexes imposées par certaines classes de machines. Espérons que d'autres personnes continueront les travaux dans ce domaine après ce mémoire.

1.7 : Intégration complète de l'algorithme de visualisation dans celui du calcul :

        Tout comme pour le strip mining, intégrer les calculs de visualisation dans la boucle de calcul permet de bénéficier des propriétés de localité spatio-temporelles des modèles FHP et similaires. Les microprocesseurs comme le Pentium, limités en bande passante vers la mémoire externe, bénéficient de ce type de programme car ils évitent ainsi de saturer inutilement le bus externe par des cache miss à répétition. L'intégration des différents algorithmes, lorsqu'ils portent sur une même donnée, permet de rééquilibrer les caractères memory bound et CPU bound et d'entrelacer la latence des mémoires avec les lourds calculs. Les coeurs OOO (comme le PII ou l'Alpha 21264) ont principalement été conçus dans cette optique, l'exécution spéculative de dizaines d'instructions permet de continuer le programme en attendant une donnée venant de la mémoire centrale.

        Lors de la conception d'un programme de calcul intensif, il est donc important de bien faire cohabiter la partie de calcul brut et les parties qui ne participent pas au calcul proprement dit : bien que le résultat ne change pas, il devient inutile s'il n'est pas présenté correctement à l'utilisateur. Le problème se complique avec l'introduction du strip mining, mais heureusement il se résoud naturellement dans notre cas précis, ce qui peut ne pas être vrai pour d'autres cas.

        Le travail n'a pas abordé la collecte de statistiques, se limitant à la mesure de la densité en particules, mais la remarque reste valable : il faut autant que possible réunir toutes les parties du code autour des m�mes données, afin de bénéficier des mémoires caches de données. Le code exécutable tient souvent largement dans la cache des instructions, la taille du code n'est pas actuellement un problème.

1.8 : Mesurer son code :

        C'est un détail qui est tout le temps oublié, mais même sans le langage assembleur, même si peu de moyens sont disponibles, il est toujours important de connaitre à tout moment le temps d'exécution de chaque bloc d'instructions. Il faut toujours "profiler" son code et comparer ses performances entre chaque modifications.

        Tous les détails qui rendent ces éléments cohérents sont très importants, car comme noté précédemment, l'ajout d'un détail peut remettre en cause le programme entier. La plus grande prudence est donc conseillée. Un de mes mots d'ordres est "coder proprement paie toujours" (proprement dans le sens de l'efficacité) car la mentalité actuelle dans le milieu informatique est que "les ordinateurs sont bien assez rapides comme ça, pourquoi se compliquer la tâche ?". Cela conduit à des situations où des programmes tournent � 10% de la puissance nominale de l'ordinateur, sans que personne ne s'inquiète. Or selon le cas et avec du logiciel optimisé, on pourrait économiser 90% du prix d'achat sur le matériel, à performance égale. Cela semble inaproprié pour les calculs de ce mémoire, mais devient très important pour un serveur central ou institutionnel co�tant des millions de francs.

2.1 : FHP-3 est memory bound et CPU bound sur x86 :

        Les mesures ont montré que l'on ne peut plus améliorer radicalement les performances de FHP-3, on se heurte toujours à une partie qui dépend trop des calculs ou de la mémoire. Le cas du x86 est dramatique car ni la mémoire ni les instructions ne sont conçues pour supporter et maintenir leur puissance théorique. Le cadre des applications bureautiques et ludiques, permettant l'expansion du marché de masse des particuliers, ne justifie pas une refonte ou un abandon de la vieille architecture du x86, même si Intel essaie de promouvoir l'IA64 pour les domaines où la recompilation est déjà une nécessité. L'IA64 réunit probablement le pire du VLIW et du SPARC en essayant de faire mieux que le x86 : nous ne sommes pas prêts de voir une architecture "sympatique" pour FHP dans le commerce avant longtemps.

2.2 : Il est possible de descendre à 3 cycles par octets :

        Malgré une substancielle augmentation de la taille du kernel de calcul, par rapport au code multisite 32 bits étudié dans le passé, il est possible d'aller encore plus vite et d'être plus flexible grâce au code multispin entrelacé. Il est peu probable que l'on puisse descendre beaucoup plus bas, par exemple 1 cycle par cellule, car la mémoire freine le processeur de plus en plus si de très longs mots (128 bits ou plus) sont utilisés. De nombreuses améliorations peuvent toutefois êtres effectuées et ne manqueront pas de voir le jour. Les codes multisites sont donc bien morts.

2.3 : A condition d'effectuer de nombreux efforts, il est possible d'utiliser des PC pour des calculs lourds :

        Comme ce mémoire le prouve, il est possible de compter sur le rapport performance/prix/disponibilité intéressant du PC si de nombreux points sont pris en compte et traités :

• il faut pouvoir remettre en question de nombreux points acquis : algorithmes, structures de données, techniques de développement...
  
• il est de plus en plus difficile de contr�ler TOUS les aspects du calcul car les architectures logicielles et matérielles sont de plus en plus complexes
  
• le résultat sera souvent proportionnel à l'effort, il faur donc pouvoir accorder beaucoup de temps et user sa patience sur des problèmes parfois incompréhensibles
  

2.4 : Un PC équivaut presque à une CAM8 :

        Les mesures effectuées au MIT en janvier 2000 montrent que les PC de bureau de dernière génération sont presque aussi rapides qu'un bloc CAM8 (8 cartes à 25MHz). Les tout derniers microprocesseurs généralistes permettent de rivaliser avec des ASIC créés il y a plusieurs années. En terme de génération équivalente, si l'on considère la règle de Moore, l'optimisation poussée du code a permis probablement de gagner trois ou quatre années par rapport à un code non optimisé. Le code est 4 fois plus rapide que le plus rapide des codes testés, ce qui permet d'affirmer que l'effort a permis de gagner 3 ans. Ce gain permet d'utiliser une machine plus vieille à vitesse égale (donc moins ch�re) ou bien de gagner 3 ans sur la machine la plus récente. Cet aspect d'économie est valable si le code original était "bâclé", mais reste dans le cadre de la démonstration du fait qu'un codage consciencieux n'est pas une perte de temps à longue échéance.

2.5 : La loi de Moore est trompeuse :

        La règle de Moore, découverte par le co-fondateur d'Intel dans les années 70, ne signifie que ce qu'elle dit : que le nombre de transistors est multiplié approximativement par quatre tous les trois ans.

        Il est donc abusif de considérer qu'un programme fonctionnera quatre fois plus vite dans quatre ans. Les grands constructeurs pr�parent de nouveaux coeurs de plus en plus étranges et de moins en moins adaptés aux algorithmes actuels. Les détails architecturaux deviennent de plus en plus complexes et il est de plus en plus difficile de tous les prendre en compte lors de la conception d'un programme. Cela veut aussi dire en filigrane que pour tourner quatre fois plus vite dans quatre ans, il faudra complétement reprendre la conception du programme à partir de zéro. Les efforts à fournir afin d'accélérer un programme deviennent de plus en plus grands, à la mesure des efforts fournis pour augmenter le nombre de transistors sur les puces. Les compilateurs ne seront plus assez sophistiqu�s et il faut d'autres moyens pour programmer. Enfin, comme de nombreux exemples le prouvent, l'histoire des ordinateurs ne suit pas une courbe monotone sur du papier à échelle semi-logarithmique : de nombreuses révolutions nous attendent et personne ne pourra plus utiliser ses sources en C écrits il y a dix ans. Si l'avenir de l'informatique est garanti pour les vingt prochaines années, les chemins empruntés sont inconnus et il faut se pr�parer maintenant à la concrétisation de projets incroyables aujourd'hui.

2.6 : Le sujet des gaz sur réseaux booléens est loin d'être tari :

        Alors que les études actuelles portent sur des modèles à virgule flottante sur des approximations de Bolzmann, la présente étude du modèle FHP est loin d'épuiser toutes les ficelles des informaticiens et des physiciens. FHP reste un "terrain de jeux" incontournable pour comprendre de nombreux problèmes et phénomènes de mécanique des fluides, de mécanique statistique, d'algorithmique, il est donc important de continuer les études dans ce domaine, non pour faire de FHP un outil, mais comme modèle à part entière, pour être étudié et approfondi (le manque d'études dans ce domaine étant le plus souvent d� à un manque de connaissance du sujet, décourageant les étudiants). Un autre piège est la simplicité apparente du modèle, qui séduit les débutants mais qui les perd ensuite. Il est certain que d'autres approches sont nécessaires.

        Le programme, même s'il fonctionne, nécessite de nombreuses améliorations plus ou moins complexes, décrites ici :

3.1 : Mesurer la pression :

        Le problème "annexe" le plus important, celui de la mesure de la densité en particules (donc de la pression), a été résolu dans ce mémoire (popcount parallèle réutilisé dans la phase de détection, puis addition semi-parallèle). Pourtant, un aspect intéressant des mesures en soufflerie n'est pas la mesure dans le fluide, mais sur les parois. Il faut donc mettre en place un dispositif "comptant les bits" qui sont réfléchies par une paroi. Cela peut nécessiter une simple "adaptation" des algorithmes des listes de modification, ou un nouveau type d'algorithme, mais le sujet est sérieux : il permettra d'implémenter un jour des "intéractions fluide/solide", comme la mise en mouvement d'un objet dur par le fluide qui l'entoure (balle dans un courant d'air, vibration d'une surface en fonction de turbulences...). C'est un domaine où les LGA n'ont pas été appliqués, malgré les potentiels certains de ce modèle. L'objection de Norman Margolus à propose des parois mobiles concerne la non-réversibilité du mécanisme, ce qui ne nous intéresse pas dans le cas des écoulements dissipatifs.

3.2 : Emission de particules :

        Le présent programme manque cruellement d'un "générateur de particules" associé à son contraire ("avaleur" ?). Le problème est caractéristique des veines de soufleries : il faut créer un vent uniforme, qui puisse imposer une vitesse contr�lable du fluide. Il faut donc "créer" et "supprimer" des particules à certains endroits, à une vitesse particulière, tout en conservant la pression totale (le nombre de particules dans la veine, voir la partie III).

        Il n'y a pas de problème conceptuel notable, mais il faut tout de même le programmer. En attendant, l'utilisateur est obligé d'utiliser un nombre fixe de particules qu'il doit programmer en assembleur : des efforts de programmation subtanciels et nombreux sont encore à fournir. De plus, il faudrait mettre au point un algorithme de mesure du "vent" afin de rétrocontr�ler ce mécanisme.

3.3 : Multi-CPU :

        Le "stub" de bi-processing symétrique a été écrit mais n'est pas utilisé car le reste de l'algorithme n'est pas terminé. Quand ce sera le cas, il faudra adapter légèrement le code (principalement : copier/coller/renommer) puis enlever le code en commentaire (en faisant attention aux bugs qui ont déjà été trouvés).

        Le véritable problème est dans le cadre d'un beowulf car de nombreux autres problèmes devront être résolus. Les prochaines années verront probablement apparaitre des codes d'exemples. Le programme devient encore plus compliqu� lorsque l'espace d'adressage n'est plus partag�.

3.4 : Extension du code et du modèle :

        Le code développé ici est un modèle simple mais suffisament représentatif des problèmes à résoudre dans des cas réels, même avec des modèles différents. Il peut donc être adapté à des modèles "thermiques" (plusieurs vitesses de particules), multiphases (plusieurs "couleurs" de particules), compressibles, avec des géométries différentes (D2Q9 ou autres) et avec plus ou moins de particules immobiles (pour résoudre le problème de l'invariance galliléenne). Un gros effort de codage devra être effectué pour chaque cas. Cependant, cet effort est d'abord de l'ordre informatique, il faut donc que le recodage du projet tienne compte des besoins d'adaptation du programme pour anticiper la résolution des problèmes futurs. Le code actuel, étude de cas pour un sujet particulier, devra encore beaucoup m�rir pour devenir une plateforme d'expérimentation encore plus générale.

3.5 : Langage script :

        La définition et la programmation d'un langage script universel permettra l'automatisation des mesures et la description des objets à l'intérieur du domaine d'étude. Actuellement, la définition de la g�om�trie des objets est inclue dans le source en assembleur, ce qui impose un réassemblage à chaque changement de géométrie. Les opérations interactives (changement de taille du tunnel, pas à pas...) ne sont pas absolument précis et nécessite de programmer des mécanismes pour contr�ler plus finement tous les paramètres.

        Un format de fichiers et des mots-clés sont actuellement en stade de réflexion mais leur implémentation nécessite un effort trop important pour être effectué actuellement : il faut que d'autres problèmes plus importants soient résolus (tout d'abord l'exactitude des collisions).

3.6 : Evolution de la plateforme :

        Les nouvelles instructions "SSE" introduites dans le Pentium III permettent d'envisager un doublement de la performance brute du programme (à fréquence d'horloge égale) avec peu de changements notables du code. Il faudrait tenir compte des tailles doublées des registres (128 bits au lieu de 64), de la configuration des MTRR et des instructions de gestion de la mémoire et des flux. Ces travaux sont déjà possibles mais seront encore plus faciles lorsque tout sera recodé avec GNL, et quand la plateforme (PIII) sera plus largement répandue (et moins chère). C'est donc une échéance de 1 à 2 ans.

Avec des adaptations plus ou moins profodes, le programme FHP peut être utilisé dans de nombreux domaines :

4.1 : Acoustique :

        Déjà utilisés dans ce domaine, les modèles de gaz sur réseaux (FHP, Bolzmann) ne sont pas encore très répandus malgré leur fort potentiel. Cela peut être résumé comme un problème d'"école", les chercheurs en acoustique étudiant les équations (directes) d'acoustique et non de mécanique statistique.

        Les LGA sont très bien adaptés pour visualiser les ondes sonores et permettent des géométries arbitrairement complexes : elles surpassent les méthodes spectrales pour les cas non triviaux. Par exemple, ils sont utilisés pour étudier les turbulences autour des voitures, et donc déterminer et améliorer leur niveau de bruit aérodynamique (leur silence).

        Les LGA peuvent permettre aussi d'effectuer des simulations d'instruments de musique. En l'absence d'interactions fluide/solide, les anche simples, doubles et lipales ne peuvent pas être simulées. Il reste cependant les vibrateurs de type "fl�te" (à biseau) qui peuvent être étudiés en attendant de meilleurs programmes. Ainsi, il sera possible d'étudier un instrument (pour l'instant en 2D) avec un PC suffisamment rapide, et de générer un son de fl�te véritablement réaliste, sans utiliser la moindre méthode de reproduction classique (analyse spectrale, analyser d'impédance du corps résonant, échantillonnage). Les applications en musique et en synthèse directe sont alléchantes. Le niveau de bruit des LGA bool�ens comme FHP est toutefois trop important pour ce type d'applications, les mod�les en virgule flottante sont donc de rigueur pour ce type de probl�mes.

4.2 : Cryptographie :

        Les LGA de type FHP ont des propriétés macroscopiques très intrigantes : contrairement àla majorit� des Automates Cellulaires classiques, ils ont la capacité de "décorr�ler" un état initial et de le transformer en bruit "brownien" de manière fondamentalement non-linéaire. Le nombre de pas de temps dépend des propriétés du "signal" initial (configuration initiale des particules). Les gaz sur réseaux booléens sont donc très intéressants dans le domaine de la cryptographie, car ils sont potentiellement réversibles, ils sont peu compliqués à calculer, leur état initial est très difficile à retrouver en l'absence des paramètres du calcul et les possibilités d'utilisation sont très nombreuses. Ils peuvent être utilisés comme générateurs de nombres aléatoires ou comme "bruiteurs", avant ou après le calcul. Le brevet des cryptographie par automates cellulaires réversibles ne s'applique pas directement et le MIT travaille sur le sujet. L'efficacité pratique de cette technique reste à démontrer mais elle dépend principalement d'une utilisation judicieuse des LGA, donc l'efficacité potentielle se situe probablement entre DES et RSA. Même si ce n'est probablement pas une révolution pour la cryptographie, c'est un élément important parmi l'arsenal des techniques déjà disponible (courbes elliptiques, nombres premiers, transposition, substitution, permutation...) avec lequel il faudra dorénavant compter. Dans ce cadre, la vitesse de calcul est un élément absolument critique.

        Il faut toutefois rester tr�s prudent dans ce domaine car aucun exemple n'a encore �t� cryptanalys�. Malgr� leurs caract�ristiques non lin�aires, leurs propri�t�s sont �tudi�es depuis longtemps et ils disposent toujours d'une entropie caract�ristique qui permet une attaque par les �quations de Bolzmann.

4.3 : Réalité Virtuelle / jeux vidéo :

        L'accélération de l'algorithme permet d'envisager la simulation des phénomènes turbulents en "temps réel" ou même plus vite, ce qui est toutefois très subjectif selon le cas. Pour les jeux vidéo, où l'exactitude des calculs et la réalité scientifique des résultats importent peu, les gaz sur réseaux offrent une opportunité incomparable pour les mondes simulés. Les techniques courament utilisées sont des simplifications peu rigoureuses et parfois irréalistes (vent dans l'herbe, nuages de "plasma"...). Avec la montée en puissance des consoles de jeu vidéo, l'utilisation de LGA n'est plus qu'une formalité actuellement et pour le futur. Non seulement des phénomènes réalistes peuvent être visualisés, mais en plus les éléments (acteurs, objets) du jeu peuvent intéragir avec le phénomène. Actuellement, aucune application n'est envisagée dans ce domaine. Espérons que ce n'est qu'une "frilosité passagère" : l'effet de mode peut favoriser une utilisation et un développement très actifs dans ce domaine. Il "suffit" juste de lancer correctement l'idée.

4.4 : Education :

        La simulation interactive de phénomènes de mécanique des fluides peut être aussi d'intérêt éducatif. En classes de physique et chimie, la disponibilité d'un "fluide virtuel" non nocif et parfaitement contr�lable peut diminuer les risques de certaines expériences et permet donc aux élèves de manipuler eux-mêmes les produits. En classe de Technologie Industrielle, cela permet aussi de simuler des vérins hydrauliques ou des pompes avec un matériel moins encombrant, afin de tester des algorithmes de contr�le de montée en charge de pompe ou de régulation par exemple.

4.5 : Industrie :

        L'industrie commence timidement à adopter les gaz sur réseaux dans quelques domaines, pour des applications pratiques plus ou moins attendues. Le calcul "industriel" pour les carrosseries et les tuyauteries permet à la société EXA de s'imposer progressivement. De nombreux domaines sont encore en attente d'un miracle que les LGA peuvent produire dans la chimie (réactions d'advection/diffusion), le pétrole (stockage et infiltration dans les sables), la climatisation (où placer la bouche d'aération dans une pièce sans déranger les usagers ?), les peintures (viscosité et écoulements), les risques industriels (souffles d'explosions et écarts entre bâtiments), l'urbanisme (advection des gaz d'échappement), la propagation des flammes et de nombreux problèmes de tous les jours qui ne sont pas encore résolus de manière définitive par les méthodes classiques.

4.6 : Derniers mots sur les applications des LGA :

        Aujourd'hui, les LGA ont fait l'objet de recherches scientifiques poussées dans des domaines où ils apportent un réel avantage par rapport à la dynamique moléculaire, aux techniques spectrales et aux équations classiques : le facteur d�terminant est l'ind�pendance totale entre la complexit� du calcul et la complexit� des g�om�tries �tudi�es.

        Le premier domaine d'adoption concerne les fluides multiphases et la séparation des phases, la tension de surface, les fluides immiscibles, les interfaces solide/liquide/gaz. Ce domaine reste proche de la m�canique statistique. L'image suivante repr�sente une s�paration de phase calcul�e par Rothman et Keller :

Un autre domaine qui b�n�ficie des avantages de LGA est l'infiltration et le mouillage en milieu poreux. Nous voyons sur l'image ci-contre l'infiltration de p�trole (hydrophobe, orange, opaque) rempla�ant de l'eau (bleue) dans un grain de sable de Fontainebleau, calcul�e en 3D par John Olson

        Les LGA se sont montrés inadaptés pour les calculs aéronautiques à cause de la faible vitesse permise (< Mach 0,3) et du faible nombre de Reynolds. La simulation des all�es de Von Karman est facile, comme nous avons pu le voir, mais l'explosion quadratique du temps de calcul rend les gaz sur r�seaux classiques inadapt�s pour des nombres de Reynolds sup�rieurs � 10000. L'exemple d'EXA montre que des ordinateurs massivement parall�les sont utilis�s avec des temps de calcul similaires par les autres techniques. Les seuls avantages sont la plus gande finesse, la r�solution dans le domaine temporel et la meilleure fiabilit� des r�sultat, le temps de calcul n'est plus un crit�re d�terminant.

        Le modèle FHP3 a été plus difficile à implémenter que prévu, sa programmation nécessite des moyens sophistiquée et des connaissances solides qui sont très diffuses dans la bilbiographie classique. La technique prévue au départ pour l'équation booléenne n'a pas fonctionné (compilation VHDL puis GCC puis NASM) pour des raisons d'efficacité. Le travail a dû être effectué entièrement à la main, une nouvelle équation a été déduite à partir de zéro : cela a consommé la moitié du temps du projet qui s'est étalé sur deux ans. Le manque d'outils logiciels adaptés a aussi ralenti la progression.

        Dans ces conditions, l'ambition initiale de "clore le sujet" n'a pu être atteinte, au contraire : ce travail a levé une grande quantité de questions de toutes sortes. Quelques-une ont été résolues mais la plupart est hors du domaine d'étude initial, elles sont plus théoriques ou trop complexes.

        La limite du projet est surtout au niveau du temps d'étude car le sujet est inimaginablement vaste et peu de personnes travaillent actuellement dans ce domaine. L'autre limite, pratique, est très importante car le "budget" (matériel, financier) est quasiment inexistant (ce fut déjà une chance de disposer d'une plateforme biprocesseur). Enfin, les logiciels sont actuellement inadaptés à la tâche de programmation efficace en langage machine : beaucoup de temps a été investi dans le développement de techniques de programmation.

        La première chose à faire est d'améliorer l'environnement de développement. Le projet GNL est l'aboutissement d'années de réflexion sur la pratique du codage en assembleur sur machines superscalaires. Les premiers fichiers sont écrits (l'interface et l'importation de fichiers C sont commencés) pour l'environnement Unix (le logiciel GNL est complétement portable). Il faut concevoir des API afin que de nombreux modules d'assistance au codage puissent communiquer avec l'utilisateur. Des modules d'entrée et de sortie doivent aussi être con�us pour chaque langage et chaque processeur supporté.

        Une fois qu'une version suffisamment complète de GNL sera prête, le programme de gaz sur réseaux de ce mémoire sera recodé, analysé, amélioré, porté et pourra accueillir les suggestions du début de ce chapitre. GNL ne se limite pas à ce projet de LGA, le but avoué est de remplacer une grande partie du codage textuel classique. Si on peut recoder entièrement le programme avec GNL et l'améliorer facilement, ce sera une sorte de première mise à l'épreuve du concept et permettra de l'appliquer à d'autres domaines sensibles : programmation de kernels de systèmes d'exploitations, de routines d'interruptions, de kernels de calculs, de moteurs d'IA ou de 3D pour les jeux vidéos, et même pourquoi pas pour le prototypage rapide d'algorithmes.

        Enfin, à mesure que les processeurs superscalaires exécutent de plus en plus d'instructions en parallèle pour chaque cycle, il devient de plus en plus difficile de programmer des compilateurs pouvant extraire l'ILP d'un programme écrit en C. Les moyennes actuelles sont de l'ordre de trois instructions par cycle dans des conditions "idéales", ce qui est à la fois insuffisant pour augmenter substanciellement la performance des processeurs et trop compliqué pour les compilateurs. GNL expose le parallélisme d'un programme à de nombreux niveaux et permet donc de mieux profiter des processeurs futurs. GNL sera ainsi capable d'aider à l'adaptation du programme vers le Pentium II (les règles de pairage sont différentes).

        Pour ce qui est des applications pratiques, le modèle FHP n'est pas aussi efficace pour des simulations à grands nombres de Reynolds. Je vais donc essayer d'appliquer mes connaissances et mes techniques dans le domaine des gaz sur réseaux non booléens comme le modèle Lattice Bolzmann ou BGK. Une voie de recherche à explorer est l'utilisation du modèle ILG (Integer Lattice Gas, sorte de compromis entre Bolzmann et FHP) avec un système de numération logarithmique (nombres entiers différent de Base-2, plus simples que les nombres en virgule flottante). La parallélisation de ce travail permettrait de simuler de larges systèmes efficacement, c'est l'objet du projet beo-kragen (beowulf de Kragen Sitaker).

[1]	Bernard Ourghanlian "Les microprocesseurs Alpha" InterEditions, 1995, ISBN 2 7296 0565 7         Excellente introduction à ce processeur révolutionnaire, écrite en français par le directeur de développement de Digital. Ce n'est pas une documentation technique pure, car elle explique aussi le pourquoi du comment de chaque aspect et de chaque décision de la conception du modèle de programmation. Bon livre sur l'architecture RISC du futur.
[2]	Henri Lilen - René-Véran Honorat "Microprocesseurs PowerPC" ed. Dunod, 1995, ISBN 2 10 002464 7         Ecrit en francais, ce livre présente le modèle de programmation de la famille PowerPC et s'attarde sur les membres commercialisés avant la sortie du livre. Plusieurs techniques architecturales sont expliquées comme l'exécution dans le désordre, mais la suite est plut�t une traduction des documents anglais. On comprend donc le PPC sans vraiment en savoir plus...
[3]	David A Patterson - John L Hennessy "Computer Organisation & Design: the hardware/software interface", Morgan Kaufmann, 1994, ISBN 1 55860 282 8         C'est LE "Patterson & Hennessy", qui explique clairement les fondements de l'architecture des presque tous les ordinateurs. Il raconte leurs aléas et leurs avancées et permet de mieux comprendre l'importance de la relation entre le logiciel et le matériel qui le fait tourner, en fonction de leur rapport performance/prix. Il ne faut pas confondre ce livre avec leur précédent ouvrage "A quantitative approach" (le "QA") qui le précède.
[4]	Gilles Deghilage "Architectures et programmation parallèles", Addison-Wesley, 1992, ISBN 2087908 023 1         "Approche pratique en environnement scientifique sur multiprocesseurs Silicon Graphics", ou comment tirer le meilleur de architectures MIMD, des compilateurs, des algorithmes... Les aspects théoriques et pratiques sont abordés et les résultats sont comparés aux autres plateformes existantes. Il introduit dans leur contexte les techniques de parallélisation de code comme le strip mining.
[5]	Hans-Peter Messmer "Pentium et compagnie", Addison-Wesley, 1994, ISBN 2 87908 074 6         Bien qu'entaché de quelques erreurs et plein de détails inutiles, l'intérêt de ce gros livre est de présenter sous toutes ses coutures le monstre CISC de 3 millions de transistors, ainsi que ses relations avec ses voisins directs: contr�leur de cache, de mémoire, de bus PCI... Cette synthèse met le doigt sur énormément d'aspects, électroniques, architecturaux ou de programmation, qui sont nécessaires pour développer efficacement, mais il vaut mieux se référer aux constructeurs et aux sites x86.org et PCguide.com pour avoir une information plus fiable.
[6]	Alfred Aho - Ravi Sethi - Jeffrey Ullman et al. "Compilateurs : principes, techniques et outils", InterEditions, 1989, ISBN 2 7296 0295 X         Le "dragon book", c'est la raison pour laquelle je n'utilise pas de compilateur. C'est aussi la raison pour laquelle j'en utilise un. Seulement, je sais maintenant quand je peux m'y fier. Accessoirement, permettrait de faire un compilateur optimiseur pour le PII en MMX si cela en valait la peine...
[7]	Michael Abrash "Le Zen de l'optimisation du code", Sybex, ISBN 2-7361-2128-7         Devrait être lu par toute personne considérant coder "correctement", bien que l'aspect "assembleur" puisse repousser les habitués du "tout C". Quand on va à la chasse à la performance, "on ne fait pas d'omelette sans casser d'oeufs", et l'auteur nous apprend à faire enfin fonctionner le merveilleux compilateur qui est entre nos deux oreilles. En plus, ce livre est facile à lire.
[8]	Ronald J. Tocci "Circuits Numériques : théorie et applications", Dunod, 1992, ISBN 2-10-001576-1         Présente entre autres les techniques de base de simplification manuelle d'équations booléennes. Malheureusement insuffisant avec plus de 5 ou 6 variables d'entrée et plus d'une variable en sortie. De plus, l'analyse se réduit aux sommes de produits alors qu'on dispose souvent du XOR.
[9]	Daniel H. Rothman et Stéphane Zaleski "Lattice-gas cellular automata", Collection Aléa Saclay, Cambridge University Press, 1997, ISBN 0-521-55201-X         Sous-titré "Simple models of complex hydrodynamics", ce livre anglais de physique statistique est écrit par deux spécialistes et fournit des bases importantes sur la théorie des gaz sur réseaux. Après un exposé complet, les domaines de recherche des auteurs transparaissent dans l'étude des infiltrations des milieux poreux et de la transition de phase. Niveau doctorat.
[10]	Valérie Guimet "Analyse numérique et simulation de problèmes d'interaction fluide-structure en régime incompressible", thèse soutenue le 20 octobre 1998         Cette thèse pour le doctorat de Mathématiques Appliquées de l'université Paris VI explique des méthodes classiques de simulation de couplage entre des objets déformables (en 1D et 2D) et des fluides turbulents (en 2D et 3D). Bien que le domaine du calcul intensif soit abordé (ressources informatiques de l'ONERA), aucune explication pratique n'est donnée.
[11]	Jean-Pierre Rivet "Hydrodynamique par la méthode des gaz sur réseaux", thèse #88NICE4215 en Mécanique fondamentale et appliquée         Dirigée par Uriel Frisch, cette thèse communique une grande somme de connaissances théoriques (par l'explication des modèles) et pratiques (par l'explication des algorithmes) sur les gaz sur réseaux bi et tridimensionnels. Le RAP-1 est décrit mais surtout des simulations tridimensionnelles sur CRAY-2 sont effectuées et expliquées. Pas de code source pourtant mais un intérêt certain pour l'efficacité.
[12]	Pierre Audibert "Méthodes de grille pour le traitement des problèmes de mécanique des fluides", mémoire de maitrise de l'université Paris 8         Cet exposé présente clairement 3 méthodes différentes avec explications et programmes sources commentés. Bien que le modèle FHP soit utilisé en dehors de ses limites, le mémoire montre clairement les avantages et inconvénients pratiques et théoriques de chaque approche (les équations différentielles classiques et la méthode des tourbillons sont aussi présentées). L'efficacité est implicitement limitée par la plateforme (compatible PC sous MS-DOS, programmé en Turbo C).
[13]	Yue Hong Qian "Gaz sur réseaux et théorie cinétique sur réseaux appliquée à l'équation de Navier-Stockes", thèse #90PA06 658 de mécanique et physique statistique         Dirigée par D. d'Humières et P. Lallemand, cette thèse commence par l'étude spectrale d'un gaz sur réseau unidimensionnel pour extrapoler vers 2 dimensions (FHP, D2Q9 etc) puis 3 dimensions (FCHC et similaires). L'intérêt est principalement de consolider et d'explorer les travaux théoriques effectués 5 ans plus t�t.
[14]	Umberto d'Ortona "Hydrodynamique et Gaz sur réseau", thèse         Explication claire, malgré quelques détails ambigus, des gaz sur réseau classiques dans les premiers chapitres. Les néophytes aprécieront de pouvoir comprendre certains détails et les solutions apportées. La thèse étudie entre autres les phénomènes de capillarité (interfaces entre une goutte d'eau et un autre fluide).
[15]	Valérie Pot "Etude Microscopique du transport et du changement de phase en milieu poreux par la méthode des gaz sur réseaux", thèse #94PA06 233         Ces travaux font l'objet d'un chapitre entier dans le livre de Stéphane Zaleski. Beaucoup de formules et pas de code.
[16]	Stéphane Zaleski "Les transitions de phase calculées", Pour la science #183, janvier 1993         Article de présentation (2 pages) annonçant les résultats encourageants de recherches sur la transition de phase en 3D, améliorant les techniques présentées par Jean-Pierre Rivet et pour des règles d'interaction non locales.
[17]	Bernard Derrida "Dynamique d'un gaz sur réseau", Pour la science #184, février 1993         Cette "récréation informatique" présente un algorithme 2D pour simuler l'évaporation ou la condensation d'une goutte d'eau. Rien de commun avec FHP mais il simule une transition de phase (liquide/gaz) sans utiliser les intéractions non locales utilisées par Zaleski, ce qui fait de cet automate cellulaire une alternative intéressante.
[18]	Pierre Lallemand "Les gaz sur réseau : un nouveau médium pour le calcul des écoulements", Revue du Palais de la Découverte, avril 1987, vol. 15, numéro 147         C'est cet article qui m'a fait découvrir le domaine des LGA. La théorie y est exposée ludiquement et le modèle FHPII est succintement expliqué, suggestivement mais sans aider pour la programmation pratique. En particulier les règles de collisions sont expliquées mais la mise en pratique n'est pas évidente (création de la table des collision).
[19]	U. Frisch, B. Hasslacher, Y. Pomeau "Lattice gas automata for the Navier-Stockes equation", Physical Review Letters, 7 avril 1986, vol.56, pp 1505-1508.         C'est l'article de référence cité par tout papier ou document portant sur les gaz sur réseaux. Appelé dans cet article "HLG" pour "Hexagonal Lattice Gas", le modèle présenté retiendra les initiales des auteurs pour la postérité. Le passage du modèle HPP (carré) au voisinage hexagonal est expliqué, trois règles de collisions sont introduites (frontale, triangulaire et avec particule immobile) et il est prouvé que cela est suffisant pour correspondre aux équations de Navier-Stockes. Pas de code mais l'existence de calculateurs dédiés est évoquée (RAP-1 ?)
[20]	Gary Doolen (éditeur) "Lattice gas methods for partial differential equations", Santa Fe Institute, Studies in the sciences of complexity, Addison-Wesley 1990, ISBN 0-201-15679-2 ou 0-201-13232-X         Probablement l'un des ouvrages les plus intéressants qui existe. Tous les acteurs du domaine de la première décade sont présents dans cette collection de papiers, traitant de tous les domaines, pratiques et théoriques. Par exemple on y trouve le modèle FHP4, sur 8 bits, qui restaure l'invariance galliléenne (en ralentissant le fluide). La variété des points de vue et des sujets abordés rend sa consultation impérative, mais pourtant il n'y a pas de code.
[21]	J. A. Somers & P. C. Rem "Obtaining numerical results from the 3D FCHC lattice gas" dans Lecture Notes in Physics, T. M. M. Verheggen (ed.), "Numerical methods for the simulation of multi-phase and complex flows" (1990) Springer-Verlag, ISBN 0-387-55278-2         Les auteurs proposent une amélioration de la technique de J. P. Rivet pour réduire la taille de la table des collisions du modèle FCHC par une analyse poussée des 1152 isométries des collisions : la table contient seulement 106496 entrées mais le code Pascal est assez complexe. Performance de 2Mc/s sur un Transputer à 400 noeuds.
[22]	Jean-Christophe Culioli "Introduction � l'optimisation" (1994) Ed. Marketing/Ellipses, ISBN 2-7298-9428-4, ref 519.8 CUL         Sous-titre : "Analyse num�rique de syst�mes complexes". Pr�sente les techniques de Newton-Raphson, de dichotomie, de programmation lin�aire et dynamique.

[23]	Christopher Adler, Bruce M. Boghosian, Eirik G. Flekkoy, Norman Margolus, Daniel H. Rothman : "Simulating Three-Dimensional Hydrodynamics on a Cellular-Automata Machine"             à paraître dans Journal of Statistical Physics (1995)             r�f�rence preprint : chaos-dyn/9508001 ou comp-gas/9508001
[24]	Quian et.al. : Lattice BGK Models for Navier-Stokes Equation,             Europhys.Lett., 17 (6), pp. 479-484 (1992)
[25]	Norman Margolus : "Integer Lattice Gases",             Phys. Rev. E 55 (April, 1997) 4137-4147.
[26]	Ales Alajbegovic, Chris Teixeira, David Hill, Andrew Anagnost, Sudheer Nayani and Ram Iyer : "The study of benchmark laminar flows using DIGITAL PHYSICS"             1997 ASME Fluids Engineering Division Summer Meeting, FEDSM'97,             June 22-26, 1997, FEDSM97-3648
[27]	Jean-Pierre Boon, Uriel Frisch et Dominique d'Humi�res : "L'hydrodynamiqe mod�lis�e sur r�seau"             dans La Recherche n� 253, avril 1993, volume 54, pages 390-399.
[28]	Dominique d'Humi�res, Pierre Lallemand : "Numerical simulations of hydrodynamics with Lattice Gas Automata in two dimensions"             dans Complex Systems 1 (1987) 599-632
[29]	Dominique d'Humi�res, Pierre Lallemand, Geoffrey Searby : "Numerical experiments on Lattice Gases : mixtures and galilean invariance"             dans Complex Systems 1 (1987) 633-647
[30]	Norman Margolus (+?) : "An embedded DRAM architecture for large-scale spatial-lattice computations" (1999)             � para�tre (?)
[31]	Kristian Lindgren, Cristopher Moore et Mats G. Nordahl "Predicting Lattice Gases is P-complete"             Santa Fe Institute Working Paper 97-04-034.
[32]	D. d' Humières and P. Lallemand : "Numerical simulations of hydrodynamics with lattice gas automata in two dimensions."             dans Complex Systems, 1:599, 1987.
[33]	"GA-586ATE users manual, 1st edition"             manuel de carte m�re pour PC, 1995.
[34]	Intel : "Pentium II Processor Developper's Manual"             r�f�rence 243502-001, octobre 1997
[35]	Bruce M. Boghosian, Jeffrey Yepez, Francis J. Alexander, Norman H. Margolus : "Integer Lattice Gases"             comp-gas 9602001 (15 fevrier 1996)             mise � jour le 22 novembre 1998
[36]	"Benchmark for the simulation of the incompressible flow around a cylinder", 1996             http://gaia.iwr.uni-heidelberg.de/~ture/papers.html             ou in Flow Simulation with High-Performance Computers II,             Notes on Numerical Fluid Mechanics, Vol 52, Vieweg 1996
[37]	J. Hardy, Y. Pomeau & O. de Pazzis : "Time evolution of two-dimensional model system. I. Invariant states and time correlation functions", 1973             J. Math. Phys. 14, pp. 1746-1759.
[38]	Fung Funh Lee : "A Scalable Computer Architecture for Lattice Gas Simulations", 1993             Technical Reports : CSL-TR-93-576             d�crit "ALGE", un "superordinateur SIMD"             M�moire de Doctorat de Philosphie, Universit� de Stanford, d�partements EE/CS
[39]	Intel : "Intel Architecture Software Developer�s Manual Volume 2: Instruction Set Reference", 1999             24319102.pdf
[40]	Kohring, G.A. "Parallelization of short- and long-range cellular automata on scalar, vector, SIMD and MIMD machines" 1991             Inter. Jour. Modern Phys. 2, 755-772.
[41]	Brosa, U. and Stauffer, D. "Vectorized multisite coding for hydrodynamic cellular automata" 1989             Jour. Stat. Phys. 57, 399-403.

                Pour des raisons évidentes, les URL collectionnées ici ne sont ni exhaustives ni garanties pour leur fraîcheur, puisque le travail présenté dans le mémoire a commencé vers 1995, au début de l'émergence d'Internet comme média pour le grand public. Internet était cependant déjà utilisé et très utile pour la communauté des chercheurs sur les LGA et on peut remarquer que les sites comme les pratiques ont peu changé. Les URL ont été vérifiées en mai 2000 et on peut penser que la plus grande partie sera valable encore pendant plusieurs années.

Usenet :
news:comp.lang.asm.x86 (ou "clax")
C'est là que les "hackeurs", les "demo makers" et autres "code gurus" apprennent à écrire leurs premières lignes de code en assembleur, configurer les registres de la carte vidéo, lire et contr�ler les périphériques. C'est ce newsgroup qui a engendré le travail coopératif sur NASM. Le "rapport signal/bruit" et le niveau peu avancé fait qu'on se tourne ensuite rapidement vers d'autres ressources plus pointues comme news:comp.arch.

[email protected] :
la "pmode-l", mailing list dédiée à la programmation avancée en mode protégé. On y apprend à passer en mode protégé "à la main" ou à profiter des trous de sécurité d'un certain système propriétaire grand public. Plus sérieusement, c'est un point de rencontre de tous les développeurs de kernel.

http://www.cryogen.com/Nasm/
Le site officiel de NASM, l'assembleur fait par et pour les gens qui programment les x86 en assembleur.

Contributeurs de NASM et autres gurus de la programmation en assembleur (ne sont ici que les meilleurs) :

http://www.x86.org/,
http://www.sandpile.org/ :
ressources indépendantes sur les architectures x86 (à consulter absolument avant de programmer !)

http://www.ddj.com/ftp/ :
Dr. Dobb's Source Code archive : plein de morceaux croustillants !

http://www.simtel.net/simtel.net/msdos/index-msdos.html et http://oak.oakland.edu/simtel.net/ :
Simtel.net FTP shareware archive (certains avec source commentés)

http://www.agner.org/assem :
Les précieux conseils d'Agner Fog

http://www.cs.virginia.edu/stream/ :
John McCalpin ([email protected]) développe depuis 1991 le benchmark STREAM, permettant de comparer des architectures de manière totalement objective : en mesurant la bande passante de la machine, disponible à partir de sources en C non optimisés. On dispose ainsi non pas de la puissance de "crête" mais de la puissance "brute" du système, celle qu'un utilisateur moyen accède dans la pratique.

http://www.senet.com.au/~cpeacock :
page de Craig Peacock sur la programmation des interfaces et périphériques (PDF et HTML pour contr�ler le port parallèle, le clavier, la souris, le port USB, les interruptions...)

http://www.cs.wisc.edu/~glew/ :
Andy Glew ([email protected]), impliqué dans l'acceptation du MMX par Intel et dans le coeur du P6, est une des figures de comp.arch. Ses opinions sur la programmation et les architectures des ordinateurs le placent en opposition de la culture "Patterson & Hennessy". Oui, il y a quelqu'un de compétent travaillant pour Intel... Si au moins il était écouté !

http://research.microsoft.com/~gbell/Computer_Structures__Readings_and_Examples/contents.html :
Gordon Bell (qui a créé le prix du même nom) a mis en ligne un de ses livres, introuvable actuellement, qui est une sorte de pierre de Rosette pour la paléoinformatique comparée. L'histoire des ordinateurs et l'étude de nombreux cas devrait être un module obligatoire dans une formation d'informaticien !

http://www.cs.cmu.edu/afs/cs/user/ralf/pub/WWW/files.html :
Ralph Brown a créé et entretenu LA référence sur l'architecture logicielle des PC, en réunissant et en documentant des milliers d'appels systèmes de MSDOS et d'autres logiciels similaires (Windows, TSR, Virus, utilitaires, DOS-extenders...). Indispensable.

Intel:
http://developer.intel.com/drg/pentiumii/appnotes/index.htm : Notes d'applications pour le PentiumII.
http://developer.intel.com/design/pentiumII/manuals/243502.htm : Manuel du développeur pour le PentiumII
http://developer.intel.com/design/mmx/manuals/ : Manuels pour développer avec les instructions MMX

VESA :
http://www.monstersoft.com/tutorial1/VESA_info.html
http://mirriwinni.cse.rmit.edu.au/~steve/vbe/vbe20.htm :
"VESA BIOS EXTENSION(VBE) Core Functions with DJGPP Code" ou comment accéder aux modes vidéo haute résolution en assembleur...

http://www.gdsoft.com/swag/downloads.html :
"SWAG" pour "SourceWare Archive Group", plein d'algorithmes et d'utilitaires pour Turbo Pascal/PC. Un des derniers vestiges de la culture "shareware" et "BBS". Au lieu de réinventer la roue, il suffit de la copier :-)

http://www.phoenix.gb.net/x86/ :
Page de Win32NASM, pour programmer en assembleur sous Windows, mais en mieux !

Usenet :
news:comp.theory.cell-automata et news:comp.theory.dynamic-sys

http://www.cs.runet.edu/~dana/ca/cellular.html :
Page dédiée aux automates cellulaires et à Cellang

Page de Bruce Boghossian, travaillant pour Thinking Machines Corp (avant sa fermeture) et responsable de la CA-list
http://physics.bu.edu/~bruceb/
http://physics.bu.edu/~bruceb/MolSim/ : "Mesoscale Modeling of Amphiphilic Fluid Dynamics" (transition de phase en 3D)

Homepage d'Oleh Baran, travaillant sur les LGA :
http://www.physics.mcgill.ca/WWW/oleh/Welcome.html

pre-print archive :
Avant d'être envoyés à des publications périodiques sur papier, les chercheurs soumettent leurs papiers à ces sites permettant d'effectuer des recherches documentaires sur des sujets pointus. Les serveurs sont souvent chargés car les documents comme les pages sont générés à la volée à partir d'une base de donnée et transformés en une grande variété de formats (PS, PDF, DOC etc) avant d'être transférés.
http://www.arc.umn.edu/publications/preprints/, http://xyz.lanl.gov, http://xxx.lanl.gov/archive/cond-mat
miroir à Jussieu :
http://xxx.lpthe.jussieu.fr, http://fr.arXiv.org/, http://fr.arXiv.org/find/cond-mat/

http://www.ph.ed.ac.uk/~jmb/thesis/tot.html :
thèse de James M. Buick, calculant FHP-3 sur une CM200 à Edinburgh.

http://poseidon.ulb.ac.be/lga.html :
"CNLPCS: Unité Automates de Gaz sur Réseau". L'Université Libre de Bruxelles a un département actif dans le domaine des gaz sur réseaux, les systèmes dynamiques, non linéaires et adaptatifs.

Site web de la CAM8 :
http://www-im.lcs.mit.edu/
http://www.im.lcs.mit.edu/broch/ (quelques exemples hauts en couleurs)

http://zanzibar.mit.edu/ (serveur sporadiquement éteint)
dédié aux applications géophysiques des LGA (infiltration dans les roches etc)

http://www.wizard.com/~hwstock/saltfing.htm
Exposé des travaux de Haarlan Stockman. C'est probablement le seul scientifique qui soit aussi impliqué dans la vitesse de son code, car c'est aussi un codeur système engagé (il a écrit pour la première fois un code de calcul de l'ensemble de Mandelbrot en virgule fixe pour 386 en 1986, paru au Doctor Dobb's Journal)

http://www.fuji-ric.co.jp/complex/complex/LGA/result/flatplate.html : (en japonais)
"Behind Flat Plate Flow", travaux des laboratoires Fuji au Japon : même les Japonais utilisent les LGA !

http://www.tele.unit.no/akustikk/person/kristiansen/sudosparrow.html
"The Sudo/Sparrow lattice gas model" : application des LGA à l'acoustique.

http://www.obs-nice.fr/cassini/HTML_FR/rivet.html
"NON-LINEAR DYNAMICS AND TURBULENCE APPLIED TO FLUIDS IN ASTRO- AND GEO-PHYSICS"
L'observatoire de Nice a été parmi les premiers laboratoires, avec l'Ecole de Physique-Chimie de Paris, à travailler sur les modèles "FHP" et "FCHC" vers 1980. Ils ont activement participé à l'élaboration et à la validation de la théorie FHP (voir thèse 1982 à Jussieu).

http://www.tu-bs.de/institute/WiR/weimar/ZAscript/ :
"Simulation with Cellular Automata" par J�rg Weimar, une étude générale des possibilités et des recherches sur les automates cellulaires.

http://www.exa.com :
Spinoff du MIT, exploite une version contestée de Lattice Bolzman en 3D pour des applications industrielles. En pratique les résultats sont à la hauteurs des milliers d'heures-CPU qui sont souvent nécessaires à une simulation à grand nombre de Reynolds...

http://amber.aae.uiuc.edu/~m-selig/ads/coord_database.html
"UIUC Airfoil Coordinates Database" : base de donnée de profils vectoriels d'ailes d'avions. J'y ai contribué avec un repackaging des données avec visualisation sous X et MSDOS.

Divers :

• http://www.santafe.edu/~moore/ :
  page de l'auteur de : ``Predicting Lattice Gases is P-complete'' (avec Mats Nordahl) Santa Fe Institute Working Paper 97-04-034.
• http://umunhum.stanford.edu/~morf/lattice.gas/lattice.gas.html
  Lien vers la thèse de Fung Fung Lee ([email protected]), "A Scalable Computer Architecture for Lattice Gas Simulations", Ph.D. dissertation, Dept. of Electrical Engineering, Stanford University, Mai 1993.
• http://www.cfd-online.com/Resources/docs.html
  Nombreux liens vers des ressources en ligne classiques et professionnelles.
• http://www.aoe.vt.edu/aoe/courses/webteach.html
  Tout pour la mécanique classique appliquée à l'aéronautique !
• http://www.aoe.vt.edu/aoe3114/calc.html
  Calculatrice en Java pour étudier les phénomènes de mécanique des fluides "classique".
• http://www5.informatik.tu-muenchen.de/forschung/visualisierung/praktikum.html (allemand)
  L'université de M�nich a étudié les LGA au département de calcul intensif et de visualisation scientifique, cette page propose des films MPEG des simulations qu'ils ont effectué. Merveilleux.
• http://raphael.mit.edu:80/Java/ : "the Java Virtual Wind Tunnel" ou la première expérience de mécanique des fluides intéractive sur Internet. Idéal aussi comme jeu ou pour tester la stabilité de Netscape : les codes (Euler, classiques) "explosent" parfois.