N-Body with Magnum and EnTT

README.md

N-Body

With Magnum and EnTT.

An adaptation of ecs_nbody which is example material for Flecs. Made as a learning exercise of Magnum and EnTT, and data-oriented design in general.

Related

• Pong
• Circles and Boxes pt. 1
• Circles and Boxes pt. 2

What's changed?

Compared to Pong, systems are methods of the application this time, meaning no more global variables. Other than that, it's pretty much the same. Too early to tell whether it's any better (or worse), but it does mean there's the possibility of state in the systems, and that the systems share state. That isn't intentional and could lead to trouble down the line. Ideally, they would carry no state, which is already the case here; they merely share a few constants like overall speed and initial sizes of things that I didn't feel comfortable exposing as global variables like in the Pong example.

Usage

This example depends on Magnum and EnTT.

• Magnum
• EnTT

Build

git clone https://gist.github.com/86d33be269f6e721e847f26e4cb299c4.git nbody
cd nbody
mkdir externals
wget -O externals/entt.hpp https://raw.githubusercontent.com/skypjack/entt/03c4267b84fafca32be264892c81fb0d17d7c2f7/single_include/entt/entt.hpp 
mkdir build
cd build
cmake ..
msbuild nbody.sln
start debug\nbody.exe

Todo

This example illustrates a few interesting bits in Flecs that isn't implemented here, primarily threading and having one system called in response to another system. In this case, that call is made using just the function call. My suspicion is that the reason this is made more complex in the Flecs example is to facilitate threading.

1. Bulk Create Entities Entities are currently being created in a loop, and startup with 20,000 entities takes a few good seconds; that's no good. There's a way of bulk creating instances, but the real hurdle is bulk initialising entities, as they are each given a unique position and velocity on startup which is (I think) where the real bottleneck is at.
2. Call forceSystem in response to gravitySystem This is what I mentioned about one system directly calling another. I think EnTT has some mechanism of hooking systems up via events; it might not be relevant here, but look into that and figure it out.
3. Threading With 2,000 entities, this example just about maxes out 1 core on my machine, with about 10% GPU utilisation. At 20,000 fps drops to about 0.5 frames/sec.
4. More Efficient Neighbor Lookup At the moment, every entity is looking up every other entity to calculate distance that is later added to as attraction force; the closer it is, the stronger the attraction. That lookup is O(n^2). Surely there must be a better method of finding only the closest ones within a given distance? Find it, implement it.
5. Trails Particle-looking circles is cool, but what would be really cool is if they could draw something like a trail after it.

CMakeLists.txt

cmake_minimum_required(VERSION 3.1)

project(nbody)

