Differentiable MPM Simulator

• Differentiable MPM is an advanced computational framework that propagates exact gradients through particle-grid interactions, supporting nonlinear materials and contact handling.
• It combines Lagrangian particles with an Eulerian grid to perform large-deformation simulation using integrated automatic differentiation and sparse Jacobian assembly.
• Key applications include inverse modeling, optimal control, and design optimization in domains like geomechanics, soft robotics, and system identification.

A differentiable Material Point Method (MPM) simulator is an advanced computational framework that enables the propagation of exact gradients through the full MPM pipeline, supporting inverse modeling, learning-augmented simulation, design optimization, and control involving complex, history-dependent, large-deformation material behavior. The core technical objective is to construct an MPM architecture where the mapping from model parameters, initial conditions, or control sequences to simulation outputs is fully differentiable—even in the presence of nonlinear constitutive models, implicit or explicit time integration, and advanced features such as contact, friction, or topology change.

1. Mathematical and Algorithmic Foundations

MPM discretizes continuum mechanics problems using both Lagrangian material points (particles) and an Eulerian mesh (grid). The discrete update cycle typically advances in three primary steps per time step:

1. Particle-to-Grid (P2G) Transfer: Particles carry mass, momentum, deformation gradient, and internal states, which are mapped onto the supporting grid nodes by interpolation with kernel functions.
2. Grid Operations: Forces—including internal (stress-based) and external—are assembled by combining nodal contributions and evolving grid momenta or velocities according to dynamics (explicit or implicit time integration). Boundary and contact conditions are imposed here if applicable.
3. Grid-to-Particle (G2P) Transfer: Updated velocities, accelerations, and other relevant quantities are interpolated back onto the particles, updating their velocities, positions, deformation gradients, and history.

Typical notations include:

• Reference configuration: $X \in \Omega_0$
• Current configuration: $x = X + u(X,t)$
• Deformation gradient: $F = \frac{\partial x}{\partial X} = I + \nabla_X u$
• Jacobian determinant: $J = \det F$

The governing (weak) form for dynamic or quasi-static momentum balance after MPM discretization yields a nonlinear residual equation,

$R(u) = f_{\text{int}}(u) - f_{\text{ext}} = 0\ ,$

where $u$ is the vector of unknown nodal displacements, $f_{\text{int}}$ is the assembled internal force vector (from particle stresses), and $f_{\text{ext}}$ includes gravity, tractions, and other external sources (Zhao et al., 13 Jul 2025).

For implicit solvers, the consistent tangent (Jacobian)

$J = \frac{\partial R}{\partial u}$

captures the non-trivial chain between grid DOFs and particle states, aggregating all $\partial \sigma/\partial \epsilon$ for all active particles.

Constitutive updates span a range of models, including:

• Large deformation Hencky elasticity: $\epsilon = \log U$, $$\sigma = (2\mu/J)\,\mathrm{dev}[\epsilon] + (\kappa/J)[\mathrm{tr}\,\epsilon]I$$
• $J_2$ elastoplasticity with return mapping and consistent tangent
• Biot poromechanics coupling fluid and solid phases (Zhao et al., 13 Jul 2025)

Explicit time-integration (FLIP, PIC, RK4) is widespread for dynamics, but implicit (Newton–Raphson–type) solvers are crucial for quasi-static, stiff, or highly nonlinear mechanics.

2. Differentiability and Automatic Differentiation Techniques

Differentiability is achieved by expressing all MPM kernels—P2G, stress update, force assembly, grid evolution, G2P, and auxiliary routines—as differentiable elemental operations assembled into a computational graph.

• Reverse-Mode Automatic Differentiation (AD):
  • Forward pass records a computation tape, with local subgraphs corresponding to each particle, including reads/writes to support nodes and constitutive kernels.
  • Backward pass propagates gradients from outputs (e.g., costs, loss, or residuals) back to any upstream parameter: grid state, particle variable, or physical parameter (Zhao et al., 13 Jul 2025, Du et al., 6 Jul 2025, Hu et al., 2018).
• Kernel-Level Differentiation:
  • In frameworks such as NVIDIA Warp or JAX, each kernel (e.g., P2G, stress update) is either traced or written with analytical Jacobians so that AD can traverse all necessary code paths (Zhao et al., 13 Jul 2025, Du et al., 6 Jul 2025).
  • In ChainQueen, all key update and transfer steps are differentiated analytically with hand-derived Jacobians to maximize efficiency and minimize memory overhead (Hu et al., 2018).
