🔒 Locked

Runge Kutta 4 For cloth simulation

Started by Shadx Jan 10, 2006 at 7:52 AM 12 replies 6k views

Original Post

Jan 10, 2006 at 7:52 AM

Hi all, This is a bit of a follow up of the following post : http://www.gamedev.net/community/forums/topic.asp?topic_id=367668 After seeing this post, I went back to my application and changed my Euler integrator for a RK4 integrator. I'm doing cloth simulation with point mass and when using the same Force and dampening for the spring as my Euler simulation, my RK4 becomes so unstable that I eventually lose my cloth. I also use the same time step for both my RK4 and Euler implementation. To make sure that my integrator was working properly, I compiled jjd's code(provided in the post referenced on top) and compared result in a free fall. Beside double vs float, the result of the free fall are similar. On a free fall, the simulation seems to fail on the introduction of small errors on x and z, caused by draging force acting on the particles. I would have thought that the integrator and calculating Hooke force with damping would easily take care of those errors. But they quickly get exagerated and the simulation gets out of control. Locking a point mass in position (or two for simulating a flag), makes the simulation much harder to debug, which is why I'm just currently dealing with free fall. I have to say that I can make it work, even with RK4, by bringing down the force of the spring and by bring up the dampening. Since RK4 should be more stable then Euler, I wouldn't expect that I had to do this though. My Euler simulation is the same except for the integration, how can it be more stable then RK4? On a free fall, the acceleration is constant. Using the integrator provided by jjd as an example, the result for the position in K2 and K3 is always the same, is this the expected result? All the example of derive function always seems to discard the time and use a constant acceleration. When I first started the implementation of RK4 on my cloth, I would recalculate the Hooke force between each node after a step (So calculate RK1 for all pointmass, then recalculate Hooke force, then calculate RK2...). I'm not sure that was the right approach, it just seemed like a very expensive operation to do. So when I changed this, I changed it so that I would perform the full RK4 calculation on each point mass. Before, my force would be different between each step (r1 to rk2 to rk3), so was my acceleration of course. But now, I calculate my Hooke force and then run the integrator, which means the acceleration is constant for every step and K2.pos == k3.pos, which seems pointless. Can anyone point me in the right direction for applying RK4 to cloth...anything beside going with verlet (that's my next project). Thanks, Shadx

jjd

Jan 10, 2006 at 8:22 AM

Hi Shadx,

the short answer is that RK4 is more stable than Euler, so there is something going wrong in the code (either mine or yours). It's not clear to me exactly what is going on in your code. For instance, I'm not sure what the following sentence means

Quote:
Before, my force would be different between each step (r1 to rk2 to rk3)

what are rk1, rk2, and rk3?

I'm going to have a look at my code and let you know if I find any problems. And there are a couple of points that you might want to consider: (1) cloth simulations often suffer from something called stiffness. This will screw up any explicit integrator like RK4 or Euler. However, I'm not sure if this is actually the problem you are having. (2) What timestep are you using? RK4 is not unconditionally stable. If the timestep is too large it will become unstable. If k is the string stiffness and m is the mass of a node in your mesh, then a rule of thumb is that your timestep should be less than sqrt(m/k).

Eelco

Jan 10, 2006 at 8:24 AM

Quote:
Original post by Shadx
Hi all,

This is a bit of a follow up of the following post : http://www.gamedev.net/community/forums/topic.asp?topic_id=367668

After seeing this post, I went back to my application and changed my Euler integrator for a RK4 integrator. I'm doing cloth simulation with point mass and when using the same Force and dampening for the spring as my Euler simulation, my RK4 becomes so unstable that I eventually lose my cloth. I also use the same time step for both my RK4 and Euler implementation.

To make sure that my integrator was working properly, I compiled jjd's code(provided in the post referenced on top) and compared result in a free fall. Beside double vs float, the result of the free fall are similar. On a free fall, the simulation seems to fail on the introduction of small errors on x and z, caused by draging force acting on the particles. I would have thought that the integrator and calculating Hooke force with damping would easily take care of those errors. But they quickly get exagerated and the simulation gets out of control.

Locking a point mass in position (or two for simulating a flag), makes the simulation much harder to debug, which is why I'm just currently dealing with free fall. I have to say that I can make it work, even with RK4, by bringing down the force of the spring and by bring up the dampening. Since RK4 should be more stable then Euler, I wouldn't expect that I had to do this though.

