The best engine to generate close packings of spheres?

[juliusc]

Hi everyone,

I'm trying to generate random close packings of spheres inside a box for a physics project I'm working on. Accuracy and stability are very important to me but speed is crucial, because I need to use 15000 spheres, maybe more, and I have to generate thousands of configurations.

The idea is to put the spheres inside an enormous box and then to reduce gradually its volume until eventually I could obtain a close packing. The collisions must have friction and the spheres must be hard spheres.

Could anyone recommend me which would be the best engine for this project?

Thanks in advance,
Julius

[Erwin Coumans]

I'm trying to generate random close packings of spheres inside a box for a physics project I'm working on. Accuracy and stability are very important to me but speed is crucial, because I need to use 15000 spheres, maybe more, and I have to generate thousands of configurations.

What are random close packings of spheres exactly? Is that one compound/composite object containing 15000 spheres? Or just 15000 independent moving spheres?

The idea is to put the spheres inside an enormous box and then to reduce gradually its volume until eventually I could obtain a close packing. The collisions must have friction and the spheres must be hard spheres.

How do you want to determine when to stop 'shrinking' the box? Based on constraint solver feedback?

Could anyone recommend me which would be the best engine for this project?

If this can lead to some demo for Bullet, I'm happy to collaborate on this,
Thanks,
Erwin

[juliusc]

Hi Eric,

Thanks for your answer

What are random close packings of spheres exactly? Is that one compound/composite object containing 15000 spheres? Or just 15000 independent moving spheres?

They are independent spheres, moving and colliding with each other. The simplest picture of a random close packing would be the configuration that you obtain when you fill a box with lots of golf balls.

How do you want to determine when to stop 'shrinking' the box? Based on constraint solver feedback?

The idea is to stop when the fraction of the total volume of spheres with respect to the volume of the box reaches some predefined value, for example 40%. When this value is reached, the spheres stop their movement and their positions and radii are saved to a file.

Since the the spheres can't penetrate between them, their total volume is just the sum of their individual volumes. At each time step the length of the box is reduced in a proportion that depends on how close we are to the desired fraction.

If this can lead to some demo for Bullet, I'm happy to collaborate on this,
Thanks,
Erwin

Thank you very much for your offer Erwin Please let me know if you need more information

[tasora]

Hi,
I developed the chrono::engine physics engine, it is available at

http://www.deltaknowledge.com/chronoengine

Here we used chrono::engine to simulate dense granular flow of spheres up to
150'000 spheres, for simulating a PBR nuclear reactor. I suppose that poses
a problem that it is similar to yours.
I also put the 100'000 spheres video on youtube, but I don't remember the url.

It is free for academic & research purposes. Feel free to give it a try, it is not 100%
completed but people is starting using it in some projects.

Interesting note: chrono::engine uses Bullet as a collision engine (with minor
modifications). Thank you Erwin for this great CD library!

Now we are working together with two US universities to port the solver on GPU,
because we want to efficiently solve problems with >500'000 spheres.

regards

Alessandro Tasora

[juliusc]

Hi Alessandro,

I've been checking the chrono::engine site and I have to tell you: your engine is a pretty impressive work of programming!!

I think it would be the perfect engine for my project if it were not for one important thing: it's not available for linux. Yep, I need to run my simulations in linux clusters. I hope you release a linux version very soon. I have the feeling I'm not the only one waiting for it

I think I'll try Bullet by now, but I've been checking ODE and AGEIA too. Do you recommend them? Are they worth the try? Or are they too game oriented?

Thanks,
Julius

[ngbinh]

Hi,
I developed the chrono::engine physics engine, it is available at

http://www.deltaknowledge.com/chronoengine

Here we used chrono::engine to simulate dense granular flow of spheres up to
150'000 spheres, for simulating a PBR nuclear reactor. I suppose that poses
a problem that it is similar to yours.
I also put the 100'000 spheres video on youtube, but I don't remember the url.

It is free for academic & research purposes. Feel free to give it a try, it is not 100%
completed but people is starting using it in some projects.

Interesting note: chrono::engine uses Bullet as a collision engine (with minor
modifications). Thank you Erwin for this great CD library!

Now we are working together with two US universities to port the solver on GPU,
because we want to efficiently solve problems with >500'000 spheres.

regards

Alessandro Tasora

Hi Alessandro!

Can you elaborate what kind of changes to Bullet you've made? Also, I know that the Chrono engine is using a contact model/solver developed by Mihai Anitescu. Is it your only solver or you have something else?

Thanks

[tasora]

Hallo,

Can you elaborate what kind of changes to Bullet you've made?

Just few changes.
First of all, I use only the btCollisionWorld (I do not need to use also the Bullet
dynamics, I have my own solver)
Then, I modified the structure with the contact info (for the persistent
manifold) so that it contains also the cached multipliers for warm starting the
solver (six float values: 3 for tangential and normal impulses in speed LCP and 3 for the
position post-stabilization LCP). Note that such structure already has a pointer for user
data, but I prefer to embed it into the structure rather than allocating / deallocating those
6 floats, for performance reasons.

