So, just for fun I'm writing a quick liquid simulation.  Because I like 
liquids, and I'm learning C#, that's why :)

Anyway, I'm running into some problems with the simulation.  Namely, a 
certain instability.  Here's what I'm doing:

For each particle, I loop through all the other particles (i=0..count-1, 
j=i+1..count), and check the distance between them (I use some 
additional bounding to speed things up, but that's not important to my 
problem here).  I compute the forces the two particles place on each 
other, and then add them to each particles' total "force" value.

Once I'm done computing interactions, I loop through all the particles, 
add the force to the velocity, and reset the force to zero.

The way I see it (and granted, I've never written a liquid sim before), 
there are two forces I need to deal with, pertaining to two stages that 
occur when liquid particles collide:

1) When the particles are farther apart, they attract each other.  This 
weak attraction would be the driving force behind surface tension.

2) When the particles are closer together, they push against each other. 
  This force must be much stronger than the first in order to maintain 
some volume.

Now, so far I'm only actually implementing the second force above. 
However, what I see is that the particles tend to settle pretty quickly 
to the floor, and then they clump together.  Apparently, the particles 
on the outside "force" the particles in the middle of the clump closer 
and closer, until they suddenly explode and things go flying everywhere.

I *think* the problem has to do with conservation of momentum.  I 
naively dealt with that by halving the force applied to each particle in 
a pair (one half of the force gets applied to each), but that's not 
quite right, is it?  The momentum of each particle shouldn't simply be 
split evenly with every individual particle it collides with... rather, 
I'll need to collect a list of interactions, and split this particle's 
momentum between *all* of those other particles...

If anyone else here has dealt with similar issues, I'd appreciate 
knowing if I'm on the right track...

-- 
...Ben Chambers
www.pacificwebguy.com

Post a reply to this message

> Once I'm done computing interactions, I loop through all the particles, 
> add the force to the velocity, and reset the force to zero.

Firstly, I wouldn't use integration like this, it's called Euler integration 
and is basically rubbish and well known for being unstable.  A far better 
method, and the de-facto standard in real-time applications, is Runge-Kutta 
4, or RK4 for short.  It is far more stable and much more accurate per CPU 
time than the Euler method.

There are numerous tutorials and explanations online of how to use RK4 in 
exactly the sort of simulation you are doing, but one of my favourites is 
here, it's very clear.  It's got small code snippets in C++ but they are 
very simple and should be easy enough to convert to any language.

http://www.gaffer.org/game-physics/integration-basics

> The way I see it (and granted, I've never written a liquid sim before), 
> there are two forces I need to deal with, pertaining to two stages that 
> occur when liquid particles collide:
>
> 1) When the particles are farther apart, they attract each other.  This 
> weak attraction would be the driving force behind surface tension.
>
> 2) When the particles are closer together, they push against each other. 
> This force must be much stronger than the first in order to maintain some 
> volume.
>
> Now, so far I'm only actually implementing the second force above.

When I did this sort of thing before, I did both forces as one using a 
single function.  Something like:

Force = A * (1/d^4) - B * (1/d^2)

Where d is the distance between two particles, and A,B are constants.

> However, what I see is that the particles tend to settle pretty quickly to 
> the floor, and then they clump together.  Apparently, the particles on the 
> outside "force" the particles in the middle of the clump closer and 
> closer, until they suddenly explode and things go flying everywhere.

Your physics time-step is too big.  WHat is happening is that from one frame 
to the next two particles are moving from far apart to *very* close 
together.  IRL (and with a smaller time-step or better integrator) this 
wouldn't happen, because there would be inbetween steps that prevent the 
particles getting this close together.

There are several ways to fix this (I would do all of them if I were you):

1) Use RK4 instead of Euler (this will let you get away with much bigger 
time-steps without explosion)
2) Use smaller time-steps, eg update the physics 10 times for every graphics 
update
3) Limit the maximum force artificially in the code to avoid explosions

So long as your time step and integration method are accurate enough to 
avoid explosions, it doesn't matter what the formula for Force is, so long 
as you apply it equally and opposite to each pair of particles momentum will 
be conserved.

I would also add some damping force, proportional to speed or speed squared, 
this will also help calm things down.

Post a reply to this message

scott wrote:

> Firstly, I wouldn't use integration like this, it's called Euler 
> integration and is basically rubbish and well known for being unstable.  
> A far better method, and the de-facto standard in real-time 
> applications, is Runge-Kutta 4, or RK4 for short.  It is far more stable 
> and much more accurate per CPU time than the Euler method.

True indeed. Now, if I could just figure out how the hell to implement 
RK4...

