MGPBD: A Multigrid Accelerated Global XPBD Solver

TL;DR

MGPBD tackles the instability and stagnation of XPBD in high-resolution, high-stiffness deformable simulations by introducing a dual-space UA-AMG preconditioned CG solver. The approach leverages lazy setup to reuse prolongators, and a lightweight near-kernel construction to bootstrap the multigrid hierarchy, achieving robust convergence with sparser coarse grids. Empirical results across diverse scenes show substantial convergence gains, improved numerical stability, and near-linear per-iteration time scaling with mesh resolution, outperforming XPBD and primal-space MG methods, especially at extreme stiffness. The work enables efficient, stable, high-fidelity deformable simulations and provides a practical, GPU-accelerated framework for integrating AMG into XPBD pipelines.

Abstract

We introduce a novel Unsmoothed Aggregation (UA) Algebraic Multigrid (AMG) method combined with Preconditioned Conjugate Gradient (PCG) to overcome the limitations of Extended Position-Based Dynamics (XPBD) in high-resolution and high-stiffness simulations. While XPBD excels in simulating deformable objects due to its speed and simplicity, its nonlinear Gauss-Seidel (GS) solver often struggles with low-frequency errors, leading to instability and stalling issues, especially in high-resolution, high-stiffness simulations. Our multigrid approach addresses these issues efficiently by leveraging AMG. To reduce the computational overhead of traditional AMG, where prolongator construction can consume up to two-thirds of the runtime, we propose a lazy setup strategy that reuses prolongators across iterations based on matrix structure and physical significance. Furthermore, we introduce a simplified method for constructing near-kernel components by applying a few sweeps of iterative methods to the homogeneous equation, achieving convergence rates comparable to adaptive smoothed aggregation (adaptive-SA) at a lower computational cost. Experimental results demonstrate that our method significantly improves convergence rates and numerical stability, enabling efficient and stable high-resolution simulations of deformable objects.

Figures (16)

• Figure 1: Power density graph after 2 iterations of MGPBD and 300 iterations of XPBD. From left to right, the red subgraph shows dual residuals before the solving, while the other green ones are after solving. Two iterations of MGPBD eliminate the two spikes at LF and HF, while 300 iterations of XPBD solving still keep the spikes. This explains the superiority of MG in solving all-frequency errors over the iterative method. The data is generated by a cloth hanging case consisting of $1025\times1025$ particles at the 11th frame.
• Figure 2: Pipeline of the setup phase of UAAMG. Input: $A$ at the finest level; Output: $P$s at all levels. It includes 3 steps: 1) Filter: filter out the weak connection with a threshold $\theta_s$ and get the strength of connection matrix $\mathbf{S}$. Here $\mathbf{A}$ is symmetric, and off-diags $|a_{ij}|$ are used as strength. 2) Aggregate: aggregate the nodes into separate groups using $\mathbf{S}$. Here is an example of cloth with $4\times4$ particles. Its dual graph (constraint node graph) is aggregated into four groups marked in yellow, green, blue, and red, respectively. 3) Inject: inject a near-kernel component $\mathbf{B}$ into the matrix $\mathbf{P}$. Each column of $\mathbf{P}$ corresponds to one aggregate, and its nonzero value corresponds to the weight within that aggregate. For every level, use QR decomposition to form the next level, where $\mathbf{R}$ of QR decomposition serves as $\mathbf{B}$ at the next level. With $\mathbf{P_l}$ at this level, the next level $A$ will be formed with Eq. \ref{['eq:Galerkin']} and used as the input to construct the next level $\mathbf{P_{l+1}}$. This procedure is repeated for all levels until the matrix size is under $400$. We finally get a set of $\mathbf{P}$s as the output of the AMG setup phase.
• Figure 3: Sparsity patterns of SA and UA. Number of nonzeros ($nnz$): 58M/587K/2K (SA) vs. 58M/136K/4K (UA). Sparsity: 8e-5/0.8%/74.2% (SA) vs. 8e-5/0.2%/4.6% (UA) VS. Data is collected from the first frame of bunny squash case with 850K tetrahedrons. The key difference of $nnz$ in level 1 makes UA faster than SA.
• Figure 4: Residual curves for different setup intervals. There is no significant difference in convergence when updating the prolongators aross $1/10/20/50/100$ frames. Data are collected from the $99$-th frame of the $270$K-resolution bunny squash scenario (Fig. \ref{['fig:bunny']}).
• Figure 5: Comparison of linear solvers for different AMG methods and different smoothers in convergence rate. a) Comparison of Different setup phases, the smoother uses weighted-Jacobi method. It shows that the nullspace method (ours) is superior. b) Different smoothers, the setup phase uses the nullspace method. It shows that the weighted-Jacobi method is the best smoother for softbody cases. The data is collected from the bunny squashing case with 12K tetrahedrons (Fig.\ref{['fig:bunny']}).