Also, I know that the Chrono engine is using a contact model/solver developed by Mihai Anitescu. Is it your only solver or you have something else?

I am developing a new solver together with Mihai Anitescu. The complementarity solver
that I implemented in chrono::engine is inspired to a fixed-point scheme originally
introduced by Mangasarian et al. some years ago, but I modified it so that it does not
simply handle LCPs, but also CCP (Cone Complementarity Problem) and NCP, more in
general, so it fits well into the multibody environment. This solver is documented
in http://www.optimization-online.org/DB_H ... /1687.html and a
more detailed paper is under reviewing process.
Also, the time stepper (and the overall model) is similar to the original one from Anitescu,
but we modified it a bit to allow a more reliable constraint stabilization. It requires a CCP
per time step. (however in chrono::engine you can optionally use another time-stepper
which performs position-stabilization of constraints as a second-pass CCP).
We hope to prepare a detailed paper with algorithmic details as soon as possible.

best regards

A.Tasora


Alessandro Tasora

[Erin Catto]

Alessandro,

I scanned over this paper. Some of the math is a bit deep and I didn't spend the time to digest it all.

How would you compare this method to Bullet's sequential impulse solver?

Does the CCP provide you with isotropic friction? Is that the main advantage?

The paper states that the solver is like Gauss-Seidel with overrelaxation. In my experience over-relaxation can be energetic when the iterations are limited. Does your solver have the same issue?

You've presented examples with large numbers of bodies. How does the solver behave with smaller numbers of bodies in stacking scenarios?

Can this solver also handle bilateral constraints?

[tasora]

Hallo Erin,

How would you compare this method to Bullet's sequential impulse solver?

As already said in other threads, a PGS-like solver is not much different
from a sequential impulse solver, under some circumstances (i.e. immediate
update of variables as in GS or SOR schemes, etc.). Anyway, our iterative
solver does not need the knowledge of the TOI, because it is the CCP
solution which takes care of this, and the order of variable solution can
be free (although symmetric or unidirectional proved to be the best for
cache-coherency). Also, we can tweak some additional factors and we can
play with the convergence assumptions.
Also, I don't have time to look how friction cone projection is done in Bullet;
ours is not strictly horizontal to contact plane (which looks counterintuitive
but it guarantees convergence even for high friction values).

Does the CCP provide you with isotropic friction? Is that the main advantage?

Yes, friction is isotropic and does not suffer from the 'polygonal' artifact
as in LCP approaches. Also, it is easy to turn the friction cone into anisothropic
'elliptical' friction as in tire-ground friction models - but I had no time to play
with this.
This is not the only advantage of the solution method... Others are the robustness,
the easy implementation, the fact that it can solve LCPs and linear systems as
sub-cases, etc. etc.

In my experience over-relaxation can be energetic when the iterations are limited. Does your solver have the same issue?

Mhhh.. If I understand what you mean, you experienced that high relaxation factor
omega can lead to 'jumpy' behaviour, right? This is because your omega exceeded
the stability interval (which, in these multibody problems is not the usual 0...2 interval
as in SOR theory for diagonaldominant systems, etc.)
There are some tricks to adjust the omega relaxation factor on the fly, but I don't
have time to discuss them today. ..

You've presented examples with large numbers of bodies. How does the solver behave with smaller numbers of bodies in stacking scenarios?

Yes I can simulate stacking scenarios, but I am not yet 100% satisfied with the results
(I mean, the solver behaves well, but to get stable stacking with few iterations
I need warm-starting as in Bullet solver: however warm starting uses the Bullet
contact persistent manifold, and I still must exploit contact persistency at its
best)

Can this solver also handle bilateral constraints?

Yes, during these years Chrono::Engine has been used to succesfully design
mechanical divices (robots, cars), and recently evolved in the current
real-time library. So bilateral constraints are supported, as well as rheonomic
constraints (either exact -without the use of PIDs- and force-based)

bye

Alessandro Tasora

[Dirk Gregorius]

Hi Alessandro,

I scanned your paper today and I have two questions:

A) As Erin pointed out using a SOR parameter other than one might add energy to the system. You state that the value may not be in the interval [0, 2]. In the past I used values in the range from 1.1 - 1.3. But these values were more a guessing. To my knowledge the SOR parameter depends on the eigenvalues of the system. For a multibody system with several hundred bodies it is not feasible in a real-time scenario to compute the eigenvalues . So what heuristics do you use to get an estimate for omega? You also mention that you can adjust omega on the fly during the iterations. I have seen these methods and some are founded on Aitken or Steffenson acceleration. I wonder if these methods can blindly be used for LCP or CCP, since the series are not strictly converging when the active set changes. Maybe you can elaborate a little bit more on this? What SOR value do you typically use? How many iterations did you save in comparison to Gauss-Seidel? How many iterations do you typically use in general? What do you suggest to tweak omega on the fly? What is the convergence number of your solver - linear, super-linear or quadratic?

B) In the sketch of your solver you suggest a forward and backward sweep over the constraints. This is know as SSOR to my knowledge and has known worse convergence (e.g. s. Golub). What improvements did you discover or is this only feasible for specific scenarios?


