From 327da119d8f1804b5a0b49e5be1dc6c5c8ad1945 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 09:43:52 +1000 Subject: [PATCH 1/7] Note: testing a solver against exact solutions UWTN 2026-012, draft. The uw.analytic suite and, more to the point, the decision to check each solution by differentiating it into the momentum balance rather than by comparing it against the kernel it was transcribed from. Four defects, two of them in vendored published sources. The measurement table's generating script ships in examples/. Underworld development team with AI support from Claude Code --- .../examples/convention_audit.py | 175 ++++++++ .../metadata.yml | 33 ++ ...esting-a-solver-against-exact-solutions.md | 388 ++++++++++++++++++ 3 files changed, 596 insertions(+) create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py create mode 100644 articles/testing-a-solver-against-exact-solutions/metadata.yml create mode 100644 articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md diff --git a/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py b/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py new file mode 100644 index 0000000..6810eaa --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py @@ -0,0 +1,175 @@ +"""Measure which convention every exact solution in `uw.analytic` obeys. + +This generates the table in the note. Nothing here reads a published paper or a +vendored C kernel: every number is formed from the solution's own symbolic +fields by differentiation, so a solution that agrees with a wrong source still +fails. That is the whole point of the exercise. + +The negative-control column is the one that makes the table mean anything. It +re-measures the momentum residual with the body force negated. A gate that +cannot fail is not a gate, and a table of small numbers with no such column is +an assertion that small numbers are small. + + pixi run -e amr-dev python convention_audit.py + +Underworld3 0.0.0 (development, 2026-08). +""" + +import numpy as np +import sympy + +import underworld3 as uw +from underworld3.analytic import _validation as V + + +def build_meshes(): + """One coarse box per dimension. Resolution is irrelevant here. + + The residuals are symbolic identities sampled at points; they do not + converge with resolution, they either hold or they do not. The mesh exists + only to carry the coordinate system the solutions are written on. + """ + + return { + 2: uw.meshing.StructuredQuadBox( + elementRes=(4, 4), + minCoords=(0.0, 0.0), + maxCoords=(1.0, 1.0), + qdegree=3, + ), + 3: uw.meshing.StructuredQuadBox( + elementRes=(2, 2, 2), + minCoords=(0.0, 0.0, 0.0), + maxCoords=(1.0, 1.0, 1.0), + qdegree=2, + ), + } + + +def stokes_solutions(): + """Every registered Stokes solution that exposes symbolic fields. + + `symbolic = False` marks a solution reached through a compiled third-party + package rather than a SymPy expression tree; there is nothing to + differentiate, so it cannot be measured this way. + """ + + for name in sorted(uw.analytic.available()): + cls = getattr(uw.analytic, name) + if not getattr(cls, "symbolic", False): + continue + if not uw.analytic.is_available(name): + continue + if not hasattr(cls, "dim"): + continue + yield name, cls + + +def trace_check(solution, points): + r"""Relative size of $\mathrm{tr}\,\sigma + d\,p$. + + Zero says the exposed stress is the TOTAL stress under + $\sigma = \tau - p\mathbf I$ with a traceless deviator. It is scaled by the + pressure because that is the term being cancelled. + """ + + dim = solution.dim + trace = sum(solution.fn_stress[i, i] for i in range(dim)) + residual = V.sample(solution, trace + dim * solution.fn_pressure, points) + pressure = V.sample(solution, solution.fn_pressure, points) + + return float(np.max(np.abs(residual))) / max(float(np.max(np.abs(pressure))), 1e-300) + + +def constitutive_check(solution, points): + r"""Relative size of $\sigma + p\mathbf I - 2\eta\dot\varepsilon$. + + Note what this can and cannot say. For a solution that published only one + of stress and strain rate, `set_fields` derived the other from exactly this + identity, so the check is structural rather than evidential -- and the + signature is visible in the output: an EXACT zero where it is structural, + order 1e-16 where two separately derived quantities happen to agree. + """ + + dim = solution.dim + worst = 0.0 + scale = 0.0 + + for i in range(dim): + for j in range(dim): + deviator = solution.fn_stress[i, j] + (solution.fn_pressure if i == j else 0) + constitutive = 2 * solution.fn_viscosity * solution.fn_strainrate[i, j] + + difference = V.sample(solution, deviator - constitutive, points) + magnitude = V.sample(solution, constitutive, points) + + worst = max(worst, float(np.max(np.abs(difference)))) + scale = max(scale, float(np.max(np.abs(magnitude)))) + + return worst / max(scale, 1e-300) + + +def momentum_with_flipped_bodyforce(solution, points): + r"""The negative control: $\nabla\cdot\sigma - \mathbf f$. + + A solution driven entirely through its boundary has no body force to flip, + so this returns the same number as the un-flipped residual. That is a + property of the problem rather than a gap in the check, and the note says + so where it happens. + """ + + original = solution.fn_bodyforce + solution.fn_bodyforce = sympy.Matrix([[-component for component in original]]) + + try: + return V.momentum_residual(solution, points) + finally: + solution.fn_bodyforce = original + + +def main(): + meshes = build_meshes() + rows = [] + + for name, cls in stokes_solutions(): + solution = cls(meshes[cls.dim]) + + if not hasattr(solution, "fn_stress"): + continue + + points = solution.sample_points(10) + + rows.append( + ( + name, + "deviatoric" if cls.stress_is_deviatoric else "total", + trace_check(solution, points), + constitutive_check(solution, points), + V.momentum_residual(solution, points), + momentum_with_flipped_bodyforce(solution, points), + V.incompressibility_residual(solution, points), + V.strainrate_consistency(solution, points), + ) + ) + + header = ( + "solution", + "source stress", + "tr(sigma) + d p", + "sigma + pI - 2 eta edot", + "momentum + f", + "momentum - f", + "div u", + "edot vs grad u", + ) + + print("| " + " | ".join(header) + " |") + print("|" + "---|" * len(header)) + + for row in rows: + cells = [row[0], row[1]] + [f"{value:.1e}" for value in row[2:]] + print("| " + " | ".join(cells) + " |") + + +if __name__ == "__main__": + main() diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml new file mode 100644 index 0000000..6e1d03c --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -0,0 +1,33 @@ +# Validated in CI against schemas/article-metadata.schema.json. +# `pixi run validate` checks this and the cross-file invariants a schema cannot +# express -- that the article file is named .md, that canonical_path +# matches the slug, and that no legacy DOI is ever paired with a new registrant. +id: UWTN 2026-012 +slug: testing-a-solver-against-exact-solutions +title: Testing a solver against exact solutions +article_type: technical-note +status: draft +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliation: Australian National University +publication_date: null +version: 1.0.0 +# The deposit writes archive_doi and repository_record_id when the note is +# published; leave them out until then. `doi` and `doi_registrant` were here +# once and are not fields the schema knows -- every note made from this +# template failed `pixi run validate` on all three of them. +license: CC-BY-4.0 +canonical_path: /testing-a-solver-against-exact-solutions/ +legacy_paths: [] +# Facets, from vocabulary.yml. Both keys must be present even when empty: a +# note with no subject is normal -- many are purely about method. +subjects: +methods: + - benchmarks-validation + - solvers + - symbolic-codegen +ghost_tags: + - Underworld Code +figures: 0 +source: native diff --git a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md new file mode 100644 index 0000000..f009cfd --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -0,0 +1,388 @@ +--- +title: Testing a solver against exact solutions +description: >- + Underworld3 ships thirteen exact Stokes solutions and checks them by + differentiating each one against the momentum balance rather than by + comparing it to the kernel it came from. Doing it that way found four + defects, two of them in published sources that have been vendored and + reused for twenty years. +date: 2026-08-16 +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliations: + - Australian National University +license: CC-BY-4.0 +keywords: + - Underworld Code + - benchmarks + - verification +exports: + - format: typst + logo: ../../static/uwtn-logo.png + series: "Underworld Technical Notes" + origin_url: https://www.underworldcode.org/testing-a-solver-against-exact-solutions/ + template: ../../templates/pdf + output: testing-a-solver-against-exact-solutions.pdf + article_id: UWTN 2026-012 + article_version: 1.0.0 + software_version: underworld3 0.0.0 +--- +A geodynamics solver is hard to test because the thing it computes is the thing +you do not know. Grid refinement tells you a calculation is converging, but not +what it is converging to. Comparing two codes tells you they agree, which is +worth something and is not correctness. The one check that settles the question +is a problem whose answer can be written down. + +Underworld3 ships thirteen of those, as `uw.analytic`: the Velic family of +manufactured and semi-analytic Stokes solutions that the geodynamics community +has used for two decades, plus an elliptical inclusion, a cylindrical annulus +flow, and a set of scalar transport solutions. They cover a viscosity jump, an +exponentially varying viscosity, a laterally oscillating one, a power-law +rheology, a dense block in three dimensions. + +This note is about the part that turned out to matter more than the solutions +themselves: how you check that an exact solution is exact. We wrote the checks +so that they consult no reference — not the paper, not the C kernel the +transcription came from, not the solver being tested. Each one is formed from +the solution's own symbolic fields by differentiation. That decision found four +defects. Two of them are in sources that have been vendored, copied and reused +since the 1990s, and one of those is invisible at the parameter value everyone +runs. + +## An exact solution is only as good as the check you apply to it + +There are three ways to test a transcribed solution, and they are not equally +strong. + +**Compare it against the kernel it was transcribed from.** This is the obvious +one and it is a good test of the transcription. It is not a test of the +mathematics at all: if the kernel is wrong, agreement is the wrong answer. Every +defect below in a published source would have passed this check with a residual +at machine precision. + +**Compare a solve against the solution.** This is what the solutions are for, +and it is how you find solver bugs. It cannot find a bug in the solution, +because a wrong exact answer and a wrong solver produce the same symptom: a +number that does not go to zero. + +**Differentiate the solution and substitute it into the equation it claims to +solve.** This consults nothing. If $\mathbf u$, $p$, $\sigma$, $\eta$ and +$\mathbf f$ are what the solution says they are, then + +$$ +\nabla\cdot\sigma + \mathbf f = 0 +\qquad\text{and}\qquad +\nabla\cdot\mathbf u = 0 +\qquad\text{and}\qquad +\sigma + p\mathbf I = 2\eta\dot\varepsilon +$$ + +hold pointwise, and there is nowhere for an error to hide. A solution that +faithfully reproduces a defective source fails this check, which is the entire +reason for preferring it. + +We call these the oracle-free gates. They are cheap — symbolic differentiation +of expressions that already exist, sampled at a handful of points — and they run +for every solution in the family on every commit. + +## One convention, and it belongs to the solver + +Before any of that means anything, the family has to agree on what the symbols +mean. Much of the benchmark literature writes the momentum balance with the body +force on the other side, or reports pressure positive in tension. A suite that +silently adopted a paper's convention would report a solver bug that was really +a sign disagreement. + +Underworld3's own Stokes solver fixes the conventions, and the suite exists to +validate that solver, so the solver wins wherever a source disagrees: + +| quantity | convention | +|---|---| +| total stress | $\sigma = \tau - p\mathbf I$ | +| pressure | positive in compression | +| momentum balance | $\nabla\cdot\sigma + \mathbf f = 0$ | +| strain rate | $\dot\varepsilon = \tfrac12(\nabla\mathbf u + \nabla\mathbf u^{T})$ | +| boundary normal traction | $\sigma_{nn}$ along the domain's outward normal | + +The published sources do not agree with each other on all of this, and the +disagreement is absorbed at exactly one declared boundary. Each solution carries +a `stress_is_deviatoric` flag, honoured only inside the base class's `set_fields` +method, which builds the exposed stress. Four of the thirteen kernels publish the +deviator $\tau$ and the rest publish the total $\sigma$; downstream of +`set_fields` there is one convention. + +Two things *are* uniform across every source we vendored, and neither is stated +in most of the files. Both had to be measured, by finite-differencing each +kernel's own published stress against its own published body force: the momentum +sign, $\nabla\cdot\sigma + \mathbf f = 0$ with $\mathbf f = -\rho\hat z$ under +unit gravity; and the pressure sign, positive in compression. + +## The measurement, and the column that makes it mean something + +Every registered symbolic Stokes solution, on a coarse box, sampled at the +solution's own sample points. All quantities are relative, normalised by the +largest term being cancelled. The generating script is +[`examples/convention_audit.py`](examples/convention_audit.py). + +| solution | source stress | $\mathrm{tr}\,\sigma + d\,p$ | $\sigma + p\mathbf I - 2\eta\dot\varepsilon$ | momentum $+\mathbf f$ | momentum $-\mathbf f$ | $\nabla\!\cdot\!\mathbf u$ | $\dot\varepsilon$ vs $\nabla\mathbf u$ | +|---|---|---|---|---|---|---|---| +| EllipticalInclusion | total | 1.1e-15 | 0 | 3.6e-15 | 3.6e-15 † | 1.7e-16 | 0 | +| SolA | total | 0 | 0 | 6.0e-17 | **2.0** | 0 | 0 | +| SolB | total | 2.1e-16 | 0 | 8.7e-15 | **2.0** | 4.8e-14 | 1.7e-14 | +| SolC | total | 0 | 0 | 4.1e-16 | **1.8** | 2.8e-17 | 5.9e-16 | +| SolCx | total | 3.3e-16 | 1.0e-16 | 2.3e-16 | **1.7** | 3.5e-17 | 4.0e-15 | +| SolDA | total | 8.6e-17 | 1.1e-16 | 1.4e-15 | **2.0** | 2.4e-17 | 2.9e-15 | +| SolDB2d | deviatoric | 0 | 0 | 9.7e-17 | **0.5** | 8.9e-16 | 0 | +| SolDB3d | deviatoric | 3.9e-15 | 0 | 5.3e-17 | **2.0** | 8.9e-16 | 0 | +| SolH | total | 3.3e-16 | 0 | 3.8e-16 | **2.0** | 2.1e-17 | 2.3e-16 | +| SolKx | total | 7.4e-16 | 0 | 3.3e-16 | **2.0** | 6.9e-18 | 1.5e-15 | +| SolKz | deviatoric | 0 | 0 | 2.4e-16 | **2.0** | 5.2e-18 | 3.4e-16 | +| SolM | deviatoric | 0 | 0 | 1.7e-16 | **2.0** | 0 | 0 | +| SolNL | deviatoric | 0 | 0 | 4.7e-16 | **2.0** | 0 | 0 | + +† EllipticalInclusion is driven entirely through its boundary and has no body +force, so negating the body force is a no-op and the control cannot fire. That +is a property of the problem rather than a gap in the check — there is no +body-force sign to certify. + +The momentum $-\mathbf f$ column is the load-bearing one. It re-measures the +momentum residual with the body force negated, and it moves from $10^{-16}$ to +order unity for every solution that has a body force. Without that column the +table would be an assertion that small numbers are small. With it, the gate is +demonstrably capable of failing, and the sign it certifies is the one the solver +assembles. + +One further signature is worth reading off the table. The +$\sigma + p\mathbf I - 2\eta\dot\varepsilon$ column contains exact zeros for +most solutions and $10^{-16}$ for a few. The exact zeros are not better +agreement — they are the solutions that published only one of stress and strain +rate, so `set_fields` derived the other from precisely this identity. For those, +the check is structural rather than evidential. Only SolNL, SolDB2d and SolDB3d +publish both and are genuinely tested by it, and the independent counterpart for +everyone else is the last column, which compares +$\tfrac12(\nabla\mathbf u + \nabla\mathbf u^{T})$ built from the *velocity* +against the strain rate built from the *stress*. + +## Four defects + +Three of these are corrections to a source and one is ours. Each was found by a +gate that consults no reference, and each is recorded rather than silently +applied: the vendored kernels stay verbatim, the correction lives in the +transcription, and the defect stays visible to anyone auditing the provenance. + +### SolA's vertical normal stress is missing its viscosity + +This is the one to take away from the note, and it was live in our own suite +until we audited it. + +The kernel computes the two normal components of the total stress on adjacent +lines: + +```c + u3 = 2.0*kn*ss_z - pp; /* zz total stress */ + txx = -2.0*Z*kn*ss_z - pp; /* xx total stress */ +``` + +The $xx$ component carries the viscosity $Z$ and the $zz$ component does not. +The sibling kernel `solB.c` writes the same line with the $Z$ present, and +solA's own deviatoric $\tau_{xx}$ carries it, so this is a defect rather than a +convention. + +The error is exactly $\tau_{zz}(1-Z)/Z$, which vanishes identically at $Z = 1$. +Unit viscosity is the only case the file's own driver exercises, and it was the +default value of the viscosity parameter in our transcription — which is why the +defect had survived being vendored, transcribed and run for years. + +Measured on the transcription as it stood: + +| viscosity $Z$ | momentum residual | deviator trace | $\dot\varepsilon$ consistency | +|---|---|---|---| +| 1.0 | 0 | 0 | 0 | +| 3.0 | 2.8e-1 | 6.7e-1 | 6.7e-1 | +| 0.25 | 6.4e-1 | 3.0e+0 | 7.5e-1 | + +$6.7\times10^{-1} = |1-3|/3$ and $3.0 = |1-0.25|/0.25$: the predicted $(1-Z)/Z$, +exactly. SolB is clean at every value tested, which is the control. + +Three independent gates fire, and they are independent in a way worth spelling +out. The deviator is no longer traceless, which is a statement about +incompressibility and has nothing to do with the body force. The strain rate no +longer matches the velocity gradient. And the momentum balance fails. All three +are silent at the default. + +We restore the missing factor on the term that lost it rather than editing the +vendored source, and we declare the correction per solution. There is an easier +repair available — tracelessness gives $\sigma_{zz}$ directly from the kernel's +correct $xx$ component — and we rejected it, because it would make tracelessness +true by construction and so retire one of the three gates that caught the defect. +The test asserts both halves: that solA's published deviator is *not* traceless +at $Z = 3$, and that solB's is. + +Anyone using this kernel at a viscosity other than 1 has a wrong $\sigma_{zz}$ +and no reason to suspect it. + +### SolM's published stress uses the wrong viscosity + +The kernel declares its viscosity as $(1 + \cos(r\pi x))\eta_0 + 1$ and then +computes its stress as $2(\eta - 1)\dot\varepsilon$. The constant part of the +viscosity is missing from the stress. + +| quantity | value | +|---|---| +| momentum residual, published stress | 2.1e-1 | +| momentum residual, stress from the kernel's own $\dot\varepsilon$ | 1.7e-16 | +| published $\tau$ against $2\eta\dot\varepsilon$ | 3.3e-1 | +| published $\tau$ against $2(\eta-1)\dot\varepsilon$ | **0** | + +The last row is the argument. A transcription slip would leave a residue; an +exact zero against $2(\eta-1)\dot\varepsilon$ says the kernel computed a +self-consistent stress for the wrong viscosity, which is a defect in the source. +Everything else SolM publishes — velocity, pressure, strain rate, viscosity, body +force — is mutually consistent, so the transcription passes the strain rate to +`set_fields` and lets the stress be derived. + +This is the clearest case for a check that consults no reference. Comparing our +SolM against its own kernel reproduces the error faithfully and reports +agreement. + +### SolC publishes a density where its siblings publish a force + +Most kernels in the family negate internally: they compute $\rho$ and hand back +$+\sigma\sin\cos$ as the force. SolC accumulates the density itself, so the +transcription has to negate. + +Measured: as summed, the momentum residual is 1.8; negated, 1.6e-16. + +Incompressibility and the free-slip boundary conditions hold either way. This +sign is invisible to everything except the momentum balance, which is the +argument for having that gate at all. + +### The elliptical inclusion ignored its matrix viscosity + +This one is ours. With a matrix viscosity of 3, the momentum residual is 6.3e-1 +while the deviator trace and the strain-rate consistency stay at machine +precision. That combination localises the fault immediately: the velocity, stress +and strain rate are mutually consistent, but they are not consistent with the +momentum balance. + +The Muskhelishvili potentials are normalised to unit matrix viscosity. Under +$\eta \to \lambda\eta$ at fixed boundary velocity, Stokes flow leaves the +velocity and strain rate alone and scales the stress *and the pressure* by +$\lambda$. The construction scaled the viscosity, and so the viscous part of +$\sigma = 2\eta\dot\varepsilon - p\mathbf I$, but left the pressure at its +unit-viscosity value. The two parts of the stress were then in different units, +which no gate looking at only one of them can see. Scaling the pressure takes the +residual to 3.9e-15. + +Worth recording why this one escaped for as long as it did: the elliptical +inclusion is the only solution in the family that assigns its stress, pressure +and strain rate directly instead of going through `set_fields`, which is the one +place the stress–pressure relationship is applied. + +## A transcription can be right about the source and wrong about the array + +SolKz is not an erratum — the kernel is correct — but it is the sharpest trap in +the family for anyone transcribing afresh. + +`solKz.c` computes a deviator, converts it to the total stress, and says so: + +```c + sum5 += u5*cos(n*M_PI*x); /* pressure */ + u6 -= u5; /* get total stress */ +``` + +The array the function returns really is the total stress. But our transcriber +reads the per-mode straight-line block and stops at the first accumulation, +because the series solutions sum over modes with `+=` and the summation happens +in SymPy rather than in C. The `sum5 +=` line precedes the `u6 -= u5` line, so +what we capture is the value *before* the conversion — the deviator. The +`stress_is_deviatoric` declaration on SolKz is correct, and it describes the +transcription's cut point rather than the kernel's output. + +Measured on the transcribed fields: read as the deviator, the momentum residual +is 1.7e-16; read as the total, 6.0e-1. The deviator we capture is exactly +traceless, which is the independent confirmation. + +Two signatures settle this kind of question cheaply on any new kernel. A deviator +is traceless, so its normal components are exact negatives of each other. And +$\tau = 2\eta\dot\varepsilon$, where the strain rate is a *different output of +the same kernel*. On SolKz the shear component agreed with $2\eta\dot\varepsilon$ +to machine precision while the normal components agreed with nothing, which +located the problem in one step. + +The related trap is component order. Several kernels in this family label the +vertical velocity `u1` and the horizontal `u2`, the opposite of what a Cartesian +$(u_1, u_2)$ invites, and `solKx.c` orders its stress components `[xx, xz, zz]` +where every sibling uses `[xx, zz, xz]`. Neither is documented in the file that +does it. A swap is caught by the momentum residual only because the two +components have different functional forms; in a symmetric problem it would be +invisible. + +## What we would tell someone building the same thing + +**Write the check so it cannot consult the answer.** Every defect above was found +by differentiating the solution and substituting it into the equation. None would +have been found by comparing against the source, and the two published defects +would have been actively concealed by it. + +**Give every gate a negative control, and run it.** Not as an argument that the +gate is sound, but as a measurement in the same table as the result. The +body-force flip costs one extra evaluation and converts a column of small numbers +into evidence. + +**Test away from the defaults.** SolA's defect is identically zero at the one +viscosity its driver exercises. A parameter sweep is now part of the suite, with +a table every registered solution must appear in, precisely because the defaults +are where defects go to hide. + +**Watch for a gate that is true by construction.** If the framework derived the +strain rate from the stress, then checking the stress against the strain rate +checks the framework's arithmetic and nothing else. It is worth knowing which of +your checks are structural, and saying so, rather than counting them all as +evidence. + +We had a small version of the same mistake while writing this up. The first +guard that decided which solutions publish both quantities did it by walking the +class hierarchy for a call to `set_fields` and reading its parameter names. It +reported the elliptical inclusion as publishing both, because that solution never +calls `set_fields` at all and the walk fell through to the base class and matched +the parameter names in *its* signature. The check looked like it was inspecting +the solution and was inspecting the thing doing the inspecting. `set_fields` now +records what it was handed. + +## Using them + +The solutions are constructed on a mesh and expose SymPy expressions, so they +compose with the rest of Underworld3 directly: + +```python +import underworld3 as uw + +mesh = uw.meshing.StructuredQuadBox(elementRes=(32, 32)) +solution = uw.analytic.SolCx(mesh, eta_B=1.0e6) + +stokes = uw.systems.Stokes(mesh) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = solution.fn_viscosity +stokes.bodyforce = solution.fn_bodyforce + +solution.apply_boundary_conditions(stokes) +stokes.solve() + +print(solution.error("velocity", stokes.u)) +``` + +Each solution states its own boundary conditions, because they are part of the +problem it answers: free slip on the walls for the Velic family, the exact +velocity for the manufactured solutions, and both radii of the annulus for the +cylindrical case. Every enclosed solution removes the pressure nullspace +explicitly. Leaving that out is a failure mode these solutions exist to catch — +a direct solve on a singular saddle returns a quiet, wrong answer with an +arbitrary pressure offset, and only an exact answer exposes it. + +The full family, the transcription machinery and the validation gates are in +`src/underworld3/analytic/`, and the developer documentation for the subsystem +covers adding a solution of your own. + +
Comments
Discussion of these notes happens in GitHub Discussions, so it stays with the source and is searchable alongside it.
Read the discussionStart one
From 9487e0cef72abb0db16f3d0f32baa2f1bf80e1c4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 10:14:27 +1000 Subject: [PATCH 2/7] Classify the exact-solutions note Underworld development team with AI support from Claude Code --- classification.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/classification.yml b/classification.yml index 68d9c8e..8ceff6b 100644 --- a/classification.yml +++ b/classification.yml @@ -317,3 +317,8 @@ introducing-the-technical-notes: article_type: news subjects: [] methods: [] + +testing-a-solver-against-exact-solutions: + article_type: technical-note + subjects: [] + methods: [benchmarks-validation, solvers, symbolic-codegen] From 6c822bdeb76c4660d995f5ee8e68217a550b752d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 17:55:56 +1000 Subject: [PATCH 3/7] Renumber the exact-solutions note to UWTN 2026-015, and stop the allocator reissuing a number UWTN 2026-012 was allocated twice: to the annulus and spherical-shell benchmarks (#11) and to this note. The benchmarks note claimed it first and keeps it. The allocator read only the working tree, so a number claimed by a note in review on its own branch was invisible to it and was offered again. It now also searches every local and remote-tracking ref, which covers open pull requests as far as they have been fetched. Underworld development team with AI support from Claude Code --- .../metadata.yml | 2 +- ...esting-a-solver-against-exact-solutions.md | 2 +- scripts/new_article.py | 34 ++++++++++++++++++- 3 files changed, 35 insertions(+), 3 deletions(-) diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml index 6e1d03c..eabaaf7 100644 --- a/articles/testing-a-solver-against-exact-solutions/metadata.yml +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -2,7 +2,7 @@ # `pixi run validate` checks this and the cross-file invariants a schema cannot # express -- that the article file is named .md, that canonical_path # matches the slug, and that no legacy DOI is ever paired with a new registrant. -id: UWTN 2026-012 +id: UWTN 2026-015 slug: testing-a-solver-against-exact-solutions title: Testing a solver against exact solutions article_type: technical-note diff --git a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md index f009cfd..1e9b790 100644 --- a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -24,7 +24,7 @@ exports: origin_url: https://www.underworldcode.org/testing-a-solver-against-exact-solutions/ template: ../../templates/pdf output: testing-a-solver-against-exact-solutions.pdf - article_id: UWTN 2026-012 + article_id: UWTN 2026-015 article_version: 1.0.0 software_version: underworld3 0.0.0 --- diff --git a/scripts/new_article.py b/scripts/new_article.py index a010c81..c60324f 100644 --- a/scripts/new_article.py +++ b/scripts/new_article.py @@ -17,6 +17,7 @@ import datetime import pathlib import re +import subprocess import sys ROOT = pathlib.Path(__file__).resolve().parent.parent @@ -44,6 +45,37 @@ def load_authors(): return registry +def ids_on_branches(year): + """Article numbers claimed on a branch, merged or not. + + The rest of the allocator reads the working tree, which cannot see a number + claimed by a note still in review on its own branch. Two notes drafted in + parallel were therefore both offered the same number, and UWTN 2026-012 was + claimed twice before anyone noticed. + + Every local and remote-tracking ref is searched, so this covers open pull + requests as far as they have been fetched. A note on a branch that has never + been pushed to a remote this checkout tracks is still invisible, which is why + `pixi run validate` checks for duplicates as well. + """ + def git(*args): + return subprocess.run(("git",) + args, cwd=str(ROOT), + capture_output=True, text=True) + + refs = git("for-each-ref", "--format=%(refname)", "refs/heads", "refs/remotes") + if refs.returncode != 0: + return set() # no git, or not a checkout: the tree is all we have + + used = set() + for ref in refs.stdout.split(): + found = git("grep", "-h", "-E", r"^id: +UWTN +[0-9]{4}-[0-9]{3}", + ref, "--", "articles/*/metadata.yml") + for match in re.finditer(r"UWTN\s+(\d{4})-(\d{3})", found.stdout): + if match.group(1) == str(year): + used.add(int(match.group(2))) + return used + + def next_article_id(year): """Allocate an ID that no existing article uses. @@ -52,7 +84,7 @@ def next_article_id(year): number already present in that year keeps a new note clear of anything the backfill will produce, and an ID that has been published never moves. """ - used = set() + used = ids_on_branches(year) for meta in ARTICLES.glob("*/metadata.yml"): match = re.search(r"^id:\s*(UWTN\s+(\d{4})-(\d{3}))\s*$", meta.read_text(encoding="utf-8"), re.M) From 37aaeea9700ce82ee75f26894b8acad9786d238a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 11:16:27 +1000 Subject: [PATCH 4/7] Banner, a provenance stamp that identifies something, and mark published MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The note stamped software_version as 'underworld3 0.0.0', and the audit script said the same. That identifies nothing: uw.__version__ reports 0.0.0 for every build, which is exactly the point the sibling note's scripts make. Both now carry the commit the numbers came from, 0addec15 — verified by re-running convention_audit.py against that build, which reproduces the note's table cell for cell. The banner is the note's own sentence about what the family covers — a viscosity jump, an exponentially varying viscosity, a laterally oscillating one — drawn as SolCx, SolKz and SolM's own fn_viscosity evaluated onto the render points. Generated from the solutions rather than a stock photograph, so there is nobody to credit. Status published, dated 2026-08-17, following UWTN 2026-011 and 2026-014. Underworld development team with AI support from Claude Code --- ...dQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh | 125 +++++++++++++ ...adBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh.h5 | Bin 0 -> 52288 bytes ...inC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh | 169 ++++++++++++++++++ ...(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh.h5 | Bin 0 -> 66040 bytes .../examples/banner.py | 61 +++++++ .../examples/convention_audit.py | 5 +- .../figures/banner.png | Bin 0 -> 8247 bytes .../metadata.yml | 9 +- ...esting-a-solver-against-exact-solutions.md | 5 +- 9 files changed, 370 insertions(+), 4 deletions(-) create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh.h5 create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh.h5 create mode 100644 articles/testing-a-solver-against-exact-solutions/examples/banner.py create mode 100644 articles/testing-a-solver-against-exact-solutions/figures/banner.png diff --git a/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh new file mode 100644 index 0000000..231fa41 --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh @@ -0,0 +1,125 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 11 "Bottom" +1 12 "Top" +1 13 "Right" +1 14 "Left" +2 99999 "Elements" +$EndPhysicalNames +$Entities +4 4 1 0 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +4 1 1 0 0 +11 0 0 0 1 0 0 1 11 2 1 -2 +12 0 1 0 1 1 0 1 12 2 4 -3 +13 1 0 0 1 1 0 1 13 2 2 -4 +14 0 0 0 0 1 0 1 14 2 3 -1 +1 0 0 0 1 1 0 1 99999 4 11 13 12 14 +$EndEntities +$Nodes +9 25 1 25 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +0 1 0 +0 4 0 1 +4 +1 1 0 +1 11 0 3 +5 +6 +7 +0.2499999999994109 0 0 +0.4999999999986921 0 0 +0.7499999999993406 0 0 +1 12 0 3 +8 +9 +10 +0.7500000000003471 1 0 +0.5000000000020595 1 0 +0.2500000000010405 1 0 +1 13 0 3 +11 +12 +13 +1 0.2499999999994109 0 +1 0.4999999999986921 0 +1 0.7499999999993406 0 +1 14 0 3 +14 +15 +16 +0 0.7500000000003471 0 +0 0.5000000000020595 0 +0 0.2500000000010405 0 +2 1 0 9 +17 +18 +19 +20 +21 +22 +23 +24 +25 +0.2499999999998184 0.2500000000006332 0 +0.2500000000002257 0.5000000000012177 0 +0.2500000000006332 0.7500000000000954 0 +0.4999999999995339 0.2500000000002257 0 +0.5000000000003759 0.5000000000003757 0 +0.5000000000012175 0.7499999999998438 0 +0.7499999999995922 0.2499999999998184 0 +0.749999999999844 0.4999999999995339 0 +0.7500000000000953 0.7499999999995921 0 +$EndNodes +$Elements +5 32 1 32 +1 11 1 4 +1 1 5 +2 5 6 +3 6 7 +4 7 2 +1 12 1 4 +5 4 8 +6 8 9 +7 9 10 +8 10 3 +1 13 1 4 +9 2 11 +10 11 12 +11 12 13 +12 13 4 +1 14 1 4 +13 3 14 +14 14 15 +15 15 16 +16 16 1 +2 1 3 16 +17 1 5 17 16 +18 16 17 18 15 +19 15 18 19 14 +20 14 19 10 3 +21 5 6 20 17 +22 17 20 21 18 +23 18 21 22 19 +24 19 22 9 10 +25 6 7 23 20 +26 20 23 24 21 +27 21 24 25 22 +28 22 25 8 9 +29 7 2 11 23 +30 23 11 12 24 +31 24 12 13 25 +32 25 13 4 8 +$EndElements diff --git a/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh.h5 b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh.h5 new file mode 100644 index 0000000000000000000000000000000000000000..2416546b68bb4cc371d1682d713d3c9db380791c GIT binary patch literal 52288 zcmeHQPi$Pp8K3p8olW9^18EZ|&<%0HZE5f(o77NXH^EL^8mI|wk!U5?*bjTr+Pip{ zXww2xfkXHskyC*iQUaWiOq)YbKqU?wQY8meaX`w71LVNLl`0`2B=~+a^S!rk?6u3r zUhl5Ik$!LHdo%OR%aTQZNy2zJrZOWWCci86oKXIW( z9wC1bs;k)g4%KH$(U`r%9^KqjZ&7jG-!@j=YCsnwYA4mi)#OYe=at8brTpCZls9(5 zn=5AXvr1h{t-1&p8!5TBe{VnKLH}!X|AJR4PH?xesccENHvQvN=Kd%4@7;HRJoH$8 zG99L&U5)L9_A%Y3V($d&Wl!!p*)|>A^>(R*u81&QDo$%GWhcDiw4H(-sbKTQX3R&w z{n-rr!-KsbMd)d01}qH``#Xqqm>bvGro1NCj;dX{|Mo}_t>+KD>bd{MqNga4-G|Vbbk7rpSJwen#cUMaE;58 zSyi}=Si*oXAPfit!hkR!31GFQNSSQ^r%-2WnEulxJy!?z;nmajCEnQd_R&K&`QA+=8Oi^D~pSNAc zudkCkC%jy)IK7<*EpuUA@`&23wWhg1clMe2TK&tqPw}eH~4$VHu~F0T!FN)iTvCw{UHV2-SPMt2lF&Doxc&Sp~^X#=fRen*V4bt zkK3N@`5U-oa+fEbSGbi}!hkR!3CYai%~xZ)29=m(&@lYo8z%Gmv)AAs^E z0(@{x5&>8A15n;dKtBLw$ZqrlP=@S;k6^#PIoPjHW@o5Mr&{}UtQ{#rUa*I+z@QYw zCwb=#{kv2$?H9)C=uRfDg!^xGg4P-8v?{r8rj(&fqj83j$GBn$`x z!hkR!3CT;>MCrTatJy?%H)|kq>paRVlctsYR|- zn<~{`*W<0FMKUbuYIyzTV7S8NNwi3qAF+f1VL%uV2801&Kp1Fs27XElz@MWLiopIB z`7a223AfuEJJKX!pUwAEo+2C|VEgz>!T} zLpVf$ALNG#_Y&X-`5W(&kf0My~$bru{uLmkH8SlAva0|}e=QVzf; z@}F417s`nG0PqVR95cIx$^d*JztsXhQ0}vU50p8!24C>8-Qs%60r)`vBMbOI8LyiG z@PV9Te(<+aOz1SVL5i46i>GS44jpz3zbQ&^4Yx5E%^nk4Nryc@h*bp1~!YblAtx;FRF1Wpvm|@8)$M*K6`GE>EJvd}^d13eak;yxE#h%$LWXdNaF}o}VURFH zI6}aFf0)4WHjcA#e2wF3wZC@>29a1{y1A)!&sXvJrRlS#9&<1WPa{adIVLz8~Pfw!&64;I-D<+@-qvkjA^fn&s#-D^M%ykLXFVdKbOVhMo(n_oWajQ zr(it%qL8|BbhBv{%+}6>A6Q(zu_?7u%+T%Udun>rh4d%fI9DQuB@^cRadO zJ@XVveK{M{n`m#f{r$?T|NWDMoCS4s- zMXp_kMlU~XT*4LH;-9mokCjehlFOv|@o=TflU`)FiCDsbFdz&F1HynXAPfitEzCd< z&3!vzEj4!Iv`K6t{3x1JoAlX9vNmz;I<(0SyU&NKmeD3}gdf%->sUx1n|Ef>h(Dh{+Lgjp z%jlX6NLVKIRtHe2oy(Kxnq|n7Dq%nv5C((+VL%uV23nK>Y-aH$7~V$1e((Fd7rFZJ z6ue1E8*G2$*?9-j=u`SOUCLZomwffto6zeB*2_9;@agaZt@;SxTytw3-fbIw#;mHX z-feq1{0gf6<_ng|`qclqulW;|0fYfzKo}4PgaKhd7!U>;#lV&5tOQMRH{oW&9zrj{ zJwLYk+Qijj9m!iiQgFA(k*83r@0y#LsyY-`pSOe_ zVez|#Z@j51{L6;XV)JjAHMf-(yFBkdNL>&u7Un}NVL%uV2801&Ko}4P+KqwgZ_eWz z?)ZK?zTJ-Rwy(aX*hOB+um6%%M>3WAyS`^G>#N0BW0l`GvwFJ2^g=^@^v$o9Wf(FQG?WtjYfQfv)f`&8f+H?ELtgmuuH? z&g-Lp8J}>~GUvQ9mWP_z`EZTP6aPBHt;7-rgaKhd7!U@80bxKGXjKMqwyRo;;XOs@ XF`W0p`-(X8b(DZ}UvAuNgzkL@iPc$Z literal 0 HcmV?d00001 diff --git a/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh new file mode 100644 index 0000000..f075e4f --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh @@ -0,0 +1,169 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +7 +2 11 "Bottom" +2 12 "Top" +2 13 "Right" +2 14 "Left" +2 15 "Front" +2 16 "Back" +3 99999 "Elements" +$EndPhysicalNames +$Entities +8 12 6 1 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +4 1 1 0 0 +5 0 0 1 0 +6 1 0 1 0 +7 0 1 1 0 +8 1 1 1 0 +1 0 0 0 1 0 0 0 2 1 -2 +2 1 0 0 1 1 0 0 2 2 -4 +3 0 1 0 1 1 0 0 2 4 -3 +4 0 0 0 0 1 0 0 2 3 -1 +5 0 0 1 1 0 1 0 2 5 -6 +6 1 0 1 1 1 1 0 2 6 -8 +7 0 1 1 1 1 1 0 2 8 -7 +8 0 0 1 0 1 1 0 2 7 -5 +9 0 0 0 0 0 1 0 2 5 -1 +10 1 0 0 1 0 1 0 2 2 -6 +11 0 1 0 0 1 1 0 2 7 -3 +12 1 1 0 1 1 1 0 2 4 -8 +11 0 0 0 1 1 0 1 11 4 1 2 3 4 +12 0 0 1 1 1 1 1 12 4 5 6 7 8 +13 1 0 0 1 1 1 1 13 4 10 6 -12 -2 +14 0 0 0 0 1 1 1 14 4 9 -4 -11 8 +15 0 0 0 1 0 1 1 15 4 1 10 -5 9 +16 0 1 0 1 1 1 1 16 4 -3 12 7 11 +1 0 0 0 1 1 1 1 99999 6 15 13 16 12 14 11 +$EndEntities +$Nodes +27 27 1 27 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +0 1 0 +0 4 0 1 +4 +1 1 0 +0 5 0 1 +5 +0 0 1 +0 6 0 1 +6 +1 0 1 +0 7 0 1 +7 +0 1 1 +0 8 0 1 +8 +1 1 1 +1 1 0 1 +9 +0.4999999999986921 0 0 +1 2 0 1 +10 +1 0.4999999999986921 0 +1 3 0 1 +11 +0.5000000000020595 1 0 +1 4 0 1 +12 +0 0.5000000000020595 0 +1 5 0 1 +13 +0.4999999999986921 0 1 +1 6 0 1 +14 +1 0.4999999999986921 1 +1 7 0 1 +15 +0.5000000000020595 1 1 +1 8 0 1 +16 +0 0.5000000000020595 1 +1 9 0 1 +17 +0 0 0.5000000000020595 +1 10 0 1 +18 +1 0 0.4999999999986921 +1 11 0 1 +19 +0 1 0.5000000000020595 +1 12 0 1 +20 +1 1 0.4999999999986921 +2 11 0 1 +21 +0.5000000000003758 0.5000000000003758 0 +2 12 0 1 +22 +0.5000000000003758 0.5000000000003758 1 +2 13 0 1 +23 +1 0.499999999998692 0.4999999999986922 +2 14 0 1 +24 +0 0.5000000000020595 0.5000000000020595 +2 15 0 1 +25 +0.4999999999986921 0 0.5000000000003758 +2 16 0 1 +26 +0.5000000000020595 1 0.5000000000003756 +3 1 0 1 +27 +0.5000000000003757 0.5000000000003756 0.5000000000003759 +$EndNodes +$Elements +7 32 1 32 +2 11 3 4 +1 1 9 21 12 +2 12 21 11 3 +3 9 2 10 21 +4 21 10 4 11 +2 12 3 4 +5 5 13 22 16 +6 16 22 15 7 +7 13 6 14 22 +8 22 14 8 15 +2 13 3 4 +9 2 18 23 10 +10 10 23 20 4 +11 18 6 14 23 +12 23 14 8 20 +2 14 3 4 +13 5 17 24 16 +14 16 24 19 7 +15 17 1 12 24 +16 24 12 3 19 +2 15 3 4 +17 1 9 25 17 +18 17 25 13 5 +19 9 2 18 25 +20 25 18 6 13 +2 16 3 4 +21 3 11 26 19 +22 19 26 15 7 +23 11 4 20 26 +24 26 20 8 15 +3 1 5 8 +25 1 9 21 12 17 25 27 24 +26 17 25 27 24 5 13 22 16 +27 12 21 11 3 24 27 26 19 +28 24 27 26 19 16 22 15 7 +29 9 2 10 21 25 18 23 27 +30 25 18 23 27 13 6 14 22 +31 21 10 4 11 27 23 20 26 +32 27 23 20 26 22 14 8 15 +$EndElements diff --git a/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh.h5 b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh.h5 new file mode 100644 index 0000000000000000000000000000000000000000..091721c2ab1ad95090f9225980d1efa236652b91 GIT binary patch literal 66040 zcmeHQTWnm#8J=UW9p_3CE-?uq8%V%F2zKI-gb-pJ+ZYHC5_1pOj(3S2Y_G95#tEPT z1q4#1UepI3poLx}+VYSG)Kb(3s3=v{s;!_^QKbsE>H`myiwXs#0R6r*^Y7U+_S&s6 zzO4U|zMV7k|MSnxcaGQRJ9Ez2$D13P7B5)2z%W%*n1opxZpz2!E`1yXVU?7{-XNioEuNn>TLSf(T0; zKXoP6(apxSVtj56*^5$uy%93tXb9YK5R}-8RnE24<$k%;REF)BC3j)JG1Zycg{Jm%*Pslwj8x(R zRyhCClF}jP!jRkyn{)3ovCX2^hYlOl8J4H&qh()XV}Rae?@Rtze5Z~W&hf?ZDxW7i ztGJJrlmTTx8Bhk40cAiLPzH*Xf%kTeJu}}fWyuS*7Q5u)SEwxYBlC~G&okXO3ykx< zUA}#=nSBjAe~Bmk8xyx1zaQQu zx%4Du$^ChK=y?6a<5&8V+~-H)iCyAub>rCPaw%Wp%2S|ovzfk3cjnO1J|`_571N`0 zN73Pj(wUxgw*Tl}Cuzg+i5`AqjU8@y__~oZn|$6X!|2iVfwGQtclSWsa_MwPHl`&H z%=}Pe+Yipa%8u9Re$lIi@=|yGXek`*C6Hv;ehj_e$7?@G3`vB|e7zrC2U?74LKTLj zwRU9sdvO7Yjs5;)AImspKRl*+rukalBB7nX; zQ>_6u0j%2$(6<)Y0$dDi1$c_N9k>M80k99&T?*_3E(3TID*LDdE(hpa59|im2kROD z`ZfYh0Q`m_X3>4CKU82ihuuq{CUri8T(i z(KiWjTVl~pKgL=FoZ_$=IkCclwj@A1w;6+NX{R4!aBK@4mLn$^m$tPIE0GgyOJ8o& zMqkD`&0z_0g7N9YZQAI|n9CgIBPV2Tr(>Hq(}8`{E_H0jIV^TK6*BD-K*`8xs z=)iWYTjjuUaGQ4epX9(Xa+|jC)=Yi<>X#>=`%{;pU!*a0|4dsm9Z&Z^*=#2zpFU3s zF=(y~Cg1 zt8gA;<+K$_qs3$${qt@4<^MI7_>C4v;)#!xrSH+|M~jp1FLr4wHmUb>-%rYtOM59Z z&ySH`7N7h!;wv`E{d&Kg?LC3e2Q)QNqxviNI0^BEXV_uucP5V+hgPeAt{xYovE zD25FsvHz^=Ow;}2-7Ed=n;k(e^LbL;ALmC)%78MU3@8K2fHI&A6axc1nJocI0iN%b z0p-9vfNPjs+v7SS*WtJ>$@M#~6><%ZYl2+2<2oMKCb{m%wK}fRab1yXm|XMY`Xbls zxURSw;5y$LAPKAmxQ2K(unt%coCBN-Q~~D!=K~i28-NRejle}fHBbX=0yYD+z!u;NtWb^=_(+Xd7Cmjm^{ZlD2Z1e$;=04|GN30wtS4eSA~0jLG81NH*f z1N(seKnrjKa3gRNz|-qnfLnpvfCIqoz#Tv<&<3;v2Z0VC1;~0l*Qq;!F5ob57tjs# z0KGs4=mYKs`hfu;3mgFkfqQ@<;3#kmxEHt&2(YIE@x14riBN? z2z!aE32IZ^@h(N}WnRJ74PO;2^UD0fi0#RLe^xT=BYKW|URQ65uh)K>e-e3#NA=XF z@vbTG`-@5ThF_k*HHG8+9zbF!h7Av-nSM=y-$7%gxA>)%mwcY&uO-h#HATr-$DjnY5~o|k5(M&Wf2CebK)@xh1VUB7;4sNlsCd5 z7ykGfFOuXXpC{EQdF1G(GN2471ImChpbRJjMZiGuuI+uKLx6F^P`x%M;Ks8y>MC8EV$(n4}$Z!Ozs!<0hXc)iH5? zw4@9u1ImChpbRJj%0MwNz>`@%5yPimc)nK#$n$KIuPf{`G|@9HW!ssuxeM-YwE$z7 zjx_o6O@`Bm<}2$#sc**TH#2pVJmDrf$}dkqM`gcdC(cl_Mn_TWkCmQ#GH&&GQXLiN zM@!0pGN2471ImChpbQiP17qo^voOKQ*A~8(DygY?>8U&_=Zcm({++x@V{>L|De7t_ z(NcbS0$S>&@7jox*>Ct)2h720Xji65gUNpnRzz5NrJvCZdrNbemKm zm-nat;#an^F3-lz$>%K@_Av5Zz`mb{hTO-T(O?g}Xy;s18tjEvayv5IPz@I6Lrcnl zGN2471ImChpbQin17~7FI{^M()^32mjWzjNq7qRNH+M5sNlonuCuBm)^lGrC{relu zlJ~=rm|$LDKH=tTPd~v2&5%TU`M*tNtzQOT$uM-pmpq5{_N&sT&F75vI`L~e;iA%B zwXenZrfRP^A6ildlmTTx8Bhk40cD`r7+8Y|RRR3X<4u755uV9@$M;C>wPoma^d@O~ zwO8~zwjvJ8_H`ZV%r?7w!y-H1^1)Qr)Sv0iHoJ4@$=9mMu$PhAtMgBxA@?z7wAav| z?VO8BdmVo>w-61Lp(Vft|qQYl)HC>*5cF z(UYW^qP^toH?b#uXZCk;OIW16=D7hU>f1Z+dPPp1C!eom*w0Aq_0PWxV>X{N+UtXV z*a;Vv_S*2T_}-MGPu?{mR|4W*T2cm-0cAiLPzIC%Wk4C2tqiQiz^j2;pdM%hCSObV z+AC30U0u`mA^NgudbL;6zUIA#vZ@5P7v|FEF5w=;Fv(4c;VbW1Sa|2$d55I+|PF1zJ8zSfDdRg+k;rsvW3{}9Q#5s$b92|4$#XzWh+q3R zo9X$ZdwY41;BNphm9}L17P;4<1PM_O@pZ%?Au1!Laf3vGWyQp(k~c`mZ%!?_yt})# zK69iu)!yHg9;ovqjUCvW{cPHgW68)-Ec;e+p(SDVnjw=!(OTr$#+A-N@~eKhk=U*1Y=J!HhYE%}Vp)tvmb8)!01U^k27ZF#K*=Wp~w;ESrG;0?>@c(f|Me literal 0 HcmV?d00001 diff --git a/articles/testing-a-solver-against-exact-solutions/examples/banner.py b/articles/testing-a-solver-against-exact-solutions/examples/banner.py new file mode 100644 index 0000000..d4c0202 --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/examples/banner.py @@ -0,0 +1,61 @@ +"""Banner: three of the viscosity structures the exact solutions cover. + +Run from the repository root: + + python3 articles/testing-a-solver-against-exact-solutions/examples/banner.py + +The note says the family covers "a viscosity jump, an exponentially varying +viscosity, a laterally oscillating one". This is that sentence: SolCx's step, +SolKz's exponential gradient, SolM's lateral oscillation, each evaluated from +the solution's own symbolic `fn_viscosity` onto a mesh variable and rendered. + +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that identifies +what this came from. +""" +import numpy as np +import underworld3 as uw +import pyvista as pv + +OUT = "articles/testing-a-solver-against-exact-solutions/figures/banner.png" +RES = 64 + +pv.global_theme.allow_empty_mesh = True +pv.global_theme.background = "white" + +panels = [] +for name, build in ( + ("SolCx", lambda m: uw.analytic.SolCx(m, eta_B=1.0e4)), + ("SolKz", lambda m: uw.analytic.SolKz(m)), + ("SolM", lambda m: uw.analytic.SolM(m)), +): + mesh = uw.meshing.StructuredQuadBox(elementRes=(RES, RES)) + solution = build(mesh) + pvm = uw.visualisation.mesh_to_pv_mesh(mesh) + # Evaluate the solution's own symbolic viscosity straight onto the render + # points. No MeshVariable in between: nothing here is being solved, so + # there is nothing to project. + values = np.asarray( + uw.function.evaluate(solution.fn_viscosity, pvm.points[:, :mesh.dim]), + dtype=float).reshape(-1) + # log10, because these span four orders and a linear map shows one band. + pvm.point_data["log10_eta"] = np.log10(np.maximum(values, 1e-30)) + print("%-6s log10 eta in [%.2f, %.2f]" + % (name, pvm.point_data["log10_eta"].min(), + pvm.point_data["log10_eta"].max())) + panels.append((name, pvm)) + +pl = pv.Plotter(shape=(1, 3), window_size=(2100, 700), off_screen=True, + border=False) +for col, (name, pvm) in enumerate(panels): + pl.subplot(0, col) + pl.set_background("white") + pl.add_mesh(pvm, scalars="log10_eta", cmap="RdBu_r", show_edges=False, + show_scalar_bar=False, lighting=False) + pl.enable_parallel_projection() + pl.view_xy() + pl.reset_camera(bounds=(0.0, 1.0, 0.0, 1.0, -0.05, 0.05)) + pl.camera.zoom(1.30) +pl.screenshot(OUT) +print("wrote", OUT) diff --git a/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py b/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py index 6810eaa..94ce9da 100644 --- a/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py +++ b/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py @@ -12,7 +12,10 @@ pixi run -e amr-dev python convention_audit.py -Underworld3 0.0.0 (development, 2026-08). +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that +identifies what these numbers came from. """ import numpy as np diff --git a/articles/testing-a-solver-against-exact-solutions/figures/banner.png b/articles/testing-a-solver-against-exact-solutions/figures/banner.png new file mode 100644 index 0000000000000000000000000000000000000000..0adee8a81c49b6500393b53b5d896b99ef50dce4 GIT binary patch literal 8247 zcmds6e^66b7Jgd0VppV=+1f=QTF1fK7S`5L5J+fcOJ!P6%WM~6F z7Idhv(*@I21g32&7J-6DXaY%Fsu>#v7AjeYA(KYH5(31S@FOJId+%!!?#oLoJMGTw z<_~_n_wGCAe&;*iIrrqm&bS!w2mBr&2*P{E^U=Er!fP`@EN=5!46c;(M^g!cVcik^ zOai+?r~GH}Yd<7P_1i`gjMsv9?Mh5Nv1IKtFHAmKTwQZqU=7$i6Eq`X?`u2$=mYP3 zvu16y{?TY(<8!Z{`>fL(n?2F|^ntgQ_4li`n`ZLJ{!5+sx!eO*}89!Mcu!da3GNTi4`lP_SS~t-WKDC{HoC`st~N{I}oI%^Dk|y)AFb znt~O|mO zubQ1yk}FdhhK`6@ZKmdmbXJ&uy!K-5urQ!cb8|d>l~nXa;DjVHphe*`)gkQa$iBwu zlrdStK*Lk%tK^$f_DYM)tz(*@?Bqak`?$}l5l&m^ge_~ZOK!QW)3!>)mdU(!OHBO; zXKbIOM4LKhenhU2Tee-v77JT%T0hkuP43EiKQpYWku`TomEWYYii?WPlECrKj?1+T z8spo-3j&ewfSVs!t(M%;xQhR~#19NfmX>ioT(o zM*l#oj@>%ALi4+NnX&DC%jZ?R%*cr3Qte@L4kIMRzn*ur{TeT1*m#)BspOd(Mtpqa zvd|QhEIYYglh^*r&GBz^L0Jrj<<^acVdj*=`?RKMVq}v^&(1X8ut>ZG;fDrCKFq3; z%fwwB`Jz0HHq^g!%4QuinBHJnN3FAyja97Vh`G_Pq~b{bblaCz-SJai+_9^ixf!MK zLW3oeb$MMLw=r*QWqxZbS6G#QY@$+C`MsQ>etqHzX|HBogg~U&XP)*_u{sTGv$~XX zEJSc_Iww;c*=EqMv>sx}IZWo)A;G1kx;^Y;(xE%kHG^VKtM$!5weIA|RC=OZVDwT2 z2Uhl*W2V*L>sRLIv_(p!*x)U+%{;nd6aEu2`jIIQ6HR)<&rsvdc6)N_zE=H$z z#y3!HF4kJIX0NPkQ>n({4PGgOb60~RzNi{rX?@?ydiO9T;hC!roF$0-hnO3nM0g+F zPU^t%-9ReJ14}O62mVjm^`ma=D)dFSqslJD5mEeujdCaw1CHB3R;oX^uMxbddhO5> z^yq`Bp``XiCu)}v35VBg0hc#_VwamTUm|M#(pQ0%sO>JI)YP|(sJRvOCRnlm`LG{f z^H>m35WC`=UlBw>Ay^GRQDthou{vL3=heUd3ta5}n|-TjEB?p5!9{Pd8h+-2MQA z-1wy({`etqySqdHNeF`m0Z8j7VbN`%4P{SL65ts@IB%mx%DeLy`j1@rLka5 zRQTWu@}(c$J8PQLX$>WEJi|ajYr(y{2X=HR&ClL*6Z(4lIeXgT%szEJ2GKYmV=1|% zU5F4(+fdFq1qIie?2AAUsCv`W<#Q*?%@>#(AqqNnIQR)Z zWX0%n2`nDMR0V|>NX}hhaCDzgo?S!y6Z|8oG>g@yo~BvP=^2XnhcEg7al}m>pn)Ac zeh*!D1&L(;f8bnC(?L%fy1Nz>YAB(0p2L{(|G+2*7YlOQ0$|_>A{tN5he<7D&95P4 zBpDpnFj?DhiRw9uPZCv4u%N;gV!pUS1Lda>8hA{^4%+cC6?=87pK_HNR%d^eu;Z`lgr7h)VOqu9p! zCwc@nKIzfUw;+dr6gwP8p4C!o3a$sm6fY3asClY(pn`dpkkoP=c31`Jj?KbE*vz|3 zxr()FZ|JFt3KngE#)twpsL|~v_9YY)XWhc1Ny9Wg^3qLdkAFB&OuCWQU0}vnK(jMUKEj#HNa4bKJe+v=e{OK)9{;?jR88|I9djx(lmu!+hGSf4_#AM6;@FEp z2xE8KUpU4e?&X5*@5dxD*02U*QsDF!lPYRch9dbX#amnwn{GfhIXeeb<(bXsP~oV^ zm?UgGnT9gFd0r3srL^7OuqrM-E?5;=ui|DJKULrt06OR{SQS++iRK>ZK7Fo*N-{rrCcc;t1;?irfnReFn*>ihK zhWZV*X{|*6>Y@zgBYP9T6Z@ry7PlsvwEC^wsq^wz11GNJTKgWqA1-eDO`-3X>egVJ z@x2hH@!4dln%8L?PhHGh_v6WYc~G@Tq?y$*#OzT)NQ36UMEeiRp>Owqq1SYB z!g~9>N*p@Dy<48rX*~g7T+qepSaw^>9uUveF(Ql#ZR(uqWtGf&qM+nT#SJc$&XhR%a+E@t%!w=Khi>1!ZgMOzSn>lHOA%39?VU3&(^TO@exXO5(=A*K$-3$s zbHU9wm#TkxBWOP7cD{y~@AMh%L%xZEI=&OSQPvJr?cI&!$ZGEx;^8t$HfYYJut8hU z>yUbrJ>bt}qKA2oZnm@Bql)y>sfb5T6w)-B8?RG2kE8!7-5PJ?n47-QWx@)nw|V4< zhtpYtUk6DHysBogEwG>mAZyeUUBYSV$pzUN_GpK0;ItmS-d(cpdx1{hDoC6i?*1~E z70DSu69yb-JY0D<(Q;;mN6yemj~^9BtanK(p=bs=pP$~B>nWUYU`g#>JkCPa!*3Gi zP5x0&d9qOitG>HQ8+FJILQ?I~2tTZlF3ZM6l&jd$p~s)Z%%T2h|4;+w&Af22{&^~b zpJ#~an9Tno+|v-0_KS!A!cG=!!ZlF#ki7AatiRY+qh7F}&oM1<5Z3;Eks}=)8T0@0 eNYN*@`xkXzZO94UFbXJ0?D&0Lbp5k`dgtE)eQ;y| literal 0 HcmV?d00001 diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml index eabaaf7..6c3a0b5 100644 --- a/articles/testing-a-solver-against-exact-solutions/metadata.yml +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -6,12 +6,12 @@ id: UWTN 2026-015 slug: testing-a-solver-against-exact-solutions title: Testing a solver against exact solutions article_type: technical-note -status: draft +status: published authors: - name: Louis Moresi orcid: 0000-0003-3685-174X affiliation: Australian National University -publication_date: null +publication_date: 2026-08-17 version: 1.0.0 # The deposit writes archive_doi and repository_record_id when the note is # published; leave them out until then. `doi` and `doi_registrant` were here @@ -30,4 +30,9 @@ methods: ghost_tags: - Underworld Code figures: 0 +# Generated from the solutions themselves rather than a stock photograph, so +# there is nobody to credit. `figures` counts figures in the body; the banner +# is not one. +banner: figures/banner.png +banner_credit: null source: native diff --git a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md index 1e9b790..f606aef 100644 --- a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -13,6 +13,7 @@ authors: affiliations: - Australian National University license: CC-BY-4.0 +banner: figures/banner.png keywords: - Underworld Code - benchmarks @@ -26,8 +27,10 @@ exports: output: testing-a-solver-against-exact-solutions.pdf article_id: UWTN 2026-015 article_version: 1.0.0 - software_version: underworld3 0.0.0 + software_version: underworld3 development @ 0addec15 --- +
+ A geodynamics solver is hard to test because the thing it computes is the thing you do not know. Grid refinement tells you a calculation is converging, but not what it is converging to. Comparing two codes tells you they agree, which is From 827bd7f380a362a757e17363839a1f83f122da8a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 11:39:22 +1000 Subject: [PATCH 5/7] Add the Barr & Houseman erratum, with Thyagarajulu Gollapalli as co-author MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The note previously stopped at the Velic family. The Barr & Houseman 1996 appendix contains the same kind of problem and the method resolves it the same way, so it belongs here rather than being held back. Their half-integer sine terms of u_theta carry one sign in the boundary datum (A8b) and the opposite in the solution (A9b), with the same mismatch in the plane-stress pair. Two printed equations disagree and reading them again will not say which is right. Incompressibility does: div u = 0 forces g' = -(3/2) A f, which integrates to (A8b)'s sign, and flipping it back makes the divergence non-zero purely in the fault's own half-integer modes. That is the note's argument applied to a paper rather than to a C kernel. The transcription is in review as underworld3#550 and follows the implementation Thyagarajulu Gollapalli has been using for fault benchmarking, so the note carries their name. Also, per review: a footnote on the SolA row. It reads clean because uw.analytic carries the correction, while the published kernel still does not — the row certifies the transcription, not the source. Underworld development team with AI support from Claude Code --- .../metadata.yml | 3 + ...esting-a-solver-against-exact-solutions.md | 60 ++++++++++++++++++- 2 files changed, 60 insertions(+), 3 deletions(-) diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml index 6c3a0b5..916f94f 100644 --- a/articles/testing-a-solver-against-exact-solutions/metadata.yml +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -11,6 +11,9 @@ authors: - name: Louis Moresi orcid: 0000-0003-3685-174X affiliation: Australian National University + - name: Thyagarajulu Gollapalli + orcid: 0000-0001-9394-4104 + affiliation: Monash University publication_date: 2026-08-17 version: 1.0.0 # The deposit writes archive_doi and repository_record_id when the note is diff --git a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md index f606aef..01d6055 100644 --- a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -5,13 +5,18 @@ description: >- differentiating each one against the momentum balance rather than by comparing it to the kernel it came from. Doing it that way found four defects, two of them in published sources that have been vendored and - reused for twenty years. + reused for twenty years, and settled a contradiction between two printed + equations in a third. date: 2026-08-16 authors: - name: Louis Moresi orcid: 0000-0003-3685-174X affiliations: - Australian National University + - name: Thyagarajulu Gollapalli + orcid: 0000-0001-9394-4104 + affiliations: + - Monash University license: CC-BY-4.0 banner: figures/banner.png keywords: @@ -51,7 +56,8 @@ transcription came from, not the solver being tested. Each one is formed from the solution's own symbolic fields by differentiation. That decision found four defects. Two of them are in sources that have been vendored, copied and reused since the 1990s, and one of those is invisible at the parameter value everyone -runs. +runs. The same method settles a place where a published appendix contradicts +itself, without anyone having to decide what the printed sign was meant to be. ## An exact solution is only as good as the check you apply to it @@ -131,7 +137,7 @@ largest term being cancelled. The generating script is | solution | source stress | $\mathrm{tr}\,\sigma + d\,p$ | $\sigma + p\mathbf I - 2\eta\dot\varepsilon$ | momentum $+\mathbf f$ | momentum $-\mathbf f$ | $\nabla\!\cdot\!\mathbf u$ | $\dot\varepsilon$ vs $\nabla\mathbf u$ | |---|---|---|---|---|---|---|---| | EllipticalInclusion | total | 1.1e-15 | 0 | 3.6e-15 | 3.6e-15 † | 1.7e-16 | 0 | -| SolA | total | 0 | 0 | 6.0e-17 | **2.0** | 0 | 0 | +| SolA ‡ | total | 0 | 0 | 6.0e-17 | **2.0** | 0 | 0 | | SolB | total | 2.1e-16 | 0 | 8.7e-15 | **2.0** | 4.8e-14 | 1.7e-14 | | SolC | total | 0 | 0 | 4.1e-16 | **1.8** | 2.8e-17 | 5.9e-16 | | SolCx | total | 3.3e-16 | 1.0e-16 | 2.3e-16 | **1.7** | 3.5e-17 | 4.0e-15 | @@ -149,6 +155,11 @@ force, so negating the body force is a no-op and the control cannot fire. That is a property of the problem rather than a gap in the check — there is no body-force sign to certify. +‡ SolA's row is clean because `uw.analytic` carries the correction described +below. There is an erratum in the published kernel, and it is not repaired +there — anyone reading `solA.c` still has it. What the row certifies is the +transcription, which is consistent. + The momentum $-\mathbf f$ column is the load-bearing one. It re-measures the momentum residual with the body force negated, and it moves from $10^{-16}$ to order unity for every solution that has a body force. Without that column the @@ -283,6 +294,49 @@ inclusion is the only solution in the family that assigns its stress, pressure and strain rate directly instead of going through `set_fields`, which is the one place the stress–pressure relationship is applied. +## An erratum settled without adjudicating the paper + +The four above are in the Velic family. The same method settles a question in a +different solution, and settles it in the way this note is arguing for: by +mathematics rather than by deciding what a scanned minus sign was meant to say. + +Barr & Houseman (1996, *GJI* **125**, 473–490) give, in their Appendix, a linear +plane-strain solution for a fault terminating inside a viscous medium — an +internal boundary carrying zero shear traction, with the tip in the interior. +It is a genuine absolute standard for a fault calculation, which otherwise has +only other discretisations to be measured against. + +The half-integer sine terms of $u_\theta$ appear with **one sign in the paper's +boundary datum (A8b) and the opposite sign in its solution (A9b)**, and the same +mismatch appears in the plane-stress pair, (A12b) against (A13b). Two printed +equations disagree, and no amount of reading them more carefully will say which +is right. + +Incompressibility does. For $u_r = A\sqrt{R}\,f(\theta)$ and +$u_\theta = \sqrt{R}\,g(\theta)$, requiring $\nabla\cdot\mathbf u = 0$ forces + +$$g'(\theta) = -\tfrac{3}{2} A f(\theta),$$ + +which integrates to (A8b)'s sign. The negative control is the same shape as the +one in the table above: flip the sign back and the divergence becomes +$(0.75\cos(\theta/2) + 2.25\cos(3\theta/2))/\sqrt{r}$ — non-zero purely in the +half-integer terms, which are the fault's own modes. Nothing here reads the +paper's intent. + +Two further things about this solution are conventions rather than errors, and +both matter to anyone comparing against it. Its pressure is **extension +positive**, so its force balance is $\partial_j\tau_{ij} + \partial_i p = 0$; +negate to compare against a compression-positive solver such as ours. And only +the plane-strain solution is a Stokes benchmark — the thin-viscous-sheet +solution of (A10)–(A13) has non-zero in-plane divergence, which is a different +equation set. + +The transcription is in review as +[underworld3#550](https://github.com/underworldcode/underworld3/pull/550) and +follows the implementation Thyagarajulu Gollapalli has been using for fault +benchmarking; a second independent implementation is the strongest check +available on a solution of this kind. + ## A transcription can be right about the source and wrong about the array SolKz is not an erratum — the kernel is correct — but it is the sharpest trap in From b7f24415f84fae53503c7e04a3d0a919fab81c9e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 21:17:53 +1000 Subject: [PATCH 6/7] The faulted medium's second erratum, and the sweep that separates it from a convention MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Underworld3 #550 is merged, so FaultedMedium is registered and the audit script picks it up on its own: the measurement table gains a row, and the note no longer describes the transcription as in review. Two additions from Thyagarajulu Gollapalli's reading of the papers and his independent implementation. Printed (A9b) carries a second typo besides the sign: cos(3 theta / 2) where sin(3 theta / 2) belongs. The form implemented reproduces (A8b) at r = R0, is divergence-free and satisfies both momentum components symbolically with U0, R0 and eta free; printed (A9b) fails boundary matching and incompressibility together, so neither reading is an alternative convention. The pressure sign is where a convention and an error look alike from one run. A resolution sweep separates them: against the negated pressure the error falls with the mesh (13.1%, 6.9%, 3.3% at h = 0.20, 0.10, 0.05) and against the paper's own sign it sits near 199% at every resolution. Differentiating the transcription agrees without a solver — momentum residual 3.6e-16 negated against 1.06 as printed. FaultedMedium shares EllipticalInclusion's footnote: both are boundary-driven with no body force, so the negated-body-force control cannot fire and the momentum residual carries the whole weight. That is why the sign measurement is stated separately rather than left to the table. Version 1.1.0. Underworld development team with AI support from Claude Code --- .../metadata.yml | 2 +- ...esting-a-solver-against-exact-solutions.md | 40 ++++++++++++++----- 2 files changed, 32 insertions(+), 10 deletions(-) diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml index 916f94f..e909d83 100644 --- a/articles/testing-a-solver-against-exact-solutions/metadata.yml +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -15,7 +15,7 @@ authors: orcid: 0000-0001-9394-4104 affiliation: Monash University publication_date: 2026-08-17 -version: 1.0.0 +version: 1.1.0 # The deposit writes archive_doi and repository_record_id when the note is # published; leave them out until then. `doi` and `doi_registrant` were here # once and are not fields the schema knows -- every note made from this diff --git a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md index 01d6055..0acbd06 100644 --- a/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -137,6 +137,7 @@ largest term being cancelled. The generating script is | solution | source stress | $\mathrm{tr}\,\sigma + d\,p$ | $\sigma + p\mathbf I - 2\eta\dot\varepsilon$ | momentum $+\mathbf f$ | momentum $-\mathbf f$ | $\nabla\!\cdot\!\mathbf u$ | $\dot\varepsilon$ vs $\nabla\mathbf u$ | |---|---|---|---|---|---|---|---| | EllipticalInclusion | total | 1.1e-15 | 0 | 3.6e-15 | 3.6e-15 † | 1.7e-16 | 0 | +| FaultedMedium | total | 2.6e-16 | 0 | 3.3e-16 | 3.3e-16 † | 1.3e-15 | 5.1e-16 | | SolA ‡ | total | 0 | 0 | 6.0e-17 | **2.0** | 0 | 0 | | SolB | total | 2.1e-16 | 0 | 8.7e-15 | **2.0** | 4.8e-14 | 1.7e-14 | | SolC | total | 0 | 0 | 4.1e-16 | **1.8** | 2.8e-17 | 5.9e-16 | @@ -150,10 +151,12 @@ largest term being cancelled. The generating script is | SolM | deviatoric | 0 | 0 | 1.7e-16 | **2.0** | 0 | 0 | | SolNL | deviatoric | 0 | 0 | 4.7e-16 | **2.0** | 0 | 0 | -† EllipticalInclusion is driven entirely through its boundary and has no body -force, so negating the body force is a no-op and the control cannot fire. That -is a property of the problem rather than a gap in the check — there is no -body-force sign to certify. +† EllipticalInclusion and FaultedMedium are driven entirely through their +boundaries and have no body force, so negating the body force is a no-op and the +control cannot fire. That is a property of the problem rather than a gap in the +check — there is no body-force sign to certify. For FaultedMedium the momentum +residual carries the whole weight, which is what makes the pressure-sign +measurement below worth stating separately. ‡ SolA's row is clean because `uw.analytic` carries the correction described below. There is an erratum in the published kernel, and it is not repaired @@ -323,6 +326,13 @@ $(0.75\cos(\theta/2) + 2.25\cos(3\theta/2))/\sqrt{r}$ — non-zero purely in the half-integer terms, which are the fault's own modes. Nothing here reads the paper's intent. +Printed (A9b) carries a second defect, of a kind mathematics settles just as +cleanly. It has $\cos(3\theta/2)$ where $\sin(3\theta/2)$ belongs. The form we +implement reproduces (A8b) exactly at $r = R_0$, is divergence-free, and +satisfies both momentum components symbolically with $U_0$, $R_0$ and $\eta$ +left free; printed (A9b) fails boundary matching and incompressibility together. +Neither reading of it is an alternative convention. + Two further things about this solution are conventions rather than errors, and both matter to anyone comparing against it. Its pressure is **extension positive**, so its force balance is $\partial_j\tau_{ij} + \partial_i p = 0$; @@ -331,11 +341,23 @@ the plane-strain solution is a Stokes benchmark — the thin-viscous-sheet solution of (A10)–(A13) has non-zero in-plane divergence, which is a different equation set. -The transcription is in review as -[underworld3#550](https://github.com/underworldcode/underworld3/pull/550) and -follows the implementation Thyagarajulu Gollapalli has been using for fault -benchmarking; a second independent implementation is the strongest check -available on a solution of this kind. +The pressure sign is the one place where a convention and an error look alike +from a single run, and a resolution sweep tells them apart. Measured on a +Underworld3 Stokes solve over a Gmsh slit disc, against the negated pressure the +error falls with the mesh — 13.1%, 6.9% and 3.3% at $h = 0.20$, $0.10$ and +$0.05$ — and against the paper's own sign it sits at about 199% at every +resolution. An error that does not move under refinement is not a discretisation +error. Differentiating the transcription says the same thing without a solver at +all: the momentum residual is 3.6e-16 with the sign negated and 1.06 with it as +printed. + +The transcription is +[underworld3#550](https://github.com/underworldcode/underworld3/pull/550), +merged, and it follows the implementation Thyagarajulu Gollapalli has been using +for fault benchmarking. Both typos and the resolution sweep above are his; a +second independent implementation is the strongest check available on a solution +of this kind, and it is what turned a suspected convention into a measured +erratum. ## A transcription can be right about the source and wrong about the array From 6efdc0e9783cc79e18f082fc9793fc7411c38b30 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 21:28:08 +1000 Subject: [PATCH 7/7] =?UTF-8?q?Keep=20the=20exact-solutions=20note=20at=20?= =?UTF-8?q?1.0.0=20=E2=80=94=20it=20has=20not=20been=20deposited?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The version bump in the previous commit was wrong. Versioning exists to relate a new deposit to an earlier one, and this note has no earlier one: it is absent from deposit-queue.txt and carries neither archive_doi nor repository_record_id, so no DOI has been minted and there is nothing for a 1.1.0 to supersede. Publishing to the site is not depositing. Per PUBLISHING.md the deposit happens only when somebody merges the request that adds a slug to the queue, so a note can be `status: published` and still be at its first version. The content changes stand; only the number goes back. Underworld development team with AI support from Claude Code --- articles/testing-a-solver-against-exact-solutions/metadata.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/articles/testing-a-solver-against-exact-solutions/metadata.yml b/articles/testing-a-solver-against-exact-solutions/metadata.yml index e909d83..916f94f 100644 --- a/articles/testing-a-solver-against-exact-solutions/metadata.yml +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -15,7 +15,7 @@ authors: orcid: 0000-0001-9394-4104 affiliation: Monash University publication_date: 2026-08-17 -version: 1.1.0 +version: 1.0.0 # The deposit writes archive_doi and repository_record_id when the note is # published; leave them out until then. `doi` and `doi_registrant` were here # once and are not fields the schema knows -- every note made from this