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 0000000..2416546 Binary files /dev/null and b/articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh.h5 differ 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 0000000..091721c Binary files /dev/null and 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 differ 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 new file mode 100644 index 0000000..94ce9da --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/examples/convention_audit.py @@ -0,0 +1,178 @@ +"""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 + +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 +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/figures/banner.png b/articles/testing-a-solver-against-exact-solutions/figures/banner.png new file mode 100644 index 0000000..0adee8a Binary files /dev/null and b/articles/testing-a-solver-against-exact-solutions/figures/banner.png differ 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..916f94f --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/metadata.yml @@ -0,0 +1,41 @@ +# 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-015 +slug: testing-a-solver-against-exact-solutions +title: Testing a solver against exact solutions +article_type: technical-note +status: published +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 +# 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 +# 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 new file mode 100644 index 0000000..0acbd06 --- /dev/null +++ b/articles/testing-a-solver-against-exact-solutions/testing-a-solver-against-exact-solutions.md @@ -0,0 +1,467 @@ +--- +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, 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: + - 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-015 + article_version: 1.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 +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. 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 + +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 | +| 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 | +| 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 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 +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 +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. + +## 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. + +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$; +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 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 + +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
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] 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)