Une galaxie dans un onglet, calculée sur ton GPU avec WebGPU

Je voulais savoir jusqu’où une simulation gravitationnelle N-corps naïve pouvait aller dans un onglet de navigateur, sans serveur et sans images précalculées. La réponse : 65 536 corps à 60 fps, soit quelques milliards de calculs de force par seconde, le tout sur ton GPU via les compute shaders WebGPU.

Démo live · Code source sur GitHub (TypeScript, zéro dépendance runtime)

Voici comment ça marche, et les quelques décisions qui l’ont rendu assez rapide pour être amusant.

La forme du problème

Chaque corps subit la gravité de tous les autres. Avec N corps, c’est N² interactions, recalculées de zéro à chaque image. À 16 384 corps, ça fait environ 268 millions de paires de forces par pas ; à 65 536, c’est 4,3 milliards. Pas d’arbre, pas d’approximation, aucun cache entre les images. C’est la somme en force brute, et la seule raison pour laquelle ça tourne en temps réel, c’est qu’un GPU adore faire le même petit calcul quelques milliards de fois en parallèle.

Ce cadrage compte : c’est d’abord un problème de calcul, et seulement ensuite un problème de rendu. WebGL peut simuler du calcul avec des astuces de rendu vers texture, mais WebGPU offre de vrais compute shaders et de la mémoire partagée on-chip, ce qui est exactement ce dont cette charge a besoin. Donc ce sera WebGPU, et la page affiche un message de repli sur les navigateurs qui ne le supportent pas encore.

Tiler la somme des forces en mémoire partagée

Le kernel évident fait boucler chaque thread sur les N corps, en lisant chaque position directement depuis la mémoire globale. Ça marche, mais c’est limité par la bande passante, parce que chaque thread martèle la VRAM pour les mêmes données dont ses voisins ont aussi besoin.

La solution est le layout N-corps classique de GPU Gems. Chaque workgroup charge un bloc de corps (une tuile) depuis la mémoire globale vers la mémoire partagée une seule fois, les threads se synchronisent, puis chacun fait sa boucle interne sur cette copie on-chip. On parcourt la liste des corps tuile par tuile, en payant la lecture en mémoire globale une fois par tuile au lieu d’une fois par paire. La mémoire partagée est environ un ordre de grandeur plus rapide d’accès que la VRAM, et comme chaque thread du workgroup réutilise la même tuile, c’est de là que vient l’essentiel du débit.

Garder le disque en disque : un intégrateur stable

La première version utilisait un simple Euler explicite, et la galaxie se défaisait lentement en bouillie. Euler perd un peu d’énergie à chaque pas, et sur des millions de pas cette erreur s’accumule jusqu’à faire disparaître la structure.

Passer à un pas de leapfrog a réglé ça. Leapfrog est l’intégrateur standard pour ce genre de simulation parce qu’il reste stable en énergie sur de longues durées : l’énergie totale oscille un peu mais ne dérive pas, donc le disque garde sa forme quasi indéfiniment. Ça coûte un petit surcroît de comptabilité, en gardant position et vitesse décalées d’un demi-pas.

C’est le seul changement qui a transformé une jolie explosion en quelque chose qui ressemble à une galaxie.

Le softening, pour que les rencontres proches n’explosent pas

La gravité newtonienne varie en 1/r², ce qui veut dire que deux corps qui se frôlent subissent une force énorme et partent à l’infini. Les vrais codes N-corps gèrent ça avec du softening, et celui-ci aussi : un petit ε² est ajouté à chaque distance au carré avant la division. La force reste finie quelle que soit la proximité des deux corps.

Le softening a un second bénéfice, gratuit. Un corps qui calcule sa propre attraction sur lui-même a r = 0, et avec le terme ε² cette auto-interaction vaut zéro au lieu d’une division par zéro. La boucle interne n’a donc pas besoin de cas particulier pour ignorer le terme propre, ce qui garde le kernel sans branchement.

Buffers ping-pong et un détail d’ordonnancement WebGPU

Les positions vivent dans deux buffers. Chaque pas lit l’un et écrit l’autre, puis on les échange. Ça évite le hazard lecture-écriture qui surviendrait en mettant à jour les positions sur place pendant que d’autres threads les lisent encore.

Le détail à connaître : à l’intérieur d’une même passe de calcul WebGPU, les dispatches sont ordonnés. Donc quand je lance plusieurs sous-pas à la suite dans une passe, le read-after-write entre eux est déjà sûr et je n’ai pas besoin d’insérer de barrières manuelles entre les dispatches. Les vitesses, elles, vivent dans un seul buffer, parce que chaque thread ne touche jamais qu’à sa propre vitesse : aucun hazard.

Dessiner 65k points qui ne font pas 1 pixel

La primitive point de WebGPU ne dessine jamais qu’un point d’1 pixel, ce qui ressemble à de la neige. Chaque corps est donc rendu comme un quad instancié, deux triangles agrandis à une taille fixe en pixels, avec une atténuation radiale douce dans le fragment shader et un blending additif pour que les corps qui se chevauchent rayonnent. La couleur est indexée sur la vitesse, donc le disque interne rapide chauffe et le bord lent reste sombre. Cette seule correspondance fait l’essentiel du travail visuel pour que ça se lise comme une galaxie et non comme un nuage de particules.

Partir de quelque chose qui a déjà l’air juste

L’état initial est un disque en rotation autour d’une masse centrale lourde. La vitesse orbitale de chaque corps est fixée à partir de la masse contenue à l’intérieur de son rayon, donc le disque démarre en équilibre rotationnel approximatif au lieu de s’effondrer vers le centre. Les rayons sont biaisés vers le centre pour que le cœur soit plus dense que le bord, une petite épaisseur verticale vient d’une normale de Box-Muller, et une légère dispersion des vitesses l’empêche de ressembler à un anneau d’horlogerie. Le PRNG est initialisé avec une graine, donc la même galaxie revient à chaque rechargement.

Où ça atterrit

Sur un GPU dédié récent, ça tient 60 fps à 65k corps, ce que le HUD rapporte comme quelques milliards d’interactions de paires par seconde. L’ensemble tient en une poignée de fichiers TypeScript et deux shaders WGSL, construit avec Vite, sans dépendance runtime, et se déploie sur GitHub Pages.

L’étape suivante évidente est de casser le plafond en N² avec Barnes-Hut ou un FMM pour viser les millions de corps, puis de semer deux disques et de les faire entrer en collision. Mais il y a quelque chose de satisfaisant à voir jusqu’où va la bête version O(N²) quand on pointe le matériel droit dessus.

Si tu l’essaies et que quelque chose semble physiquement faux, ou si tu as un tiling plus rapide, ça m’intéresse.