This is a 2D gravitational N-body simulation using the Particle-Mesh (PM) method, implemented entirely in WebGPU compute shaders. The PM method approximates gravitational interactions by depositing particle masses onto a regular grid, solving for the gravitational potential on that grid, then interpolating the resulting forces back to the particles.
The gravitational potential is related to the mass density through the Poisson equation,
which in Fourier space becomes a simple algebraic relation, This allows the potential to be computed efficiently via FFT. The gravitational acceleration is then , which can also be computed in Fourier space by multiplying by .
The PM method is considered somewhat obsolete because it cannot accurately model close interactions between particles. The grid resolution sets a minimum length scale below which gravitational forces are artificially softened. Modern cosmological simulations typically use the Particle-Particle Particle-Mesh (PM) method, which augments the PM calculation with direct particle-particle summation for nearby pairs, or tree-based methods that adaptively refine force calculations.
Nonetheless, PM lends itself well to GPU architecture and remains a useful pedagogical tool and performs well for problems where small-scale structure is unimportant or where the softening is acceptable.
Initial conditions
The simulation begins with particles distributed uniformly at random, each with a small random velocity drawn from a Gaussian distribution. Two physical scales govern how the system evolves.
The first is the Jeans length , which sets the scale above which gravitational collapse overcomes thermal pressure (here represented by the initial velocity dispersion). Density perturbations larger than will grow and collapse, while smaller perturbations oscillate as acoustic waves. With and mean density , the Jeans length is approximately .
The second consideration is shot noise from the finite number of particles. The initial density field is not perfectly uniform but has Poisson fluctuations from the random particle positions. With particles distributed over grid cells, each cell contains on average particles, and the relative density fluctuation is . These fluctuations seed the initial perturbations that gravity amplifies.
Implementation Details
The simulation begins each timestep by depositing particle masses onto the grid. Each particle contributes to the four nearest grid cells using bilinear interpolation weights, which ensures smooth density fields and conserves total mass. This “cloud-in-cell” scheme is standard in PM codes.
With the density field assembled, we transform it to Fourier space via a 2D FFT. The Poisson equation becomes , which we solve by simple division (setting the mode to zero to remove the arbitrary constant in the potential). The gradient is computed simultaneously in Fourier space by multiplying by , yielding both components of the force field.
After inverse FFTs return the force field to real space, particles sample their accelerations using the same bilinear interpolation used for mass deposition. This symmetry is important for momentum conservation. The particles are then advanced using a second-order midpoint (leapfrog) integrator.
The FFT naturally imposes periodic boundary conditions, so particles that exit one edge of the domain reappear on the opposite side. The overall complexity is for the grid operations (where is the grid resolution) plus for the particle operations, making the method efficient even for millions of particles.