hardly so. in terms of gained stability, its certainly not worth the extra cycles you need to throw at it. i believe it allows for only about twice as big steps compared to euler, but its atleast 4 times as expensive to evaluate.

Quote:

On a free fall, the acceleration is constant. Using the integrator provided by jjd as an example, the result for the position in K2 and K3 is always the same, is this the expected result? All the example of derive function always seems to discard the time and use a constant acceleration. When I first started the implementation of RK4 on my cloth, I would recalculate the Hooke force between each node after a step (So calculate RK1 for all pointmass, then recalculate Hooke force, then calculate RK2...). I'm not sure that was the right approach, it just seemed like a very expensive operation to do.

So when I changed this, I changed it so that I would perform the full RK4 calculation on each point mass. Before, my force would be different between each step (r1 to rk2 to rk3), so was my acceleration of course. But now, I calculate my Hooke force and then run the integrator, which means the acceleration is constant for every step and K2.pos == k3.pos, which seems pointless.

the results should be the same for freefall for any > 1 order method, since freefall is itself a fenomena that is contant in its 2nd order derivative, and thus can be integrated exactly.

however, the hookean force is exactly which gives rise to a not-constant in any derivative at all situation. the point of higher order methods is to better approximate this non-constant force over the timestep with more samples, in order to reach higher accuracy. so if you dont resample your not-constant forces at every intermediate step, using a higher order method serves no point at all.

Quote:

