METHOD.md

METHOD — the hanging bar, the material, and the diagnostics

This is the methodological core of the comparison: the one physical problem, how the two codes are made to describe exactly it, the three Newton solvers, the FEM reference, and the pure-numpy diagnostics that turn "looks soft" into numbers. Everything here is defined once in common/params.py — the single source of truth both sides import.

1. The problem: a soft bar hanging under its own weight

A soft block is clamped on its top face (max-z) and hangs downward; gravity (−z) stretches it. We start with this case because its closed-form deformation (the 1-D self-weight bar) scores every solver against an answer for the displacement field itself, not only against each other. (The other scenarios carry analytic anchors too — Hertz, the Coulomb μ·W plateau, the Neo-Hookean stress law — but for a force or stress, not the whole deformed body.)

quantity	symbol	value	source
block size	BLOCK_LX × LY × LZ	0.6 × 0.6 × 1.6 m	GRID_DIM = 6×6×16, CELL = 0.1 m
long (hanging) axis	z	clamped face = top (max-z)	FIXED_FACE = "top"
density	ρ	100 kg/m³	DENSITY
→ mass	ρV	57.6 kg (tet-lumped, both codes)	mesh_io.lumped_masses
Lamé parameters	μ, λ	1.0×10⁴, 5.0×10⁴ Pa	K_MU, K_LAMBDA
→ Young's modulus, Poisson	E, ν	≈ 2.83×10⁴ Pa, ≈ 0.417	lame_to_E_nu()
gravity	g	9.80665 m/s² (−z)	GRAVITY

Gravity is set explicitly on the Newton builder (the vector (0, 0, −GRAVITY), via newton_run._solver.make_builder) because Newton's default is −9.81; this matches FEM to the same constant. The resulting nominal (average) axial strain is ≈ 2.8 % (tip elongation / length); the local axial strain is zero at the free tip and peaks at ≈ 5.5 % at the clamp (ρgL/E). Small enough that linear bar theory is a meaningful anchor, large enough that the nonlinear terms are exercised. (μ is the shear modulus — resistance to change of shape; λ sets the volumetric stiffness — resistance to change of volume; the row above converts them to the equivalent Young's modulus E and Poisson ratio ν.)

The closed-form anchor

For a 1-D bar of length L hanging from a clamp, the downward displacement of a material point originally at depth s = z_top − z below the clamp is

u(s) = (ρ g / E) · (L·s − s²/2),     tip elongation u(L) = ρ g L² / (2E).

With the numbers above the tip elongation is ≈ 44.3 mm (params.analytic_hanging_displacement, printed by params.summary()). The matching 1-D strain energy is U = ρ²g²A L³ / (6E) (energies.analytic_hanging_strain_energy).

Honest scope: this 1-D solution ignores Poisson contraction and 3-D effects, so it is a sanity anchor, not the ground truth — the FEM solve is the reference the fast solvers are scored against. On the shared mesh the tet FEM tip is 43.2 mm, 2.5 % below the analytic 44.3 mm. Refining the mesh raises it, by a shrinking step: 43.75 and 43.96 mm on the 2× and 4× refined tet meshes, 43.68 mm on the finest hex mesh of the convergence study. Neither sequence has reached its limit, so most of the 2.5 % is discretisation error. How much of the remaining gap (below 1 %) is the 3-D/Poisson effect the 1-D formula omits is not resolved by these meshes (EXPERIMENTS.md §3).

2. Same problem, two codes

The point of the project is that the only difference between the two simulations is the solver. Three mechanisms enforce that:

