 |
|
|
|
|
|
| |
| |
|
|
|
|
| |
| |
|
|
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
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> 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.
That does indeed seem like at least *some* improvement...
(In C, almost every single time anybody touched printf(), the compiler
would utter something about "suspicious pointer conversion". To this day
nobody knows what that actually means... Even the tutors didn't seem to
know.)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Assuming that you're writing this simulation for graphics purposes
rather than physical simulation where accuracy is a major concern
there's some other methods that you might consider.
First, a bit of terminology. There's two main paradigms I'm aware of
which are used to simulate fluids: Eulerian and Lagrangian. In the
Eulerian approach you consider what's happening in a small fixed piece
of space and look at how that changes based on the pieces of space
around it. In the Lagrangian approach, on the other hand, you consider
what happens to pieces of the fluid as they move around. The approach
you describe is then a Lagrangian approach.
As you've noticed, stability can be a problem when you integrate out
these sorts of equations to determine how the fluid should move over
time. Changing the integration mechanism (such as to a Runge-Kutta
integrator) or decreasing the time-step can help this to a degree, but
there exist fluid simulation methods which are unconditionally stable.
That is, the simulation will never `blow up' no matter what you do, even
with very large time steps in the simulation.
Probably the simplest such technique is the "Stable Fluids" algorithm
invented by Jos Stam. It uses some concepts from both the Eulerian and
Lagrangian viewpoints, and is this referred to as a "semi-Lagrangian"
method. There's a very easy to ready paper on the algorithm that even
gives some code for how to implement it:
http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf
The algorithm itself should be also pretty easy to implement, and is
also useful in a educational sense as it uses some interesting
techniques (such as, iirc, explicitly projecting the flow onto a
divergence-free space to achieve an incompressible flow).
There are also some fully particle based fluid simulation methods which
are more similar to yours in their apploach, but I'm not aware of any
which are as simple to implement as the stable fluids algorithm. You
can check out some early examples:
http://graphics.ethz.ch/Downloads/Publications/Papers/2003/mue03b/p_Mue03b.pdf
http://www.cs.utah.edu/vissim/papers/ParticleFluid/ParticleFluidsLowRes.pdf
Or see a more recent paper on such a method at:
http://graphics.stanford.edu/~barta/research/adams-sig07.pdf
I suspect implementing one of these methods (or one of the newer
Eulerian / semi-Lagrangian ones) would be more work than you care to
spend, but I think stable fluids should be pretty reasonable to do
without a crazy amount of work.
Note, however, that all the stuff I mentioned is based on simulating the
Navier-Stokes equations. If you're only looking for something that
looks sortof-kinda like a fluid rather than behaving as a proper fluid
then there's probably all sorts of shortcuts and hacks you can do, and
in this case I agree with scott's suggestions so I'd start playing
around with using a better integrator and/or adding some sort of
artificial dampening.
Chambers wrote:
> 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...
>
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> 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)
...
> 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.
Here's a question: The derivative is based on the forces applied to the
current particle. The only way these forces change is when a particle
collides with other particles (or ceases to collide).
So, for each of the above, am I to run a full collision detection?
--
...Ben Chambers
www.pacificwebguy.com
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> Ah, grasshopper, you are still thinking in imperative terms. ;-)
BTW, that grasshopper thing... does that come from the old TV show "Kung
Fu"? I could have sworn I saw it in a movie, though. In fact, I asked
at a couple local video stores what movie it was, and nobody knew... :(
--
...Ben Chambers
www.pacificwebguy.com
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> 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?
AFAIK, yes, since they're defined as immutable.
In other words, Haskell doesn't have variables, only constants.
--
...Ben Chambers
www.pacificwebguy.com
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> 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.
What's that like on performance though? What would it be like if your state
block consisted of hundreds of 4D vectors, and the system was being updated
at 1000 Hz? Would the GC be noticeable?
> (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".)
Ah ok I see, so it's like you every variable is defined as "const" in the
C++ sense, you can't modify it.
> animate state = do
> dispaly state
> let new_state = integrate state
> animate new_state
I assume that the compiler is clever enough to not keep on the stack any
data that is not needed, ie old versions of state that are no longer needed.
If you attempted that method in C++ you would fill up the RAM very quickly
;-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Here's a question: The derivative is based on the forces applied to the
> current particle. The only way these forces change is when a particle
> collides with other particles (or ceases to collide).
>
> So, for each of the above, am I to run a full collision detection?
Yes, inside the "acceleration" function, which is called from each
"evaluate" you put your code for calculating the acceleration. This, in
your case, would presumably call the routine to check for collisions and
work out the forces for all the particles.
Note that you need to put all the positions and velocities for every
particle inside some state block and evaluate them in parallel using this
algorithm - you can't do one at a time it doesn't work like that.
You may say that this is 4x as much work per integration step, but it will
allow your time-step to be far more than 4x bigger and still be completely
stable.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> Note that you need to put all the positions and velocities for every
> particle inside some state block and evaluate them in parallel using
> this algorithm - you can't do one at a time it doesn't work like that.
That part, at least, seems obvious to me.
It certainly seems like I've got my work cut out for me :) I almost
feel that, if I'm going to implement something this complex, I may as
well do something useful with it!
--
...Ben Chambers
www.pacificwebguy.com
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> 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.
>
> What's that like on performance though? What would it be like if your
> state block consisted of hundreds of 4D vectors, and the system was
> being updated at 1000 Hz? Would the GC be noticeable?
Erm... without trying it, I couldn't say off the top of my head. Notice
though that since this is "the normal way" to do stuff in Haskell, the
compiler designers will be basing all their design decisions on the fact
that people are going to be writing code this way.
I would say that for *vast* numbers of vectors, you might need to move
over to using mutable arrays (uses less space and less GC time). But
without benchmarking it I couldn't tell you exactly how much of a
difference it makes. (I would also suspect that if you have "hundreds"
of vectors you probably need random access to them too...)
>> (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".)
>
> Ah ok I see, so it's like you every variable is defined as "const" in
> the C++ sense, you can't modify it.
Well, kinda, yeah.
It's like, a function can only have one body. Nobody would find that
surprising. Like in C, you can't write a function, and then later change
that function to something else. (While we're on the subject, in Haskell
a "function name" is just a normal variable who's value happens to be a
function...)
A basic property of Haskell code [and any "pure" functional language] is
that if you take any variable and replace all occurrances of that
variable with the expression that defines its value, the meaning of the
program is unchanged. [Although possibly efficiency is reduced.]
Notationally this is how you "execute" such a program.
>> animate state = do
>> dispaly state
>> let new_state = integrate state
>> animate new_state
>
> I assume that the compiler is clever enough to not keep on the stack any
> data that is not needed, ie old versions of state that are no longer
> needed. If you attempted that method in C++ you would fill up the RAM
> very quickly ;-)
Er, yes, in C++ (or most "normal" programming languages) the above loop
would be an *absurdly* bad idea. In Haskell, this is the normal way to
achieve looping, and hence the compiler goes to extreme lengths to
optimise it well.
(If "state" is small enough, it may even get permanently allocated into
CPU registers and updated in-place. It depends on exactly what it
contains... If it's two floats like in the examples, then probably yes.
If it's a set of 200 vectors, obviously no...)
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> Are you
>> basically saying that variables are write-once in Haskell?
>
> AFAIK, yes, since they're defined as immutable.
>
> In other words, Haskell doesn't have variables, only constants.
Well, depending on how you define "variable" and "constant".
Take this for example: double x = 2*x
Every time you call this function, "x" can have a different value. It is
a "variable". However, DURING EACH INVOCATION, "x" can only have one
value. It is a "constant".
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> Yes, inside the "acceleration" function, which is called from each
> "evaluate" you put your code for calculating the acceleration. This, in
> your case, would presumably call the routine to check for collisions and
> work out the forces for all the particles.
>
> Note that you need to put all the positions and velocities for every
> particle inside some state block and evaluate them in parallel using
> this algorithm - you can't do one at a time it doesn't work like that.
Question: Would it be worth using some kind of space partitioning to
enable you to process only "nearby" particles?
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Question: Would it be worth using some kind of space partitioning to
> enable you to process only "nearby" particles?
I think he said he was using something like that already.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> Question: Would it be worth using some kind of space partitioning to
>> enable you to process only "nearby" particles?
>
> I think he said he was using something like that already.
Ah, OK.
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> It's like, a function can only have one body. Nobody would find that
> surprising.
Unless you're using Erlang, I guess. :-)
"This is beyond dynamic binding. This is
make-up-your-mind-after-it's-already-running binding."
--
Darren New / San Diego, CA, USA (PST)
"That's pretty. Where's that?"
"It's the Age of Channelwood."
"We should go there on vacation some time."
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Chambers wrote:
> BTW, that grasshopper thing... does that come from the old TV show "Kung
> Fu"? I could have sworn I saw it in a movie, though. In fact, I asked
> at a couple local video stores what movie it was, and nobody knew... :(
I have no idea. My boss used to call me grasshopper. No idea why...
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
On Tue, 01 Apr 2008 21:04:23 +0100, Orchid XP v8 <voi### [at] dev null>
wrote:
>Chambers wrote:
>
>> BTW, that grasshopper thing... does that come from the old TV show "Kung
>> Fu"? I could have sworn I saw it in a movie, though. In fact, I asked
>> at a couple local video stores what movie it was, and nobody knew... :(
>
>I have no idea. My boss used to call me grasshopper. No idea why...
Walk the rice paper and leave no mark, Grasshopper.
--
Regards
Stephen
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
While we're on the subject... I've managed to get particle systems like
this to work before. (E.g., my chaos pendulum example, and various other
things.) What I haven't ever managed to do is figure out how to build
myself my own wave tank.
I know I've asked this before, but what's the algorithm for this?
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Right. Well maybe it would be helpful to fix a few typos...
> evaluate (v,p) t dt (dv,dp) =
> let
> v' = v + dt*dv
> p' = p + dt*dp
> in (v', a (t+dt) v' p')
Err... that's incorrect. The acceleration function ("a") needs to be a
parameter somewhere. Or it needs to be a top-level function. Let's make
it top-level, like in the C++ version:
acceleration t v p = ???
evaluate (v,p) t dt Nothing = (v, acceleration t)
evaluate (v,p) t dt (Just (dv,dp)) =
let
v' = v + dt*dv
p' = p + dt*dp
in (v', acceleration (t+dt) v' p')
[I made the derivatives optional too. Did I get it right?]
> 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)
...and if we're using a top-level acceleration function, you don't need
that "a" parameter there (which incidentally causes a name clash).
integrate t dt (v,p) =
let
a = evaluate (v,p) t dt Nothing
b = evaluate (v,p) t (dt/2) (Just a)
c = evaluate (v,p) t (dt/2) (Just b)
d = evaluate (v,p) t dt (Just 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)
I'll have to try this out later and see if it actually works...
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
You know what? Screw it. Let's start again...
data State = State {velocity, position :: !Double}
data Derivative = Derivative {dvelocity, dposition :: !Double}
evaluate :: State -> Double -> Double -> Maybe Derivative -> Derivative
evaluate s t dt Nothing = Derivative {dvelocity = acceleration s t,
dposition = velocity s}
evaluate s t dt (Just d) =
let
s' = State {velocity = velocity s + dvelocity d * dt, position =
position s + dposition d * dt}
in Derivative {dvelocity = acceleration s' (t+dt), dposition =
velocity s'}
integrate :: State -> Double -> Double -> State
integrate s t dt =
let
a = evaluate s t dt Nothing
b = evaluate s t (dt/2) (Just a)
c = evaluate s t (dt/2) (Just b)
d = evaluate s t dt (Just c)
dv = (dvelocity a + 2 * (dvelocity b + dvelocity c) + dvelocity d) / 6
dp = (dposition a + 2 * (dposition b + dposition c) + dposition d) / 6
in State {velocity = velocity s + dv * dt, position = position s +
dp * dt}
acceleration :: State -> Double -> Double
acceleration = error "not implemented yet"
Does that look correct?
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> Does that look correct?
Well, I just ran it on a simple harmonic oscilator and got a sine wave,
so that's always a good start...
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> I know I've asked this before, but what's the algorithm for this?
The wave equation.
http://en.wikipedia.org/wiki/Wave_equation
2nd time differential of "wave value" is proportional to the 2nd spatial
differential of "wave value"
2nd time differential is just acceleration (use it in an integrator like RK4
to get the new positions).
2nd spatial differential in 2D is just d2u/d2x + d2u/d2y
See attached POV code.
Post a reply to this message
Attachments:
Download 'us-ascii' (3 KB)
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> I know I've asked this before, but what's the algorithm for this?
>
> The wave equation.
>
> http://en.wikipedia.org/wiki/Wave_equation
>
> 2nd time differential of "wave value" is proportional to the 2nd spatial
> differential of "wave value"
>
> 2nd time differential is just acceleration (use it in an integrator like
> RK4 to get the new positions).
>
> 2nd spatial differential in 2D is just d2u/d2x + d2u/d2y
>
> See attached POV code.
Right. So for each pixel, measure the difference against the four
neighboring pixels, total that up, multiply by the wave propogation
constant, and that's the total acceleration for this pixel. Seems simple
enough.
What happens if the wave propogation medium propogates different
frequencies at different velocities?
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Right. So for each pixel, measure the difference against the four
> neighboring pixels, total that up, multiply by the wave propogation
> constant, and that's the total acceleration for this pixel. Seems simple
> enough.
>
> What happens if the wave propogation medium propogates different
> frequencies at different velocities?
The wave propagation constant is not constant, but a function of
frequency (dispersion relation).
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> What happens if the wave propogation medium propogates different
> frequencies at different velocities?
Well then your c is a function of frequency, it's also probably a function
of "wave value" and "wave steepness" too in real life. (higher and steeper
waves might propogate faster than lower waves). Good luck calculating the
frequency at a point though :-) I'm thinking you might be able to do
something clever with fourier though...
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> What happens if the wave propogation medium propogates different
>> frequencies at different velocities?
>
> Well then your c is a function of frequency, it's also probably a
> function of "wave value" and "wave steepness" too in real life. (higher
> and steeper waves might propogate faster than lower waves). Good luck
> calculating the frequency at a point though :-) I'm thinking you might
> be able to do something clever with fourier though...
Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a
simple all-pass filter in the feedback loop. OTOH, it's only trying to
simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could do
some sort of interesting convulation operation?
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> Well then your c is a function of frequency, it's also probably a
>> function of "wave value" and "wave steepness" too in real life. (higher
>> and steeper waves might propogate faster than lower waves). Good luck
>> calculating the frequency at a point though :-) I'm thinking you might
>> be able to do something clever with fourier though...
>
> Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a
> simple all-pass filter in the feedback loop. OTOH, it's only trying to
> simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could do
> some sort of interesting convulation operation?
I think if you are simulating a musical instrument, it would be better to
get the resonant parts working first before you start trying to simulate
different propagation speeds at different frequencies. I would imagine that
is a much smaller secondary effect.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>>> Well then your c is a function of frequency, it's also probably a
>>> function of "wave value" and "wave steepness" too in real life.
>>> (higher and steeper waves might propogate faster than lower waves).
>>> Good luck calculating the frequency at a point though :-) I'm
>>> thinking you might be able to do something clever with fourier though...
>>
>> Well, Reaktor's "Steam Pipe" instrument simplates dispersion using a
>> simple all-pass filter in the feedback loop. OTOH, it's only trying to
>> simulate a 1D wave. At exactly 1 point in space. Hmm, maybe you could
>> do some sort of interesting convulation operation?
>
> I think if you are simulating a musical instrument, it would be better
> to get the resonant parts working first before you start trying to
> simulate different propagation speeds at different frequencies. I would
> imagine that is a much smaller secondary effect.
Oh yeah, sure. I'm just going off at a tangent now. ;-)
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
|
|