BTW:
I just noticed in an early post in this thread, that you also warmstart the post-stabilization. I had huge problems with this and disabled it. Did you have also any problems?

Cheers,
-Dirk

[tasora]

Hallo Dirk,

Some informations about the choiche of the omega overrelaxation factor.
Currently I don't have an optimal algorithm for the best choiche of
that factor - in most cases I prefer to keep omega=1 (although this is not
guaranteed to produce a contractive mapping under all circumstances, it
proved to be a 'reliable guess' for most situations which I tested, on average).

Of course we also developed a basic algorithm to adjust omega (as cited
in the article) but it is somewhat experimental: we start with an optimistic
omega value (say, omega= 2), then if a projective iteration produces an
increase in the residual, the iteration is repeated with a smaller omega.
Anyway: this schema 'wastes' some cpu resources because of the
repetition of the iteration cycle, though it seems that few adjustements are
enough (and the next time step can start with the last 'good' omega, with
some slow increase to capture the changing nature of the scene).
In pratice, I never use this scheme because I don't like to 'undo' an
iteration cycle, either because it wastes cpu time, either because I
need some additional RAM for undoing the multipliers. So I was thinking of
reducing the omega 'on the fly', not on a iteration-cycle basis, but on
a multiplier-by-multiplier basis (assuming that residual is computed
as infinite-norm, an increase in the residual can be monitored even
before terminating all the sweep on the multipliers); this would
lead to a fast 'local' undo, but I have not yet looked in detail
what happens to the convergence theorems etc.
I think it's worth while exploring this overrelaxation topic, but
currently I have other priorities (parallel execution on GPU).

Also, as you point out, the convergence can be slow in systems with
strange eigenvectors, but I don't know how easy it would be to accelerate
the convergence with Aitken methods etc.
Formerly, I was rather thinking of alternating this fixed-point method with
some Krylov steps (ex. GS, etc.) as suggested in Kocvara.. in case you
have simple bilateral contraints, the GS phase would converge at least
as fast as for linear systems.. But I still haven't found a proper way
to make a decent Krylov iteration without needing to waste RAM or
introducing heavy computations (nor I am too optimistic on the
robustness of Krylov steps with overconstrained systems..)

By the way, about overconstraining: the iterative approach
can give a 'noisy' solution to problems with too many constraints.
I mean, think to this example: a door with two hinges... one
vertical reaction is unneeded. For this example, the solver still
gets a solution (in fact the residual will be really low or zero), but
in some timestep the reaction is into one hinge, then in other
timesteps is in the other hinge, or in mixed cases, etc., with a 'random'
pattern, of course. Not wrong solutions, I mean, but I bet that
the user would prefer to see that the solution is always the same set,
i.e. it would be more intuitive if the the multipliers were
the 'least square result', so if one plots the reactions, he
does not get a 'noisy' zigzag pattern.
I still don't see an easy way to fix this problem (operating
the solver on the normal equations will give numerical and
implementation difficulties) and I don't know if playing with
constraint-force-mixing terms could help a bit.. Have you
an idea?

About warm-starting: I had no major problems with warmstarting
the solver for position stabilization... The only problems wiht the
warmstarting, on my side, come from the fact that the algorithm
for the manifold persistency can be a bit unpredictable for ill cases
such as cube on cube, etc.

As for SSOR vs SOR: in fact I noticed few differences in convergence
speed, but for some 'lucky' reordering of constraints , the SOR is very
good in one direction, and very bad in the other, so SSOR is just a
way to be 'on the safe side' in case one meets these cases As for
a random ordering, maybe that it can affect negatively RAM cache
performance, so I never took it into consideration.

Well, I hope that this 'solver' thread can open some debate - I think
that this is an interesting research topic, which still need much
improvements and ideas.

regards,

Alessandro Tasora

[Dirk Gregorius]

Hi Alessandro,

sorry for the late reply.

1) Over-relaxation
I think Erin brought it to the point. An omega greater one might result in a gain of energy if you stop iterating before convergence. Since in games we usually break after a fixed number of iterations most solvers use the Gauss-Seidel approach. For an academic problem this might be a good solution if you can wait for convergence.

2) Convergence acceleration is an intersting topic. Erin came up with an idea in the thread below which is similar Aitkens acceleration. If I start looking into this problem again I will most probably start there: http://www.bulletphysics.com/Bullet/php ... f=4&t=1325. There are some postings of Erin about issues with a Projected Conjugate Gradient implementation (I assume that CG is a Krylov space method out of the back of my head).

3) The noisy problem you describe is also mentioned in some paper of Brian Mirtich. In constructional engineering we used "Elasticity-Theory" (sorry I don't know the english word) to solve these problems. I have no clue if this can be used in the context of constrained dynamics. My favored example is that you might have e.g. a spider crawling over a pillow and because of the wrong contact forces you get wrong deformations in the pillow. The GS usually seems to spread the contact forces equally. So in your example I would have expected that the hinges share the vertical load instead of this zigzag pattern.

4) SSOR
Do you have an example for such a case?


I agree, maybe we should make a sticky solver thread where we collect ideas for the solver.



Cheers,
-Dirk