• Contact, Collision, and Conditional Logic:
  • Differentiable collision handling is realized via smoothed velocity projections, sign-based contact splitting, and compatible node-particle transfers (Vasile et al., 10 Nov 2025).
  • All non-smooth operations are implemented as smooth functions or as subgradient-producing primitives (Bolliger et al., 15 Dec 2025, Vasile et al., 10 Nov 2025, Du et al., 6 Jul 2025).

AD supports gradients with respect to:

• Initial or boundary conditions
• Control variables
• Material constitutive parameters
• Embedded neural-network closures

3. Sparse Jacobian Construction and Solver Implementation

Efficient construction of the global tangent (Jacobian) is critical for scalable implicit simulation. Differentiable MPM implementations exploit highly localized particle-grid interactions to dramatically reduce AD passes:

• In, e.g., GIMP or CPDI (grid-interpolated particle or convected particle domain interpolation), each particle influences a $b^3$ block of grid nodes; only these blocks must be considered for assembling the global sparse Jacobian (Zhao et al., 13 Jul 2025).
• A block-wise seeding algorithm is used, constructing one seed per local support offset (total $b^3$), so that each AD backward pass yields all rows of $J$ for a single block across the domain.
• The resulting Jacobian is written into block-CSR or regular CSR matrix formats, compatible with GPU-accelerated sparse linear solvers (Zhao et al., 13 Jul 2025).

For solving the resulting nonlinear systems:

• Newton–Raphson iterations update $u$ using $J \Delta u = -R$, with each inner linear solve employing iterative solvers with (block) preconditioning (AMG, block BiCGStab, or direct sparse factorization as problem size dictates) (Zhao et al., 13 Jul 2025).
• Accelerator support enables end-to-end GPU computation, with negligible Jacobian assembly time compared to total runtime at large scale.

4. Constitutive Models and Inverse Modeling Capabilities

Differentiable MPM supports a range of constitutive behaviors:

• Hyperelasticity (e.g., Saint Venant–Kirchhoff, Neo-Hookean, Hencky)
• Elastoplasticity (e.g., $J_2$, Drucker–Prager, von Mises, with tension cutoff and return mapping)
• Poromechanics (Biot theory coupling solid and fluid mass balance)
• Newtonian fluids or non-Newtonian rheologies
• Frictional and cohesive contact
• Learned closures using neural networks (MLP-parameterized spatial fields) (Du et al., 6 Jul 2025)

Inverse modeling, system identification, and optimal control are made feasible by treating simulation as a mapping $y = f(\theta)$ from model parameters $\theta$ to physical observables. Differentiable MPM enables:

• Gradient-based optimization of control trajectories, as demonstrated for active damping of hyperelastic ropes, yielding faster convergence and lower final energies than sampling-based baselines (Bolliger et al., 15 Dec 2025).
• Parameter identification from target observations, e.g., adjusting Young’s modulus to match indentation responses or reconstructing spatially varying friction from data (Du et al., 6 Jul 2025, Vasile et al., 10 Nov 2025).

The optimization loop consists of forward simulation (unroll steps; or implicit solve), loss evaluation, AD-based backward pass to compute gradients, and an optimizer step (e.g., Adam, L-BFGS) updating the control or parameter variables.

5. GPU-Accelerated Implementation and Performance

All modern differentiable MPM simulators employ GPU-based acceleration to handle the large computational requirements:

• Data layout leverages structure-of-arrays, coalesced access, and shared memory tiling to maximize effective bandwidth and occupancy on modern GPUs (e.g., NVIDIA A100, 1080Ti) (Zhao et al., 13 Jul 2025, Hu et al., 2018).
• Kernels for P2G/G2P, force assembly, and reductions are written in CUDA (or JAX lower-level XLA) to enable full device parallelization.
• AD backward passes are similarly kernelized, with per-time-step caching (memoization) or checkpointing to balance memory and computation (Hu et al., 2018, Du et al., 6 Jul 2025).
• Multi-GPU scaling, with domain decomposition and halo exchange for both primal and adjoint variables, allows nearly linear strong scaling at the millions-of-particles scale (Zhao et al., 13 Jul 2025).

