FeatherPGS grasp stability

The contact model is sound. Its defaults cannot express a real contact set, and one row family cannot express a squeeze.

Given enough constraint rows and enough iterations, FeatherPGS holds a shaken object as tightly as SolverMuJoCo and tolerates a far larger time step before the grasp itself fails — with no contact tuning at all. Getting there means changing defaults that silently discard most of the scene's contacts, and avoiding two options that make things worse.

Start here

This configuration produced every measurement on this page. Two of the three settings are structural — they decide whether the scene's contacts can be represented at all — and the third is required for CUDA graph capture.

solver = newton.solvers.SolverFeatherPGS(
    model,
    pgs_mode="matrix_free",       # "split" cannot hold this scene's rows cheaply
    dense_max_constraints=256,    # scene peak is 201; the default 32 drops most of it
    mf_max_constraints=4096,
    pgs_iterations=12,            # 20 halves the creep at proportional cost
    double_buffer=False,          # required for CUDA graph capture
    # friction_mode="bisection_desaxce",   # 0.72 mm -> 0.40 mm, matrix-free only
    # pgs_beta, pgs_cfm, compliance and friction settings left at their defaults
)

Then supply the arm's gravity torque through control.joint_f, because the solver will not. And check what your scene actually needs:

solver = newton.solvers.SolverFeatherPGS(model, row_watermark=True, ...)
# ... run the scene ...
print(solver.constraint_row_watermarks())
# {'dense_high_water': 201, 'mf_high_water': 12, 'contact_high_water': 71, ...}

One prerequisite before any of this is measurable

A position-controlled arm needs a gravity term, and no Newton solver computes one yet — it is a platform gap being closed across all of them rather than a FeatherPGS defect. It matters here only because without it the test arm never reaches its target pose, so nothing downstream can be measured at all. Every arm-scene number on this page supplies it as a feed-forward torque through control.joint_f, leaving every gain and solver setting untouched.

Why these and not others

With the arm in position and the rows in place, the solver options can be judged on their own. Only two move the number much.

The two parameters people usually reach for do nothing useful. The grasp survives every value of the position-correction factor from zero to 0.9 and every value of the regularisation term across four orders of magnitude — but the residual creep under the position-correction factor bounces between 0.33 mm and 4.50 mm with no pattern, so there is nothing to tune toward.

The contact-row budget is 32. The arm scene needs 201, and going over is silent.

Every contact between the arm and anything else becomes one normal row plus two friction rows in a fixed-size per-world buffer. Drive rows and joint-limit rows are allocated from the same buffer first, so contacts are the first thing sacrificed when it fills.

The solver's own high-water telemetry gives the real requirement over a full grasp and shake: 201 dense rows at peak, from up to 71 simultaneous contacts. The constructor default is dense_max_constraints=32, which is ten contacts. The gap is not a matter of degree — below 96 rows the gripper cannot pick the cube up at all.

What overflow does

Slots are claimed with an atomic counter. A contact that does not fit rolls the counter back and marks itself dead:

    else:
        # Dense path
        slot = wp.atomic_add(world_slot_counter, world, slots_needed)
        if slot + slots_needed > max_constraints:
            # Roll back the counter so finalize sees only filled slots
            wp.atomic_add(world_slot_counter, world, -slots_needed)
            contact_slot[c] = -1
            contact_path[c] = -1
            return

• The dropped contact contributes nothing at all — no normal row, no friction, no force. A finger simply stops touching the cube for that step.
• There is no warning. The entire 15,760-line solver contains no warnings.warn and no logging. The only way to find out is to opt into row_watermark=True.
• Which contacts survive is decided by an atomic race, which also fixes the order the solve sweeps the rows in. Nothing orders them by penetration depth or by pair, and no contact keeps a stable row index between steps — which is why the dense path has no identity-matched warm start.

And the budget is expensive to raise on the default solve mode

The default pgs_kernel="tiled_row" stages the whole constraint matrix in static shared memory:

// Packed LOWER triangle of C in row-major (i*(i+1)/2 + j), j<=i
__shared__ float s_Ctri[TILE_TRI];      // TILE_TRI = M*(M+1)/2

__shared__ float s_lam[TILE_M];
__shared__ float s_rhs[TILE_M];
__shared__ float s_diag[TILE_M];
__shared__ int   s_rtype[TILE_M];
__shared__ int   s_parent[TILE_M];
__shared__ float s_mu[TILE_M];