• One shared mesh. A single generator, common/mesh_io.py block_tet_mesh, builds the tet mesh of the block: nodes on a regular grid, each cubic cell split into five tets, alternating between two mirror-image splits so neighbouring cells share their face diagonals. Newton builds every soft body from it (newton_run._solver.add_soft_block hands the nodes and tets to Newton's add_soft_mesh); the FEM tet run builds its dolfinx mesh from the same arrays (hanging_bar_mesh, create_mesh(comm, cells, element, x)) and evaluates its solution at those nodes, so the two results are comparable node-for-node (the point evaluation uses dolfinx's point-location API) without either side reading the other's output; compare.hanging_bar checks that the nodes agree. The grid, node numbering and split are exactly those of Newton's add_soft_grid; tests/test_mesh.py checks that add_soft_grid builds this mesh and that, once its masses are corrected (next point), it gives a model identical to add_soft_mesh's in every array (Newton leaves particle_colors uninitialised unless the model is coloured, so for XPBD and SemiImplicit only its shape is compared). Tets are defensively re-oriented to positive volume for dolfinx (orient_tets_positive; the generated ones already are).

• One mass. Newton's add_soft_grid gives every node a full cell of mass (cell volume × ρ) instead of sharing each cell's mass among its nodes, so the body comes out (nx+1)(ny+1)(nz+1)/(nx·ny·nz) times too heavy — 83.3 kg instead of 57.6 kg here, and 1.45× on the load that stretches the bar, because the bottom node layer alone carries 2.7× its share. Every Newton model therefore gets consistent tet-lumped masses, ρ·V_tet/4 per tet node (newton_run._solver.lump_tet_masses, before the clamp zeroes the fixed face). They sum to ρV, they are exactly the P1 gravity load the FEM side integrates, and they are what Newton's own add_soft_mesh assigns, to float32 rounding (tests/test_mesh.py); lump_tet_masses pins them to the one definition in mesh_io.lumped_masses. tests/test_masses.py checks the factors on the shared mesh.

