newton.solvers.SolverImplicitMPM#

class newton.solvers.SolverImplicitMPM(model, config, *, temporary_store=None, verbose=None, enable_timers=False)[source]#

Bases: SolverBase, CouplingInterface

Implicit MPM solver for granular and elasto-plastic materials.

Implements an implicit Material Point Method (MPM) algorithm roughly following [1], extended with a GPU-friendly rheology solver supporting pressure-dependent yield (Drucker-Prager), viscosity, dilatancy, and isotropic hardening/softening.

This variant is particularly well-suited for very stiff materials and the fully inelastic limit. It is less versatile than traditional explicit MPM but offers unconditional stability with respect to the time step.

Call register_custom_attributes() on your ModelBuilder before building the model to enable the MPM-specific per-particle material parameters and state variables (e.g. mpm:young_modulus, mpm:friction, mpm:particle_elastic_strain).

Multi-world models retain the shared FEM topology by default. Set Config.separate_worlds to use an independent FEM environment for each world. In isolated mode, every MPM particle must belong to a local world and particles must be stored contiguously by world. World-local colliders affect only that world; global static colliders and global colliders backed by kinematic bodies affect every world, while global colliders backed by dynamic bodies are rejected.

A sparse grid is rebuildable when Config.max_active_cell_count is positive, Config.grid_padding is zero, the velocity basis is "Q1", and the strain and collider bases support rebuilding. Cell and node capacities are totals across all FEM environments, and resolved capacities must satisfy upper <= lower <= leaf <= active.

Outer graph capture requires CUDA, an enabled memory pool, conditional graph support, enable_timers=False, positive active-cell capacity, and either a fixed or rebuildable sparse grid. Construction creates persistent capture resources internally. Run an uncaptured warm-up step for rheology modes that build an inner graph lazily, and keep model topology, solver configuration, captured buffers, and step arguments fixed across replays. Keep reset() outside capture and call check_sparse_grid_rebuild_status() after replay to report sparse-grid capacity failures.

[1] https://doi.org/10.1145/2897824.2925877

Parameters:

  • model (newton.Model) – The model to simulate.

  • config (Config) – Solver configuration. See SolverImplicitMPM.Config.

  • temporary_store (fem.TemporaryStore | None) – Optional Warp FEM temporary store for reusing scratch allocations across steps.

  • verbose (bool | None) – If True, enable verbose solver output. If False, suppress details. If None, enable verbose output when wp.config.log_level is configured for debug logging.

  • enable_timers (bool) – Enable per-section wall-clock timings. This must be False when step() is recorded in an outer graph because section timers synchronize the device.

Raises:

ValueError – If an isolated multi-world model contains global MPM particles, particle world IDs outside [-1, model.world_count), or particles that are not stored contiguously by world.

class Config(max_iterations=250, tolerance=0.0001, solver='auto', warmstart_mode='auto', collider_velocity_mode='forward', voxel_size=0.1, grid_type='sparse', grid_padding=0, max_active_cell_count=-1, max_leaf_node_count=-1, max_lower_node_count=-1, max_upper_node_count=-1, transfer_scheme='apic', integration_scheme='pic', critical_fraction=0.0, air_drag=1.0, collider_normal_from_sdf_gradient=False, collider_basis='S2', strain_basis='P0', velocity_basis='Q1', separate_worlds=False)#

Bases: object

Configuration for SolverImplicitMPM.

Per-particle properties can be configured using custom attributes on the Model. See SolverImplicitMPM.register_custom_attributes() for details.

__init__(max_iterations=250, tolerance=0.0001, solver='auto', warmstart_mode='auto', collider_velocity_mode='forward', voxel_size=0.1, grid_type='sparse', grid_padding=0, max_active_cell_count=-1, max_leaf_node_count=-1, max_lower_node_count=-1, max_upper_node_count=-1, transfer_scheme='apic', integration_scheme='pic', critical_fraction=0.0, air_drag=1.0, collider_normal_from_sdf_gradient=False, collider_basis='S2', strain_basis='P0', velocity_basis='Q1', separate_worlds=False)#
air_drag: float = 1.0#

Numerical drag for the background air.

collider_basis: Literal['Q1', 'S2', 'pic', 'pic8', 'pic27'] | str = 'S2'#

Collider basis function. Defaults to "S2"; pass "Q1" to restore the previous trilinear collider basis. Common values are "Q1" (trilinear), "S2" (quadratic serendipity), or "pic", "pic8", "pic27" (particle-based with optional max points per cell). Any "picN" form with integer N is accepted.

collider_normal_from_sdf_gradient: bool = False#

Compute collider normals from sdf gradient rather than closest point

