NVIDIA Newton vs. FEniCSx (FEM) — a soft-body accuracy benchmark

A quantitative, apples-to-apples comparison of one deformable soft body simulated two ways:

• NVIDIA Newton (on Warp, CUDA) — three solvers: XPBD (Extended Position-Based Dynamics: fast positional projection), VBD (implicit), SemiImplicit (explicit, differentiable);
• FEniCSx / dolfinx — implicit FEM, used as the reference solve.

Same mesh (one generator builds it for both codes), same mass (ρV), same material, same gravity — only the solver differs. The material is one Lamé pair (μ, λ), handed to each code in the parametrisation it expects: FEM uses a compressible Neo-Hookean law; Newton's VBD and SemiImplicit a stable Neo-Hookean (Smith et al. 2018) that equals it at small strain; Newton's XPBD a two-constraint form with no small-strain shear stiffness (see docs/METHOD.md §2). The goal is to make it measurable how far Newton's fast solvers deviate from an accurate FEM solve, and exactly why.

Terms used below (details in docs/METHOD.md):

• frame — 1/60 s of simulated time; a substep is one time-integration step inside a frame; an iteration is one solver sweep over the mesh within a substep; the budget is substeps × iterations.
• settle — run the dynamic solver until no node moves by more than 0.1 % of the deflection over a 40-frame window; a run that hits the 600-frame limit first is reported as cap.
• tip — the bottom-face node on the bar's centre axis.
• equilibrium residual — the out-of-balance nodal force of a settled state, evaluated with the FEM's material law (≈ 0 for a true static equilibrium).
• θ* — the stiffness scale that best fits SemiImplicit's result to FEM (a differentiable fit; it characterises SemiImplicit, not XPBD).

Key result (hanging bar, Newton 1.6.0 on an A100): with the set-up matched, the implicit VBD and the explicit SemiImplicit land on the FEM tip (−0.2 % and −0.4 %). FEM on this mesh is 2.5 % below the 1-D analytic bar; refining the mesh closes most of that gap (0.8 % below on the 4× refined mesh, still rising), so most of it is discretisation error. XPBD settles ≈ 3.1× soft at its 10-iteration budget. Its tet law is two constraints with no small-strain shear stiffness, and Newton does not carry the constraint multipliers across iterations, so its effective stiffness grows with the iteration count. It leaves a finite force-balance residual, and its answer keeps moving with the iteration count (1.18× FEM at 256 iterations): it is a budget-dependent projection, not a converging solve. Numbers with provenance: docs/STATUS.md; the mechanism: 10_hanging_bar, 30_convergence.

What an answer costs (same bar, 833 to 40 625 nodes): time until the tip stays within 1 % of FEM on the same mesh (figures/performance_report.txt, logs/fem_performance.log):

mesh (nodes)	SemiImplicit [s]	VBD, fastest budget [s]	XPBD	FEM warm static solve [s]
833	0.07	2.7	none of the budgets tried	0.12
40 625	1.1	23	none of the budgets tried	15

SemiImplicit reaches the answer soonest on both meshes. The comparison sets dynamic relaxations on the GPU against a static solve on the CPU, so it is not a per-step comparison. The caveats are in docs/STATUS.md → Accuracy vs compute time; the walkthrough is 35_performance.

Friction: read on a tilted floor, XPBD's friction is exactly Coulomb for a lone particle. On its own soft block it is not: the block slides below μ, and at 64 iterations almost without friction. VBD and SemiImplicit creep under any load below μ (40_friction §4).

