My particle system is simulating a Cartesian frame of reference when it should be simulating gravity

[pliant crescent]

New update: i tried changing the data i send to the GPU to vec4, giving w the value 1.0f, and i send it like this:

// Position
    glBindBuffer(GL_ARRAY_BUFFER, this->positionBuffer);
    glEnableVertexAttribArray(0);
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, (void*)0);
    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, this->positionBuffer);
    glBufferData(GL_ARRAY_BUFFER, this->nParticles * sizeof(glm::vec4), &this->positions[0], GL_DYNAMIC_DRAW);```

but the result is even weirder. Now It's creating groups of particles like in the picture.

changed the code to this:

vec3 Fe = vec3(0.0);
    vec3 velocity = vec3(0.0);

    // TODO: Experimentar fazer o for apenas para os n bodies mais próximos
    for(int i = 0; i < vertexPosition.length(); i++)
    {
        if (id == i)
        {
            continue;
        }

        vec3 r = vertexPosition[i].xyz - vertexPosition[id].xyz;

        vec3 Gforce = getGravitationalForce(r, vertexMass[id], vertexMass[i]);

        // Calculate the gravitational force
        Fe += Gforce;
    }

    // Calculate the new velocity
    velocity.xyz += Fe * deltaTime * 100.0;


    // Calculate the new position
    vertexPosition[id].xyz += velocity;```

what might be the error?

[pliant crescent]

I fixed it! I should take the w value in consideration!

It's kind of a bummer but we simply can't leave w as it is, we also need to make the calculations to it

[lime atlas]

Can you show me what kind of calculations you are doing? I am kind of curious.

[pliant crescent]

vec4 getGravitationalForce(vec4 r, float m1, float m2)
{
    float d = length(r);
    r = normalize(r);

    // Calculate the gravitational force
    vec4 Gforce = G * m1 * m2 / (d * d) * r;

    return Gforce;
}```

It's Newton's formula for gravitational force

the end result is way wackier than i anticipated

But I guess that's physics for you

[lime atlas]

Ohhh can you try this instead:

[pliant crescent]

(And i know i should calculate d after, but the result was way funnier this way)

[lime atlas]

vec3 getGravitationalForce(vec3 r, float m1, float m2)
{
    float d = length(r);
    float inv_d = 1.0 / d;
    r = r * inv_d ;

    // Calculate the gravitational force
    vec3 Gforce = (G * m1 * m2) * ( inv_d * inv_d ) * r;

    return Gforce;
}

Try this one 👀

Also just optimized it a tiny bit to avoid extra division.

If that works you can turn vec4 into vec3.

And drop the .xyz from r.xyz.

[pliant crescent]

Really interesting but I had to make a few changes, but yeah, the result is basically the same. this is because there are a lot of particles clustered in the center, so all of them get pulled into the center

[lime atlas]

Does it work?

If it does then use it, not having to work on a full vec4 and avoiding divisions can help a tiniest bit with performance lol.

[pliant crescent]

vec3 getGravitationalForce(vec4 r, float m1, float m2)
{
    float d = length(r.xyz);
    float inv_d = 1.0 / d;
    r = r * inv_d ;

    // Calculate the gravitational force
    vec3 Gforce = (G * m1 * m2) * ( inv_d * inv_d ) * r.xyz;

    return Gforce;
}```

had to change the third line

[lime atlas]

Oh okay nice.

[pliant crescent]

since r is vec4 but let me make an adjustment

[lime atlas]

vec3 getGravitationalForce(vec3 r, float m1, float m2)
{
    float d = length(r);
    float inv_d = 1.0 / d;
    r *= inv_d ;

    // Calculate the gravitational force
    vec3 Gforce = (G * m1 * m2) * ( inv_d * inv_d ) * r;

    return Gforce;
}

[pliant crescent]

Bam now it receives a vec3

[lime atlas]

Like this one.

Yay.

I am obsessed with microoptimizing everything lol.

[pliant crescent]

also r should be normalized before calculating d

So they don't jump around

[lime atlas]

That is included.

Look at what is happening to r.

Normalization formula is vec/length(vec)

What is happening to r is:

r * (1.0/length(r))

That is the same as normalize(r).

[pliant crescent]

oh that's true

that's honestly genius

[lime atlas]

The reason I did that is because you need to divide by d^2 afterwards.

Oh no, it is not.

