DiffGIPC: Differentiable Mixed-Coordinate Contact Dynamics on GPU

Kemeng Huang; Bowen Ren; Xinyu Lu; Minchen Li; Taku komura

cs.GROct 13, 2026

DiffGIPC: Differentiable Mixed-Coordinate Contact Dynamics on GPU

Authors: Kemeng Huang, Bowen Ren, Xinyu Lu, Minchen Li, Taku komura

Organizations: The University of Hong Kong, Hong Kong SAR, China · TransGP, Hong Kong SAR, China · Carnegie Mellon University, USA

Abstract

Contact-rich inverse problems require efficient gradients and reliable enforcement of non-penetration. Existing differentiable methods provide fast simulation and rigid-deformable coupling, but contact treatments based on penalties or velocity-level constraints can permit penetration or positional drift. Incremental Potential Contact (IPC) maintains non-penetration through geometric barriers and {continuous collision detection (CCD)}. However, repeated forward and backward evaluations remain costly in full-space finite-element (FEM) implementations, particularly when nearly rigid components retain many vertex {degrees of freedom (DoFs)}. We present DiffGIPC, a GPU-accelerated differentiable IPC framework for optimizing interacting deformable and nearly rigid bodies. Our affine-body derivatives enable {rest-shape} design, initial-state optimization, and articulated control while retaining only twelve motion DoFs per affine body. These gradients account for how changes in geometry, placement, and actuation affect other objects through collision and friction. We preserve these interactions in one coupled FEM-affine adjoint system per frame, reducing the backward-system size without discarding the dependencies that drive the inverse problem. GPU-optimized backward computation further reduces the cost of repeated gradient evaluations. A reusable PyTorch interface connects supported physical parameters to user-defined parameter maps and losses across design, identification, motion, placement, and articulated control. Experiments demonstrate peak speedups over DiffIPC of approximately 150.0x with full-space FEM and 4133.3x with affine bodies for combined forward simulation and backward differentiation.

Figures & tables

Panel of Figure 1 — Figure 1. Optimizing the placement of seven affine bodies. Keeping their shapes fixed, DiffGIPC adjusts only the horizontal initial positions to form a stable stack. The panels show the unoptimized configuration and the simulated results after 10, 50, and 78 iterations. Four panels show the initial stone placement and optimization updates 10, 50, and 78 at selected rendered states.

Figure 2. DiffGIPC’s backward pipeline. Each reverse frame assembles one coupled FEM-ABD adjoint system, retaining collision, friction, and joint interactions with twelve DoFs per affine body. The GPU solve and local accumulation propagate gradients to earlier states and physical parameters for optimization through PyTorch. Saved FEM and affine states feed local derivative evaluation, matrix assembly, one coupled adjoint solve, and gradient accumulation. The PyTorch loss supplies gradients, and the parameter map returns physical-parameter gradients to the optimizer.

Figure 3. Reverse-time accumulation in the frame-wise adjoint. Dashed arrows indicate forward state dependencies. Blue arrows carry contributions from the solved

\boldsymbol{\lambda}_{t}

to the earlier right-hand sides

\mathbf{b}_{t-1}

and

\mathbf{b}_{t-2}

. For the dynamic predictor shown, inertia and friction contribute to

\mathbf{b}_{t-1}

, while inertia contributes

-\mathbf{M}\boldsymbol{\lambda}_{t}

\mathbf{b}_{t-2}

. Each right-hand side also includes its direct loss gradient. Three time levels are stacked vertically in a single-column figure. Saved states occupy the left lane, and the coupled solve and earlier adjoint right-hand sides occupy the right lane. Dashed arrows point toward newer states. Blue arrows carry the current adjoint's contributions to the two earlier frames.

Panel of Figure 4 — Figure 4. Garment Hanger shape design. Two shared Bézier control scales modify both shoulders of the affine rest surface under frictional cloth contact. The initial and optimized configurations are shown at frame 33. The optimized hanger retains both garment straps. Initial and optimized garment hangers shown side by side with their supported cloth at frame 33.

