Propagation options

61 · Solver and runtime

A twelve-link chain hung from one end swings down, lands and folds onto itself, so it has ground contacts and contacts between its own links. With the colored sweep and the propagation response, contacts between two links of one articulation stay dense rows by default; propagation_same_articulation_rows moves them to propagation rows. propagation_cached_response applies cached body responses instead of walking the tree each iteration, and propagation_color_order picks the order of the colored batches.

drag to orbit · scroll to zoom

0.00 s

Cost

Time per step, CUDA graph replay, median of 7 interleaved rounds on one NVIDIA RTX PRO 6000 Blackwell Server Edition; ratios are against the first row. Solver allocation is the GPU memory the solver allocates at construction, before warm-up and graph capture. Scene: the folding twelve-link chain in each world.

Option	1 world	64 worlds	1024 worlds	solver allocation, 1024 worlds	result in this tutorial
serial + immediate	0.98 ms	1.14 ms	1.65 ms	211 MB	reference
colored + propagation (defaults)	1.02 ms 1.04×	1.17 ms 1.03×	1.98 ms 1.21×	374 MB	96 dense rows
propagation_same_articulation_rows=True	2.09 ms 2.13×	5.90 ms 5.19×	12.87 ms 7.82×	374 MB	0 dense rows
propagation_cached_response=False	1.12 ms 1.14×	1.25 ms 1.11×	1.96 ms 1.19×	362 MB	same result
propagation_color_order="rotating"	1.02 ms 1.04×	1.17 ms 1.03×	1.99 ms 1.21×	374 MB	same result
propagation_constraint_max=64, free_body_constraint_max=8	1.02 ms 1.03×	1.17 ms 1.03×	1.94 ms 1.18×	254 MB	fits, less memory

What to look at

The chain lands the same way in every run; same_articulation_rows turns the 96 dense self-contact rows into propagation rows.
The 96 dense rows do not fit the default articulation_constraint_max of 32. As propagation rows they do: the last run uses the default capacity and lands like the others. Separate runs that kept them dense at 32 dropped rows, and some went non-finite.
The divergence comes from which contacts keep their rows: by default that set changes every step, and contacts that penetrated while dropped are pushed out when they return, which adds energy. In 6 runs each it diverged in 1 with this colored propagation solve and in 4 with the default immediate response; with deterministic solver and pipeline order it repeated exactly and stayed finite.

Key parameter

propagation_same_articulation_rows / propagation_cached_response / propagation_color_order

Options of articulated_contact_response="propagation"; propagation_constraint_max sizes its rows.

Try changing

Set propagation_constraint_max=64 and shrink free_body_constraint_max.
Try propagation_color_order="longest_first".

Concept

Propagation rows

What it is

The propagation response solves each contact as a small row on the touched bodies and pushes its impulse through the tree each iteration. Self-contacts of one articulation normally stay dense rows with their own capacity. Cached responses replace the tree walk with small precomputed matrices, and the color order changes which contacts Gauss-Seidel favours.

In this scene

Moving self-contacts to propagation rows takes 0 dense rows instead of 96, so the default dense capacity of 32 is enough.

Considerations

Cost: same_articulation_rows was 2× slower in one world and 6× at 1024 worlds (12.9 vs 2.0 ms); cached_response=False was 10% slower in one world; setting propagation_constraint_max=64 together with free_body_constraint_max=8 reduced solver-owned memory from 374 to 254 MB at 1024 worlds.
Interactions: rotating color order changed nothing here; it targets symmetric grasps on sorted contact streams.
Use it when: same_articulation_rows when dense capacity overflows; size propagation_constraint_max from row_watermark to save memory.

The script

cd tutorials && python 61_propagation_options.py prints each variant's final metric; --record out.json also writes the replay. Every tutorial imports replay.py from the same folder.

61_propagation_options.py
replay.py

# SPDX-FileCopyrightText: Copyright (c) 2026 Skild AI and the Steel contributors
# SPDX-License-Identifier: Apache-2.0

"""Propagation options: same-articulation rows, cached responses and the colored contact order."""

import newton
import numpy as np
import warp as wp
from newton.solvers import SolverSteel

from replay import Recorder

DT = 1.0 / 240.0
SUBSTEPS = 8
LINKS = 12
BASE = {"solve_mode": "colored", "articulated_contact_response": "propagation"}
OPTIONS = {
    "colored + propagation (defaults)": {},
    "propagation_same_articulation_rows=True": {"propagation_same_articulation_rows": True},
    "propagation_cached_response=False": {"propagation_cached_response": False},
    'propagation_color_order="rotating"': {"propagation_color_order": "rotating"},
    # the 96 dense self-contact rows would not fit the default 32; as propagation rows they do
    "same_articulation_rows, articulation_constraint_max=32": {
        "articulation_constraint_max": 32,
        "propagation_same_articulation_rows": True,
    },
}


