n-Body Galaxy Simulation using Compute Shaders on GPGPU via Unity 3D

Following on from my prior article about GPGPU, I thought I try my hand at n-Body simulation.

Here’s my first attempt at n-Body simulation of a galaxy where n equals 10,000 stars with the algorithm running as a compute shader on the GPU.  Compute shaders are small HLSL-syntax-like kernels but can run in massive parallel compared to a regular CPU equivalent.

It’s important to remember that conventional graphics including those in games, do not use GPGPU for star fields.  Generally this is accomplished by either manipulating the star coordinates CPU-side and uploading to the GPU every frame or so; or utilise a small shader that renders a pre-uploaded buffer of stars that may or may not move over time.  In any event, it would be unlikely to perform n-Body simulation using conventional non-GPGPU means due to the slow performance.

My simulation entails n² gravity calculations per frame or 10,000² = 100,000,000 calculations per frame at 60 frames per second – something that is quite impossible if I attempted to do it from the CPU-side even if the GPU was still rendering the stars.

Here’s the compute kernel (note I’m hosting the Gist on GitHub as a “c” file in order to display syntax hi-lighting.  Readers should be aware the correct file extension is .compute):

	// Copyright 2014 Michael A. R. Duncan
	// You are free to do whatever you want with this source
	//
	// File: Galaxy1Compute.compute
	// Each #kernel tells which function to compile; you can have many kernels
	#pragma kernel UpdateStars
	#include "Galaxy.cginc"
	// blackmagic
	#define BLOCKSIZE 128
	RWStructuredBuffer<Star> stars;
	Texture2D HueTexture;
	// refer to http://forum.unity3d.com/threads/163591-Compute-Shader-SamplerState-confusion
	SamplerState samplerHueTexture;
	// time ellapsed since last frame
	float deltaTime;
	const float Softening=3e4f;
	#define Softening2 Softening * Softening
	static float G = 6.67300e-11f;
	static float DefaultMass = 1000000.0f;
	// Do a pre-calculation assuming all the stars have the same mass
	static float GMM = G*DefaultMass*DefaultMass;
	[numthreads(BLOCKSIZE,1,1)]
	void UpdateStars (uint3 id : SV_DispatchThreadID)
	{
	uint i = id.x;
	uint numStars, stride;
	stars.GetDimensions(numStars, stride);
	float3 position = stars[i].position;
	float3 velocity = stars[i].velocity;
	float3 A=float3(0,0,0);
	[loop]
	for (uint j = 0; j < numStars; j++)
	{
	if (i != j)
	{
	float3 D = stars[j].position – stars[i].position;
	float r = length(D);
	float f = GMM / (r * r + Softening2);
	A += f * normalize(D);
	}
	}
	velocity += A * deltaTime;
	position += velocity * deltaTime;
	if (i < numStars)
	{
	stars[i].velocity = velocity;
	stars[i].position = position;
	stars[i].accelMagnitude = length(A);
	}
	}

view raw

GalaxyCompute.c

hosted with ❤ by GitHub

Here’s my Unity MonoBehaviour controller that initialises and manages the simulation:

	#region Copyright
	// Copyright 2014 Michael A. R. Duncan
	// You are free to do whatever you want with this source
	// File: Galaxy1Controller.cs
	#endregion
	#region
	using Assets.MickyD.Scripts;
	using UnityEngine;
	#endregion
	public class Galaxy1Controller : MonoBehaviour
	{
	private const int GroupSize = 128;
	private const int QuadStride = 12;
	#region Fields
	/// <summary>
	/// The galaxy radius
	/// </summary>
	/// <remarks>This will appear as a Property Drawer in Unity 4</remarks>
	[Range(10, 1000)]
	public float GalaxyRadius = 200;
	public Texture2D HueTexture;
	public int NumStars = 10000;
	public ComputeShader StarCompute;
	public Material StarMaterial;
	private GameManager _manager;
	private ComputeBuffer _quadPoints;
	private Star[] _stars;
	private ComputeBuffer _starsBuffer;
	private int _updateParticlesKernel;
	#endregion
	// Use this for initialization
	#region Properties
	private Vector3 StartPointA
	{
	get { return new Vector3(GalaxyRadius, 0, 0); }
	}
	private Vector3 StartPointB
	{
	get { return new Vector3(-GalaxyRadius, 0, 0); }
	}
	#endregion
	#region Methods
	private void CreateStars(int offset, int count, Vector3 T, Vector3 V)
	{
	for (var i = offset; i < offset + count; i++)
	{
	var star = _stars[i];
	star.color = Vector3.one; // white
	star.position = Random.insideUnitSphere*GalaxyRadius + T;
	star.velocity = V;
	_stars[i] = star;
	}
	}
	private void OnDestroy()
	{
	// must deallocate here
	_starsBuffer.Release();
	_quadPoints.Release();
	}
	private void OnDrawGizmos()
	{
	Gizmos.color = Color.yellow;
	Gizmos.DrawWireSphere(transform.position, GalaxyRadius);
	Gizmos.DrawWireSphere(transform.position + StartPointA, GalaxyRadius);
	Gizmos.DrawWireSphere(transform.position + StartPointB, GalaxyRadius);
	}
	private void OnRenderObject()
	{
	if (!SystemInfo.supportsComputeShaders)
	{
	return;
	}
	// bind resources to material
	StarMaterial.SetBuffer("stars", _starsBuffer);
	StarMaterial.SetBuffer("quadPoints", _quadPoints);
	// set the pass
	StarMaterial.SetPass(0);
	// draw
	Graphics.DrawProcedural(MeshTopology.Triangles, 6, NumStars);
	}
	private void Start()
	{
	_updateParticlesKernel = StarCompute.FindKernel("UpdateStars");
	if (_updateParticlesKernel == -1)
	{
	Debug.LogError("Failed to find UpdateStars kernel");
	Application.Quit();
	}
	_starsBuffer = new ComputeBuffer(NumStars, Constants.StarsStride);
	_stars = new Star[NumStars];
	var n = NumStars/2;
	var offset = 0;
	CreateStars(offset, n, StartPointA, new Vector3(-10, 5, 0));
	offset += n;
	CreateStars(offset, n, StartPointB, new Vector3(10, -5, 0));
	_starsBuffer.SetData(_stars);
	_quadPoints = new ComputeBuffer(6, QuadStride);
	_quadPoints.SetData(new[]
	{
	new Vector3(-0.5f, 0.5f),
	new Vector3(0.5f, 0.5f),
	new Vector3(0.5f, -0.5f),
	new Vector3(0.5f, -0.5f),
	new Vector3(-0.5f, -0.5f),
	new Vector3(-0.5f, 0.5f),
	});
	_manager = FindObjectOfType<GameManager>();
	}
	private void Update()
	{
	// bind resources to compute shader
	StarCompute.SetBuffer(_updateParticlesKernel, "stars", _starsBuffer);
	StarCompute.SetFloat("deltaTime", Time.deltaTime*_manager.MasterSpeed);
	StarCompute.SetTexture(_updateParticlesKernel, "hueTexture", HueTexture);
	// dispatch, launch threads on GPU
	var numberOfGroups = Mathf.CeilToInt((float) NumStars/GroupSize);
	StarCompute.Dispatch(_updateParticlesKernel, numberOfGroups, 1, 1);
	}
	#endregion
	}

view raw

Galaxy1Controller.cs

hosted with ❤ by GitHub