[Sure, I have a dictionary of mathematics, and there's Wikipedia and 
there's MathWorld, but all involve pretty abstract mathematics that is 
difficult to figure out how to apply to a really complicated 
differential euqation.]

> There are numerous tutorials and explanations online of how to use RK4 
> in exactly the sort of simulation you are doing, but one of my 
> favourites is here, it's very clear.  It's got small code snippets in 
> C++ but they are very simple and should be easy enough to convert to any 
> language.
> 
> http://www.gaffer.org/game-physics/integration-basics

Mmm, that's a nice article. I'm still having trouble following it 
though. Perhaps a Haskell example?

   integrate t dt a (v,p) = (v + dt*(a t v p), p + dt*v)

OK, well that seems trivial enough. So how about RK4?

   evaluate (v,p) t dt (dv,dp) =
     let
       v' = v + dt*dv
       p' = p + dt*dp
     in  (v', a (t+dt) v' p')

   integrate t dt a (v,p) =
     let
       a = evaluate (v,p) t dt -- er... where's the rest?
       b = evaluate (v,p) t (dt/2) a
       c = evaluate (v,p) t (dt/2) b
       d = evaluate (v,p) t dt     c
       dv = (fst a + 2 * (fst b + fst c) + fst d) / 6
       dp = (snd a + 2 * (snd b + snd c) + snd d) / 6
     in (v + dt*dv, p + dt*dp)

Er... is that actually correct? And what's with the formula for A? I'm 
not really understanding how this is supposed to work...

> Your physics time-step is too big.

Yep, I've seen this happen lots before.

> There are several ways to fix this (I would do all of them if I were you):
> 
> 1) Use RK4 instead of Euler (this will let you get away with much bigger 
> time-steps without explosion)
> 2) Use smaller time-steps, eg update the physics 10 times for every 
> graphics update
> 3) Limit the maximum force artificially in the code to avoid explosions

Indeed.

For my chaos pendulum simulations, I added a constant term to the 
denominator in the "inverse square law" force calculations to prevent 
division by zero and absurdly large forces.

I intended (but never got round to) implementing a system where the 
computer integrates a time step, then integrates it again with twice the 
resolution, and if the results differ too much, it subdivides further...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message

>> http://www.gaffer.org/game-physics/integration-basics
>
> Mmm, that's a nice article. I'm still having trouble following it though. 
> Perhaps a Haskell example?
>
>   integrate t dt a (v,p) = (v + dt*(a t v p), p + dt*v)
>
> OK, well that seems trivial enough. So how about RK4?
>
>   evaluate (v,p) t dt (dv,dp) =
>     let
>       v' = v + dt*dv
>       p' = p + dt*dp
>     in  (v', a (t+dt) v' p')
>
>   integrate t dt a (v,p) =
>     let
>       a = evaluate (v,p) t dt -- er... where's the rest?

The only difference is that in the evaluate function without the supplied 
derivative, you don't need to update the state, simply return the derivative 
at time t using the initial state.

>       b = evaluate (v,p) t (dt/2) a
>       c = evaluate (v,p) t (dt/2) b
>       d = evaluate (v,p) t dt     c
>       dv = (fst a + 2 * (fst b + fst c) + fst d) / 6
>       dp = (snd a + 2 * (snd b + snd c) + snd d) / 6
>     in (v + dt*dv, p + dt*dp)
>
> Er... is that actually correct? And what's with the formula for A? I'm not 
> really understanding how this is supposed to work...

What it's doing is this:

1) Work out the derivative at time t.  This is simply a matter of 
calculating the new acceleration at time t (as the derivative of velocity) 
based on the current state, and returning the current velocity (as the 
derivative of position).

2) Using the derivative you just calculated, a, calculate a new state at 
time t+dt/2 using normal Euler integration.  Return the derivative of the 
new state.

3) Same as step 2, but use b instead of a, and store the result as 
derivative c.

4) Same as 2, but this time use the derivative c to do the whole step dt 
from the first point.

5) Take a weighted average of all the derivatives you just caluclated, a-d, 
this gives the most accurate derivative.

6) Use this weighted average to actually advance your state by dt using 
normal Euler integration.

If you want to actually understand the maths behind it, well that's another 
story (and another Google search, or two :-).

> I intended (but never got round to) implementing a system where the 
> computer integrates a time step, then integrates it again with twice the 
> resolution, and if the results differ too much, it subdivides further...

A good idea, it's called "adaptive time step".

Post a reply to this message

>>   integrate t dt a (v,p) =
>>     let
>>       a = evaluate (v,p) t dt -- er... where's the rest?
> 
> The only difference is that in the evaluate function without the 
> supplied derivative, you don't need to update the state, simply return 
> the derivative at time t using the initial state.

