# How to Set Up MPM Simulations for Granular Materials in Newton

> Learn to set up MPM simulations for granular materials in Newton. Register custom attributes, configure physics, and run stable substeps for effective collision handling.

- Repository: [Newton Physics/newton](https://github.com/newton-physics/newton)
- Tags: how-to-guide
- Published: 2026-03-19

---

**Register custom MPM attributes via `SolverImplicitMPM.register_custom_attributes()`, configure granular physics in `SolverImplicitMPM.Config`, and run double-buffered substeps with `solver.step()` and `_project_outside()` for robust collision handling.**

The Newton physics engine provides a high-performance implicit Material Point Method (MPM) solver designed for simulating granular materials like sand, gravel, and soil. By leveraging CUDA-accelerated kernels and a sparse grid architecture, Newton's `SolverImplicitMPM` enables stable simulations of stiff, frictional particles at large scales. This guide explains the exact implementation details found in the Newton source code, covering attribute registration, solver configuration, and the simulation loop.

## Registering MPM Custom Attributes

Before adding any particles or geometry, you must register MPM-specific custom attributes to the `ModelBuilder`. This step creates the necessary namespaces and data structures for per-particle material parameters and state fields.

In [`newton/_src/solvers/implicit_mpm/solver_implicit_mpm.py`](https://github.com/newton-physics/newton/blob/main/newton/_src/solvers/implicit_mpm/solver_implicit_mpm.py), the `SolverImplicitMPM.register_custom_attributes` method injects a `mpm:` namespace into the model:

```python
builder = newton.ModelBuilder()
SolverImplicitMPM.register_custom_attributes(builder)

```

This registration adds **material attributes** (attached to the model) such as `young_modulus`, `friction`, `hardening`, and `yield_stress`:

```python
builder.add_custom_attribute(
    newton.ModelBuilder.CustomAttribute(
        name="young_modulus", frequency=newton.Model.AttributeFrequency.PARTICLE,
        assignment=newton.Model.AttributeAssignment.MODEL,
        dtype=wp.float32, default=1.0e15, namespace="mpm"
    )
)

# ... similar calls for poisson_ratio, damping, hardening, friction, yield_pressure

```

It also registers **state attributes** (attached to each particle's `State` object) including deformation tracking fields:

```python
builder.add_custom_attribute(
    newton.ModelBuilder.CustomAttribute(
        name="particle_elastic_strain", frequency=newton.Model.AttributeFrequency.PARTICLE,
        assignment=newton.Model.AttributeAssignment.STATE,
        dtype=wp.mat33, default=identity, namespace="mpm"
    )
)

# ... particle_Jp, particle_qd_grad, particle_transform

```

After registration, these fields are accessible via `model.mpm.young_modulus` and `state.mpm.particle_elastic_strain`.

## Configuring SolverImplicitMPM for Granular Media

Granular materials require specific rheological parameters to simulate correctly. The `SolverImplicitMPM.Config` class, defined in [`solver_implicit_mpm.py`](https://github.com/newton-physics/newton/blob/main/solver_implicit_mpm.py), exposes all user-tunable knobs for the implicit MPM solve.

Typical configurations for sand-like materials include:

- **Grid and Transfer**: Set `grid_type="sparse"` for memory efficiency and `transfer_scheme="apic"` to preserve angular momentum during particle-grid transfers.
- **Solver Settings**: Use `solver="gauss-seidel"` (default) for the rheology linear system, with `max_iterations=250` and `tolerance=1e-6` for convergence.
- **Material Parameters**: Set `young_modulus=1e15` for stiff grains, `friction=0.68` for inter-grain friction, and `hardening=0.0` for perfectly plastic behavior.

Configuration example:

```python
options = SolverImplicitMPM.Config()
options.voxel_size = 0.1
options.grid_type = "sparse"
options.transfer_scheme = "apic"
options.solver = "gauss-seidel"
options.max_iterations = 250
options.tolerance = 1e-6
options.friction = 0.68
options.young_modulus = 1e15

```

After finalizing the model, propagate these parameters to per-particle arrays:

```python
model.mpm.friction.fill_(options.friction)
model.mpm.young_modulus.fill_(options.young_modulus)

solver = SolverImplicitMPM(model, options)

```

## Constructing the Particle Geometry

With attributes registered, emit particles using `builder.add_particle_grid`. The particle resolution should align with the MPM voxel size to ensure proper grid coverage. The implementation in [`newton/examples/mpm/example_mpm_granular.py`](https://github.com/newton-physics/newton/blob/main/newton/examples/mpm/example_mpm_granular.py) demonstrates emitting a block of particles and adding collision geometry:

```python
import numpy as np
import warp as wp

# Define particle grid based on voxel size

voxel = 0.1
lo, hi = np.array([-1, -1, 1.5]), np.array([1, 1, 3.5])
particles_per_cell = 3
res = np.ceil(particles_per_cell * (hi - lo) / voxel).astype(int)
cell_sz = (hi - lo) / res
mass = np.prod(cell_sz) * 1000.0  # density 1000 kg/m³

radius = np.max(cell_sz) * 0.5

builder.add_particle_grid(
    pos=wp.vec3(lo),
    rot=wp.quat_identity(),
    vel=wp.vec3(0.0),
    dim_x=res[0] + 1,
    dim_y=res[1] + 1,
    dim_z=res[2] + 1,
    cell_x=cell_sz[0],
    cell_y=cell_sz[1],
    cell_z=cell_sz[2],
    mass=mass,
    jitter=2.0 * radius,
    radius_mean=radius,
)

# Add ground plane and optional colliders

builder.add_ground_plane(cfg=newton.ModelBuilder.ShapeConfig(mu=0.5))
model = builder.finalize()
model.set_gravity([0, 0, -10])

```

## Executing the MPM Time Step

The simulation loop follows a double-buffered pattern where `state_0` and `state_1` alternate as source and target. According to [`example_mpm_granular.py`](https://github.com/newton-physics/newton/blob/main/example_mpm_granular.py), each frame consists of multiple substeps that call `SolverImplicitMPM.step` followed by collision projection:

```python
state_a = model.state()
state_b = model.state()
dt = 1.0 / 60.0 / 2  # two substeps per frame

for _ in range(self.sim_substeps):
    solver.step(state_a, state_b, None, None, self.sim_dt)
    solver._project_outside(state_b, state_b, self.sim_dt)
    state_a, state_b = state_b, state_a

```

The `step` method in [`solver_implicit_mpm.py`](https://github.com/newton-physics/newton/blob/main/solver_implicit_mpm.py) performs:

- **Binning and rasterization**: Particles are binned into grid cells and velocity fields are populated via APIC/PIC transfer.
- **Matrix assembly**: A sparse linear system couples strain degrees of freedom to velocity degrees of freedom.
- **Implicit solve**: Gauss-Seidel or Jacobi iterations solve the rheology system, using warm-start data from `LastStepData` when available.
- **Advection**: Particle positions, velocities, and deformation gradients (`particle_transform`) are updated using the solved grid velocities.

The `_project_outside` method enforces collision constraints by pushing particles out of collider signed-distance fields (SDFs), ensuring no interpenetration with the ground plane or complex meshes.

## Accelerating Simulations with CUDA Graphs

For maximum performance on CUDA devices, Newton supports CUDA Graph capture when using a fixed grid. The `ImplicitMPMScratchpad` and `LastStepData` structures are allocated once in `SolverImplicitMPM.__init__` and reused every substep to minimize memory overhead.

If `grid_type="fixed"`, wrap the substep loop in `wp.ScopedCapture` to record the sequence of kernels, then replay with `wp.capture_launch(self.graph)` for near-zero CPU launch overhead in subsequent frames.

## Summary

- **Register attributes first**: Always call `SolverImplicitMPM.register_custom_attributes(builder)` before adding particles to initialize the `mpm:` namespace for material and state fields.
- **Configure granular physics**: Use `SolverImplicitMPM.Config` to set stiff elastic moduli (`young_modulus=1e15`), inter-grain friction (`friction=0.68`), and APIC transfer for angular momentum preservation.
- **Align particles with voxels**: Emit particles via `builder.add_particle_grid` using a resolution that matches the MPM `voxel_size` to ensure proper grid coverage.
- **Run double-buffered substeps**: Alternate between two state objects in the simulation loop, calling `solver.step()` for the MPM integration and `_project_outside()` for collision handling.
- **Leverage CUDA graphs**: Capture the step loop as a CUDA graph when using fixed grids to eliminate CPU overhead during playback.

## Frequently Asked Questions

### What is the difference between APIC and PIC transfer schemes in Newton MPM?

The **APIC** (Affine Particle-in-Cell) scheme preserves angular momentum and reduces numerical dissipation by transferring affine velocity gradients between particles and the grid, making it ideal for granular materials that exhibit rotational behavior. The **PIC** (Particle-in-Cell) scheme transfers only velocity, which is cheaper computationally but introduces more numerical damping. Set `transfer_scheme="apic"` in `SolverImplicitMPM.Config` for accurate granular flow, as implemented in [`solver_implicit_mpm.py`](https://github.com/newton-physics/newton/blob/main/solver_implicit_mpm.py).

### How do I set material properties like friction and stiffness for individual particles?

After calling `model = builder.finalize()`, material properties are exposed as Warp arrays under the `model.mpm` namespace. Fill these arrays uniformly or per-particle using Warp indexing:

```python
model.mpm.friction.fill_(0.68)
model.mpm.young_modulus.fill_(1e15)

```

These arrays map directly to the custom attributes registered earlier and are read by the rheology solver during the implicit integration step.

### Why does the MPM solver require `register_custom_attributes` before adding geometry?

The `ModelBuilder` must know the memory layout for per-particle data before allocating particle buffers. `register_custom_attributes` declares the `mpm:` namespace fields—such as `particle_elastic_strain` and `particle_Jp`—so that `builder.finalize()` can allocate the correct amount of state memory. Adding particles before registration would result in missing MPM data structures, causing runtime errors when `SolverImplicitMPM` attempts to access strain or deformation data during the rasterization phase.

### How does Newton handle collisions between granular particles and complex colliders?

Newton uses a **signed-distance field (SDF)** approach via `solver._project_outside()`. After each MPM substep, this method pushes particles out of collider geometry using the SDF of shapes added via `builder.add_shape_box()` or similar methods. This ensures non-penetration constraints are satisfied independently of the MPM grid resolution, allowing granular materials to interact robustly with concave meshes, wedges, and ground planes as demonstrated in [`example_mpm_granular.py`](https://github.com/newton-physics/newton/blob/main/example_mpm_granular.py).