• One material, in each code's own parametrisation. FEniCSx uses a compressible Neo-Hookean strain-energy density — a hyperelastic (rubber-like) law whose stored energy ψ depends on each element's stretch via the deformation gradient F = I + ∇u:

  ψ = μ/2 (I_C − 3) − μ ln J + λ/2 (ln J)²,    P = ∂ψ/∂F   (1st Piola–Kirchhoff stress: force per unit *rest* area)
  

  Here C = FᵀF, I_C = tr C is the sum of the squared principal stretches (3 at rest), and J = det F is the local volume ratio (1 at rest).

  add_soft_mesh documents its k_mu, k_lambda as the Lamé pair, but Newton's three solvers read them differently (Newton 1.6.0 source; newton_run._solver.lame_for_solver hands each solver the pair it needs):

  solver	tet law	fed (μ, λ) it models, at small strain	handed
  VBD	stable Neo-Hookean, converts internally (λ_NH = λ + μ)	(μ, λ)	(μ, λ)
  SemiImplicit	stable Neo-Hookean of Smith et al. 2018, ψ = k_μ/2 (I_C−3) − k_μ/2 log(I_C+1) + k_λ/2 (J−α)² (α a constant offset that keeps the rest state stress-free), coefficients used as given	(¾μ, λ − ⅝μ) — E 1.32× too small	Smith's map (4/3·μ, λ + 5/6·μ) → (μ, λ)
  XPBD	two constraints, I_C − 3 and J − 1, compliances 1/(V₀k_μ), 1/(V₀k_λ) (V₀ the tet's rest volume)	no shear stiffness, volumetric 4μ + λ	(μ, λ) — no map exists

  So VBD and SemiImplicit model the FEM's material exactly at small strain (they differ from it only in the large-strain shape of the law), while XPBD's constraint energy, μ/2 (I_C−3)² + λ/2 (J−1)², has no small-strain shear stiffness: I_C − 3 is second order in a pure shear. That is a property of the method, not a parameter choice. The maps come from linearising each kernel's energy about the rest state F = I: Smith's ψ fed (k_μ, k_λ) has shear modulus ¾k_μ and Lamé λ = k_λ − ⅝k_μ, and inverting that gives k_μ = 4/3·μ, k_λ = λ + 5/6·μ (lame_for_solver). The 4μ + λ volumetric entry of XPBD is the same linearisation of its constraint energy under a pure dilation. tests/test_material_mapping.py linearises transcriptions of the three kernels and checks every entry of this table.

  The pure-numpy diagnostics (§4) evaluate every result with the FEM's ψ, so Newton's settled state and FEM's solution are scored on identical footing.

The two FEM element variants

fenics_run/run_hanging_bar.py --element {tet|hex}:

• tet (default) — built from the shared mesh, genuine simplex linear (P1) tetrahedra (basix Lagrange/tetrahedron/1, full Galerkin quadrature), node-for-node with Newton. These are real 4-node constant-strain simplices — not the collapsed/degenerate-hexahedron tets with reduced integration + hourglass control that explicit codes (e.g. LS-DYNA) often use for "tets". Real P1 simplex tets are known to be over-stiff in bending and near-incompressibility (locking) because a constant-strain element is too stiff to represent those modes, so this is the honest like-for-like element, not the most accurate one. → data/fem_result.npz
• hex — an independent structured hexahedral mesh of the same block geometry (create_box, trilinear Hex8, full 2×2×2 integration — no reduced-integration hourglass modes). Hex8 is far less prone to shear/volumetric locking, so it is a second, more accurate FEM reference. It uses the same 6×6×16 cell grid and the same 833 nodes / 2499 DOFs as the tet run (only the cell type and split differ), so it is the same fineness — it just costs more per element (8 Gauss points, denser LU). It is built independently, so its node ordering differs from Newton's and it is compared via the displacement profile, not node-for-node. → data/fem_result_hex.npz

Reporting both tells the reader how much of any Newton–FEM gap is solver and how much is element choice — the hex result brackets the tet result from the accurate side.

The FEM solve itself

Compressible Neo-Hookean elastostatics, solved with dolfinx 0.11's SNES-based NonlinearProblem (Newton line-search, direct LU; these meshes are small). Gravity is ramped over 5 load increments for nonlinear robustness; the converged solution is independent of that count (the convergence study checks exactly this). Per-step SNES iteration counts and the converged reason are printed.

3. The three Newton solvers

A frame is 1/60 s of simulated time. Each frame is split into substeps (time-integration steps), and within a substep XPBD and VBD sweep the mesh a number of times (iterations); substeps × iterations is a solver's budget.

newton_run/run_hanging_bar.py --solver {xpbd|vbd|semi_implicit}. All three settle the same grid, with the same masses and the material in their own parametrisation (§2), toward a static state; they differ in how they get there and in what that state is. With the set-up of §2 the observed tips are VBD 43.1 mm and SemiImplicit 43.0 mm against FEM 43.2 mm (analytic 1-D bar 44.3 mm), and XPBD 134.9 mm — numbers and provenance in STATUS.md.

solver	class	kind	budget	writes
XPBD	SolverXPBD	positional projection (fast)	32 substeps × 10 iters	newton_result.npz
VBD	SolverVBD	implicit (Vertex Block Descent)	10 substeps × 50 iters	newton_result_vbd.npz
SemiImplicit	SolverSemiImplicit	explicit, force-based, differentiable	32 substeps	newton_result_semi.npz

• XPBD is the canonical run. It enforces its two tet constraints by projecting positions, not by solving a force balance, so its settled state leaves a finite equilibrium residual (§4). Two properties of Newton's XPBD tet constraint decide where it settles. First, the constraint energy has no small-strain shear stiffness (§2). Second, the Lagrange multiplier is not carried across iterations: the kernel computes a fresh multiplier each iteration, so each iteration applies another compliant correction. The effective stiffness therefore grows with the iteration count instead of converging to an elastic equilibrium. At 10 iterations it settles ~3.1× soft. It never fully comes to rest: its tip is stationary to 0.1 mm, but a sub-millimetre node jitter persists and the kinetic energy plateaus instead of decaying (the settle rule reports "cap"). The convergence study sweeps the iteration count to 256 to show where the tip goes (EXPERIMENTS.md §3), and compare/xpbd_model.py predicts that curve without running XPBD: the static equilibrium of the constraint energy scaled by 0.9·N (relaxation × iterations) under gravity, with no free parameter. It matches the measured tip to 1.5 % up to 64 iterations; beyond that XPBD settles softer than the model. This is the behaviour of the method as Newton implements it — and the residual is how it is measured — not a defect of the set-up.
• VBD minimises the backward-Euler objective by block coordinate descent over a coloured vertex graph (it requires builder.color(), which build_model runs only for this solver). Being implicit and converting the Lamé pair itself, it lands on the FEM answer (−0.2 % at the tip, node-for-node RMS error 0.07 mm), and at its 10-substep budget its iteration count is inert on this bar: 5 to 100 iterations give the same tip, 43.13 mm. That is not general — at 4 substeps, or on a finer mesh, the settled answer depends on the iteration count (EXPERIMENTS.md §5).
• SemiImplicit is the explicit, force-based integrator. It is the one Warp can differentiate, so it is what the θ* stiffness fit (diffsim.py) and the material test use. Handed Smith's map of the Lamé pair it lands on FEM too (−0.4 %); fed the raw pair it would be 1.32× too compliant. It is not used to claim XPBD's accuracy.

Settling to a static answer

A dynamic solver hit with full gravity at t=0 rings (under-damped). To measure the static equilibrium rather than a snapshot of an oscillation, the settle loop applies a per-frame velocity damping SETTLE_VEL_DAMP = 0.97 that drains transient kinetic energy. At rest v ≈ 0, so the factor does nothing to the equilibrium itself — it only removes the transient.

When to stop is decided on positions, not on the kinetic energy of one frame (common/settle.py): a ringing body has zero velocity at the turning point of every swing, so a kinetic-energy threshold can fire at the extreme of a swing instead of at rest (tests/test_settle.py reproduces this on a synthetic trace). A run counts as settled once no node coordinate moved by more than 0.1 % of the largest displacement over a trailing window of 40 frames — longer than one period of the bar's first axial mode (~23 frames at 60 fps) — after at least MIN_SETTLE_FRAMES; otherwise it stops at MAX_FRAMES = 600 and is reported as "cap" with the window's spread. tests/test_settle.py checks that a turning point and a slow creep are not called rest.

On a GPU the substep loop of a frame is captured once into a CUDA graph and replayed (newton_run._solver.FrameStepper) — the same kernels in the same order, without the per-launch overhead that otherwise dominates this small mesh.

4. The diagnostics — and what actually backs them

compare/energies.py is pure numpy and computes the same quantities for either solver from (rest_q, final_q, tets):

• Strain energy U — the Neo-Hookean ψ integrated over the rest tets.

• Jacobian J = det F per tet and the total volume change — how much the solver compresses/inflates the material.

• Internal nodal forces f = −∂U/∂x from the 1st Piola stress P = μ(F − F⁻ᵀ) + λ ln(J) F⁻ᵀ, assembled to nodes by H = −V·P·Dm⁻ᵀ (V the tet's rest volume, Dm its rest edge matrix, whose columns are three edge vectors).

• Equilibrium residual r = f_internal + f_gravity per node. This is the sharpest "solve vs. project" measure: for a true static equilibrium (FEM) the free-node residual ≈ 0; for XPBD it is finite and measures how far the projected state is from the material's force balance. At the clamped nodes, −Σ r_z is the support reaction, which must equal the total weight ρgV (≈ 565 N for this block) — a built-in consistency check, but only at a true force balance. For a Newton solver the reaction departs from the weight by as much as its settled state departs from the FEM law's balance. Observed values (campaign 20261004T162836Z, recomputed from the data/ hanging-bar results with equilibrium_residual):

  solver	support reaction / weight	free-node residual RMS [N]
  FEM tet	1.000	≈ 0
  VBD	0.980	0.07
  SemiImplicit	0.977	0.07
  XPBD	0.340	9.4

  The XPBD mismatch is its non-equilibrium signature, not a failing check. For Newton, r is computed with the FEM's ψ above, not the law the solver integrated, so it folds a constitutive-mismatch term into the solve-vs-project gap. For VBD and SemiImplicit (stable Neo-Hookean, equal at small strain) that term is consistent with the observed 2 % reaction deficit; it grows with strain, and the local strain peaks at ≈ 5.5 % at the clamp. For XPBD, whose constraint energy differs from ψ even at small strain, it is part of what the residual reads.

The honesty backbone: the tests

The correctness claims the comparison rests on are tested, not asserted (pytest tests/, or each file as a script — pure numpy, no GPU/FEM needed). tests/test_energies.py:

1. test_nodal_forces_match_energy_finite_difference — the analytic f = −∂U/∂x matches a central finite difference of the strain energy (relative error < 10⁻⁶).
2. test_uniaxial_stress_matches_closed_form — the volume-averaged axial 1st Piola stress for F = diag(1,1,λ) reproduces the closed-form compressible Neo-Hookean law to machine precision (≤ 10⁻⁹) across λ ∈ [0.7, 1.5].

Plus two structural invariants: zero strain energy at the rest state, and internal forces that sum to zero (no net self-force). The set-up of §2 and §3 is backed the same way:

• tests/test_material_mapping.py — transcriptions of the three Newton 1.6.0 tet kernels, linearised about the rest state: SemiImplicit fed the raw pair is (¾μ, λ − ⅝μ) and fed Smith's map is exactly (μ, λ); VBD is (μ, λ); XPBD's constraint energy has zero shear modulus. (These check the transcriptions; the kernels are quoted in lame_for_solver.)
• tests/test_mesh.py — the generated mesh is conforming, positively oriented and fills the block; Newton's add_soft_grid builds exactly these nodes and tets (the recorded runs are on this mesh); add_soft_mesh fed the generator gives the lumped masses and, once lump_tet_masses has corrected add_soft_grid's, the same model in every array (nodes, tets, rest poses, materials, surface triangles and edges, radii, VBD's colouring).
• tests/test_masses.py — on the shared mesh, the lumped masses sum to ρV and the one-cell-per-node masses give the 83.3 kg body, 1.405× the free mass (78.4 vs 55.8 kg) and 1.446× the stretching load quoted above.
• tests/test_settle.py — the settle rule does not call a turning point or a creep "rest".
• tests/test_xpbd_model.py — the XPBD model's gradient matches a central finite difference of its energy, its minimiser is a static equilibrium (the supports carry the hanging weight), and the relaxation it assumes is SolverXPBD's default.
• tests/test_friction_law.py — each solver's particle-ground friction against Coulomb.
• tests/test_pipeline.py — the pipeline stages equal the notebook's, and the engine pins are identical everywhere.

These tests are the only things in the repo described as "validated"; anything depending on a GPU or a FEM install is reported as an observation with provenance (see STATUS.md). tests/test_newton_smoke.py additionally runs every Newton entry point at tiny budgets on Warp's CPU backend — a check that the code runs on the pinned engines, not of any number.

5. Reading the result

python -m compare.hanging_bar overlays whichever Newton results exist (XPBD / VBD / explicit) against FEM tet/hex and the analytic bar, writes the figures, and prints a text report: per solver the tip drop (the bottom-face node on the centre axis — the same node every run logs), the largest drop over all free nodes (a separate quantity: under an uneven load it sits at a corner), and the settle verdict (settled at frame N, or "cap"); then, per Newton solver, the node-for-node error vs. FEM tet (RMS / max [mm], mean relative [%]) and the tip ratio to FEM. The strain-energy, Jacobian and equilibrium-residual diagnostics from compare/energies.py are surfaced in the analysis notebook rather than this CLI report. The 3-D scene (compare/scene.py, matplotlib, headless-safe) renders the deformed bar coloured by displacement so the reader sees the stretch directly. The notebook 10_hanging_bar.ipynb walks through it verdict-first.

See EXPERIMENTS.md for the convergence study, which quantifies XPBD's budget dependence, and the differentiable θ* fit, which characterises SemiImplicit against FEM, and CONTACT.md for the indentation and drop scenarios.