collider_velocity_mode: Literal['forward', 'backward'] = 'forward'#

Collider velocity computation mode. 'forward' uses the current velocity, 'backward' uses the previous timestep position.

critical_fraction: float = 0.0#

Fraction for particles under which the yield surface collapses.

grid_padding: int = 0#

Number of empty cells to add around particles when allocating the grid.

grid_type: Literal['sparse', 'dense', 'fixed'] = 'sparse'#

Type of grid to use.

A capacity-bounded "sparse" grid is rebuilt in place. Dense grids read dynamic bounds on the host.

integration_scheme: Literal['pic', 'gimp'] = 'pic'#

Integration scheme controlling shape-function support.

max_active_cell_count: int = -1#

Maximum number of active grid cells across all worlds.

A positive value reserves persistent sparse-grid capacity and bounds active subsets of dense and fixed grids. -1 retains per-step sparse-grid allocation and uses exact active counts elsewhere. Call check_sparse_grid_rebuild_status() after graph replay to detect sparse-grid overflow.

max_iterations: int = 250#

Maximum number of iterations for the rheology solver.

max_leaf_node_count: int = -1#

Maximum NanoVDB leaf-node count across all worlds.

This independently bounds leaf topology for a rebuildable sparse grid. -1 reserves one leaf per max_active_cell_count, which is the worst case for arbitrarily scattered active cells. Set an explicit value only when application-level spatial bounds guarantee a tighter limit. All node-capacity settings are validated at construction but used only for rebuildable sparse grids.

max_lower_node_count: int = -1#

Maximum NanoVDB lower internal-node count across all worlds.

-1 estimates the initial packed topology and reserves 16 times its lower-node count, capped by the resolved leaf-node capacity. An explicit value budgets spatial-spread headroom independently from active cells. Only used for rebuildable sparse grids.

max_upper_node_count: int = -1#

Maximum NanoVDB upper internal-node count across all worlds.

-1 estimates the initial packed topology and reserves 16 times its upper-node count, capped by the resolved lower-node capacity. Upper nodes cover regions 4096 voxels wide and are substantially larger than lower or leaf nodes, so applications with known spatial bounds can use this field to budget them explicitly. Only used for rebuildable sparse grids.

separate_worlds: bool = False#

Use independent FEM environments for each world in a multi-world model.

The default False retains the legacy shared-grid behavior. Set to True to isolate grid mass, momentum, stress, and collider response by world. Isolated multi-world models require every MPM particle to belong to a local world and particles to be stored contiguously by world. See Implicit MPM world isolation for how this mode interprets world assignment and collider ownership.

Experimental

Isolated multi-world MPM configuration and behavior may change without prior notice.

solver: Literal['auto', 'gs', 'gauss-seidel', 'gs-soa', 'gauss-seidel-soa', 'gs-batched', 'gauss-seidel-batched', 'jacobi', 'cg', 'cr', 'gmres'] | Sequence[Literal['auto', 'gs', 'gauss-seidel', 'gs-soa', 'gauss-seidel-soa', 'gs-batched', 'gauss-seidel-batched', 'jacobi', 'cg', 'cr', 'gmres']] = 'auto'#

Solver to use for the rheology solver. "auto" selects "gs" for Q1 velocity basis and "gs-batched" for higher-order bases (B2, B3). Accepted values: "auto", "gs" (or "gauss-seidel"), "gs-soa" (or "gauss-seidel-soa"), "gs-batched" (or "gauss-seidel-batched"), "jacobi", "cg", "cr", "gmres". Pass an ordered sequence to warmstart solvers left-to-right, e.g. ("cr", "gs") or ("cg", "jacobi", "gs").

strain_basis: Literal['P0', 'P1d', 'Q1', 'Q1d', 'pic', 'pic8', 'pic27'] | str = 'P0'#

Strain basis function. Common values are "P0", "P1d", "Q1", "Q1d", or particle-based "pic", "pic8", "pic27". Any "picN" form with integer N is accepted.

tolerance: float = 0.0001#

Tolerance for the rheology solver.

transfer_scheme: Literal['apic', 'pic'] = 'apic'#

Transfer scheme to use for particle-grid transfers.

velocity_basis: Literal['Q1', 'B2', 'B3'] = 'Q1'#

Velocity basis function. Common values are "Q1", "B2", or "B3".

voxel_size: float = 0.1#

Size of the grid voxels.

warmstart_mode: Literal['none', 'auto', 'particles', 'grid', 'smoothed'] = 'auto'#

Warmstart mode to use for the rheology solver.

