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

GPU General Purpose Computing with Unity Compute Shaders

Compute Shaders in Unity3D are a way to perform general purpose computing on a GPU card.  This is known as General Purpose computing on Graphics Processing Units or GPGPU.  This means you can use the GPU to perform operations that may not necessarily have anything to do with graphics.  e.g. protein folding simulations (ya I don’t know what that is either).  Not everything is suitable, only those operations that can run concurrently and make use of the ridiculous amount of parallelism on a GPU.

Unity’s Compute Shaders are currently limited to supporting DirectX 11 only but I suspect this will change in the future when you consider that Unity’s aim is to be cross-platform.  DX11 compute shaders is also known as DirectCompute.

Compute shaders syntax is just as per HLSL so you should feel right at home.

Compute shaders are great.  Normally from what I understand, if you had a ton of stuff to manipulate well you had to do it CPU side then upload it to the GPU – a very expensive operation not to be done frequently.

With compute shaders you merely allocate the data CPU side and send it once to the GPU.  From then on a little program (or kernel) updates the data for you and the output is passed to a shader for rendering.  All without being touched by the CPU.

Here’s my first attempt at making a GPGPU galaxy modelled with 1,000,000 stars where the stars are moved on the GPU side of things.  It’s looks more like a pizza at the moment, oh well.  I can manage 60 FPS on 1920×1080.

 

 

I’ll follow up with a non-pizza looking galaxy soon and some source code as documentation on anything GPGPU is quite rare.

Maybe if we get enough of these GPUs we can find out what “42” means?