Correction (2026-10). Results before 2026-10 carried set-up errors (node masses, SemiImplicit's material map, contact geometry); docs/STATUS.md → Corrections records what changed and why.

The references are layered: analytic → FEM → Newton

Where a closed-form solution exists, it anchors both sides — so a skeptic can follow the whole trust chain, not just take FEM on faith:

• the 1-D self-weight bar (hanging bar) — tip elongation ρgL²/2E (an approximate anchor; it omits Poisson contraction and 3-D effects);
• the confined uniaxial Neo-Hookean stress law (material test) — the diagnostics' stress evaluation reproduces it to machine precision (tests/test_energies.py); the FEM and SemiImplicit material curves are observations against it;
• the Coulomb μ·W plateau and N = W (friction);
• the Hertz sphere-on-half-space force F = (4/3)·E*·√R·δ^(3/2) (indentation, an approximate anchor for a finite soft slab; the full contact method — penalty / AL law, Hertz derivation, Newmark drop — is in docs/CONTACT.md).

What the FEM reference is (for non-FEM readers): a finite-element solve in FEniCSx/dolfinx that discretises the body into small elements and solves the discretised force balance to convergence — a genuine static equilibrium (net nodal force ≈ 0), unlike a fast positional solver that only projects positions. Its material is a compressible Neo-Hookean law (hyperelastic, i.e. rubber-like) at the Lamé parameters μ (shear modulus — resistance to change of shape) and λ (volumetric stiffness). The FEM is itself checked against these analytic solutions in the simple cases; the fast Newton solvers are then scored against it on the same mesh and material. So the comparison is analytic → FEM → Newton, with each link measured; the diagnostics behind the measurements are tested (tests/).

The scenarios (named for what they do)

scenario	what it does	the point
hanging bar	a soft bar stretches under self-weight	a closed-form deformation (1-D tip elongation) to score every solver against node-for-node — with the FEM solve as the reference
indentation	a rigid sphere is pressed into a soft slab	FEM gives a calibrated contact-force curve; the fast XPBD gives deformation, not a force (VBD/explicit selectable via --solver; method in docs/CONTACT.md)
drop	a sphere is dropped onto a block (dynamic impact)	transient impact; FEM Newmark + contact vs. Newton solvers (implicit VBD is the natural match; see docs/CONTACT.md)
friction	a block is dragged on a rigid floor	FEM friction force + dissipated work vs. analytic μ·W; Newton slip only — and each solver's friction law on a tilted floor, apart from the block
material test	confined uniaxial squeeze/stretch	stress vs. stretch into large strain (constitutive fidelity)
convergence	swept budgets / meshes	discretisation error (FEM) vs. solver-budget error (XPBD iterations/substeps, VBD iterations)
accuracy vs. time	the hanging bar, timed, on three mesh sizes	what each solver's answer costs: solver time per frame, frames to an answer, error against FEM

Each scenario has newton_run/run_<x>.py, fenics_run/run_<x>.py and compare/<x>.py (the overlay). The analysis lives in the 10/15/20/25/30/35/40_* notebooks, each written as a 10-minute read for a skeptical expert (verdict first, then the mechanism), with a rendered 3D scene of what is being simulated. The contact scenarios (indentation / drop / friction) run all three Newton solvers via --solver on the same soft-body contact path (Newton's particle–shape contact), so the implicit VBD — not just XPBD — is the implicit counterpart to the implicit FEM; docs/CONTACT.md covers the version caveats and what "no calibrated contact force" does and does not mean.

The three Newton solvers

• XPBD — positional projection of two constraints per tet (I_C − 3, J − 1); fast; leaves a finite equilibrium residual. Its constraint energy has no small-strain shear stiffness and Newton does not carry the constraint multipliers across iterations, so its effective stiffness depends on the iteration budget (≈ 3.1× soft at 10 iterations). compare/xpbd_model.py predicts where it settles from that constraint energy alone, without running it — to 1.5 % up to 64 iterations; beyond that XPBD settles softer than the model (1.18× FEM at 256 iterations). The canonical run.
• VBD — Vertex Block Descent; implicit (minimises the backward-Euler objective). Converts the Lamé pair internally and lands on FEM; on this bar at its 10-substep budget the iteration count is inert, on finer meshes it is not (the finest needs 20 iterations to get within 1 %).
• SemiImplicit — explicit, force-based; the differentiable one → used by the diffsim θ* fit and the material test. Lands on FEM when handed Smith's map of the Lamé pair.

All three are set up through one module (newton_run/_solver.py): material per solver, consistent tet-lumped masses, the clamp, the collision pipeline, and a frame stepper that is CUDA-graph-captured on a GPU.

Platform

Newton needs a CUDA GPU for the real runs, so the repo runs on any CUDA-capable machine; Warp's CPU backend runs every Newton script too, just slowly (CI uses that for a smoke test). The easiest GPU is Google Colab — interactively via 00_setup.ipynb (File → Open notebook → GitHub), or from the terminal (below). FEniCSx is CPU (fem-on-colab on the same Colab VM, or conda-forge locally). The engines are pinned exactly — newton 1.6.0, warp-lang 1.17.0 (pyproject.toml → uv.lock) — and every result records the versions it ran on.

If the two large stacks clash in one runtime (numpy/petsc/mpi), run the Newton and FEM sides in separate notebooks — they only exchange data/*.npz.

Quickstart

uv sync --extra dev --extra newton                 # locked env: Newton/Warp + ruff, mypy, pytest
uv run pytest                                      # tests (no GPU; the Newton smoke runs on Warp's CPU)

python -m newton_run.run_hanging_bar               # XPBD (CUDA; --device cpu works, slowly); add --solver vbd | semi_implicit
python -m fenics_run.run_hanging_bar --element tet  # FEM on the shared mesh (node-for-node)
python -m fenics_run.run_hanging_bar --element hex  # FEM, independent hex mesh
python -m compare.hanging_bar                       # overlay + figures/ + a text report

(Without uv: pip install -r requirements.txt installs the package with the same pinned engines.) Then the other scenarios: run_indentation, run_drop, run_friction, run_stress_strain, convergence, performance (each on both newton_run / fenics_run, then compare.<scenario>), newton_run.friction_law and newton_run.diffsim — or all of them in order with python -m common.pipeline (--list shows the stages; --only 7.a 8 --update re-runs some and keeps the rest). Each stage appends an OK/ERR line + key results to logs/summary.txt.

On a Colab GPU, from the terminal

tools/setup-colab-auth.sh                           # once per machine, in a real terminal (OAuth code paste)
tools/colab-campaign.sh --dry-run "common.pipeline --clean"      # the plan + cost disclosure, nothing allocated
tools/colab-campaign.sh --fem --timeout 18000 \
    --fetch data --fetch figures --fetch logs "common.pipeline --clean"   # ≈ 3.6 h; `colab usage` before/after shows the cost

The campaign packages the working tree, allocates an A100 through the colab CLI, installs the locked engines (and FEniCSx via fem-on-colab with --fem), runs the targets, brings the requested result directories home (replacing the local ones) and releases the VM in a trap — Ctrl-C and errors included. It allocates billable hardware, so it asks once with the whole plan (--yes for unattended runs). Driver logs and the session notebook land in out/campaigns/ (git-ignored).

Layout

src/common/params.py        single source of truth (geometry, material, gravity, paths, per-scenario params)
src/common/mesh_io.py       the one mesh source for Newton and FEM (tet orientation, lumped masses, the axis tip node)
src/common/settle.py        when a relaxation run is at rest (positions over a window, not a KE threshold)
src/common/arrays.py        numpy array type aliases used in annotations
src/common/provenance.py    versions / device / commit stamped into every result
src/common/pipeline.py      the whole comparison as ordered stages (== 00_setup.ipynb)
src/common/runlog.py        stage runner (live stream + logs/summary.txt)
src/newton_run/   run_hanging_bar · run_indentation · run_drop · run_friction · friction_law · run_stress_strain · convergence · performance · diffsim · _solver (shared Newton set-up)
src/fenics_run/   run_hanging_bar · run_indentation · run_drop · run_friction · run_stress_strain · convergence · performance
src/compare/      hanging_bar · indentation · drop · friction · stress_strain · convergence · performance · energies · scene · xpbd_model
tests/            mesh · energies · material_mapping · masses · settle · xpbd_model · friction_law · pipeline · newton_smoke
tools/            setup-colab-auth.sh · colab-campaign.sh   (Colab GPU from the terminal)
10_hanging_bar · 15_material · 20_contact · 25_dynamic · 30_convergence · 35_performance · 40_friction   (analysis notebooks)
00_setup.ipynb     install + run the whole pipeline interactively (any CUDA GPU; Colab = easy path)

compare/energies.py (pure numpy) computes the diagnostics identically for both solvers — strain energy, Jacobian (det F), internal nodal forces f = −∂U/∂x, and the equilibrium residual f_internal + f_gravity (the sharpest "solve vs. project" measure). Its correctness claims are backed by tests/test_energies.py; the Newton set-up (the shared mesh, each solver's material law, the masses, the settle rule) is backed by the other tests. compare/scene.py renders the deformed body in 3D (matplotlib, headless-safe), solver-agnostic.

Documentation

• docs/METHOD.md — problem setup, material matching, element variants (tet/hex locking), the three solvers, the energy/residual diagnostics.
• docs/CONTACT.md — the indentation contact law (penalty / Augmented Lagrangian, no C++) and the dynamic drop.
• docs/EXPERIMENTS.md — friction (and the friction law on a tilted floor), the material test, the convergence study, the differentiable θ* fit, and accuracy vs. compute time.
• docs/STATUS.md — what is tested, what is an observation with provenance, what is still open, with the result numbers and their provenance.
• CLAUDE.md — concise project guide / conventions.

Development

uv sync --extra dev --extra newton    # the locked environment (.venv)
uv run ruff check . && uv run ruff format --check .
uv run mypy                           # type check of src/ and tests/: every function fully annotated (config in pyproject.toml)
uv run pytest                         # all tests; NEWTON_VS_FEM_SKIP_SMOKE=1 skips the Newton CPU smoke
pre-commit install                    # optional: ruff + hygiene hooks on each commit

CI (.github/workflows/ci.yml) runs on uv with --locked (a resolution-drift gate): ruff lint + format check, mypy, the tests — including a CPU smoke run of every Newton entry point at tiny budgets that fails on an API break or a Newton/Warp deprecation — and a package build. It produces no results: GitHub Actions has no CUDA GPU and no fem-on-colab, so the numbers come from GPU runs (the campaign above or 00_setup.ipynb). The src/ layout means the packages must be installed (uv sync / pip install -e .) before python -m … resolves them.