diff --git a/10_hanging_bar.ipynb b/10_hanging_bar.ipynb index 5630bdd..f80b8c5 100644 --- a/10_hanging_bar.ipynb +++ b/10_hanging_bar.ipynb @@ -414,7 +414,7 @@ "## Takeaway (the 10-minute summary)\n", "\n", "* On this test, whose **deformation has a closed form**, FEM lands within a few percent of the analytic bar — but all three Newton solvers settle markedly softer at this run's budgets (XPBD ≈ 3.7×, VBD ≈ 3.2×, explicit ≈ 2.0× the FEM tip); the **explicit** solver is the closest of the three and VBD barely improves on XPBD (§1).\n", - "* The reasons differ by solver: **XPBD never reaches a force balance** — it leaves a finite equilibrium residual (§3) — because it *projects* positions instead of *solving* R(u) = 0 (mechanical, not a tuning artefact). **VBD** *is* implicit and would approach FEM with enough iterations, but at the budget used its block Gauss-Seidel has not converged on this 17-layer bar, so it stays ~3.2× soft — a convergence-budget effect, not a force-balance one.\n", + "* The reasons differ by solver: at any fixed budget **XPBD leaves a finite equilibrium residual** (§3) because it *projects* positions rather than *solving* R(u) = 0 — more iterations do drive that residual down (see `30_convergence`), but the cost to close it grows, so at a real-time budget it reads soft. **VBD** *is* implicit and would approach FEM with enough iterations, but at the budget used its block Gauss-Seidel has not converged on this 17-layer bar, so it stays ~3.2× soft — a convergence-budget effect, not a force-balance one.\n", "* That is the trade you are buying: **XPBD is fast on the GPU (§6)** at the cost of an unforced equilibrium; **VBD** is implicit but, at the fast budget here, still settles far soft (closing the gap costs many more iterations — see `30_convergence`); **FEM** is the force-balanced reference. Every diagnostic is computed identically for both sides; the underlying force/energy model is **unit-tested** (`pytest tests/`)." ] } diff --git a/30_convergence.ipynb b/30_convergence.ipynb index 8b7b9c8..7c4ac9c 100644 --- a/30_convergence.ipynb +++ b/30_convergence.ipynb @@ -9,7 +9,7 @@ "\n", "The Newton-vs-FEM gap has two different sources, one per solver family:\n", "\n", - "* **Newton (XPBD)** — *Extended Position-Based Dynamics* — is an iterative positional projection — its *effective stiffness* depends on the number of solver iterations and substeps. We sweep both and watch the equilibrium-residual RMS fall and the tip drop move toward (but not reach) the FEM / analytic value: more iterations lower the residual, yet the XPBD tip plateaus well above FEM (~2.4×) — more compute narrows the gap without closing it.\n", + "* **Newton (XPBD)** — *Extended Position-Based Dynamics* — is an iterative positional projection — its *effective stiffness* depends on the number of solver iterations and substeps. We sweep both and watch the equilibrium-residual RMS fall and the tip drop move toward the FEM / analytic value: more iterations lower the residual **monotonically** and the tip keeps descending — at the largest budget tested it is still ~2.4× FEM and falling, and the per-doubling gains shrink, so closing the residual gets progressively more expensive.\n", "* **FEM (FEniCSx)** converges under mesh **h-refinement**. We watch the tip drop and strain energy reach a mesh-independent limit, compare it to the analytic 1-D bar (FEM credibility), and confirm the *converged* answer is independent of the number of gravity load increments (a check on the nonlinear solve).\n", "\n", "This separation — *solver* error (XPBD) vs *discretisation* error (FEM) — is the heart of the Newton-vs-FEM difference. Produce the data with `python -m newton_run.convergence` and `python -m fenics_run.convergence`." @@ -51,10 +51,16 @@ "More **iterations** tighten the positional projection (substeps mainly buy stability — the\n", "substep sweep here is essentially flat). The **equilibrium residual** is the leftover net\n", "nodal force at the settled state — how far the body is from a true force balance (zero for an\n", - "exact solver, finite for XPBD); we report its RMS over free nodes, in newtons. It is the\n", - "cleanest \"have we converged\" measure: with more iterations it should fall and the tip drop\n", - "should move toward the FEM / analytic line — though even at the largest budget the XPBD tip\n", - "**plateaus well above** FEM (here ~2.4×); more compute narrows the gap without closing it." + "exact solver, finite for XPBD); we report its RMS over free nodes, in newtons. It is the cleanest \"have we converged\" measure: with more iterations it falls\n", + "**monotonically** and the tip drop moves steadily toward the FEM / analytic line. At the\n", + "largest budget tested (32 iters) the tip is at ~2.4× FEM and **still descending** — the\n", + "residual has not bottomed out. The per-doubling gains shrink, so closing the residual gets\n", + "progressively more expensive: that escalating cost is the price of a fast positional solver,\n", + "and the residual is the gauge for where a given real-time budget lands.\n", + "\n", + "**Note on the substep sweep:** the `substeps = 4` point is below the explicit stability floor\n", + "for this dt — its energy and residual blow up — so it is a stability artefact, not a point on\n", + "the convergence trend; read the substep sweep from `substeps ≥ 8`, where it is flat." ] }, { @@ -201,11 +207,13 @@ "source": [ "## Takeaway\n", "\n", - "* XPBD's answer is **iteration-budget dependent**: too few iterations leave a large equilibrium\n", - " residual (it under-resolves the elastic balance). More iterations lower the residual and move\n", - " the tip toward FEM, but it **plateaus well above** the FEM/analytic value (here ~2.4× even at\n", - " the largest budget), and the substep sweep is essentially flat (substepping mainly buys\n", - " stability). More compute narrows the gap without closing it — the price of a fast positional solver.\n", + "* XPBD's answer is **iteration-budget dependent**: too few iterations leave a large\n", + " equilibrium residual (it under-resolves the elastic balance). More iterations lower the\n", + " residual **monotonically** and move the tip toward FEM — at the largest budget tested\n", + " (32 iters) it is at ~2.4× FEM and **still descending**, not at a floor; the per-doubling\n", + " gains shrink, so closing the residual gets progressively more expensive. The substep sweep,\n", + " by contrast, is essentially flat (substepping mainly buys stability). That escalating cost —\n", + " not an intrinsic floor — is the price of a fast positional solver.\n", "* FEM's answer is **discretisation dependent** but converges monotonically under h-refinement to a\n", " budget-independent equilibrium, and is independent of the load-step schedule. That separation —\n", " *solver* error vs *discretisation* error — is the essence of the Newton-vs-FEM difference." diff --git a/README.md b/README.md index 7bd76ea..258b7cc 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ -# How accurate is NVIDIA Newton, really? — Newton vs. FEniCSx (FEM) + +# NVIDIA Newton vs. FEniCSx (FEM) — a soft-body accuracy benchmark A quantitative, apples-to-apples comparison of one **deformable soft body** simulated two ways: diff --git a/docs/METHOD.md b/docs/METHOD.md index f61673c..19bb966 100644 --- a/docs/METHOD.md +++ b/docs/METHOD.md @@ -90,7 +90,7 @@ is the solver. Two mechanisms enforce that: 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 poor to represent those modes, so this is the honest + 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 diff --git a/docs/STATUS.md b/docs/STATUS.md index 42b80b4..df30e56 100644 --- a/docs/STATUS.md +++ b/docs/STATUS.md @@ -104,7 +104,7 @@ What is **still numerically open** (quality, not wiring): | item | what is unsettled | |---|---| -| **SemiImplicit drop** | numerically unstable at this dt/substep budget — the energy blows up (peak strain energy 7536 J vs FEM 3.4 J); the only clearly-broken numeric this run | +| **SemiImplicit drop** | numerically unstable at this dt/substep budget — the energy blows up (peak strain energy 7536 J vs FEM 3.4 J); the only clearly-unstable numeric this run | | **FEM drop** (`fenics_run/run_drop.py`) | dt / damping are not tuned to convergence; the transient numbers are observations, not a settled benchmark | | **VBD / SemiImplicit indentation contact** | the `soft_contact` penalty is too soft — the sphere sinks ~33 mm through the 40 mm indent (strain energy ~0.1 J vs XPBD's ~13 J), so **XPBD** is the only Newton solver that geometrically resolves the indentation; the soft-contact stiffness for the implicit/explicit solvers is unsettled | | **θ\* fit** (`newton_run/diffsim.py`) | one converged fit characterising the **SemiImplicit** solver vs FEM (θ\* = 1.79); not cross-validated across budgets/scenarios |