Conjugate Gradient based MLCP solver source
Posted: Fri Jan 05, 2018 5:56 am
Hello!
Several years ago we already shared this kind of solver, but last time our demonstration framework was not very compelling (for example, it was lacking collision detection and hence only very basic scenes with constraint force limits could be tested) - so I decided to revisit it now as I've got some free time, and quickly wrote a testbed that compares PGS with CG-based solver. PGS comes in two flavors - classic and shuffled - the last one allows to use various schemes to permute constraint rows each iteration, for example random shuffling, even-odd shuffling and so on.
Our solver is based on the MPRGP (Modified Proportioning with Reduced Gradient Projections) by Z. Dostál; we modified the original MPRGP to support box limits (low/high, instead of just lower limit), and incorporated preconditioner that is fairly cheap - double (left+right) Jacobi preconditioner. Solver also benefits from the J*M^-1*J^T decomposition that significantly speeds up matrix-vector multiplications.
Additionally, I added a hack to make joints more stiff, which is similar to shock propagation in a way - constraint masses are scaled locally for each joint, allowing to setup open loop systems where scales are set so that impulse propagates better from root to nodes. It is still a hack with artifacts, for example stack is still becoming unnaturally resistant to tumbling with large mass scaling parameters - but quite tolerable for mass scaling of, say, 2. The aim of that testbed was to see what solver matches which situation better, and using hacks seems alright for real-time game simulation.
The old speculations were somewhat confirmed - for fast real-time systems where quality/accuracy is not the goal - PGS still works better, as it guarantees lower error each iteration, and can behave reasonably well with even very low iteration count, another plus is that each iteration is very lightweight. Local mass scaling can help pretty noticeably even with relatively low mass scaling parameters.
In the systems where accuracy/joint stiffness is required, CG-based solver fits much better. And although each iteration is heavier than one iteration of PGS, and the method requires spectral radius estimation (additional overhead), ultimately it converges faster with larger iteration counts (>50), and for stiff systems provides better behavior within the same frametime. The method is also naturally parallel, contrary to PGS which requires coloring or splitting, so it will benefit from GPU acceleration (which we also experimented with several years ago and got significant speedup even on GTX260-era HW). Additional note, however, is that the solver calculates allowed stepsize each iteration, so it should definitely benefit from the island detection (rather than solving all the islands in a single system).
The project also incorporates FEM-based deformables (corotational FEM constraints that fit into the very same system with the regular rigid body joints, so coupling is natural), which comes in handy when testing stiff systems.
For anybody interested, source code is available on github - but the render uses Vulkan API, so LunarG Vulkan SDK installation is required (details in the readme.md file). Our old paper, describing both our approach to formulating FEM-based constraints (which actually very similar to the one by M. Servin et al., from "Interactive Simulation of Elastic Deformable Materials") and the CG-based MLCP solver, is also available on github.
Several years ago we already shared this kind of solver, but last time our demonstration framework was not very compelling (for example, it was lacking collision detection and hence only very basic scenes with constraint force limits could be tested) - so I decided to revisit it now as I've got some free time, and quickly wrote a testbed that compares PGS with CG-based solver. PGS comes in two flavors - classic and shuffled - the last one allows to use various schemes to permute constraint rows each iteration, for example random shuffling, even-odd shuffling and so on.
Our solver is based on the MPRGP (Modified Proportioning with Reduced Gradient Projections) by Z. Dostál; we modified the original MPRGP to support box limits (low/high, instead of just lower limit), and incorporated preconditioner that is fairly cheap - double (left+right) Jacobi preconditioner. Solver also benefits from the J*M^-1*J^T decomposition that significantly speeds up matrix-vector multiplications.
Additionally, I added a hack to make joints more stiff, which is similar to shock propagation in a way - constraint masses are scaled locally for each joint, allowing to setup open loop systems where scales are set so that impulse propagates better from root to nodes. It is still a hack with artifacts, for example stack is still becoming unnaturally resistant to tumbling with large mass scaling parameters - but quite tolerable for mass scaling of, say, 2. The aim of that testbed was to see what solver matches which situation better, and using hacks seems alright for real-time game simulation.
The old speculations were somewhat confirmed - for fast real-time systems where quality/accuracy is not the goal - PGS still works better, as it guarantees lower error each iteration, and can behave reasonably well with even very low iteration count, another plus is that each iteration is very lightweight. Local mass scaling can help pretty noticeably even with relatively low mass scaling parameters.
In the systems where accuracy/joint stiffness is required, CG-based solver fits much better. And although each iteration is heavier than one iteration of PGS, and the method requires spectral radius estimation (additional overhead), ultimately it converges faster with larger iteration counts (>50), and for stiff systems provides better behavior within the same frametime. The method is also naturally parallel, contrary to PGS which requires coloring or splitting, so it will benefit from GPU acceleration (which we also experimented with several years ago and got significant speedup even on GTX260-era HW). Additional note, however, is that the solver calculates allowed stepsize each iteration, so it should definitely benefit from the island detection (rather than solving all the islands in a single system).
The project also incorporates FEM-based deformables (corotational FEM constraints that fit into the very same system with the regular rigid body joints, so coupling is natural), which comes in handy when testing stiff systems.
For anybody interested, source code is available on github - but the render uses Vulkan API, so LunarG Vulkan SDK installation is required (details in the readme.md file). Our old paper, describing both our approach to formulating FEM-based constraints (which actually very similar to the one by M. Servin et al., from "Interactive Simulation of Elastic Deformable Materials") and the CG-based MLCP solver, is also available on github.