"auto" uses particle-backed stress for rebuildable sparse grids and P1d/Q1d strain bases, and grid-backed stress otherwise. Grid-backed "grid" and "smoothed" modes are not supported for rebuildable sparse grids because their topology changes in place.

classmethod register_custom_attributes(builder)#

Register MPM-specific custom attributes in the ‘mpm’ namespace.

This method registers per-particle material parameters and state variables for the implicit MPM solver.

Attributes registered on Model (per-particle):

  • mpm:young_modulus: Young’s modulus in Pa

  • mpm:poisson_ratio: Poisson’s ratio for elasticity

  • mpm:damping: Elastic damping relaxation time in seconds

  • mpm:friction: Friction coefficient

  • mpm:yield_pressure: Yield pressure in Pa

  • mpm:tensile_yield_ratio: Tensile yield ratio

  • mpm:yield_stress: Deviatoric yield stress in Pa

  • mpm:hardening: Hardening factor for plasticity

  • mpm:hardening_rate: Hardening rate for plasticity

  • mpm:softening_rate: Softening rate for plasticity

  • mpm:dilatancy: Dilatancy factor for plasticity

  • mpm:viscosity: Viscosity for plasticity [Pa·s]

Attributes registered on State (per-particle):

  • mpm:particle_qd_grad: Velocity gradient for APIC transfer

  • mpm:particle_elastic_strain: Elastic deformation gradient

  • mpm:particle_Jp: Determinant of plastic deformation gradient

  • mpm:particle_stress: Cauchy stress tensor [Pa]

  • mpm:particle_transform: Overall deformation gradient for rendering

__init__(model, config, *, temporary_store=None, verbose=None, enable_timers=False)#
check_sparse_grid_rebuild_status()#

Raise if a rebuildable sparse grid exceeded its reserved capacity.

Rebuildable sparse-grid failures accumulate until a valid reset(). Call this method outside graph capture after replay; the check synchronizes the solver device. It is a no-op when the solver has no asynchronous sparse-grid status buffer.

Raises:

RuntimeError – If rebuild status is inspected during graph capture, or if a rebuild reported a capacity or topology failure.

collect_collider_impulses(state)#

Collect current collider impulses and their application positions.

Returns a tuple of 3 arrays:

  • Impulse values in world units.

  • Collider positions in world units.

  • Collider id, that can be mapped back to the model’s body ids using the collider_body_index property.

coupling_eval_gravity_acceleration(out_body_acceleration, out_particle_acceleration)#

Evaluate gravity acceleration applied internally by the MPM solver.

coupling_harvest_proxy_particle_forces(particle_local_to_proxy_global, out_particle_f, *, particle_qd_before, state, state_out, contacts, dt)#

Convert MPM proxy momentum changes and collider impulses to forces.

coupling_harvest_proxy_wrenches(body_local_to_proxy_global, out_body_f, *, body_qd_before, state, state_out, contacts, dt)#

Convert MPM collider grid impulses to proxy-body wrenches.

coupling_notify_input_state_update(state, flags, *, iteration_restart=False, dt=0.0)#

Synchronize deformable collider meshes after particle input-state updates.

coupling_rewind_proxy_body(body_local_to_proxy_global, state, coupling_forces, body_gravity_acceleration, dt)#

Remove lagged velocity-level proxy wrenches from collider velocities.

coupling_rewind_proxy_particle(particle_local_to_proxy_global, state, coupling_forces, particle_gravity_acceleration, dt)#

Remove lagged velocity-level proxy forces from proxy particle velocities.

notify_model_changed(flags)#
project_outside(state_in, state_out, dt, gap=None)#

Project particles outside of colliders, and adjust their velocity and velocity gradients

Parameters:

  • state_in (State) – The input state.

  • state_out (State) – The output state. Only particle_q, particle_qd, and particle_qd_grad are written.

  • dt (float) – The time step, for extrapolating the collider end-of-step positions from its current position and velocity.

  • gap (float | None) – Maximum distance for closest-point queries. If None, the default is the voxel size times sqrt(3).

reset(state, world_mask=None, flags=None)#

Reset implicit MPM history for all or selected worlds.

Particle history is reset when flags includes either PARTICLE_Q or PARTICLE_QD. If flags is None, all particle history is reset. Particle- and grid-backed warm starts are cleared for selected worlds. Grid-backed warm starts cannot be selectively cleared on a shared multi-world grid; use Config.separate_worlds or a full reset. A full reset clears every warm-start field. Sparse-grid rebuild status is always cleared at a valid reset boundary, and the previous-collider-pose cache is refreshed from state. When present, the final mask entry selects global particle-backed history and collider poses whose world index is -1.