# Add module path in case this is project root
if(PROJECT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
    set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/../../modules/" ${CMAKE_MODULE_PATH})
endif()

find_package(Magnum REQUIRED
    GL
    MeshTools
    Primitives
    Shaders
    Sdl2Application)

set_directory_properties(PROPERTIES CORRADE_USE_PEDANTIC_FLAGS ON)

add_executable(nbody nbody.cpp)
target_link_libraries(magnum-primitives PRIVATE
    Magnum::Application
    Magnum::GL
    Magnum::Magnum
    Magnum::MeshTools
    Magnum::Primitives
    Magnum::Shaders)

install(TARGETS nbody DESTINATION ${MAGNUM_BINARY_INSTALL_DIR})

nbody.cpp

#include <vector>

#include <Magnum/Platform/Sdl2Application.h>
#include <Magnum/GL/DefaultFramebuffer.h>
#include <Magnum/GL/Mesh.h>
#include <Magnum/GL/Renderer.h>
#include <Magnum/MeshTools/Compile.h>
#include <Magnum/MeshTools/Transform.h>
#include <Magnum/Primitives/Circle.h>
#include <Magnum/Primitives/Square.h>
#include <Magnum/Shaders/Flat.h>
#include <Magnum/Trade/MeshData2D.h>
#include <Magnum/Timeline.h>

#include "externals/entt.hpp"

using namespace Magnum;
using namespace Math::Literals;  // For _rgbf

static entt::registry Registry;

class GravityParam;


class Nbody : public Platform::Application {
public:
    explicit Nbody(const Arguments& arguments);

private:
    void drawEvent() override;

    void initSystem();
    void forceSystem(GravityParam* param);
    void gravitySystem();
    void movementSystem();
    void colorSystem();
    void renderSystem();

    uint16_t _nbodies { 1000 };        /* Number of entities */
    float _centralMass { 12000.0f };
    float _initalC { 12000.0f };
    uint8_t _speed { 2 };
    float _zoom { 0.1f };
    float _drag { 10000.0f };

    Timeline _timeline;
    uint16_t _count { 0 };
};

/* Return "random" value between 0-1 */
auto random01() -> float {
    return static_cast<float>((double)rand() / (RAND_MAX + 1.0));
}

/* Components */

using Mass = float;
using Color = Color3;

struct Circle {
    float radius;
};

struct Position : public Vector2 {
    using Vector2::Vector2;
};

struct Velocity : public Vector2 {
    using Vector2::Vector2;
};

struct GravityParam {
    entt::entity me;
    const Position* position;
    Velocity force;
};


/* Systems */

void Nbody::initSystem() {

    // Parameters of this system
    const auto maxRadius = 70.0f;
    const auto massVariation = 0.8f;
    float baseMass { 0.1f };

    Registry.view<Position, Velocity, Mass, Circle>().each([=](auto& pos,
                                                               auto& vel,
                                                               auto& mass,
                                                               auto& circle) {
        pos.x() = rand() % 8000 - 4000;
        pos.y() = rand() % 200 - 100;
        mass = baseMass + random01() * massVariation;

        if (pos.x() || pos.y()) {
            float radius = pos.length();
            auto normal = pos / radius;
            auto rotation = normal.perpendicular();
            float velocity = sqrt(_initalC / radius / mass / _speed);

            vel.x() = rotation.x() * velocity;
            vel.y() = rotation.y() * velocity;
        }

        circle.radius = maxRadius
                      * (mass / (baseMass + massVariation))
                      + 1;
    });
}


void Nbody::forceSystem(GravityParam* param) {
    Registry.view<Position, Mass>().each([=](auto entity,
                                             const auto& pos,
                                             const auto& mass) {

        if (entity != param->me) {
            auto diff = (*param->position) - pos;
            auto distance = dot(diff, diff);

            if (distance < _drag) {
                distance = _drag;
            }

            float distance_sqr = sqrt(distance);
            float force = mass / distance;

            diff *= force / distance_sqr;
            param->force += diff;
        }
    });
}


void Nbody::gravitySystem() {
    Registry.view<Velocity, Position, Mass>().each([&](auto entity,
                                                       auto& vel,
                                                       const auto& pos,
                                                       const auto& mass) {
        GravityParam param;
        param.me = entity;
        param.position = &pos;
        param.force = {0, 0};

        this->forceSystem(&param);

        vel.x() += param.force.x() / mass;
        vel.y() += param.force.y() / mass;
    });
}


void Nbody::movementSystem() {
    Registry.view<Position, Velocity>().each([&](auto& pos, const auto& vel) {
        pos.x() -= _speed * vel.x();
        pos.y() -= _speed * vel.y();
    });
}


void Nbody::colorSystem() {
    Registry.view<Color, Velocity>().each([&](auto& color, const auto& vel) {
        float f = vel.length() / 8 * sqrt(_speed);
        if (f > 1.0) f = 1.0;

        float f_red = f - 0.2f;
        if (f_red < 0) f_red = 0.0f;
        f_red /= 0.8f;

        float f_green = f - 0.7f;
        if (f_green < 0) f_green = 0.0f;
        f_green /= 0.3f;

        color = Color3{ f_red, f_green, f * 0.4f + 0.6f };
    });
}



void Nbody::renderSystem() {
    GL::Renderer::setClearColor(0x000000_rgbf);
    GL::defaultFramebuffer.clear(GL::FramebufferClear::Color | GL::FramebufferClear::Depth);
    GL::defaultFramebuffer.setViewport({{}, windowSize()});

    auto shape = MeshTools::compile(Primitives::circle2DSolid(20));
    auto shader = Shaders::Flat2D{};

    auto projection = Matrix3::projection(Vector2(windowSize() / _zoom));

    Registry.view<Position, Circle, Color>().each([&](auto& pos, auto& circle, auto& color) {
        shader.setTransformationProjectionMatrix(
            projection *
            Matrix3::translation(pos) *
            Matrix3::scaling(Vector2(circle.radius / 2, circle.radius / 2))
        );

        shader.setColor(color);
        shape.draw(shader);
    });
}


Nbody::Nbody(const Arguments& arguments) :
    Platform::Application{
        arguments,
        Configuration{}.setTitle("Nbody")
                       .setSize({640, 480}),
        GLConfiguration{}.setSampleCount(8)
    }
{
    this->setSwapInterval(0);

    entt::entity entity;
    for (int ii = 0; ii < _nbodies; ii++) {
        entity = Registry.create();
        Registry.assign<Position>(entity);
        Registry.assign<Velocity>(entity);
        Registry.assign<Mass>(entity);
        Registry.assign<Color>(entity);
        Registry.assign<Circle>(entity);
    }

    // Initialize values
    this->initSystem();

    // Have all circles move about the center
    Registry.replace<Position>(entity, 0.0f, 0.0f);
    Registry.replace<Velocity>(entity, 0.0f, 0.0f);
    Registry.replace<Mass>(entity, _centralMass);
    Registry.replace<Circle>(entity, 2.0f);
}

void Nbody::drawEvent() {
    this->gravitySystem();
    this->movementSystem();
    this->colorSystem();
    this->renderSystem();

    swapBuffers();
    redraw();
}

// Call cross-platform application loop
MAGNUM_APPLICATION_MAIN(Nbody)

Argh! Got confused by the forks and replied on the wrong one. Sorry @alanjfs ... that comment should have gone to your gist.

No problem, good to see @skypjack having a look at how this could be improved. Will be keeping an eye on this. :)

