search

MADEAPPS / newton-dynamics

Newton Dynamics is an integrated solution for real time simulation of physics environments.
http://www.newtondynamics.com
Other
928 stars 183 forks source link

Orbital assertion failure #319

Closed lperkin1 closed 1 year ago

lperkin1 commented 1 year ago

Got another assertion failure. I'm sure it's something I'm doing wrong, but I don't know how to figure out what. The simulation seems to be stable when it's not in debug, and I get the results of collisions I would expect, so I don't know if this is a problem or not. That said, at these scales, I don't know if I would notice incorrect results anyway.

Here's a test that causes the failure. I put the simulation time step to 100 so it doesn't take as long to have the problem, but it's the same if I have the 1/60th of a second time step (it just takes 30k iterations to get there).

#include "ndNewton.h"
#include <gtest/gtest.h>
class csBodyNotify : public ndBodyNotify {
public:
    csBodyNotify(ndBigVector const& defaultGravity) : ndBodyNotify(defaultGravity) {}

    virtual ~csBodyNotify() {}
    virtual void OnApplyExternalForce(ndInt32 threadIndex, ndFloat32 timestep) {
        ndVector dir(GetBody()->GetPosition() * -1);
        ndFloat32 dis(dir.m_x * dir.m_x + dir.m_y * dir.m_y + dir.m_z * dir.m_z);
        dir = dir.Normalize();
        ndFloat32 f(ndFloat32(6.67e-11f) * (m_mass / dis));
        SetGravity(dir * f);
        ndBodyNotify::OnApplyExternalForce(threadIndex, timestep);
    }

    ndFloat32 m_mass;

};

TEST(Extremes, OrbitalCollisions)
{
    ndWorld world;
    world.SetSubSteps(1);
    ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000)));

    ndMatrix matrix(ndGetIdentityMatrix());

    ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

    ndBodyKinematic* staticbody = new ndBodyKinematic();
    staticbody->SetCollisionShape(planetshape);
    staticbody->SetMatrix(matrix);
    staticbody->SetMassMatrix(ndFloat32(5.9722e+24), planetshape);
    ndSharedPtr<ndBody> staticPtr(staticbody);
    world.AddBody(staticPtr);

    matrix.m_posit.m_y = 6357000 + 1021140;

    for (int i = 0; i < 20; i++) {
        ndBodyDynamic* movingbody = new ndBodyDynamic();
        csBodyNotify* n = new csBodyNotify(ndVector(ndFloat32(0), ndFloat32(0), ndFloat32(0), ndFloat32(0)));
        n->m_mass = ndFloat32(5.9722e+24) * ndFloat32(1.5f);
        movingbody->SetNotifyCallback(n);
        movingbody->SetCollisionShape(shapeinst);
        movingbody->SetMatrix(matrix);
        movingbody->SetMassMatrix(ndFloat32(10), shapeinst);
        movingbody->SetVelocity(ndVector(ndFloat32(9000), ndFloat32(0), ndFloat32(0), ndFloat32(0)));
        //movingbody->SetDebugMaxLinearAndAngularIntegrationStep(ndPi, 2.0f);
        ndSharedPtr<ndBody> movingPtr(movingbody);
        world.AddBody(movingPtr);
        matrix.m_posit.m_y += 100005;
    }

    for (int i = 0; i < 48; i++)
    {
        //world.Update(1.0f / 60.0f);
        world.Update(100.0f);
        world.Sync();
    }

    world.CleanUp();
}

I don't really want to pull you away from the controller for my being stupid, but hopefully it's something easy to spot.

JulioJerez commented 1 year ago

ndAssert(m_separationDistance < ndFloat32(1.0e6f)); that assert of to find out when the distance between two object is too far away, and do no ne to keep tha pair, but take so long of a time step, I am no sure what can happens

I suppose that //world.Update(1.0f / 60.0f); world.Update(100.0f);

is just a test, that will break the integrators, you say that I happen at normal time step, in 30k step, I will try that and see what happens.