It is just microoptimizations and I love doing them.

But people like @undone hazel are going to start bullying you for thinking about them too much.

So do not focus on microoptimizations on this server.

[pliant crescent]

When making big big projects it's mandatory

[lime atlas]

You are going to get bullied by people here if you start microoptimizing ALUs in public.

Be careful.

[pliant crescent]

Hahaha xD

Thanks!

How much less efficient is using vec4 than vec3?

[lime atlas]

Well it makes everything up to 1.333333333333333333333333333333333333333333... times more expensive.

In that case it is not just ALUs that get worse, it is also the memory bandwidth.

Here is a small example, btw, if you want to see:

vec4 length:
sqrt(x * x + y * y + z * z + w * w)
That is around 5 cheap operations (mul, mad, mad, mad, mul) and one not so each operation (inversesqrt).
vec3 length:
sqrt(x * x + y * y + z * z)
That is around 3 cheap operations (mul, mad, mad, mul) and one not so each operation (inversesqrt).

So in this case the difference of ALU performance is negligible.

But you still need more memory to convey vec4s.

And those memory slots could be taken up by something else.

The answer is: it depends.

@pliant crescent Speaking of, I just noticed another microoptimization.

They make me so excited.

[pliant crescent]

I understand, noted!

[lime atlas]

vec3 getGravitationalForce(vec3 r, float m1, float m2)
{
    float inv_d = inversesqrt(dot(r, r));
    r *= inv_d ;

    // Calculate the gravitational force
    vec3 Gforce = (G * m1 * m2) * ( inv_d * inv_d ) * r;

    return Gforce;
}

Try that.

That is a tiny miserable unimportant microoptimization.

But look how pretty it makes your code.

A dot function is just a dot product, in this case it would do:
r.x * r.x + r.y * r.y + r.z * r.z

[pliant crescent]

yeah i thought of that, i thought you wanted to make the code more readable

[lime atlas]

Well.

I just...

Get my dopamine injections from doing this.

I am mad and inversesqrt'd.

[pliant crescent]

Sometimes making the code smaller makes it less readable

[lime atlas]

inversesqrt is (usually) a tiny bit faster than sqrt btw.

And yeah.

Sqrt is usually implemented as 1.0/inversesqrt.

That is why inversesqrt is faster.

[pliant crescent]

I just need to find a way to make a random position generator where the concetration of particles is the same in all the sphere

there's a high concentration of particles in the center

[lime atlas]

Yeah

Yeah because you are attracting them hard.

Do you want something that just distributes particles all over a sphere?

Within its volume?

[pliant crescent]

yes

[lime atlas]

Without any attraction?

Oh okay that can be done with a simple function.

Wait...

[pliant crescent]

atraction aside

[lime atlas]

Well if you keep the attraction they are going to drop.

Unless you have full inertia.

For that you need to store their velocities.

Are you storing the particle velocities?

[lime atlas]

Oh wait I think I get it now.

Show me your initial particle position generator code.

[pliant crescent]

i have a buffer for it but i don't really use it for this case

[lime atlas]

Yeah okay.

Just show me the generation code 👀

[pliant crescent]

