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

Abstract: 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.

• The paper introduces DiffPhD, a differentiable simulation framework that robustly solves elastodynamics with material heterogeneity and contact, achieving strict gradient consistency.
• It leverages stiffness-aware heterogeneity routing and adaptive trust-region filtering to overcome challenges in non-convex, hyperelastic models.
• The approach delivers order-of-magnitude speedups via persistent, GPU-resident operator factorization, enabling efficient inverse design and real-to-sim robotic manipulation.

DiffPhD: Differentiable Projective Dynamics for Heterogeneous, Hyperelastic, and Contact-Rich Elastodynamics

Overview

The paper "DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration" (2605.14526) introduces DiffPhD, a framework for differentiable simulation of elastodynamics that robustly addresses material heterogeneity, large-strain hyperelasticities (notably Neo-Hookean), and large-scale frictional contact, all with GPU-resident acceleration and a persistent, amortized factorization-based computational core. The solver supports end-to-end optimization for system identification, inverse design, trajectory optimization, and real-to-sim transfer in robotics, achieving strict gradient consistency and substantial speedups and robustness gains compared to prior differentiable solvers, especially in regimes with extreme stiffness contrast and contact density.

Context and Problem Formulation

Recent advances in differentiable physics-based simulation have enabled inverse problems in graphics and robotics, but existing differentiable Projective Dynamics (PD) pipelines break down in practical scenarios with (i) strong spatial material heterogeneity, (ii) hyperelastic deformations with non-convex and indefinite element energy landscapes, and (iii) frictional contact with large and active contact sets. PD-based methods depend on pre-factored linear global operators, but high stiffness contrast and non-convexity destroy conditioning and gradient reliability.

The canonical energy models, particularly compressible Neo-Hookean, introduce element-level Hessian indefiniteness under large strain or near inversion. Meanwhile, heterogeneous Young’s modulus and Poisson ratio variation, especially with $\mu_{\max}/\mu_{\min} \geq 100\times$, produce global operators with unbounded condition numbers, making previous differentiation and optimization schemes numerically unstable or prohibitively slow.

Methodological Innovations

Stiffness-Aware Heterogeneity Routing

The central architecture of DiffPhD routes all material heterogeneity through projective weights, making the global PD stiffness operator $\bm{A}$ explicitly encode the contrast, while maintaining a uniform, mesh-average local prox-map in the per-element Newton iteration. This strategy ensures that the fixed-point contraction assumption required for Anderson Acceleration (AA) holds, amortizing the cost of the pre-factorized global operator across iterations and passes (forward, backward, contact).

Figure 1: Two consequences of material heterogeneity in Neo-Hookean PD—the spectrum of the global tangent stiffness and the element-localized non-convexity.

Figure 2: Stiffness-aware projective assembly: heterogeneity routed into $\bm{A}$, not local prox, avoiding mesh-level oscillations.

State-Adaptive Trust-Region Filtering for Hyperelasticity

To address breakdowns in the backward gradient due to indefinite prox-map Hessians in the Neo-Hookean model, DiffPhD uses an adaptive trust-region eigenvalue filtering, selecting (per element and per iteration) among unfiltered, clamped, or absolute-value-filtered Hessian surrogates in the backward pass according to a model fidelity criterion. This mechanism guarantees strict gradient descent directions and preserves differentiability where prior methods fail due to indefinite or poorly conditioned local solves.

Figure 3: Trust-region filtering on the prox-map Hessian, applied in the backward pass only, stabilizes gradients in high-strain regimes.

Unified Persistent Factorization and GPU Pipeline

DiffPhD employs a single, METIS-reordered, nested-dissection-based sparse Cholesky factorization for the global operator $\bm{A}$, which is uploaded once to the GPU and shared across forward solves, adjoint backward passes, and dense frictional contact Delassus computations. This persistent factor enables highly efficient, exact two-SpMV applications for every operator inversion and SpMM for batches of contact basis vectors.

Figure 4: The persistent SPD factor $\bm{A}$ is shared verbatim across forward, backward, and contact passes.

Results

Robustness to Extreme Material Heterogeneity

DiffPhD remains stable and strictly gradient-consistent across cases with up to $100\times$ stiffness contrast, where baseline differentiable PD and domain decomposition preconditioners diverge or produce spurious gradients. The Armadillo and Crab benchmarks demonstrate DiffPhD's capacity to partition and resolve multi-material bodies with minimal deformation artifact at high-contrast boundaries, while prior methods fail or require per-iteration costly refactorization.

Figure 5: Armadillo under varying stiffness contrast; DiffPhD remains stable as the upper torso/legs ratio increases up to $100\times$.

Figure 6: Crab pull-up—homogeneous model (bottom) vs. heterogeneous with stiff shell (top); only heterogeneity enables true load transfer.

Figure 7: Crab pull-up, differentiable solver comparison showing only DiffPhD remains stable in high-contrast, contact-rich configurations.

Performance and Speed