Figure 5. PCG computation and execution control. (a) Five stages of one iteration (Algorithm 2 ). (b) A host-controlled loop retrieves a convergence value after each GPU iteration. (c) A conditional CUDA Graph keeps iteration and stopping on the GPU. Frame restoration, Hessian assembly, and gradient accumulation remain outside the graph. Both execution modes use the same PCG recurrence and stopping criterion. The top row shows matrix multiplication, solution and residual update, preconditioning, reduction and search-direction update, and the stopping test. Two CPU and GPU lane diagrams compare a host-controlled loop with a device-controlled conditional graph. Only the host-controlled loop returns a convergence scalar to the CPU at every iteration.

Panel of Figure 6 — Figure 6. Collision-mediated initial-velocity optimization. The first bunny’s initial velocity is adjusted so that the second bunny’s horizontal centroid reaches its target. The panels compare frame 40 before optimization and after 1, 3, and 10 iterations. Four panels show the initial placement and optimization updates 1, 3, and 10 at simulation frame 40.

Panel of Figure 7 — Figure 7. Trajectory matching through contact between two deformable bodies. The first body’s initial velocity is optimized to reproduce the target motion. Opaque surfaces show the simulation and translucent surfaces the target at frame 100, before optimization and after 1, 3, and 10 iterations. Trajectory matching for two contacting deformable bodies. The shared initial velocity of the first body is optimized to reproduce a prescribed motion. Opaque surfaces show the simulated bodies and translucent surfaces show the target configuration. The panels show frame 100.

Example	System	Forward						Backward				Forward + Backward
Example	System	Linear	CCD	Hessian	Search	Other	Total	Hessian	Linear	Other	Total	Forward + Backward
Bunny collision (Fig. 6 )	DiffIPC	4.5	10	4.6	1.1	1.9	22	1.4	1.1	0.44	2.9	25
	DiffGIPC	0.65	0.044	0.037	0.017	0.002	0.75	0.0069	0.22	0.0057	0.24	0.98
	Speedup	$6.9\times$	$232.0\times$	$123.9\times$	$67.7\times$	$963.8\times$	$30.0\times$	$204.7\times$	$4.9\times$	$76.3\times$	$12.5\times$	$25.8\times$
Tentacle matching (Fig. 7 )	DiffIPC	140	14	120	66	11	350	7.4	11	2.8	21	370
	DiffGIPC	8.1	2.4	2.2	1.4	0.011	14	0.04	1	0.024	1.1	15

Table 1. Mean time per forward and backward evaluation (seconds) on the shared FEM benchmarks. Forward + Backward gives the sum of the two per-evaluation means.

Panel of Figure 8 — Figure 8. Rest-shape optimization of a loaded FEM support. Under fixed constraints and downward loading, the objective minimizes squared first-Piola stress. The panels show the deformed support at frame 20 before optimization and after 20, 50, and 200 iterations. Four front-view stress renderings compare the initial support and the designs after 20, 50, and 200 updates.

Panel of Figure 9 — Figure 9. Fitting the material response of a constrained cube to marker positions. Each column pairs side (top) and oblique (bottom) views at frame 19, before optimization and after 2, 5, and 20 iterations. Opaque surfaces denote the simulation, translucent surfaces the target, and arrows the Dirichlet constraints. A two-row, four-column image sequence for material identification. Side views form the top row and oblique views the bottom row. Each column pairs the two views of the same optimization iteration, ordered 0, 2, 5, and 20 from left to right. Solid surfaces show the simulated deformation, translucent surfaces show the target, and arrows show the constraint directions.

Panel of Figure 10 — Figure 10. Friction optimization on a curved slide. A bunny-slide friction coefficient is adjusted so that the bunny reaches its target position at frame 100. The panels show the unoptimized simulation and results after 1, 3, and 5 iterations. A deformable bunny slides along a curved surface. The bunny-slide friction coefficient is optimized so that its position at the specified time matches a target position. The panels show final frame 100 for the initial coefficient and after one, three, and five optimization updates.

Figure 11. Optimization histories for the shared FEM benchmarks. Blue and red denote DiffGIPC and DiffIPC, respectively. The horizontal axes count accepted parameter updates from the initial parameters (iteration 0). Loss scales are task-specific, and Loaded Support is remeshed between updates. Five side-by-side plots show loss against optimization iteration for bunny collision, tentacle matching, loaded-support shape design, material fitting, and friction optimization. Blue curves denote DiffGIPC, and red curves denote DiffIPC.