That is 2M² + 26M bytes — quadratic in the row budget. CUDA's default per-block limit for a static shared-memory allocation is 48 KiB, and the opt-in above that applies only to dynamic shared memory, so on a device holding to the default the budget caps at 150 rows. Some devices allow more: on the RTX PRO 6000 Blackwell used for these measurements, 201 rows loads and 256 does not, putting its limit between 84 and 134 KiB.

Nothing checks any of this at construction. The solver already has a shared-memory estimator with a helpful error message, but it is wired into a different kernel family; the PGS kernel is built straight from dense_max_constraints with no device query, and the failure arrives as a CUDA module-load error inside the first step.

The arm scene on the constructor's own defaults

On its own defaults the solver never picks the cube up, entirely because of the 32-row budget, and it does so quietly — the arm runs the full scripted sequence over an object that never moves. Raising the split path's budget until it fits costs about three times the wall-clock time, because a 201-row tile is 84 KiB of shared memory driven by a single warp. Matrix-free mode is not an optimisation here; it is the only configuration that is both correct and fast.

A hard contact model cannot express a squeeze between two immovable surfaces.

Strip the scene down to two pads and one object and this becomes visible in isolation. It matters because it is the shape most people reach for when writing a minimal grasp test, and it has no answer.

Static pads: the object slides, and the rate is arbitrary

A 40 mm cube, 0.5 mm of overlap per side, gravity down, an oscillating force shaking it sideways for ten seconds. It does not fall — it never reaches a frame with zero contacts — but it does not stay put either:

Why there is no answer to find

Push the object equally hard from both sides and it does not move, whether you push with one newton or a thousand. Nothing in its motion distinguishes them, so nothing in the solve pins the pair of forces down: both contacts only say "do not penetrate", and neither has a compliance that would let a penetration depth stand in for a force. A compliant contact model does not have this problem, which is why SolverMuJoCo is untroubled by the same scene.

What FeatherPGS produces instead comes from the position-correction term, which asks for the penetration to be removed within one step — from both sides at once, in a sweep that visits one row and then the other. Three measurements pin that down.

First, the answer moves with the iteration count, on fixed geometry:

This is also why the collision representation matters so much here and so little once the pads are driven. The box presents eight contact points and the mesh eleven, so the two land on different values of an undetermined force, and the friction capacity that follows differs with them. Nothing about triangles is at fault; each representation simply converges to a different accident.

Second, squeezing harder makes it worse, and past a threshold the object is ejected:

Third, take the shake away, take gravity away, and simply place the object between the pads. Nothing in the scene is trying to move it. It moves anyway:

Turning the position-correction factor down settles it completely:

At pgs_beta = 0 the object does not move at all — the drift is exactly zero, not small. Every millimetre it travels comes from the position-correction term, and it grows far faster than linearly in that term because the push feeds back into the next step's penetration. The two opposed rows would cancel if solved together; they are not, because the sweep visits one and then the other.

The usual remedy is a compliance term, which turns the penetration back into a force. FeatherPGS has one — dense_contact_compliance — but it does not reach here. Sweeping it from zero to 1e-4 over the ejecting coupon produces bit-identical results, because contacts between a free object and world-fixed pads land on the free-body row family, whose diagonal is seeded with pgs_cfm and nothing else. A gripper's contacts land on a different family that does have the term.

It does not repeat, either

The box is bit-identical under every combination of double buffering and parallel streams, and bit-identical when simply run again. The mesh is not: two runs with nothing changed at all agree until frame 573 and then part company, ending 1.9 mm apart. Contact rows are claimed with an atomic counter, so the sweep order is whatever the GPU happened to schedule. On a well-posed problem that is a rounding-level detail. Here it decides the answer.

Driven pads: well posed, and the geometry difference collapses

Replace the welded pads with a real finger drive and grip force becomes a defined quantity: stiffness times how far past the surface the finger is commanded. Nothing about the contact model changes.

Each cell is downward slip / sideways lag in millimetres. Downward slip is the object sliding through the grasp; sideways lag is it trailing the shake, which is inertia rather than slipping.