Er... right... I think...

So just return v and the instantenous acceleration at time t?

It's news to me that C++ supports optional function arguments.

[And while we're on the subject, this is the first example of C++ code 
I've ever seen which appears to be comprehensible. Usually it just looks 
like gibberish...]

So anyway, it computes several derivatives and takes their weighted 
average. That seems simple enough in principle. I rather suspect that in 
practice it's spine-tinglingly easy to screw up some tiny detail that 
makes the entire algorithm not work properly. There's *a lot* of 
individual steps there to get right...

>> I intended (but never got round to) implementing a system where the 
>> computer integrates a time step, then integrates it again with twice 
>> the resolution, and if the results differ too much, it subdivides 
>> further...
> 
> A good idea, it's called "adaptive time step".

Mmm, as with all adaptive algorithms, the key is to get the maximum step 
size "small enought", and to avoid infinite recursion. (I also wondered 
if maybe it would be better to just look at the magnitude of force 
changes or something to decide how "stable" the system is at this point...)

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message

> So just return v and the instantenous acceleration at time t?

Yes, exactly.

> It's news to me that C++ supports optional function arguments.

Of course it does (see printf) but you don't use that here.  You just define 
two functions with the same name but with different arguments (this is where 
the article text is a bit misleading IMO).  The compiler knows which one to 
actually call because it matches up the arguments.

> [And while we're on the subject, this is the first example of C++ code 
> I've ever seen which appears to be comprehensible. Usually it just looks 
> like gibberish...]

Maybe because it's written for the purpose of clarity rather than trying to 
make it look "clever" by being totally unreadable.

> So anyway, it computes several derivatives and takes their weighted 
> average. That seems simple enough in principle. I rather suspect that in 
> practice it's spine-tinglingly easy to screw up some tiny detail that 
> makes the entire algorithm not work properly. There's *a lot* of 
> individual steps there to get right...

Yeh, try it first with a simple model like a mass on a spring, or a 
pendulum.  You can then check the results with the actual theory of those 
models, and even compare to standard Euler at the same time.

When I first wrote my RK4 code for the thing I'm currently working on, I 
accidentally updated the initial state for each of the 4 derivative 
calculations.  It was only when I noticed that the velocities being reported 
didn't match up with the actual distances being covered (like covering 1km 
in 10 seconds at a speed of 20 m/s!) that I found the bug.

> Mmm, as with all adaptive algorithms, the key is to get the maximum step 
> size "small enought", and to avoid infinite recursion. (I also wondered if 
> maybe it would be better to just look at the magnitude of force changes or 
> something to decide how "stable" the system is at this point...)

A common method used in car simulations (where the tyres are very stiff and 
hence need a very small time step or the thing explodes easily) is to check 
that energy is (almost) conserved.  Going back to the simple mass/spring 
model, if you keep track of the kinetic energy in the mass and the spring 
energy stored, it should be fairly trivial to include a check if the total 
energy has suddenly spiked.  If so you could then do some smaller time 
steps.  Of course if you are getting input from outside the system you need 
to account for this too.

Post a reply to this message

>> It's news to me that C++ supports optional function arguments.
> 
> Of course it does (see printf) but you don't use that here.  You just 
> define two functions with the same name but with different arguments 
> (this is where the article text is a bit misleading IMO).  The compiler 
> knows which one to actually call because it matches up the arguments.

Ah, right, so it's overloaded but that's not shown in the article?

(BTW, I thought printf() was only for C, not C++?)

>> [And while we're on the subject, this is the first example of C++ code 
>> I've ever seen which appears to be comprehensible. Usually it just 
>> looks like gibberish...]
> 
> Maybe because it's written for the purpose of clarity rather than trying 
> to make it look "clever" by being totally unreadable.

Did the Haskell version look "readable"?

> When I first wrote my RK4 code for the thing I'm currently working on, I 
> accidentally updated the initial state for each of the 4 derivative 
> calculations.  It was only when I noticed that the velocities being 
> reported didn't match up with the actual distances being covered (like 
> covering 1km in 10 seconds at a speed of 20 m/s!) that I found the bug.

One of the nice things about Haskell is that you can't make that 
mistake. Since you never update anything in-place, you can't 
accidentally destroy a result you need to use again later, and by left 
wondering where this "wrong" data is coming from.

...unfortunately you *can* still accidentally return the wrong damn 
varible from your function. :-S

> A common method used in car simulations (where the tyres are very stiff 
> and hence need a very small time step or the thing explodes easily) is 
> to check that energy is (almost) conserved.  Going back to the simple 
> mass/spring model, if you keep track of the kinetic energy in the mass 
> and the spring energy stored, it should be fairly trivial to include a 
> check if the total energy has suddenly spiked.  If so you could then do 
> some smaller time steps.  Of course if you are getting input from 
> outside the system you need to account for this too.