Parameters:

  • state (State) – Simulation state whose MPM history is modified in place.

  • world_mask (wp.array[wp.bool] | None) –

    Optional one-dimensional Warp boolean mask on the model device with shape (model.world_count + 1,). The final entry selects global objects whose world index is -1. If None, reset all worlds and global objects.

    Deprecated since version 1.5: Passing a mask with shape (world_count,) is deprecated. Use shape (world_count + 1,) with a final False entry to select local worlds only.

    Experimental

    Selective per-world MPM reset behavior may change without prior notice.

  • flags (StateFlags | int | None) – Optional state bitmask. If None, reset all particle history.

sample_render_grains(state, grains_per_particle)#

Generate per-particle point samples used for high-resolution rendering.

Parameters:

  • state (State) – Current Newton state providing particle positions.

  • grains_per_particle (int) – Number of grains to sample per particle.

Returns:

A wp.array with shape (num_particles, grains_per_particle) of type wp.vec3 containing grain positions.

Return type:

array

setup_collider(*, collider_meshes=None, collider_body_ids=None, collider_margins=None, collider_friction=None, collider_adhesion=None, collider_projection_threshold=None, collider_particle_ids=None, model=None, body_com=None, body_mass=None, body_inv_inertia=None, body_q=None, collider_world_ids=None)#

Configure collider geometry and material properties.

By default, collisions are set up against all shapes in the model with newton.ShapeFlags.COLLIDE_PARTICLES. Use this method to customize collider sources, materials, or to read colliders from a different model.

Parameters:

  • collider_meshes (list[Mesh] | None) – Warp triangular meshes used as colliders.

  • collider_body_ids (list[int] | None) – For dynamic colliders, per-mesh body ids.

  • collider_margins (list[float] | None) – Per-mesh signed distance offsets (m).

  • collider_friction (list[float] | None) – Per-mesh Coulomb friction coefficients.

  • collider_adhesion (list[float] | None) – Per-mesh adhesion (Pa).

  • collider_projection_threshold (list[float] | None) – Per-mesh projection threshold (m).

  • collider_particle_ids (list[list[int] | wp.array[wp.int32] | None] | None) – For deformable mesh colliders, model particle ids corresponding to each mesh vertex.

  • model (Model | None) – The model to read collider properties from. Default to solver’s model.

  • body_com (array | None) – For dynamic colliders, per-body center of mass.

  • body_mass (array | None) – For dynamic colliders, per-body effective mass. When omitted, bodies flagged with newton.BodyFlags.KINEMATIC have zero effective mass. An explicit array is authoritative.

  • body_inv_inertia (array | None) – For dynamic colliders, per-body inverse inertia.

  • body_q (array | None) – For dynamic colliders, per-body initial transform.

  • collider_world_ids (list[int] | None) –

    Per-collider Newton world IDs. Custom meshes default to global (-1). In isolated mode, body-backed colliders infer their body’s world and require any supplied ID to match. Shared-mode default discovery globalizes colliders. IDs must be -1 or in [0, model.world_count).

    Experimental

    Per-world MPM collider filtering may change without prior notice.

Raises:

ValueError – If collider-aligned inputs have different lengths, a world ID is invalid, an isolated external collider model has a different world_count, or an isolated global body-backed collider is dynamic. Replicate a global dynamic collider per world, make it static or kinematic, or disable Config.separate_worlds.

step(state_in, state_out, control, contacts, dt)#

Advance the simulation by one time step.

Transfers particle data to the grid, solves the implicit rheology system, and transfers the result back to update particle positions, velocities, and stress.

Parameters:

  • state_in (State) – Input state at the start of the step.

  • state_out (State) – Output state written with updated particle data. May be the same object as state_in for in-place stepping.

  • control (Control) – Control input (unused; material parameters come from the model).

  • contacts (Contacts) – Contact information (unused; collisions are handled internally).

  • dt (float) – Time step duration [s].

update_particle_frames(state_prev, state, dt, min_stretch=0.25, max_stretch=2.0)#

Update per-particle deformation frames for rendering and projection.

Integrates the particle deformation gradient using the velocity gradient and clamps its principal stretches to the provided bounds for robustness.

update_render_grains(state_prev, state, grains, dt)#

Advect grain samples with the grid velocity and keep them inside the deformed particle.

Parameters:

  • state_prev (State) – Previous state (t_n).

  • state (State) – Current state (t_{n+1}).

  • grains (array) – 2D array of grain positions per particle to be updated in place. See sample_render_grains.

  • dt (float) – Time step duration.

property collider_body_index: array#

Array mapping collider indices to body indices.

Returns:

Per-collider body index array. Value is -1 for colliders that are not bodies.

property voxel_size: float#

Grid voxel size used by the solver.