On Tue, Jul 4, 2023 at 7:13 PM Luke Perkins @.***> wrote:

Got another assertion failure. I'm sure it's something I'm doing wrong, but I don't know how to figure out what. The simulation seems to be stable when it's not in debug, and I get the results of collisions I would expect, so I don't know if this is a problem or not. That said, at these scales, I don't know if I would notice incorrect results anyway.

Here's a test that causes the failure. I put the simulation time step to 100 so it doesn't take as long to have the problem, but it's the same if I have the 1/60th of a second time step (it just takes 30k iterations to get there).

include "ndNewton.h"

include <gtest/gtest.h>

class csBodyNotify : public ndBodyNotify { public: csBodyNotify(ndBigVector const& defaultGravity) : ndBodyNotify(defaultGravity) {}

virtual ~csBodyNotify() {} virtual void OnApplyExternalForce(ndInt32 threadIndex, ndFloat32 timestep) { ndVector dir(GetBody()->GetPosition() -1); ndFloat32 dis(dir.m_x dir.m_x + dir.m_y dir.m_y + dir.m_z dir.m_z); dir = dir.Normalize(); ndFloat32 f(ndFloat32(6.67e-11f) (m_mass / dis)); SetGravity(dir f); ndBodyNotify::OnApplyExternalForce(threadIndex, timestep); }

ndFloat32 m_mass;

};

TEST(Extremes, OrbitalCollisions) { ndWorld world; world.SetSubSteps(1); ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000)));

ndMatrix matrix(ndGetIdentityMatrix());

ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

ndBodyKinematic* staticbody = new ndBodyKinematic(); staticbody->SetCollisionShape(planetshape); staticbody->SetMatrix(matrix); staticbody->SetMassMatrix(ndFloat32(5.9722e+24), planetshape); ndSharedPtr staticPtr(staticbody); world.AddBody(staticPtr);

matrix.m_posit.m_y = 6357000 + 1021140;

for (int i = 0; i < 20; i++) { ndBodyDynamic movingbody = new ndBodyDynamic(); csBodyNotify n = new csBodyNotify(ndVector(ndFloat32(0), ndFloat32(0), ndFloat32(0), ndFloat32(0))); n->m_mass = ndFloat32(5.9722e+24) * ndFloat32(1.5f); movingbody->SetNotifyCallback(n); movingbody->SetCollisionShape(shapeinst); movingbody->SetMatrix(matrix); movingbody->SetMassMatrix(ndFloat32(10), shapeinst); movingbody->SetVelocity(ndVector(ndFloat32(9000), ndFloat32(0), ndFloat32(0), ndFloat32(0))); //movingbody->SetDebugMaxLinearAndAngularIntegrationStep(ndPi, 2.0f); ndSharedPtr movingPtr(movingbody); world.AddBody(movingPtr); matrix.m_posit.m_y += 100005; }

for (int i = 0; i < 48; i++) { //world.Update(1.0f / 60.0f); world.Update(100.0f); world.Sync(); }

world.CleanUp(); }

I don't really want to pull you away from the controller for my being stupid, but hopefully it's something easy to spot.

— Reply to this email directly, view it on GitHub https://github.com/MADEAPPS/newton-dynamics/issues/319, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6EPJA55G255DBDHLIIPWDXOTEU5ANCNFSM6AAAAAAZ6J5KVA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

JulioJerez commented 1 year ago

so you have ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000))); ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

I am guessing one is the earth, and the other something likea large asteroid, is that right? I guess that at that scale, closest distance of a 1000000 meters is not unreasonable.

I will debug it, to see if this is right, but My guess is that the check has to be float dependent. many of those hardcoded checks were set since the engine was for float32, anyway I will check it tomorrow morning.

On Tue, Jul 4, 2023 at 8:21 PM Julio Jerez @.***> wrote:

ndAssert(m_separationDistance < ndFloat32(1.0e6f)); that assert of to find out when the distance between two object is too far away, and do no ne to keep tha pair, but take so long of a time step, I am no sure what can happens

I suppose that //world.Update(1.0f / 60.0f); world.Update(100.0f);