Mmm, it has a flavour...

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message

> Ah, right, so it's overloaded but that's not shown in the article?

Yup, it is in the source code however.

> (BTW, I thought printf() was only for C, not C++?)

Well it is for C, but anything you can do in C you can do in C++.  Not that 
it's a good idea though...

> Did the Haskell version look "readable"?

Not really, I tried to add some to your code to show what I meant with the 
other evaluate function, but I couldn't work out what was going on!

> One of the nice things about Haskell is that you can't make that mistake. 
> Since you never update anything in-place, you can't accidentally destroy a 
> result you need to use again later, and by left wondering where this 
> "wrong" data is coming from.

Seems logical, but without thinking about it too much, how would you keep 
track of the state of the system while it's running?  Are you basically 
saying that variables are write-once in Haskell?  SO I couldn't do something 
like:

currentState = Integrate(currentState);

I would need to instead do:

frame++
state[frame] = Integrate( state[frame-1] )

Post a reply to this message

>> Ah, right, so it's overloaded but that's not shown in the article?
> 
> Yup, it is in the source code however.

Ah, OK. I thought it was just a typo or something...

>> (BTW, I thought printf() was only for C, not C++?)
> 
> Well it is for C, but anything you can do in C you can do in C++.  Not 
> that it's a good idea though...

Ah yes, C++ is a superset of C. But I was thinking (hoping?) that in C++ 
you have something less hellish than prinf() to play with...

>> Did the Haskell version look "readable"?
> 
> Not really, I tried to add some to your code to show what I meant with 
> the other evaluate function, but I couldn't work out what was going on!

No, Haskell doesn't support overloaded functions, so you'd have to make 
some moderately radical (oxymoron?!) alterations to make it work.

> Seems logical, but without thinking about it too much, how would you 
> keep track of the state of the system while it's running?

Ah, grasshopper, you are still thinking in imperative terms. ;-)

A true Haskeller would define a notionally infinite list which 
notionally contains the entire state of the system at every future point 
in time. Examining the elements of this list causes the integration to 
actually be performed, and assuming you let go of the start of the list, 
the GC will delete each state after it has been calculated.

In other words, define an infinite list, and then map some monadic 
function over it which displays the state of the screen or dumps it to a 
file or otherwise "does" something with it. And the result will be more 
or less like this C++ example.

Now, if you want to use *user input* to alter the integration as it 
proceeds, then you must take an somewhat different approach. (See below.)

> Are you 
> basically saying that variables are write-once in Haskell?

You *could* put it that way. It would be more correct to say that a 
variable is *defined* as having a specific value in Haskell. Consider this:

   let node1 = Node {next = node2, prev = node4}
       node2 = Node {next = node3, prev = node1}
       node3 = Node {next = node4, prev = node2}
       node4 = Node {next = node1, prev = node3}
   in ...

As you can see, it's not so much that you can only assign once to a 
variable, it's that you are *defining* what a variable's value *is*. The 
compiler will chase dependencies - even circular ones - and resolve them 
for you, figuring out what order to initialise everything and so forth.

(You can't really speak of assigning several values to the same variable 
"one after the other" because [pure] Haskell function don't have any 
notion of "time".)

> So I couldn't do something like:
> 
> currentState = Integrate(currentState);

No, you couldn't.

You could, however, write

   animate state = do
     dispaly state
     let new_state = integrate state
     animate new_state

which will display the initial state, call your integrate function to 
compute the next state, display that state, and so forth. (This method 
*does* allow you to accept user input and do something different based 
on what the user input is.)

A third way would be to use explicitly mutable state - something like this:

   state <- takeMVar my_state
   writeMVar (integrate state)

But there's no particular advantage in doing that here, and it just 
takes more typing to write out.

-- 
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*

Post a reply to this message

Invisible escribió:
>>> (BTW, I thought printf() was only for C, not C++?)
>>
>> Well it is for C, but anything you can do in C you can do in C++.  Not 
>> that it's a good idea though...
> 
> Ah yes, C++ is a superset of C. But I was thinking (hoping?) that in C++ 
> you have something less hellish than prinf() to play with...

You do.

C:
printf("Number is: %d\n", somenumber);

...and be careful types match between format string and arguments 
passed, and particularly that you passed *enough* arguments

C++:
cout << "Number is: " << somenumber << "\n";

...and there is no duplicated information that could get out-of-sync 
(format and arguments).

It's just you can still use the C way if you want to.

Post a reply to this message