def tf(x, y, z):
    return wp.transform((x, y, z), wp.quat_identity())


def build():
    """A twelve-link chain hung from one end; it swings down and folds onto itself."""
    builder = newton.ModelBuilder()
    builder.add_ground_plane()
    parent, joints = -1, []
    for i in range(LINKS):
        link = builder.add_link(xform=tf(0.08 * i + 0.04, 0.0, 0.5))
        rot = wp.quat_from_axis_angle(wp.vec3(0.0, 1.0, 0.0), np.pi / 2)
        builder.add_shape_capsule(
            link, xform=wp.transform((0.0, 0.0, 0.0), rot), radius=0.025, half_height=0.02
        )
        anchor = tf(0.0, 0.0, 0.5) if parent < 0 else tf(0.04, 0.0, 0.0)
        axis = (0.0, 1.0, 0.0) if i % 3 else (0.0, 0.3, 1.0)
        joints.append(
            builder.add_joint_revolute(
                parent,
                link,
                axis=axis,
                parent_xform=anchor,
                child_xform=tf(-0.04, 0.0, 0.0),
                armature=0.002,
                damping=0.002,
            )
        )
        parent = link
    builder.add_articulation(joints)
    return builder.finalize()


def simulate(rec, label):
    model = build()
    options = dict(
        BASE, articulation_constraint_max=256, free_body_constraint_max=256, row_watermark=True
    )
    solver = SolverSteel(model, **dict(options, **OPTIONS[label]))
    pipeline = newton.CollisionPipeline(model)
    contacts = pipeline.contacts()
    state_0, state_1, control = model.state(), model.state(), model.control()

    run = rec.variant(label, model)
    for _ in range(rec.frames(4.0)):
        marks = solver.constraint_row_watermarks()
        run.frame(
            state_0,
            dense_rows=marks["articulation_raw_high_water"],
            propagation_rows=marks["propagation_raw_high_water"],
            tip_z=state_0.body_q.numpy()[-1, 2],
        )
        for _ in range(SUBSTEPS):
            pipeline.collide(state_0, contacts)
            solver.step(state_0, state_1, control, contacts, DT)
            state_0, state_1 = state_1, state_0
    run.check(solver)


if __name__ == "__main__":
    rec = Recorder()
    for label in OPTIONS:
        simulate(rec, label)
    rec.save()

# SPDX-FileCopyrightText: Copyright (c) 2026 Skild AI and the Steel contributors
# SPDX-License-Identifier: Apache-2.0

"""Record Steel tutorial runs as replay JSON for the tutorials page, or print a summary.

Every tutorial imports this module. ``python NN_name.py`` runs the variants and prints
the final metric of each; ``--record out.json`` also writes body poses, shapes and metric
series that the page replays in 3D.
"""

import argparse
import base64
import json
import os
import subprocess
import sys
import time

import newton
import numpy as np
import warp as wp

from source_hash import source_sha256

_RENDERED_TYPES = {
    newton.GeoType.PLANE,
    newton.GeoType.SPHERE,
    newton.GeoType.CAPSULE,
    newton.GeoType.ELLIPSOID,
    newton.GeoType.CYLINDER,
    newton.GeoType.BOX,
    newton.GeoType.MESH,
    newton.GeoType.CONE,
    newton.GeoType.CONVEX_MESH,
}


class Recorder:
    """Collect the variants of one tutorial; ``fps`` is the replay frame rate."""

    def __init__(self, fps: int = 30):
        parser = argparse.ArgumentParser()
        parser.add_argument(
            "--record", metavar="OUT.json", help="write a replay file for the tutorials page"
        )
        parser.add_argument("--device", default=None, help="Warp device, for example cuda:0")
        self.args = parser.parse_args()
        if self.args.device:
            wp.set_device(self.args.device)
        self.fps = fps
        self.script = os.path.abspath(sys.argv[0])
        self.variants: list[Variant] = []

    def frames(self, seconds: float) -> int:
        """Number of replay frames in ``seconds``, including the first."""
        return int(round(seconds * self.fps)) + 1

    def variant(
        self,
        label: str,
        model: newton.Model,
        *,
        ghost: bool = False,
        slot: int | None = None,
        scene: str | None = None,
    ) -> "Variant":
        """Start one run; ``ghost`` overlays run ``slot``, ``scene`` groups runs for charts."""
        variant = Variant(label, model, self.fps, ghost, slot, scene)
        self.variants.append(variant)
        return variant

    def save(self) -> None:
        """Print a summary and, with ``--record``, write the replay file."""
        for variant in self.variants:
            last = variant.metrics[-1] if variant.metrics else {}
            flag = "  (rows dropped)" if variant.overflow else ""
            print(f"{variant.label:<48} {json.dumps(last)}  {variant.wall:.2f} s{flag}")
        if not self.args.record:
            return
        source_sha = source_sha256(self.script)
        doc = {
            "script": os.path.basename(self.script),
            "script_sha256": source_sha,
            "newton_commit": _newton_commit(),
            "newton_version": newton.__version__,
            "warp_version": wp.config.version,
            "gpu": wp.get_device().name,
            "fps": self.fps,
            "variants": [v.to_json() for v in self.variants],
        }
        with open(self.args.record, "w") as f:
            json.dump(doc, f, separators=(",", ":"))
        print(f"wrote {self.args.record}")