glm::vec3 createPosition(float radius)
{
    float theta = 2.0 * glm::pi<float>() * rand() / (float)RAND_MAX;
    float phi = glm::pi<float>() * rand() / (float)RAND_MAX;

    float positionX = radius * glm::sin(phi) * glm::cos(theta);
    float positionY = radius * glm::sin(phi) * glm::sin(theta);
    float positionZ = radius * glm::cos(phi);

    glm::vec3 position = glm::vec3(positionX, positionY, positionZ);

    return position;
}```

[lime atlas]

There we go, as I thought, polar coordinates!

Because of that you have a very high density of possible positions near the centre.

[pliant crescent]

that's how I learned it in Calculus

[lime atlas]

Because there the sphere contracts.

Yeah you should not be using that one for this.

We need something more... cartesian.

[pliant crescent]

Let me check my calculus notes and the internet

that tip was amazing

Thanks!

[lime atlas]

I know what that works like, the problem is,

That polar spherical coordinates have very high coordinate density near the centre.

What I mean is that if you move in a sphere in the centre you cover a much tinier volume than if you move in a sphere on the outside.

And both those spheres have equal distribution in polar coordinates.

Which results in extremely dense inner spheres and very sparse outer spheres.

[pliant crescent]

I also just remembered something super important, they are all going to the center because the total force of all particles combined is always going to point to the center, hence no matter the particle they will be attracted to it

[lime atlas]

If your intertia gets reduced by something eventually - yes/

[pliant crescent]

That makes a lot of sense!

[lime atlas]

But with no friction - no.

Wait, I am going to make a more even generator...

.

[pliant crescent]

So I should store the velocity?

And keep using it?

[lime atlas]

If you don't want things to just collapse into the centre.

Because if you change only the position, then the things will just snail into the centre and stop there.

Btw, sin/cos can be expensive on the GPU (not on latest Nvidia ones).

Not important but I am a microoptimization junkie.

[pliant crescent]

you're a mathematical genius, that's what you are

I need to add you

[lime atlas]

I am not, I am a queen of microoptimization junkies who gets bullied on this server for getting dopamine rushes from comparing dots vs mads.

[lime atlas]

Btw, is this C?

[pliant crescent]

but in all seriousness I now know that a lot of things that i've learned and forgotten were important and I never knew why

[lime atlas]

I need to know to make sure I don't write it in a wrong language.

[pliant crescent]

c++

[lime atlas]

Oh good okay.

Is there a cube root in C++?

std::cbrt() found it.

[pliant crescent]

I bet glm has it

or standard

Both have mathematical functions

[lime atlas]

glm::vec3 randomPositionInsideASphere(float maximum_radius) {
    float inv_rand_max = 1.0f / (float)RAND_MAX;
    float u = rand() * inv_rand_max; //random number between 0 and 1
    float x1 = rand() * inv_rand_max * 2.0f - 1.0f; //random number between -1 and 1 for x.x
    float x2 = rand() * inv_rand_max * 2.0f - 1.0f; //random number between -1 and 1 for x.y
    float x3 = rand() * inv_rand_max * 2.0f - 1.0f; //random number between -1 and 1 for x.z
    float mag = std::max(0.0001f, std::sqrt(x1*x1 + x2*x2 + x3*x3)); //get length of x, avoid being 0
    float inv_mag = 1.0 / mag;
    x1 *= inv_mag; //turn x into a direction through normalization
    x2 *= inv_mag; //turn x into a direction through normalization
    x3 *= inv_mag; //turn x into a direction through normalization
    float c = std::cbrt(u) * maximum_radius; //rescale c with a cbrt to make the spherical density exactly uniform
    glm::vec3 position = glm::vec3(x1*c, x2*c, x3*c); //get a random point by multiplying random position from the centre with a random direction
    return position;
}

@pliant crescent Done. You can probably vectorize those x1 x2 x3 operations but I don't know C++ well so idk. This should work unless I made some weird mistake but the comments should help with knowing what is going on.

We are rescaling the random radius with a cube root because that exactly creates a uniform distribution for a random direction and a random radius chosen.

And this works entirely without sine or cosine.

The only problem with this method is that some particles might spawn in the exact centre and get stuck there but at least it will be extremely rare this time.

It is really cool how the cube root function literally pulls values away from the centre and prevents them from being overrepresented there.

[pliant crescent]

Hey! Big tip: Don't be dumb like me and don't randomize the radius 🙂

Here's the final code!

glm::vec3 createPosition(float radius) {
    float u = rand() / (float)RAND_MAX; //random number between 0 and 1
    float x1 = rand() / (float)RAND_MAX * 2.0 - 1.0; //random number between -1 and 1 for x.x
    float x2 = rand() / (float)RAND_MAX * 2.0 - 1.0; //random number between -1 and 1 for x.y
    float x3 = rand() / (float)RAND_MAX * 2.0 - 1.0; //random number between -1 and 1 for x.z
    float mag = std::max(0.0001f, std::sqrt(x1*x1 + x2*x2 + x3*x3)); //get length of x, avoid being 0
    float inv_mag = 1.0 / mag;
    x1 *= inv_mag; //turn x into a direction through normalization
    x2 *= inv_mag; //turn x into a direction through normalization
    x3 *= inv_mag; //turn x into a direction through normalization
    float c = std::cbrt(u) * radius; //rescale c with a cbrt to make the spherical density exactly uniform
    // glm::vec3 position = glm::vec3(x1*c, x2*c, x3*c); //get a random point by multiplying random position from the centre with a random direction
    return glm::vec3(x1*c, x2*c, x3*c);
}