Accelerated N-body sim with OpenGL graphics & automatic CUDA->SYCL conversion using dpct.
Forked from https://github.com/salel/nbody
This nbody simulation can be run with any of:
- CUDA
- DPC++ CUDA backend
- DPC++ OpenCL CPU backend
Source code for the CUDA version is in ./src/ while ./src_sycl/ contains the semi-automatically converted SYCL code.
The rendering components of this code are independent of the CUDA/SYCL backend, and depend on:
- GLM
- GLFW
- GLEW
These can be installed with apt:
sudo apt update
sudo apt install libglew-dev libglfw3-dev libglm-dev libxxf86vm-dev libxcursor-dev libxinerama-dev libxi-dev
The implementation relies on OpenGL 4.5.
The CUDA version of this code requires the CUDA runtime to be installed on your machine.
The DPC++ CUDA backend version also requires the CUDA runtime.
The DPC++ OpenCL backend requires an OpenCL runtime. To run specifically on the CPU, you must install the OpenCL runtime for your CPU.
Both DPC++ backends require the DPC++ compiler to compile the SYCL code.
This project uses CMake for build configuration. Build scripts for CUDA and DPC++ are located in ./scripts/. Note that these scripts include some hardcoded paths from our dev machine, and so will not work out-the-box.
The CMake option -DBACKEND allows to select which backend ("CUDA" or "DPCPP") to build. CUDA is built by default. The name of the built binary is suffixed with the backend (nbody_cuda or nbody_dpcpp).
The DPC++ backend, in turn, supports both an OpenCL & CUDA backend, both of which are built by default. If you are building on a machine without CUDA support, you can switch off the DPC++ CUDA backend with the flag -DDPCPP_CUDA_SUPPORT=off.
The script ./scripts/run_dpct.sh calls a containerized version of the Intel® DPC++ Compatibility Tool to automatically convert the CUDA components of this project into SYCL. A docker container was used because the dev machine has an incompatible version of the CUDA driver. This should be adapted based on your environment.
The Intel® DPC++ compatibility tool offers options for intercepting complex builds, but current dev environment restrictions require me to run the tool inside a docker container. This complicates things, so for now I'm just doing single source conversion on the simulator.cu file.
The script ./scripts/run_nbody.sh will run the nbody simulation, selecting a different binary based on the -b flag, where -b can be cuda or dpcpp. Subsequent positional arguments are passed on to the nbody binary. These positions args are described in the Simulation section. For example, to run on the DPC++ OpenCL host backend with 25600 (100 * 256) particles, executing 10 timesteps per rendered frame:
./scripts/run_nbody.sh -b dpcpp 100 10
Note that this script runs nbody with the default X window, as opposed to using xvfb. This makes it unsuitable for running on a remote machine.
run_nbody.sh is a simple wrapper around the nbody_* binaries with some environment variables set; the sections below describe how to launch the binaries directly.
The sycl-ls tool allows you to check for available backends on the system. For example, on a system with Intel OpenCL CPU runtime & CUDA runtime, the output is:
> sycl-ls
[opencl:cpu:0] Intel(R) OpenCL, Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz 3.0 [2021.13.11.0.23_160000]
[opencl:cpu:1] Intel(R) OpenCL, Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz 3.0 [2021.13.11.0.23_160000]
[cuda:gpu:0] NVIDIA CUDA BACKEND, NVIDIA GeForce RTX 3060 0.0 [CUDA 11.6]
[host:host:0] SYCL host platform, SYCL host device 1.2 [1.2]
By specifying the environment variable SYCL_DEVICE_FILTER, it's possible to switch between running with the CUDA backend and the OpenCL host backend. For example:
SYCL_DEVICE_FILTER=cuda ./nbody_dpcpp
will run on the CUDA backend, whereas:
SYCL_DEVICE_FILTER=opencl:cpu ./nbody_dpcpp
will run on a CPU through the OpenCL backend. Note the correspondence between options for SYCL_DEVICE_FILTER and the output of sycl-ls.
Note: Selection between DPC++ backends at runtime is possible because CMakeLists.txt specifies building the SYCL code for both CUDA (nvptx64-nvidia-cuda) & OpenCL (spir64) targets:
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fsycl -fsycl-targets=spir64,nvptx64-nvidia-cuda -fsycl-unnamed-lambda")
No changes to the code were required, but there were a couple of bugs which are worked around.
Firstly, when building for multiple targets (-fsycl-targets), there is a recent bug which causes failure to link to static libraries. The workaround for this is to switch from building imgui as a static to a shared library.
Secondly, I encountered the common CL header bug (see here and here). This turned out to be triggered for the spir64 backend because the CUDA headers were included only via -I and not via -internal-isystem. This caused them to take precedence over SYCL CL headers. The solution was to not include CUDA headers in src_sycl/CMakeLists.txt, which turned out to be unnecessary anyway.
OpenGL & CUDA are capable of interoperating to share device memory, but this will not play well with the Intel® DPC++ Compatibility Tool. Instead, computed particle positions are migrated back to the host by CUDA/SYCL, then sent back to OpenGL via mapping.
The DiskGalaxySimulator class handles the physics of n-body interaction. The computation of interparticle forces, velocity & updated particle positions are handled by the CUDA kernel particle_interaction.
The equation solved by this code is equivalent to Eq. 1 here, with the simplifying assumption that all particles have unit mass and there is no external/background force. This becomes:
The force vector on each particle (F) is the sum of gravitational forces from all other particles. For each particle interaction, the attractive force is inversely proportional to the square distance between them. This force is equal to the gravitational constant (G) multiplied by the unit vector pointing between the particles, divided by the square of this distance. The equation above has this last term slightly rearranged to avoid unnecessary computation.
Given the assumption of unit mass, the force (F) is equal to the acceleration, and so at each timestep, the force vector F is multiplied by the timestep size (dt) and added to the velocity vector. The position of each particle is then updated by the velocity multiplied by the timestep size (dt).
A drag factor (damping) is used to regulate the velocity. At each timestep, the velocity is multiplied by the drag term, slowing the particles. The maximum force between very close particles is also limited for stability; this is achieved via an epsilon term (distEps) which is added to the distance between each particle pairing.
The parameters described in this section can all be adjusted via command line arguments, as follows:
./nbody_cuda numParticles simIterationsPerFrame damping dt distEps G numFrames
Note that numParticles specifies the number of particles simulated, divided by blocksize (i.e. setting numParticles to 50 produces 50*256 particles). simIterationsPerFrame specifies how many steps of the simulation to take before rendering the next frame and numFrames specifies the total number of simulation steps before the program exits. For default values for all of these parameters, refer to sim_param.cpp.
You can get quite a wide range of 'galactic' behaviours by playing with the parameters described above.
Initial velocity of stars is a stable orbital velocity, computed with an implicit value for gravity of G = 1. The default value during the simulation, however, is G = 2. So by default the galaxy collapses inwards quite quickly, but by reducing G closer to 1, you can make a more stable, rotating galaxy.
The damping factor is a drag term. By default damping = 0.999998 but by reducing this value to e.g. 0.999, stars will tend to form local clusters before collapsing in towards the galactic centre.
distEps serves as a stabilising parameter to prevent numerical instability at larger timestep sizes. Setting this value very small (1.0e-10) will produce more 'explosive' simulations. This is unrealistic for n-body gravitational interaction, but it looks dramatic.
If you want to speed up the evolution of the galaxy, set a larger timestep size (dt) or increase the number of steps taken per frame (simIterationsPerFrame). Either change will increase the total simulation time per rendered frame. If you reach a sufficiently high timestep size that you get unstable explosive behaviour, increase the value of distEps and this should stabilise things. Note that there is a separate discussion below about altering the ratio of compute/render time to, for instance, visually highlight a performance difference between platforms.
Render targets for all passes except the last use dimensions a bit larger than the window, to prevent popping. This is used when some effects affect neighboring pixels (bloom, ssao..) and must be taken into account even when off-screen.
Each particle is rendered as a fixed-size flare, generated from a gaussian. Particle color depends on velocity, blue at low speeds and purple at high speeds. Additive blending is set, so dense regions look bright. The render target is RGBA16F, because GL_R11F_G11F_B10F looks yellow on subsequent render passes.
Bloom is applied through a separable Gaussian blur, applied once in the horizontal and then the vertical direction. The 1D Gaussian kernel is computed by RendererGL::gaussKernel and optimized to minimize texel lookups by RendererGL::optimGaussKernel following this guide. At present, a gaussian window of 49 pixels with sigma = 10.0 is used. Multiple passes are possible (ping pong between two RGBA16F frame buffers), though at present we execute only one blur in each direction.
Note that unlike typical bloom processing, there is no extraction of bright light sources prior to blurring, because the scene (bright stars on a dark background) makes this obsolete.
The appearance & performance of the blur is controlled by four variables which are not currently exposed as arguments to nbody_[backend] but which could be manually modified as desired. The two arguments to the gaussKernel function (sigma and halfwidth) effectively define the 'spread' of the blur. Higher values for sigma result in wider blurring, whereas halfwindow defines the actual width of the pixel window which is sampled. Higher values of halfwindow will decrease performance, as more texel lookups are required. As a general rule, when increasing sigma, it will likely be necessary to increase halfwindow to avoid an obvious visual cut-off at the edge of the window. Conversely, a wide halfwindow with a small sigma reduces performance unnecessarily, because texels with negligible contribution will be sampled.
Blur downscaling is a common technique to improve blur performance; the image is downsampled by the factor blur_dsc defined in renderer_gl.cpp, then the regular blur filter is performed, and finally the image is upscaled again. This is a very cheap way of enhancing the blur effect, but there is an associated artefact:
If this artefact is unacceptable, set blur_dsc = 1 to turn off downscaling. Note however that this will significantly reduce the blurriness, and compensating with wider halfwindow or more passes (see below) will cost a lot of rendering time.
Enhanced blurring can also be achieved by executing multiple passes. This is controlled by nPasses, and is set to 1 by default. Due to the dominance of blur in the render pipeline, total rendering time should scale pretty much linearly with nPasses, so increasing it is a potentially expensive option.
The average luminance of the scene is computed from the HDR target into a downscaled R16F target. Then we generate mipmaps to obtain the average luminance on the smallest mipmap (1x1). (Could also be obtained from a 2x2 texture but screen-size targets always seem to resolve down to odd dimensions)
The exposure of the final render is obtained from the average luminance, and the HDR and Bloom targets are combined and converted to LDR. Gamma correction is also applied. Tada.
If you run nbody_cuda on a remote machine with X-forwarding, sending the rendered frames across the net will be a significant bottleneck. This can be worked around by making use of Xvfb which provides a virtual X display. You can then read from the memory mapped file to write to e.g. MP4 output.
The script ./scripts/xvfb.sh runs nbody_cuda in this manner, producing a video file output.mp4. Note that this script will run the simulation until manually terminated.
We've previously discussed the desire for a simulation which is visibly slower when the physics kernel isn't well optimized. With current default settings, the rendering takes longer (~55ms) than the simulation (10ms). However, altering three of the simulation parameters provides almost complete control of the ratio of render to simulation time.
Firstly, the number of particles (numParticles above) has a large effect on the simulation time, as the computation scales with O(n2). By default, 12.8k particles (50 * 256) are rendered, but increasing this to 64k particles (250 * 256), the simulation time increases from 10ms to ~170ms.
Alternatively, simulation time can be arbitrarily raised or lowered by changing both timestep size (dt above) and simulation steps per rendered frame (simIterationsPerFrame, above). By default, a timestep size of 0.005 is used, and 4 simulation steps are taken per rendered frame (Note that scripts/xvfb.sh overrides these default values with dt = 0.001 and simIterationsPerFrame = 5).
To increase the simulation time by a factor of 5, for example, simply divide dt by 5 and multiply simIterationsPerFrame by 5. This will produce almost identical output. Take care with increasing dt to get the opposite effect; above a certain value, the simulation will become unstable & you may see this manifest as unphysical behaviour (very fast moving stars exploding out from the centre). Instability at large dt can be mitigated, to an extent, by increasing distEps or damping.
A significant portion of the rendering time is the bloom filter. The bloom section has some tips about how to control this.
This repo previously reported faster performance from SYCL than CUDA, but this was due to an erroneous translation in the Intel® DPC++ Compatibility Tool from __frsqrt_rn to sycl::rsqrt. The former has higher precision and runs slower than the latter. This has now been rectified so that the original CUDA code calls rsqrt.
With this bug rectified, and without any further modification to the CUDA code or migrated SYCL code, the SYCL code used to be considerably slower because the Intel® DPC++ Compatibility Tool used to insert a cast to double in the rsqrt call:
coords_t inv_dist_cube =
sycl::rsqrt((double)dist_sqr * dist_sqr * dist_sqr);
This was presumably because the tool was unaware of the equivalence of rsqrt and sycl::rsqrt. However, inspecting PTX reveals that the generated instructions are the same, so the cast to double is unnecessary. Removing the cast to double leaves a 40% performance gap between CUDA and SYCL. This is no longer necessary as newer versions of the Intel® DPC++ Compatibility Tool no longer insert the cast.
The root cause of this 40% performance gap appears to be different handling of the branch instruction:
if (i == id) continue;
in the main loop in simulation.dp.cpp. Whereas NVCC handles this via instruction predication, DPC++ generates branch & sync instructions. By replacing this branch instruction with an arithmetic expression:
force += r * inv_dist_cube * (i != id);
in both the CUDA & SYCL code, we get comparable performance between the two using our hardware set up (RTX 3060). For 5 steps of the physical simulation (1 rendered frame) with 12,800 particles, both CUDA and SYCL take ~5.05ms (RTX 3060).