Can anyone point me in the right direction for applying RK4 to cloth...anything beside going with verlet (that's my next project).

Thanks,

Shadx

i dont think there is much merit in simulating cloth with rk4 to begin with. its main advantage over euler is higher accuracy, but i suppose interactivity is more of a concern to you than physical perfection.

Shadx

Jan 10, 2006 at 8:54 AM

Hi jjd,

Sorry about not being more specific. rk1 = hf(x,yi) rk2 = hf(x+dt, yi + rk1/2) etc. At one point I had my simulation working in such a way that I would do:

CalculateForce();
//go through each point mass and get RK1
CalculateRK1();
//using the position and velocity in rk1, recalculate the force
CalculateForce();
CalculateRK2();
etc...

So my force was being recalculated between each steps. I'm not sure that was the right approach which is why I changed it to:
CalculateForce();
IntegrateRK4();

I do understand that I have to carefully chose my time step, my force and my damping factor. Which is why I was comparing with my Euler integrator. My feeling was that if my Euler integrator could take the parameter, so would the RK4. My mass is 0.5, my force K = 10000, damping is 0.5 and my timestep is 0.005. Your calculation would say that 0.007 would be acceptable. I'm suprised damping doesn't enter in that equation (sqrt(m/k)). I forgot to say that I also tried to constraint my springs if they get too far apart, the only thing this changes is how my piece of cloth makes it's exit off the screen (Do you prefer exploding fashion? Or simply a disapearing act?)

Thanks for the help.

Shadx

Shadx

Jan 10, 2006 at 9:02 AM

Quote:
Original post by Eelco
i dont think there is much merit in simulating cloth with rk4 to begin with. its main advantage over euler is higher accuracy, but i suppose interactivity is more of a concern to you than physical perfection.

What I was aiming at was for higher stability so that I could increase my time step a bit. You are right, I'm not really looking for physical perfection, but I'm looking for swift reaction to force or collision with objects. Should I be looking at another method for spring surface or volume?

Shadx

jjd

Jan 10, 2006 at 9:03 AM

Quote:
Original post by Shadx
Hi jjd,

Sorry about not being more specific. rk1 = hf(x,yi) rk2 = hf(x+dt, yi + rk1/2) etc. At one point I had my simulation working in such a way that I would do:

CalculateForce();
//go through each point mass and get RK1
CalculateRK1();
//using the position and velocity in rk1, recalculate the force
CalculateForce();
CalculateRK2();
etc...

This is the correct way to do the calculation.

Quote:

So my force was being recalculated between each steps. I'm not sure that was the right approach which is why I changed it to:
CalculateForce();
IntegrateRK4();

I do understand that I have to carefully chose my time step, my force and my damping factor. Which is why I was comparing with my Euler integrator. My feeling was that if my Euler integrator could take the parameter, so would the RK4. My mass is 0.5, my force K = 10000, damping is 0.5 and my timestep is 0.005. Your calculation would say that 0.007 would be acceptable. I'm suprised damping doesn't enter in that equation (sqrt(m/k)). I forgot to say that I also tried to constraint my springs if they get too far apart, the only thing this changes is how my piece of cloth makes it's exit off the screen (Do you prefer exploding fashion? Or simply a disapearing act?)

Thanks for the help.

Shadx

Don't worry too much about the sqrt(m/k), that's just a rough guideline. Try reducing your stepsize to 0.0007. If that works fine, it might just be that your stepsize was too large.

jjd

Jan 10, 2006 at 9:07 AM

Quote:
Original post by Shadx
Quote:
Original post by Eelco
i dont think there is much merit in simulating cloth with rk4 to begin with. its main advantage over euler is higher accuracy, but i suppose interactivity is more of a concern to you than physical perfection.

What I was aiming at was for higher stability so that I could increase my time step a bit. You are right, I'm not really looking for physical perfection, but I'm looking for swift reaction to force or collision with objects. Should I be looking at another method for spring surface or volume?

Shadx

If you're more concerned about stability, there is a nice article or two on gamasutra about modelling constraints using Verlet. The other option that comes to mind is to use an implicit integrator. However, implicit integration is expensive and more complicated than what you've been doing so you'd better be sure you want/need to go down that path.

ury

Jan 10, 2006 at 10:23 AM

Quote:
Original post by jjd
If you're more concerned about stability, there is a nice article or two on gamasutra about modelling constraints using Verlet. The other option that comes to mind is to use an implicit integrator. However, implicit integration is expensive and more complicated than what you've been doing so you'd better be sure you want/need to go down that path.

Using implicit integration for cloth simulation can be extremely fast. Infact, just as fast as verlet integration.

The only hard/slow/scary thing that you have to do is to precompute the inverse of an extremely sparse NxN matrix.

A very detailed article can be found here.

When you use implicit integration, you can simulate large cloth, increase your time-step and mess up your cloth pretty badly before it "explodes".

Shadx

Jan 10, 2006 at 11:24 AM

Thanks everyone,

I'll start by going back and changing it to what it was originally (recalculating the force between each step). The RK4 calculation on my point mass isn't a complete waste at least and it was a good exercise. I'll look at other methods that a better suited for my goals.

Cheers!

Shadx

jjd

Jan 10, 2006 at 12:02 PM

Quote:
Original post by ury
Quote:
Original post by jjd
If you're more concerned about stability, there is a nice article or two on gamasutra about modelling constraints using Verlet. The other option that comes to mind is to use an implicit integrator. However, implicit integration is expensive and more complicated than what you've been doing so you'd better be sure you want/need to go down that path.

Using implicit integration for cloth simulation can be extremely fast. Infact, just as fast as verlet integration.

The only hard/slow/scary thing that you have to do is to precompute the inverse of an extremely sparse NxN matrix.

A very detailed article can be found here.

When you use implicit integration, you can simulate large cloth, increase your time-step and mess up your cloth pretty badly before it "explodes".

True! I meant that it was more expensive per iteration, but you can recoup the cost by taking much larger steps.

Lutz

Jan 11, 2006 at 7:12 AM

Quote:
Original post by Shadx
Sorry about not being more specific. rk1 = hf(x,yi) rk2 = hf(x+dt, yi + rk1/2) etc. At one point I had my simulation working in such a way that I would do:

This is not quite correct! It must be

rk1 = hf(x     , yi        )rk2 = hf(x+dt/2, yi + rk1/2)rk3 = hf(x+dt/2, yi + rk2/2)rk4 = hf(x+dt  , yi + rk3  )yi = yi + (rk1+2*(rk2+rk3)+rk4)/6

I assume that the function hf multiplies with dt, otherwise rk1-rk4 must be multiplied with dt afterwards, i.e.

rk1 = dt * hf(x, yi)

and so on. Maybe this is the reason for the errors.

Lutz

Shadx

Jan 11, 2006 at 12:01 PM

Quote:
Original post by Lutz

This is not quite correct! It must be
rk1 = hf(x     , yi        )rk2 = hf(x+dt/2, yi + rk1/2)rk3 = hf(x+dt/2, yi + rk2/2)rk4 = hf(x+dt  , yi + rk3  )yi = yi + (rk1+2*(rk2+rk3)+rk4)/6
I assume that the function hf multiplies with dt, otherwise rk1-rk4 must be multiplied with dt afterwards, i.e.
rk1 = dt * hf(x, yi)
and so on. Maybe this is the reason for the errors.

Lutz

hehe, no, I wish it was that simple. I was just showing what I meant by rk and forgot to put in a half time step for my k2 function. My RK4 implementation seems to be alright, for how it applies to individual point mass anyway.

For cloth simulation, I tried it with implementing the RK4 by calculating the force once and then applying the integrator or calculating the the force between every step of the integrator, and neither gave me good result. I will try putting in verlet instead in my next coding session.

Btw, h normally is dt...I also had explained that wrong...should be hf(t+h/2,y+rk1/2)

Thanks though,

Shadx

JeffLander

Jan 12, 2006 at 9:53 AM

As was said earlier, if you don't reevaluate your force function between integration steps (in something like RK4 or Midpoint, etc) you are defeating the purpose of that integration system and it is a waste of time.

Think about it this way with Euler steps for clarity:

Straight forward Euler:
Calc Force
Integrate pos and vel with step = 1.0

If you were wanting a better fit, you might want two evaluations like with a midpoint style method (or just smaller Euler steps):
Calc Force
Integrate pos and vel with step = 0.5
Calc Force
Integrate pos and vel with step = 0.5

If you just do one force calc, it is the same as just taking the full single step with no benefit:
Calc Force
Integrate pos and vel with step = 0.5
Integrate pos and vel with step = 0.5 (second time is using the same force)

That all said, RK4 is probably overkill for your problem. It doesn't give you a lot more stability for stiff systems (like cloth) and has extra expense.

Beware about throwing around the Verlet thing. Verlet is just an integrator like Euler. In fact it is exactly equivalent to what many call semi-implicit Euler. No need to explain why unless you want to go into it. But the "velocity-less" Verlet that Jakobsen talks about in his paper is just an Euler method.

The trick of that paper that everyone loves really not the integrator it is the constraint solving through projection which has nothing to do with the integrator. This is where you just move the points to satisfy rigid constraints. But for cloth you may want rigid links but probably want some stretch just not as much as you get with springs that work in Euler methods.

So the solution (first proposed by Provot a while back) is to use springs with pretty stiff constants that give you good response but don't blow up. Figure out what you want your max over/under stretch to be. Then use projection methods to keep that stretch limit. So if your limit is 10% stretch, project to those limits if the springs go beyond them.

Works great. You get cloth that has some stretch but doesn't look like elastic. Plus it behaves in a more uniform manner in a way that straight projection doesn't.

For simulating more realistic cloth with actual stiff stretch properties, Implicit methods are probably the only solution. But to just get the look, I think Provot had it right.

Shadx

Jan 12, 2006 at 2:46 PM

Jeff,

I went back to my code and changed it so that I recalculate the the forces between each step. I must have done something wrong the first time I tried this approach because it works well now. I can now use a time step of 0.01 with a spring constant of 5000, which works well enough for the simulation I'm doing. Again, thanks for the help on this.

I have started looking at the paper that ury posted and I'm finding it quite interesting. I should have been clear about my goal from the begining instead of blindly trying to use RK4 for something that it's ill equiped to achieve. So before going forward, I'll try finding a solution that scales well enough to the problem at hand.

For those interested, I'm currently looking at the following paper "Interactive animation of structured deformable objects" and it's looking really promissing.

Quote:
So if your limit is 10% stretch, project to those limits if the springs go beyond them.

I'm currently applying constraints based on a maximum possible stretched, because it's rather expensive, I only do it at the end of the rk4 integration instead of doing it between each step.

I also have collision and response done between spheres and AABB, but I've got some work to do since when it gets into constant contact, my point mass can sometime "wiggle" its way through the bounding volume :) At this point, I use collision detection and responsed based on the indivual point mass, without taking into account that those point mass forms a cloth. Same thing for calculating drag force. I need to change all this. At least I'm having fun!

Regards,

Shadx

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.