Erwin Coumans wrote:Iterative solvers are very popular nowadays. However most methods don't copy well with large mass ratios. How do people deal with this ?
Well, there is a bunch of iterative methods, I guess you refer to matrix splitting methods? I will consider only these in the following discussion.
If you rephrase the iterative method as a matrix equation, for GS this is something like
x = T x + c
with
T = (D+L)^{-1} U
c = (D+L)^{-1} b
In order for GS to converge at all the spectral radius of T must be less than 1. General speaking the smaller spectral radius the better it is (or should I say the faster the convergence).
Large mass ratios seem to affect the magnitude of the eigenvalues of the system matrix. Thus the spectral radius of the A-matrix is affected.
So the trick is to get a small spectral radius! However, this is not easy.
General speaking preconditioning and diagonal scaling (diagonal scaling can be seen as a left-right preconditioning) can be used to ``change'' the spectrum of the eigenvalues. General speaking, the new preconditioned system looks like
A' = (C^{-1} A C^{-T}) C^{ T}
b' = C^{-1} b
I have tried using C = sqrt( 1.0 / diag(A) ) (this is diagonal scaling). If one compares a plot of the residual error of the new system with the old system then one see a small displacement of the curves. Diagonal scaling seems to lower the convergence rate by a constant but it does not change the asymptotic behavior.
To really get this preconditioning to work one needs to find a C-matrix which is a good approximation to A^{-1}) If this can be achieved one would have A' \approx I. Using diagonal scaling C is of course not a really good approximation to the inverse of A. Diagonal scaling seems to work best for diagonal dominant systems. In my experience multibody dynamics problems does not appear to be that. In genral matrix splitting methods seem to work best on diagonal dominant systems. Just consider the contrived case of using Jacobi on a diagonal system. Here T will be the identity matrix, and you will hit the solution in one iteration.
Shock-propagation can also be understood as a two stage blocked iterative scheme with a kind of preconditioner. The outer stage consist of deviding your problem into stack-layers. Each stack layer corresponds to a single block. Then one iterates over the blocks in a bottom-up fastion (this is a strange mixture of Jacobi/Gauss-Seidel like iteration. The diagonal blocking is similar to a Jacobi scheme, however layers must be processed in a certain order, which makes the processing of blocks sequential in nature, i.e. similar to Gauss-Seidel). For each block the mass ratios are changed by fixating some of the bodies, before solving the block.
Relaxation or damping also seem to work rather well. That is use
A' = A + epsilon I
For non-negative epsilon. Now as n-> infty we have epsilon -> zero. This is termed constaint force mixing in the ODE community. To really see what it means consider the equivalent QP. That is minimize
f(lambda) = lambda^T b + lambda^T A lambda + epsilon lambda^T lambda
s.t
lambda >= 0
Using KKT-conditions leads to an LCP with the A' matrix above. Looking closely at f one see that the magnitude of lambda is minimized by the last term. Also note lambda always will be dampen by this scheme. This leads to weaker constraint forces..... This can be used to model soft constraints but if you really want hard-constraints then you have lost!!!
However in my experiece relaxation seem to really change the slope of the convergenece rate. That is to say for the convergence rate p given by
lim_{n \rightarrow \infty} \frac{|e^n|}{|e^{n-1}|^p} = c
where c is called the convergence constant. Larger epsilon usually in my experience means lower c. Try looking at log-plots of the residual error and the different slopes becomes apparent. For GS type of solvers p=1 (this is really bad:-()
