Coupling Steel with XPBD cloth

64 · Coupling

Steel simulates rigid bodies and articulations but no particles, so it cannot simulate cloth on its own. SolverCoupledProxy runs two solvers on one model: here Steel owns two boxes, a six-link chain and a stack of six boxes, and XPBD owns a cloth fixed along two edges. XPBD sees each Steel body as a proxy body with Steel's effective mass, and the contact impulses it computes are fed back to Steel. The comparison run lets XPBD simulate everything.

drag to orbit · scroll to zoom

0.00 s

What to look at

Both runs drop the boxes and the chain into the cloth, which sags about 30 cm and holds them.
Steel's chain joints stay exact (largest anchor gap 0.00006 mm); XPBD's joints open by up to 0.08 mm.
The stack beside the cloth rests within 0.1 mm in both runs.

Key parameter

SolverCoupledProxy(model, entries, coupling)

Each Entry names a solver factory and the bodies, joints or particles it owns. A Proxy maps source bodies into the destination solver; Config.iterations sets the proxy relaxation passes.

Try changing

Raise Config(iterations=2) and compare the cloth sag.
Swap SolverXPBD for SolverVBD in the second entry.
Give Steel a driven arm that presses into the cloth.

Concept

Proxy coupling

What it is

Two solvers step the same model, each owning part of it. Before the destination solver steps, the source bodies are copied in as proxies with the source's effective mass and inertia. The contact impulses the destination applies to those proxies are harvested and fed back to the source as forces, lagged by one solve; Config.iterations repeats the exchange within a step.

In this scene

Steel owns all rigid bodies and joints, XPBD the 289 cloth particles; XPBD alone also handled these rigid bodies well.

Considerations

Cost: Two solves per step plus the exchange: 14.6 vs 9.7 ms per eager step with XPBD alone (one world, printed by the script).
Accuracy: Feedback lags by one solve, so stiff contact between the two sides can need smaller steps or more proxy iterations.
Interactions: At most two entries per coupled solver; set Steel options in its factory and check it via coupled.solver("steel").
Use it when: A Steel robot or articulated scene needs cloth or soft bodies; for a scene that is mostly deformables with a few free boxes, one deformable solver is simpler and faster.

The script

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

64_coupling.py
replay.py

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

"""Coupling: Steel simulates a chain, boxes and a stack; XPBD the cloth they land on."""

import time

import newton
import numpy as np
import warp as wp
from newton.solvers import SolverSteel, SolverXPBD
from newton.solvers.experimental.coupled import SolverCoupledProxy

from replay import Recorder

DT = 1.0 / 960.0
SUBSTEPS = 32  # 30 Hz frames
HX = 0.06  # half length of a chain link


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


def build():
    builder = newton.ModelBuilder()
    builder.add_ground_plane()
    cfg = newton.ModelBuilder.ShapeConfig(density=100.0)
    for x, z in ((0.2, 0.9), (-0.15, 1.1)):
        box = builder.add_body(xform=tf(x, 0.0, z))
        builder.add_shape_box(box, hx=0.08, hy=0.08, hz=0.08, cfg=cfg)
    # a free-floating six-link chain: its joints are tree joints in Steel, constraints in XPBD
    link = builder.add_link(xform=tf(-0.3, 0.0, 1.4))
    builder.add_shape_box(link, hx=HX, hy=0.03, hz=0.03, cfg=cfg)
    joints = [builder.add_joint_free(link)]
    for _ in range(5):
        child = builder.add_link()
        builder.add_shape_box(child, hx=HX, hy=0.03, hz=0.03, cfg=cfg)
        joints.append(
            builder.add_joint_revolute(
                link,
                child,
                axis=newton.Axis.Y,
                target_kd=0.5,
                parent_xform=tf(HX, 0.0, 0.0),
                child_xform=tf(-HX, 0.0, 0.0),
            )
        )
        link = child
    builder.add_articulation(joints)
    # a stack of six 6 cm boxes on the floor beside the cloth
    for i in range(6):
        box = builder.add_body(xform=tf(0.0, 0.8, 0.03 + 0.0602 * i))
        builder.add_shape_box(box, hx=0.03, hy=0.03, hz=0.03, cfg=cfg)
    rigid_bodies = list(range(builder.body_count))
    # a 1 m cloth fixed along two edges, 0.5 m above the ground
    builder.add_cloth_grid(
        pos=wp.vec3(-0.5, -0.5, 0.5),
        rot=wp.quat_identity(),
        vel=wp.vec3(0.0),
        dim_x=16,
        dim_y=16,
        cell_x=1.0 / 16,
        cell_y=1.0 / 16,
        mass=0.1,
        fix_left=True,
        fix_right=True,
        add_springs=True,
        spring_ke=1.0e3,
        spring_kd=1.0,
        tri_ke=1.0e3,
        tri_ka=1.0e3,
        tri_kd=1.0e-1,
        edge_ke=1.0e2,
        edge_kd=1.0e-1,
        particle_radius=0.01,
    )
    model = builder.finalize()
    model.soft_contact_ke = 1.0e2
    model.soft_contact_kd = 1.0
    return model, rigid_bodies