From one millimetre of squeeze upward, FeatherPGS holds the object with essentially zero downward slip, and it does not care whether the object is a box or a triangle mesh. SolverMuJoCo on the same coupon needs four millimetres of squeeze to get close, and at one millimetre it lets the object slide 21 mm. That ordering is the right way round: holding an object firmly with a modest squeeze is what a hard contact model is for.

On this grasp FeatherPGS still holds the object at a 16.7 ms step — one physics step per control frame — with 1.6 mm of slip, while SolverMuJoCo loses it beyond about 2 ms.

Friction: a biased anchor, a mismatched combination rule, and a better formulation you have to ask for.

Friction anchors trade a downward slide for an upward push

The friction anchor options stop the static coupon's mesh from sliding, and replace the slide with travel in the other direction:

"It no longer falls" is not the same as "it is held". Both options convert a downward slide into an upward push, which is what you would expect from a friction row carrying a positional bias rather than only a velocity target. A grasp that slowly extrudes the object upward out of the fingers will pass a drop test and fail everything else.

# FeatherPGS
mu = 0.0
mat_count = 0
if shape_a >= 0:
    mu += shape_material_mu[shape_a]
    mat_count += 1
if shape_b >= 0:
    mu += shape_material_mu[shape_b]
    mat_count += 1
if mat_count > 0:
    mu /= float(mat_count)

# SolverMuJoCo
resolved_friction = wp.max(geom_friction[worldid, g1], geom_friction[worldid, g2])

Extra velocity iterations make the grasp worse, because speculative contacts become hard walls.

pgs_velocity_iterations solves positions with a position-correction bias, then cleans up velocities without one. It should improve a grasp. It does the opposite, and the reason is one hard-coded argument.

A contact whose surfaces are still apart has a positive gap. The correct velocity-level constraint for it is "you may close at up to gap over dt", which is what the position pass builds. The velocity pass rebuilds every right-hand side with the speculative scale forced to zero:

self._stage4_compute_rhs_world(
    dt,
    bias_scale=0.0,
    contact_speculative_scale=0.0,          # <-- not a parameter
    joint_limit_speculative_scale=1.0,
    output=self.rhs_unbiased,
)
if self._has_free_rigid_bodies:
    self._compute_mf_rhs_bias(dt, bias_scale=0.0, speculative_scale=0.0,
                              output=self.mf_rhs_unbiased)

Dropping the position bias is correct and is the point of the pass. Dropping the speculative term is not: the constraint becomes "you may not approach at all", applied to pairs that are millimetres apart. Every one turns into an invisible wall for the duration of the velocity iterations.

The arm scene's collision gap is 5 mm per shape, so during the shake there are typically 14 touching contacts and 32 speculative ones. Shrinking the gap by 50× so almost no speculative contacts exist removes the regression and lets the pass do what it is supposed to:

How it compares to SolverMuJoCo: same accuracy for more compute, and far more headroom on the step.

What the grasp looks like when it works

The measurement throughout is the object's position in the gripper's own frame, referenced to the start of the shake. That removes the commanded motion and leaves only relative movement — slip, squeeze and rattle.

The control loop stays at 60 Hz; only the number of physics substeps changes. Speed is measured with the physics loop inside a captured CUDA graph so it reflects the solver rather than Python.

Two things fall out. FeatherPGS starts 2.6× behind on creep at the smallest step and the two converge as the step grows — at 8.3 ms they are within 5% of each other. And the cost ratio barely moves: SolverMuJoCo is 2.7–2.8× cheaper per frame at every step, because a flat per-substep cost across a 16× change in substep count means both are latency-bound rather than compute-bound at one environment. A single 9-DOF arm produces tens of kernel launches per substep, almost all with grids of one to thirty-two threads. The Cholesky factorisation of the 9×9 mass matrix runs on one thread — small_dof_threshold=12 selects the serial kernel below thirteen degrees of freedom — and the Gauss-Seidel solve is launched as one tile of thirty-two threads per world, in both solve modes. There is nothing wrong with that at scale; it is the wrong shape for one arm.

Three costs are avoidable at any scale:

• Three debug copies run on every step in every configuration. They are not behind a flag, and nothing in the repository ever reads the buffers they write.
• Twenty-one zero-fills alongside twenty-four kernel launches in the contact-row build stage, ten of the fills only mattering in a mode that is off by default.
• Thirteen scoped-timer objects per step, plus the generator chain of the parallel-stream dispatch, which short-circuits to nothing when there is only one articulation size but is still constructed at every call site.