is just a test, that will break the integrators, you say that I happen at normal time step, in 30k step, I will try that and see what happens.

On Tue, Jul 4, 2023 at 7:13 PM Luke Perkins @.***> wrote:

Got another assertion failure. I'm sure it's something I'm doing wrong, but I don't know how to figure out what. The simulation seems to be stable when it's not in debug, and I get the results of collisions I would expect, so I don't know if this is a problem or not. That said, at these scales, I don't know if I would notice incorrect results anyway.

Here's a test that causes the failure. I put the simulation time step to 100 so it doesn't take as long to have the problem, but it's the same if I have the 1/60th of a second time step (it just takes 30k iterations to get there).

include "ndNewton.h"

include <gtest/gtest.h>

class csBodyNotify : public ndBodyNotify { public: csBodyNotify(ndBigVector const& defaultGravity) : ndBodyNotify(defaultGravity) {}

virtual ~csBodyNotify() {} virtual void OnApplyExternalForce(ndInt32 threadIndex, ndFloat32 timestep) { ndVector dir(GetBody()->GetPosition() -1); ndFloat32 dis(dir.m_x dir.m_x + dir.m_y dir.m_y + dir.m_z dir.m_z); dir = dir.Normalize(); ndFloat32 f(ndFloat32(6.67e-11f) (m_mass / dis)); SetGravity(dir f); ndBodyNotify::OnApplyExternalForce(threadIndex, timestep); }

ndFloat32 m_mass;

};

TEST(Extremes, OrbitalCollisions) { ndWorld world; world.SetSubSteps(1); ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000)));

ndMatrix matrix(ndGetIdentityMatrix());

ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

ndBodyKinematic* staticbody = new ndBodyKinematic(); staticbody->SetCollisionShape(planetshape); staticbody->SetMatrix(matrix); staticbody->SetMassMatrix(ndFloat32(5.9722e+24), planetshape); ndSharedPtr staticPtr(staticbody); world.AddBody(staticPtr);

matrix.m_posit.m_y = 6357000 + 1021140;

for (int i = 0; i < 20; i++) { ndBodyDynamic movingbody = new ndBodyDynamic(); csBodyNotify n = new csBodyNotify(ndVector(ndFloat32(0), ndFloat32(0), ndFloat32(0), ndFloat32(0))); n->m_mass = ndFloat32(5.9722e+24) * ndFloat32(1.5f); movingbody->SetNotifyCallback(n); movingbody->SetCollisionShape(shapeinst); movingbody->SetMatrix(matrix); movingbody->SetMassMatrix(ndFloat32(10), shapeinst); movingbody->SetVelocity(ndVector(ndFloat32(9000), ndFloat32(0), ndFloat32(0), ndFloat32(0))); //movingbody->SetDebugMaxLinearAndAngularIntegrationStep(ndPi, 2.0f); ndSharedPtr movingPtr(movingbody); world.AddBody(movingPtr); matrix.m_posit.m_y += 100005; }

for (int i = 0; i < 48; i++) { //world.Update(1.0f / 60.0f); world.Update(100.0f); world.Sync(); }

world.CleanUp(); }

I don't really want to pull you away from the controller for my being stupid, but hopefully it's something easy to spot.

— Reply to this email directly, view it on GitHub https://github.com/MADEAPPS/newton-dynamics/issues/319, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6EPJA55G255DBDHLIIPWDXOTEU5ANCNFSM6AAAAAAZ6J5KVA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

JulioJerez commented 1 year ago

is this right?

virtual void OnApplyExternalForce(ndInt32 threadIndex, ndFloat32 timestep)

{ ndVector dir(GetBody()->GetPosition() -1); ndFloat32 dis(dir.m_x dir.m_x + dir.m_y dir.m_y + dir.m_z dir.m_z); dir = dir.Normalize(); ndFloat32 g(ndFloat32(6.67e-11f) * (m_mass / dis)); SetGravity(dir.Scale (g)); ndBodyNotify::OnApplyExternalForce(threadIndex, timestep); }

I am reading a value of G of 10, but those two objects are really far away, I would expect the gravity to be quite small at that distance. anyway I committed, see if that is better.

On Tue, Jul 4, 2023 at 8:32 PM Julio Jerez @.***> wrote:

so you have ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000))); ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