def make_solver(model, rigid_bodies, coupled):
    if not coupled:
        return SolverXPBD(model, iterations=40)
    return SolverCoupledProxy(
        model=model,
        entries=[
            SolverCoupledProxy.Entry(
                name="steel",
                solver=lambda view: SolverSteel(view, articulation_constraint_max=256),
                bodies=rigid_bodies,
                joints=list(range(model.joint_count)),
            ),
            SolverCoupledProxy.Entry(
                name="xpbd",
                solver=lambda view: SolverXPBD(view, iterations=40),
                particles=list(range(model.particle_count)),
            ),
        ],
        # XPBD sees each Steel body as a proxy with Steel's effective mass; impulses flow back
        coupling=SolverCoupledProxy.Config(
            proxies=[
                SolverCoupledProxy.Proxy(source="steel", destination="xpbd", bodies=rigid_bodies)
            ]
        ),
    )


def joint_gap_mm(body_q):
    """Largest distance between the two anchors of a chain joint [mm]."""
    gaps = []
    for parent in range(2, 7):
        a, b = body_q[parent], body_q[parent + 1]
        pa = wp.transform_point(
            wp.transform(wp.vec3(*a[:3]), wp.quat(*a[3:])), wp.vec3(HX, 0.0, 0.0)
        )
        pb = wp.transform_point(
            wp.transform(wp.vec3(*b[:3]), wp.quat(*b[3:])), wp.vec3(-HX, 0.0, 0.0)
        )
        gaps.append(np.linalg.norm(np.array(pa) - np.array(pb)))
    return 1000.0 * max(gaps)


def simulate(rec, label, coupled):
    model, rigid_bodies = build()
    solver = make_solver(model, rigid_bodies, coupled)
    pipeline = newton.CollisionPipeline(model)
    contacts = pipeline.contacts()
    state_0, state_1, control = model.state(), model.state(), model.control()
    newton.eval_fk(model, model.joint_q, model.joint_qd, state_0)

    run = rec.variant(label, model)
    elapsed = 0.0
    for frame in range(rec.frames(4.0)):
        q = state_0.body_q.numpy()
        cloth_z = state_0.particle_q.numpy()[:, 2].min()
        stack_sag_mm = 1000.0 * (0.03 + 0.06 * 5 - q[-1, 2])
        run.frame(
            state_0,
            particles=state_0.particle_q,
            joint_gap_mm=joint_gap_mm(q),
            cloth_low_z=cloth_z,
            stack_sag_mm=stack_sag_mm,
        )
        wp.synchronize()
        start = time.perf_counter()
        for _ in range(SUBSTEPS):
            state_0.clear_forces()
            pipeline.collide(state_0, contacts)
            solver.step(state_0, state_1, control, contacts, DT)
            state_0, state_1 = state_1, state_0
        wp.synchronize()
        if frame > 0:  # the first frame includes kernel compilation
            elapsed += time.perf_counter() - start
    if coupled:
        run.check(solver.solver("steel"))
    print(f"{label}: {1000.0 * elapsed / ((rec.frames(4.0) - 1) * SUBSTEPS):.3f} ms per eager step")


if __name__ == "__main__":
    rec = Recorder()
    simulate(rec, "XPBD for everything", False)
    simulate(rec, "Steel rigid + XPBD cloth, coupled", True)
    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.