Panel of Figure 12 — Figure 12. Cloth rest-scale fitting on a fixed mannequin. The blue panels show the initial drape and results after 1, 3, and 5 updates. The orange panel shows the target. Five garment-focused views with one common crop defined from the target. The four blue panels show the initial, intermediate, and optimized drapes. The orange panel shows the target. Sleeves and hems remain visible so changes in garment size can be compared at a common scale. Orange identifies the target, not stress.

Panel of Figure 13 — Figure 13. Grasp-angle control. Initial (left) and optimized (right) grasps at frame 60. Joint-target optimization moves the grasp toward the $0.45$ -rad target. Coral guides follow the fingers. Dashed blue guides indicate visual reference directions. Initial grasp beside optimized grasp, with finger and reference guides.

Panel of Figure 14 — Figure 14. Articulated cup catching through plate-trajectory optimization. With the arm held at rest (left), the cup falls from the plate. After optimization (right), the plate catches and retains it. The middle column shows the initial IK trajectory, which also fails to retain the cup. Four rows show the same saved frames (94, 111, 120, and 239). Rest-arm controls are on the left, the failing initial IK trajectory is in the middle, and optimized controls are on the right.

Panel of Figure 15 — Figure 15. Dual-arm towel grasp by optimizing one vertical offset applied to both grippers’ target positions. Before optimization (left), the arms fail to grasp the towel. After optimization (right), they grasp the towel and then lift and fold it. Rows show the motion at $0$ , $2.00$ , $2.40$ , and $3.30$ s. Four rows pair the same saved frames of the dual-arm towel grasp. Before optimization is on the left, and after optimization is on the right. The optimized arms grasp, lift, and fold the towel.

Figure 16. Optimization histories for the additional tasks. Stone placement changes its centroid tolerance after 50 updates. Garment Hanger, cup catching, and Dual-arm grasp use constant offsets that set their final plotted values to zero. Six loss plots for stone placement, cloth scale, Garment Hanger shape design, grasp angle, cup catching, and dual-arm grasp.

Example	System	Dim	Relative FD error
Stone placement (Fig. 1 )	DiffGIPC	14	$7.53\times 10^{-11}$
Garment Hanger (Fig. 4 )	DiffGIPC	2	$5.45\times 10^{-3}$
Bunny collision (Fig. 6 )	DiffGIPC	3	$1.21\times 10^{-6}$
Bunny collision (Fig. 6 )	DiffIPC	3	$1.01\times 10^{-4}$
Tentacle matching (Fig. 7 )	DiffGIPC	3	$3.33\times 10^{-5}$
Tentacle matching (Fig. 7 )	DiffIPC	3	$3.83\times 10^{-11}$

Table 2. Finite-difference checks at the initial optimization parameters.

Resolution / case	System	Forward						Backward				Forward + Backward
Resolution / case	System	Linear	CCD	Hessian	Search	Other	Total	Hessian	Linear	Other	Total	Forward + Backward
Bunny collision (Fig. 6 ): resolution study
933 vertices	DiffGIPC ABD+ABD	0.025	0.039	0.041	0.021	0.0023	0.13	0.0081	0.0042	0.0032	0.016	0.14
	DiffGIPC ABD+FEM	0.58	0.049	0.058	0.03	0.0025	0.72	0.011	0.22	0.0085	0.24	0.96
	DiffGIPC FEM+FEM	0.65	0.044	0.037	0.017	0.002	0.75	0.0069	0.22	0.0057	0.24	0.98
	DiffIPC FEM+FEM	4.5	10	4.6	1.1	1.9	22	1.4	1.1	0.44	2.9	25

Table 3. Mean forward and backward evaluation times (seconds) for the bunny-collision resolution study and Garment Hanger rest-shape design. Forward + Backward gives the sum of the two per-evaluation means.

Example	$N_{\mathrm{fwd}}$	$N_{\mathrm{bwd}}$	Forward mean							Backward mean					$T_{\mathrm{eval}}$
			Linear		CCD	Hessian	Search	Other	Total	Hessian	Linear		Other	Total
			PCG	Others	CCD	Hessian	Search	Other		Hessian	PCG	Others	Other