DiffPhD attains order-of-magnitude wall-clock improvements over previous differentiable solvers in both forward and backward passes, especially evident in contact-rich and heterogeneous scenes. The forward solve on a 2.2M-DoF heterogeneous model is accelerated by $8.7\times$ compared to the next-strongest baseline, with backward costs similarly improved or better, all while maintaining strict gradient correctness.

Large-Scale and High-Density Frictional Contact

The framework supports mesh-to-mesh and codimensional contact scenarios of escalating challenge, including soft-hard, soft-soft, and degenerately high-density contact. The unified contact pipeline leverages the persistent factor for efficient exact Delassus and contact system solution. DiffPhD substantially outperforms baseline differentiable PD frameworks as contact density and set size increase.

Figure 8: Gatorman, a mesh-to-mesh contact benchmark—with a rigid sword impacting a soft target, resolved without penetration or oscillation.

Figure 9: Heterogeneous codimensional napkin; deformations accurately reflect varying compliance subject to increasing contact ratios.

Differentiable Inverse and Optimization Applications

DiffPhD directly enables system identification (SI), initial-state optimization, and trajectory optimization (TO) in regimes previously inaccessible due to solver fragility or gradient inaccuracy. Strong claims include:

• Convergence to loss optima $23.9\times$ lower than the strongest baseline in slender, highly non-uniform SI settings (e.g., plant identification).
• Reliable SI over multi-material bodies with asymmetric, heterogeneous impacts (e.g., the three-material Bouncing Ball).
• Strict gradient correctness through contact transitions, shown in initial-state and muscle-driven optimization under heterogeneous contact.

  Figure 10: Bouncing Ball SI—joint parameter recovery for homogeneous (left) and heterogeneous (right) settings.

  

  Figure 11: Plant SI—optimized trajectories match ground truth across multiple sub-regions; per-region Young’s moduli are simultaneously identified.

  

  Figure 12: Bunny inverse design—DiffPhD brings the center of mass to the target; gradient filtering stabilizes backward passes through contact transitions.

  

  Figure 13: Bunny inverse design, parameter and loss convergence; DiffPhD finds a $3.05\times$ lower loss than DiffPD.

  

  Figure 14: Routing Tendon—muscle-driven actuation, showcasing convergence for complex fibre layouts.

  

  Figure 15: Torus locomotion—optimized peristaltic gait yields $\bm{A}$0 longer displacement.

Real-to-Sim Robotic Manipulation

DiffPhD bridges simulation and real-world manipulation tasks. Using real (PokeFlex) and simulated (Google Robot) data, it recovers contact locations and material distributions via gradient-based optimization, supporting digital twins with precise contact alignment.

Figure 16: Real2Sim alignment—DiffPhD reproduces real-world manipulator contact and deformation, supporting spatially varying compliance within a unified pipeline.

Ablation Studies and Failure Mode Analysis

Systematic ablations confirm that the observed performance and robustness gains require both (i) the adaptive trust-region eigenvalue filtering and (ii) persistent GPU-resident operator reuse. Fixed projection rules alone are insufficient to eliminate gradient stalls in system identification; naive heterogeneity routing induces forward or backward oscillations; and loss-only diagnostics fail to reveal gradient correctness, demanding per-component invariants and specialized per-frame metrics.

Limitations

Velocity-level NCP enforcement does not guarantee strict position-level nonpenetration in regimes of sustained static contact (e.g., grasp closure under gravity), occasionally resulting in centimeter-level drift over long trajectories. Material Jacobian differentiation is energy-specific and must be symbolically specified for each new constitutive law. Both axes are acknowledged as engineering limitations, not fundamental barriers to extension.

Theoretical and Practical Implications

The architecture of DiffPhD establishes that persistent operator amortization—rather than mere arithmetic throughput per iteration—is the key design axis for scalable differentiable elastodynamics, especially for heterogeneous, hyperelastic, and friction/contact-dominated regimes. The method is not specific to PD: the same per-element, per-iteration local filtering, global stiffness routing, and factor persistence pattern applies equally to Newton-type and position-based differentiable solvers.

Future Directions

The authors propose symbolic differentiation of the prox-map for generic, user-added hyperelastic energies and position-level NCP augmentation for quasi-static contact and manipulation. The broader implication is that future differentiable physics engines, especially for robotic manipulation, should route all sources of heterogeneity and non-convexity into a shared, persistent global operator, and adaptively filter only where local pathology arises, decoupling robustness from global loss of structure.

Conclusion

DiffPhD sets a new benchmark for GPU-accelerated, end-to-end-differentiable elastodynamics in graphics and robotics, robustly addressing the hardest cases of material contrast, hyperelastic non-convexity, and dense frictional contact without loss of gradient consistency. The technical innovations—stiffness-aware projective operator assembly, adaptive trust-region filtering, and unified persistent factorization—systematically resolve architectural and numerical failure modes endemic to prior pipelines. The design principles demonstrated here generalize beyond PD, suggesting a roadmap for future differentiable simulation in high-fidelity, optimization-driven applications.