I am guessing one is the earth, and the other something likea large asteroid, is that right? I guess that at that scale, closest distance of a 1000000 meters is not unreasonable.

I will debug it, to see if this is right, but My guess is that the check has to be float dependent. many of those hardcoded checks were set since the engine was for float32, anyway I will check it tomorrow morning.

On Tue, Jul 4, 2023 at 8:21 PM Julio Jerez @.***> wrote:

ndAssert(m_separationDistance < ndFloat32(1.0e6f)); that assert of to find out when the distance between two object is too far away, and do no ne to keep tha pair, but take so long of a time step, I am no sure what can happens

I suppose that //world.Update(1.0f / 60.0f); world.Update(100.0f);

is just a test, that will break the integrators, you say that I happen at normal time step, in 30k step, I will try that and see what happens.

On Tue, Jul 4, 2023 at 7:13 PM Luke Perkins @.***> wrote:

Got another assertion failure. I'm sure it's something I'm doing wrong, but I don't know how to figure out what. The simulation seems to be stable when it's not in debug, and I get the results of collisions I would expect, so I don't know if this is a problem or not. That said, at these scales, I don't know if I would notice incorrect results anyway.

Here's a test that causes the failure. I put the simulation time step to 100 so it doesn't take as long to have the problem, but it's the same if I have the 1/60th of a second time step (it just takes 30k iterations to get there).

include "ndNewton.h"

include <gtest/gtest.h>

class csBodyNotify : public ndBodyNotify { public: csBodyNotify(ndBigVector const& defaultGravity) : ndBodyNotify(defaultGravity) {}

virtual ~csBodyNotify() {}
virtual void OnApplyExternalForce(ndInt32 threadIndex, ndFloat32 timestep) {
    ndVector dir(GetBody()->GetPosition() * -1);
    ndFloat32 dis(dir.m_x * dir.m_x + dir.m_y * dir.m_y + dir.m_z * dir.m_z);
    dir = dir.Normalize();
    ndFloat32 f(ndFloat32(6.67e-11f) * (m_mass / dis));
    SetGravity(dir * f);
    ndBodyNotify::OnApplyExternalForce(threadIndex, timestep);
}

ndFloat32 m_mass;

};

TEST(Extremes, OrbitalCollisions) { ndWorld world; world.SetSubSteps(1); ndShapeInstance shapeinst(new ndShapeSphere(ndFloat32(50000)));

ndMatrix matrix(ndGetIdentityMatrix());

ndShapeInstance planetshape(new ndShapeSphere(ndFloat32(6357000)));

ndBodyKinematic* staticbody = new ndBodyKinematic();
staticbody->SetCollisionShape(planetshape);
staticbody->SetMatrix(matrix);
staticbody->SetMassMatrix(ndFloat32(5.9722e+24), planetshape);
ndSharedPtr<ndBody> staticPtr(staticbody);
world.AddBody(staticPtr);

matrix.m_posit.m_y = 6357000 + 1021140;

for (int i = 0; i < 20; i++) {
    ndBodyDynamic* movingbody = new ndBodyDynamic();
    csBodyNotify* n = new csBodyNotify(ndVector(ndFloat32(0), ndFloat32(0), ndFloat32(0), ndFloat32(0)));
    n->m_mass = ndFloat32(5.9722e+24) * ndFloat32(1.5f);
    movingbody->SetNotifyCallback(n);
    movingbody->SetCollisionShape(shapeinst);
    movingbody->SetMatrix(matrix);
    movingbody->SetMassMatrix(ndFloat32(10), shapeinst);
    movingbody->SetVelocity(ndVector(ndFloat32(9000), ndFloat32(0), ndFloat32(0), ndFloat32(0)));
    //movingbody->SetDebugMaxLinearAndAngularIntegrationStep(ndPi, 2.0f);
    ndSharedPtr<ndBody> movingPtr(movingbody);
    world.AddBody(movingPtr);
    matrix.m_posit.m_y += 100005;
}

for (int i = 0; i < 48; i++)
{
    //world.Update(1.0f / 60.0f);
    world.Update(100.0f);
    world.Sync();
}

world.CleanUp();

}