class Variant:
    """One recorded run: a pose per frame for every body, plus named metric values."""

    def __init__(self, label, model, fps, ghost, slot, scene):
        self.label, self.model, self.fps, self.ghost, self.slot = label, model, fps, ghost, slot
        self.scene = scene
        self.poses: list[np.ndarray] = []
        self.metrics: list[dict] = []
        self.awake: list[list[int]] = []
        self.particles: list[np.ndarray] = []
        self.overflow = False
        self.start = time.perf_counter()
        self.wall = 0.0

    def frame(self, state: newton.State, awake=None, particles=None, **metrics) -> None:
        """Record ``state.body_q``, optional awake flags, particle positions and metric values."""
        self.poses.append(state.body_q.numpy().astype(np.float32).copy())
        if particles is not None:
            self.particles.append(particles.numpy().astype(np.float32).copy())
        self.metrics.append({k: _round(v) for k, v in metrics.items()})
        if awake is not None:
            self.awake.append([int(x) for x in awake.numpy()])
        self.wall = time.perf_counter() - self.start

    def check(self, solver) -> None:
        """Note whether the solver dropped rows; call once after the run."""
        try:
            solver.check_constraint_capacity()
        except RuntimeError:
            self.overflow = True

    def to_json(self) -> dict:
        poses = (
            np.stack(self.poses)
            if self.poses
            else np.zeros((0, self.model.body_count, 7), np.float32)
        )
        return {
            "label": self.label,
            "fps": self.fps,
            "ghost": self.ghost,
            "slot": self.slot,
            "scene": self.scene,
            "bodies": int(self.model.body_count),
            "frames": int(poses.shape[0]),
            "poses": base64.b64encode(poses.astype(np.float32).tobytes()).decode(),
            "shapes": _shape_table(self.model),
            "metric": self.metrics,
            "awake": self.awake or None,
            "particles": self._particles_json(),
            "wall_s": round(self.wall, 3),
            "overflow": self.overflow,
        }

    def _particles_json(self):
        if not self.particles:
            return None
        frames = np.stack(self.particles)
        tris = self.model.tri_indices
        radius = self.model.particle_radius.numpy()
        return {
            "count": int(frames.shape[1]),
            "frames": base64.b64encode(frames.astype(np.float32).tobytes()).decode(),
            "radius": round(float(radius.mean()), 5) if radius.size else 0.0,
            "tris": [int(v) for v in tris.numpy().reshape(-1)]
            if tris is not None and tris.shape[0]
            else None,
        }


def _shape_table(model):
    types = model.shape_type.numpy()
    scales = model.shape_scale.numpy()
    xforms = model.shape_transform.numpy()
    bodies = model.shape_body.numpy()
    flags = model.shape_flags.numpy()
    shapes = []
    for i, t in enumerate(types):
        if int(t) not in _RENDERED_TYPES or not flags[i] & int(newton.ShapeFlags.VISIBLE):
            continue
        shape = {
            "b": int(bodies[i]),
            "t": int(t),
            "s": [round(float(v), 5) for v in scales[i]],
            "x": [round(float(v), 6) for v in xforms[i]],
        }
        if int(t) in (newton.GeoType.MESH, newton.GeoType.CONVEX_MESH):
            mesh = model.shape_source[i]
            shape["v"] = [round(float(v), 5) for v in np.asarray(mesh.vertices).reshape(-1)]
            shape["i"] = [int(v) for v in np.asarray(mesh.indices).reshape(-1)]
        shapes.append(shape)
    return shapes


def _round(value):
    value = float(value)
    return round(value, 5) if np.isfinite(value) else None


def _newton_commit():
    if os.environ.get("NEWTON_COMMIT"):
        return os.environ["NEWTON_COMMIT"]
    try:
        root = os.path.dirname(os.path.dirname(newton.__file__))
        return subprocess.check_output(["git", "-C", root, "rev-parse", "HEAD"], text=True).strip()
    except (OSError, subprocess.CalledProcessError):
        return None

Recorded with SolverSteel on Newton 1.9.0.dev0, Warp 1.19.0, NVIDIA RTX PRO 6000 Blackwell Server Edition. Replays read the state back every frame, so they are not a speed benchmark. Script · page text

Steel docs and tutorials are licensed under Apache-2.0. Third-party notices. "Skild", "Skild AI" and the Skild logo are trademarks of Skild AI and are not licensed under Apache-2.0.