Stone placement (Fig. 1 )	80	78	0.097	0.0093	0.079	0.12	0.024	0.0088	0.34	0.051	0.067	0.003	0.023	0.14	39
Bunny collision (Fig. 6 )	15	12	0.62	0.027	0.044	0.037	0.017	0.002	0.75	0.0069	0.22	0.005	0.0057	0.24	14
Tentacle matching (Fig. 7 )	11	10	6.8	1.3	2.4	2.2	1.4	0.011	14	0.04	0.99	0.024	0.024	1.1	170
Garment Hanger (Fig. 4 )	32	30	2.8	0.87	0.88	1.5	0.8	0.0088	6.9	0.033	0.2	0.01	0.022	0.27	230

Table 4. Complete optimization workloads and mean phase times (seconds) of the eleven examples.

N_{\mathrm{fwd}}

counts forward evaluations, including line-search trials and verification runs.

N_{\mathrm{bwd}}

counts complete gradient evaluations used for optimization. Linear comprises PCG and Others, and Hessian combines local and global assembly.

Example	Graph	Forward (s)						Backward (s)						Total (s)
		Linear			Other	Total		Linear			Other	Total
		PCG		Others	Other	Total		PCG		Others	Other	Total
Stone placement (Fig. 1 )	On	0.097		0.0093	0.24	0.34		0.067		0.003	0.074	0.14		0.49
Stone placement (Fig. 1 )	Off	0.18	$(1.8\times)$	0.0094	0.24	0.42	$(1.2\times)$	0.12	$(1.8\times)$	0.003	0.073	0.2	$(1.4\times)$	0.62	$(1.3\times)$
Bunny collision (Fig. 6 )	On	0.62		0.027	0.1	0.75		0.22		0.005	0.013	0.24		0.98
Bunny collision (Fig. 6 )	Off	0.97	$(1.6\times)$	0.026	0.1	1.1	$(1.5\times)$	0.35	$(1.6\times)$	0.0051	0.011	0.36	$(1.5\times)$	1.5	$(1.5\times)$

Table 5. CUDA Graph enabled (On) versus disabled (Off). Times are mean seconds per evaluation; parentheses in the Off rows give speedups from enabling the graph. Linear comprises PCG and Others; the separate Other is Total minus Linear. The rightmost Total combines forward and backward time.

Appendix figures & tables2 assets

Supplementary material from the paper’s appendix.

Appendix

Example	$h$	Horizon	$d_{\mathrm{rel}}$	$\kappa$	$\widehat{d}_{f}$	#DOF	Newton $\epsilon_{N}$	Friction	$E$	$\kappa_{\mathrm{A}}$
Stone placement (Fig. 1 )	$0.03$	$1000$	$10^{-3}$	$10^{5}$	$10^{-3}$	$84$	$10^{-3}$	$0.3$	–	$10^{7}$
Bunny collision (Fig. 6 )	$0.05$	$40$	$10^{-3}$	$10^{10}$	–	$5598$	$3\times 10^{-4}$	$0$	$4{\times}10^{7}$	–
Tentacle matching (Fig. 7 )	$0.005$	$100$	$10^{-3}$	$10^{7}$	$10^{-3}$	$20688$	$10^{-2}$	$1$	$2{\times}10^{4}$	–
Garment Hanger (Fig. 4 )	$1/120$	$600$	$10^{-3}$	$10^{5}$	$10^{-3}$	$51357$	$10^{-2}$	$0.05$	$5{\times}10^{4}$	$10^{8}$
Loaded Support (Fig. 8 )	$0.02$	$20$	$3{\times}10^{-4}$	$10^{6}$	–	$22545$	$10^{-3}$	$0$	$10^{8}$	–
Material fitting (Fig. 9 )	$0.1$	$10$	$10^{-3}$	$10^{7}$	–	$12567$	$10^{-2}$	$0$	$10^{6}$	–

Appendix

Table 6. Simulation settings for the design, identification, and control tasks in Section 7 . The columns report time step

h

, relative collision activation distance

d_{\mathrm{rel}}

, barrier coefficient

\kappa

, friction-smoothing velocity threshold

\widehat{d}_{f}

, motion DoFs, relative Newton tolerance

\epsilon_{N}

, friction coefficient, Young’s modulus

E