Empirically, runtimes of large 3D problems (e.g., 512k–2.7M particles) reach forward+Jacobian step times of order 22–450 seconds per 1000 steps (single precision), and the Jacobian assembly constitutes only a few percent of total solve time on high-end hardware (Zhao et al., 13 Jul 2025, Du et al., 6 Jul 2025). Sparse-AD approaches yield $>8\times$ speedup over dense AD at 2D/3D scale (Zhao et al., 13 Jul 2025).

6. Application Domains and Impact

Differentiable MPM simulators have been demonstrated in a spectrum of advanced applications:

• Geomechanics and Coupled Poromechanics: Forward and inverse analysis involving large-deformation elastoplasticity, Biot consolidation, triaxial tests, and soil-structure interaction (Zhao et al., 13 Jul 2025).
• Soft Robotics and Control: System identification, actuator optimization, and controller policy learning for soft robots with thousands of state and design variables, outperforming or accelerating relative to reinforcement learning and sampling-based techniques (Hu et al., 2018, Bolliger et al., 15 Dec 2025).
• Physics-Based Morphing and Shape Control: Space-time control of per-particle deformation gradients enables physically grounded morphing of elastic bodies, supporting topological transitions under automatic differentiation for loss-driven design (Xu et al., 24 Sep 2024).
• Visual System Identification with Arbitrary Colliders: Integration with differentiable rendering (e.g., Gaussian splatting, NeRF-based views), enabling end-to-end parameter identification from multi-view videos with non-planar rigid-body collision (Vasile et al., 10 Nov 2025).

The unification of forward simulation and exact differentiation supports robust, scalable, and extensible modeling pipelines for computational mechanics, scientific machine learning, and inverse problems in engineering and graphics.

7. Representative Workflows and Illustrative Pseudocode

All major pipelines consist of the following modular workflow:

1. Setup: Discretize the domain into particles and background grid, assign initial/boundary conditions and constitutive parameters.
2. Forward Simulation: For $T$ time steps (explicit/integrator or implicit/newton-backend), execute P2G, grid update, G2P, and state update for all particles and nodes.
3. Loss Evaluation: Compute a scalar-valued cost (physics error, misfit, kinetic energy, control cost, visual observation) at the final (or sequence of) state(s).
4. Backward Pass (AD): Compute $\nabla_\theta \mathcal{L}$ with respect to all upstream variables of interest (controls, parameters, initial states).
5. Optimization Loop: Update variables using a gradient-based optimizer, e.g., Adam, L-BFGS, or (for PDE-constrained design) projected gradient or SQP.

Canonical pseudocode for control optimization (Bolliger et al., 15 Dec 2025):

function RolloutAndCost(S0, U[0..T-1]): S ← S0 C ← 0 for t in 0..T-1: S ← MPM_RK4_Step(S, U[t]) C ← C + ½ * sum_p (m_p * ‖v_p(S)‖²) return C U ← initial_guess() for iter in 1..N_steps: C ← RolloutAndCost(S0, U) grads ← grad(RolloutAndCost, wrt U)(S0, U) U ← Adam_Update(U, grads)

This end-to-end differentiable pipeline enables advanced workflows in inverse geomechanical modeling (Zhao et al., 13 Jul 2025, Du et al., 6 Jul 2025), robot co-design (Hu et al., 2018, Xu et al., 24 Sep 2024), and data-driven system identification (Vasile et al., 10 Nov 2025).

---

• "GeoWarp: An automatically differentiable and GPU-accelerated implicit MPM framework for geomechanics based on NVIDIA Warp" (Zhao et al., 13 Jul 2025)
• "Differentiable Material Point Method for the Control of Deformable Objects" (Bolliger et al., 15 Dec 2025)
• "ChainQueen: A Real-Time Differentiable Physical Simulator for Soft Robotics" (Hu et al., 2018)
• "JAX-MPM: A Learning-Augmented Differentiable Meshfree Framework for GPU-Accelerated Lagrangian Simulation and Geophysical Inverse Modeling" (Du et al., 6 Jul 2025)
• "Gaussian-Augmented Physics Simulation and System Identification with Complex Colliders" (Vasile et al., 10 Nov 2025)
• "A Differentiable Material Point Method Framework for Shape Morphing" (Xu et al., 24 Sep 2024)