Real-Time Physics Simulation Forum

Energy drift

2006-08-24T08:46:31+00:00

Thanks Dirk, I'll look into that.

@Mewert : I guess what he means is that calculating this LCP does not cost a lot, as the error is small the calculation of the right small impulse, should be faster than a whole LCP at the position level.

Statistics: Posted by gee — Thu Aug 24, 2006 8:46 am

Energy drift

2006-08-21T22:48:58+00:00

> 1) The Jacobian is the same for both LCP, so when implemented correctly the second LCP comes for free

Maybe I'm missing something but I don't see a whole lot of savings here.

Here is the system you are solving:
Ax = b subject to LCP constraints
Ay = c "

When solving this system with a sequential impulses technique, you compose the blocks of A on the fly.

So all I get for free by solving both systems at once ( rather than in two passes ) is I avoid 1 maxtrix-matrix multiply, for composing the Jacobian. Which is not much, really. This is at the cost of having to keep around the ( stack, or SPU ) memory allocated for y, c and a temporary vector or two. Plus increasing my code size a bit. There is the nice fact that you get good cache coherency, etc.. so there should be some savings in that respect. When looking at SPU's limited memory, I'm not crazy out about the idea, though.

Have I missed something? Is there anything else I should be getting for free? Can I avoid an additional complemenatiry test or something?

Maybe the free bit you are refering to is the precomputed J elements ( where A = JMJ').

Statistics: Posted by mewert — Mon Aug 21, 2006 10:48 pm

Energy drift

2006-08-18T14:52:56+00:00

I made the same observation (e.g. when implementing cloth using sequentiel impulses). There seems to be a dependency between the number of constraints and the Baumgarte coefficent. The formula I used was:

tau ~ 1 / Sqrt( n )

where n is the number of constraints

Actually I read through a lot of papers recently in order to find better stabilization. The problem with Baumgarte stabilization is that it combats the drift, but not necessarily on the best trajectory. I think for contacts it is really ok, but for several joints building an articulated structure it might be a problem (see e.g. papers of Uri Asher on DAE). I am currently interessted in high quality joints (e.g. to couple animation and physics), but iterative solvers don't seem to be good enough. In the newest paper of R. Weinstein she actually moves from here iterative solver approach to a direct least square solver that minimes the residum...

-Dirk

Statistics: Posted by Dirk Gregorius — Fri Aug 18, 2006 2:52 pm

Energy drift

2006-08-18T14:04:32+00:00

Hello,

1) True.
2) Yes, that I know, I just thought from Erin answer that there was a way with only 1 LCP (trying to resolve 2 systems at the same time did not seem like a good idea so I wondered).
3) And yes I know, I was just using something simple. By the way using a baumgarte correction for now, I have to use a really small coefficient, while with joints 0.8 is fine most of the time. With contacts I have something around 1e-6 (well that is for around 200 balls in contact, for 3-6 balls a higher value is ok).

Thanks

Statistics: Posted by gee — Fri Aug 18, 2006 2:04 pm

Energy drift

2006-08-17T16:43:32+00:00

Basically you are right, you can see it as two LCPs if you want. The advantage of what I suggested is:

1) The Jacobian is the same for both LCP, so when implemented correctly the second LCP comes for free
2) Since you now know the "penalty" impulse to catch up the position error you can simple find the displacement from this without adding to the velocity.

dv = P * m => dx/dt = P * m => dx = P * m * dt

So you apply this displacement directly and dont add it to the velocity

HTH,
-Dirk

PS: Please note that in your formula here you are missing the scaling constant for the Baumgarte term - this:

J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

must be for example

J*W*JT * lambda = -0.1 * C(x) / dt - J*( v(t) + W * f_ext * dt )

See Erin's GDC 2005 presentation for a derivation on how to choose values here...

Statistics: Posted by Dirk Gregorius — Thu Aug 17, 2006 4:43 pm

Energy drift

2006-08-17T08:56:36+00:00

Hello,

I'll try that later and see how this works. But this is still solving 2 LCPs as I thought, not just one

Also another way would be to compute
J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

and then compute :
J*W*JT * lambda = -J*( v(t) + W * f_ext * dt )
using the old lambda values for warm start no ? As the errror should be small the lambda values should be close to new ones I think.