creating all buffers, setting up a mesh and compiling the shader 60 times a second will never be fast -- do that just once ;)

I suppose; but I didn't bother in this case, because the vast majority of time was spent in that loop happening after creating those meshes and shaders. I'm thinking either instancing, or merging meshes into one could help.

But even that probably won't help much, because I think rendering overall isn't the bottleneck here, but calculating the distance to each nbody, per nbody. The original example mentions this as well, and I had a chat with the author on Gitter about it. The algorithm itself has got a lot of room for optimisation. It's currently O(N^2).

to keep the circles as circles easiest would be to texture them (a line shader that would make drawing tiny circles easier is on my TODO list)

That's an interesting idea. I would have thought that sampling a texture would be more expensive than computing the few number of vertices that make up each circle. But maybe not.

With that in mind, would it make sense to draw quads this way too? :O That is, to draw a triangle large enough to fit a quad as texture, such that you only ever draw one triangle per entity?

This is a great example material, would you be willing to submit this to the magnum-examples repo (and the Pong, too)?

Sure, I'll put it on the TODO. :)

I would have thought that sampling a texture would be more expensive than computing the few number of vertices that make up each circle

Depends on where that vertex transformation happens -- if it would be on the GPU (as it is now), then the vertex transformation would probably be faster (unless you get sub-pixel triangles, which you probably are); but if you would be putting everything into a single buffer uploaded from the CPU every time, then the less you have to calculate (and upload) from the CPU, the better. Texturing a triangle this way could be even faster since it's one primitive per circle less (but you get some unnecessary overdraw, OTOH).

calculating the distance to each nbody, per nbody

The "brute force" optimization would be moving this whole calculation to the GPU and doing some transform feedback. But then you kinda don't need any ECS anymore...