Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion 10_hanging_bar.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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/`)."
]
}
Expand Down
28 changes: 18 additions & 10 deletions 30_convergence.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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`."
Expand Down Expand Up @@ -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."
]
},
{
Expand Down Expand Up @@ -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."
Expand Down
3 changes: 2 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
# How accurate is NVIDIA Newton, really? — Newton vs. FEniCSx (FEM)
<!-- headline: confirm wording -->
# NVIDIA Newton vs. FEniCSx (FEM) — a soft-body accuracy benchmark

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

Expand Down
2 changes: 1 addition & 1 deletion docs/METHOD.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion docs/STATUS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
Expand Down