I don't really want to pull you away from the controller for my being stupid, but hopefully it's something easy to spot.

— Reply to this email directly, view it on GitHub https://github.com/MADEAPPS/newton-dynamics/issues/319, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6EPJA55G255DBDHLIIPWDXOTEU5ANCNFSM6AAAAAAZ6J5KVA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

lperkin1 commented 1 year ago

The assert is working now. Thanks.

Good catch on the gravity calculation. I've got a 1.5 kg object at 7,378 km from a 5.9722e24 kg planet, so I think newtons should be N = ((6.6743e-11 1.5 5.9722e24) / 54,436,949,859,600) = 10.9834. Looks like I then just applied that as the acceleration directly instead of dividing it by the mass. I knew something was wrong (can't have more than 9.81 acceleration in orbit of an Earth sized object) but I was more concerned with the interface than the actual simulation.

I'll be sure to fix it in my examples though so I don't look like too much of an idiot :)

JulioJerez commented 1 year ago

Here is a trick that you can apply to make the sim more stable.

Since the gravity is given by an implicit formulat, You can using a simi implicitly method.

That is instead of using the gravity at position p and time t, you can actually calculate the partial derivatives respect to time and position. This gives you a linear system of 4x4

The you locally solve what would be the position at time t + dt, by solve that system. And them use the gravity at that predicted position as the gravity at time T.

The newton engine uses that method for all joint systems, by for constant gravity does nothing since the time derivative of a const is zero.

On Wed, Jul 5, 2023, 10:50 AM Luke Perkins @.***> wrote:

Closed #319 https://github.com/MADEAPPS/newton-dynamics/issues/319 as completed.

— Reply to this email directly, view it on GitHub https://github.com/MADEAPPS/newton-dynamics/issues/319#event-9734669606, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6EPJB2V5G7BUBLPRN5XQ3XOWSNHANCNFSM6AAAAAAZ6J5KVA . You are receiving this because you commented.Message ID: @.***>

JulioJerez commented 1 year ago

What that does is that introduces a drag that such up some momentum from the numerical integration.

On Wed, Jul 5, 2023, 11:36 AM Julio Jerez @.***> wrote:

Here is a trick that you can apply to make the sim more stable.

Since the gravity is given by an implicit formulat, You can using a simi implicitly method.

That is instead of using the gravity at position p and time t, you can actually calculate the partial derivatives respect to time and position. This gives you a linear system of 4x4

The you locally solve what would be the position at time t + dt, by solve that system. And them use the gravity at that predicted position as the gravity at time T.

The newton engine uses that method for all joint systems, by for constant gravity does nothing since the time derivative of a const is zero.

On Wed, Jul 5, 2023, 10:50 AM Luke Perkins @.***> wrote:

Closed #319 https://github.com/MADEAPPS/newton-dynamics/issues/319 as completed.

— Reply to this email directly, view it on GitHub https://github.com/MADEAPPS/newton-dynamics/issues/319#event-9734669606, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6EPJB2V5G7BUBLPRN5XQ3XOWSNHANCNFSM6AAAAAAZ6J5KVA . You are receiving this because you commented.Message ID: @.***>

lperkin1 commented 1 year ago

I remember seeing you discussing that with someone on the forum. I'll have to see if I can find that again and look a little closer. I'm just a humble country programmer, so while I understand the concept, it'll take a bit to figure the math. I only barely got through calculus and that was a long time ago :)

Building my wrapper, I occasionally look through your code and see what you're doing. Half the time I can figure the purpose with some headache, but not why it works, the other half just seems entirely like magic. So I do appreciate you taking the time to drop some hints for us idiots.