self._stage3_compute_v_hat(state_in, state_aug, dt)
self._clamp_rigid_velocity_limits(self.v_hat)
wp.copy(self._debug_stage3_qd_work, self.qd_work)
wp.copy(self._debug_stage3_joint_qdd, state_aug.joint_qdd)
wp.copy(self._debug_stage3_v_hat, self.v_hat)

Where the time is best spent

There are three ways to spend more time on the same second of simulation: smaller substeps, more solver iterations per substep, or — on an unmerged branch — frozen inner substeps that reuse the mass-matrix factorisation and the contact Jacobians.

Iterations are the cleanest lever, and they are what closes the accuracy gap: twenty iterations reaches 0.31 mm, exactly SolverMuJoCo's number on this scene, and thirty-two beats it at 0.22 mm. The price is 20.9 ms per frame against SolverMuJoCo's 5.2 ms — four times the compute to match. Past thirty-two there is nothing left to buy. The wrinkle is that convergence is not monotone — eight iterations is worse than four and much worse than twelve. Under-converged Gauss-Seidel does not give a slightly softer grasp; it gives a different one. That is an argument for running enough iterations to be past the noisy region rather than for tuning the number.

The frozen inner substeps do not help here

Two solver substeps with no inner loop gives 3.43 mm at 1.95 ms. Adding four frozen inner substeps — a four-times-smaller integration step — gives 3.51 mm at 5.82 ms. Refining the integration changed nothing, because the error is not in the integration; it is in the contact basis, which the frozen loop deliberately does not refresh. At the same effective step size, plain substepping gives 1.30 mm.

It does earn its keep rescuing a step too large to work at all: at one solver substep the plain solve drops the cube, and four frozen inner substeps turn that into a 7.47 mm hold at 5.9× realtime. But the behaviour is not monotone — eight loses the grasp and sixteen drops the cube again — which is what you would expect when positions are advanced far from where the Jacobians were built. As a speed lever for grasping it needs the contact basis refreshed on some cadence, not frozen for the whole step.

The default execution settings cannot be captured into a CUDA graph.

Every speed number here depends on capturing the 16-substep physics loop once and replaying it. With the constructor defaults, that capture fails.

Only double_buffer is responsible, and the reason is structural. The double-buffer memset forks onto a private stream and never joins back within the same step:

if self._memset_stream is not None:
    with wp.ScopedTimer("DB_Memset", ...):
        with wp.ScopedStream(self._memset_stream):        # forks the capture
            for size in self.size_groups:
                if self._mass_update_global_flag:
                    self._H_bufs[self._buf_idx][size].zero_()
                self._J_bufs[self._buf_idx][size].zero_()
        self._memset_done_event[self._buf_idx] = self._memset_stream.record_event()
        self._buf_idx = 1 - self._buf_idx

wp.ScopedStream synchronises on entry but not on exit. The only join is the wait_event at the top of the next step, so the final step of any captured region leaves an unjoined fork and capture fails. On the first captured step the event being waited on was also recorded outside the capture, which is separately illegal.

A method exists that would fix the second half of this, seed_double_buffer_events(). It has zero callers anywhere in the repository, and its docstring says it must be called inside graph capture, which is the opposite of what would work. No test in the repository captures a FeatherPGS step(), which is why this survived.

Smaller things that will each cost somebody an afternoon

In order of how much grasp stability each one buys

Gravity compensation is deliberately not on this list. A position-controlled arm needs one, and FeatherPGS has none — but neither does any other Newton solver, so it belongs to the platform rather than to this solver. Supply it however you like; everything below is what remains once you have.

Pinned versions and the harness

git clone --branch fpgs-main https://github.com/dylanturpin/newton-collab.git
cd newton-collab && git checkout d95efe1d1bb07497a4e51bf9563f36bc8821ab73
uv venv .venv --python 3.12 && uv sync --extra examples

git clone https://github.com/StoneT2000/newton-tests.git
cd newton-tests && git checkout 32c1544bac713bb722e9f805c26dc99413f0e359

from harness import make_env
from gravcomp import GravityCompensator

env  = make_env("fpgs")                 # original gains, nothing else changed
comp = GravityCompensator(env)

for frame in range(900):
    comp.apply()                        # writes tau_g into control.joint_f
    env.step()