I will try this for contact, and maybe keeping it as the first equation for joints (well I will try both, but as you said it does not seem to be always a good idea).

Thank you for your help.

Statistics: Posted by gee — Thu Aug 17, 2006 8:56 am

Energy drift

2006-08-16T17:12:04+00:00

You can think of it like this. Instead of solving:

J*W*JT * lambda = -C(x) / dt - J*( v(t) + W * f_ext * dt )

you solve the two MLCPs individually:

J*W*JT * lambda1 = -C(x) / dt
J*W*JT * lambda2 = -J*( v(t) + W * f_ext * dt )

Instead of superpositioning both results, you bybass the velocity for the position correction (Baumgarte term).

I think this is little different from the original post-stabilization idea like e.g. in the Cline thesis, where you really compute a second MLCP *after* the position update - with a new Jacobian and delta velocity. The idea of Erin is that the position correction should not effect the velocity state (the persistent momentum) of the rigid bodies. Note that this works actually not well for all scenarios. I think Erin mentioned in this forum that he ran into problems with the brige example. For contacts it seems to be a nice idea and let's contacts appear much more rigid...

-Dirk

Statistics: Posted by Dirk Gregorius — Wed Aug 16, 2006 5:12 pm

Energy drift

2006-08-16T08:16:51+00:00

Hello,

thanks for the reply but I cannot see how you do that. I've tried to look at box2D, but I'm not getting it completely.

If I do a 'normal' MLCP to resolve collisions I will only get one impulse array for my contacts right ? then if I want to have impulses with and without the baumgart error term I need to do this twice no ? I think this is what I'm thinking wrong, but as I'm not the doing doing the LCP I might miss some knowledge here, and my colleague does not seem to know another way as well.

By the way, on monday I tried a position correction only, and stacking was not that better also, and to do a good one I had to recalculate the constraints which is doing twice the job per iteration, for a not so good result (this without contact caching).

Thanks

Statistics: Posted by gee — Wed Aug 16, 2006 8:16 am

Energy drift

2006-08-15T02:18:37+00:00

Yes, the impulse splitting can be done in separate loops, but it may be more efficient to combine the loops. Only velocity impulses should be cached.

Statistics: Posted by Erin Catto — Tue Aug 15, 2006 2:18 am

Energy drift

2006-08-14T08:42:33+00:00

If the drift is due to penetration depth recovery, there are solutions (besides using continuous collision detection of course):
Split the constraint solver impulse, and only use the penetration depth impulse for the integration, but not for the next frames velocity.

See Erin's comment on the ODE mailing list:

I didn't have time to check the patch, but I updated my 2D physics engine
Box2D. This engine uses sequential impulses, which are mathematically
equivalent to PGS. Here's how I did it:

- Each body has a bias velocity that is only used for correcting overlap.
- A separate bias impulse is accumulated for each contact.
- The bias velocity and impulse are initialized to zero each time step (no
warm starting).

The bias impulse and velocity are completely independent of the regular
impulses and velocity, so they could be solved in parallel. However, it was
convenient to solve them together and that should be more cache friendly.

You can get the code here: http://www.gphysics.com/files/Box2D.zip

You can switch between the two methods by using the macro
BIAS_PRESERVES_MOMENTUM in Arbiter.cpp. With the split velocity I was able
to use a large bias factor (0.. This means that overlap is reduced 80% per
step. Such a large bias factor leads to serious bouciness with the merged
velocity.

I also tried using a split velocity for joints. The results were not nearly
so nice, especially for my suspension bridge test. I think that particular
arrangement is more stable with some compliance in the joints.

Erin

Hello,

I'm sorry but I cannot see the difference between doing that and doing another LCP at the position update step. Why doing it at the velocity one ?

My main problem today is that error correction applied on the velocity completely brake the contact caching thing (and of course it's not stable at all with stacking). I have other problems but which can be solved with all I can read here and there.

By the way I'm writing a simulator based mostly on kenny's thesis, and I would like to thank him, and all of you for every thing I've learned through your papers and this forum.

Statistics: Posted by gee — Mon Aug 14, 2006 8:42 am