for volumetric FEM or cloth stretching, and affine stiffness

\kappa_{\mathrm{A}}

. For Dual-arm grasp,

0.3

is the arm–cloth friction coefficient and

0.2

is the cloth–table coefficient.

Physical group	Differentiated quantities	Role
Tetrahedral material	$\mu$ , $\lambda$	Stable Neo-Hookean material response
Initial FEM motion	Selected initial vertex velocities	Initial dynamic state
Deformable rest shape	Selected rest coordinates	Volumetric and cloth rest geometry
Cloth material	Stretch, shear, and bending coefficients	Cloth constitutive response
Friction	Mesh self and ground friction coefficients	Frictional contact
Initial affine motion	$\mathbf{q}_{0,b}$ , $\dot{\mathbf{q}}_{0,b}$	Initial affine state

Appendix

Table 7. Physical parameter families exposed to task-specific design variables through the interface in Section 6.3 . Each family contributes local derivatives to the shared coupled adjoint.

Differentiable Physics

Differentiable physics makes physical simulators usable for gradient-based learning, parameter estimation, design, and control.

See all 66 papers

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration
May 14, 2026Shih-Yu Lai, Sung-Han Tien, Jui-I Huang +9Robotics Simulation Hyperelasticity
National Taiwan University · MoonShine Animation Studio · National University of Singapore +3
Novel Algorithms for Smoothly Differentiable and Efficiently Vectorizable Contact Manifold Construction
Apr 19, 2026Onur Beker, Andreas René Geist, Anselm Paulus +1Signed Distance Functions Contact-Rich Robotic Manipulation
University of Tübingen
IsaacIPC: Coupling High-Fidelity Simulation and Realistic Rendering for Contact-Rich Robotic Systems
May 23, 2026Qixin Liang, Zhongqing HanContact-Rich Robotic Manipulation Robotics Simulation
Anker Humanoid Lab · Shenzhen, China · The University of Hong Kong +1

May 14, 2026J/K move · Enter open · S save

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Shih-Yu Lai, Sung-Han Tien, Jui-I Huang, Yen-Chen Tseng+8

National Taiwan University · MoonShine Animation Studio · National University of Singapore+3

Robotics Simulation Hyperelasticity Differentiable Programming Differentiable Physics Physics-Based Simulation Deformable Object Simulation

Open paper AskPDF arXiv

Differentiable simulation of soft bodies is a foundation for system identification, trajectory optimization, and Real2Sim transfer. Yet, existing methods such as the differentiable Projective Dynamics (DiffPD) struggle when faced with heterogeneous materials with extreme stiffness contrasts, hyperelasticity under large deformations, and contact-rich interactions, which are common scenarios in the real world. We present DiffPhD, a unified GPU-accelerated differentiable Projective Dynamics framework for heterogeneous materials that tackles these intertwined challenges simultaneously. Our key insight is a careful integration of: (i) stiffness-aware projective weights to embed heterogeneity into the global system; (ii) trust-region eigenvalue filtering lifted to the backward pass for stable hyperelastic gradients and a type-II Anderson Acceleration scheme with dual-gate convergence to stabilize forward iteration under large stiffness contrasts; and (iii) a unified GPU pipeline that reuses a single sparse factor across forward, backward, and contact computations, with stiffness-amplified Rayleigh damping folded into the same factor for heterogeneity-aware dissipation at zero recurring cost. DiffPhD achieves strict gradient accuracy while delivering up to an order-of-magnitude speedup over prior differentiable solvers on heterogeneous, hyperelastic, contact-rich benchmarks. Crucially, this speedup does not come at the cost of stability: DiffPhD remains convergent on stiffness contrasts up to 100x where prior PD solvers degrade. This unlocks end-to-end gradient-based optimization on regimes previously bottlenecked by either solver fragility or per-iteration cost -- shell--joint composite creatures, soft characters wielding stiff weapons, and soft-gripper robotic manipulation -- all handled within a single forward--backward pass.

Abstract

Differentiable Physics

Explore similar work

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Novel Algorithms for Smoothly Differentiable and Efficiently Vectorizable Contact Manifold Construction

IsaacIPC: Coupling High-Fidelity Simulation and Realistic Rendering for Contact-Rich Robotic Systems