UWTN 2026-014: Setting Up Full Multigrid

articles/setting-up-full-multigrid/examples/banner.py

@@ -0,0 +1,61 @@
+"""Banner: the three levels a `refinement=2` mesh carries.
+
+Run from the repository root:
+
+    python3 articles/setting-up-full-multigrid/examples/banner.py
+
+The subject of the note is that a mesh built with `refinement` keeps its coarser
+ancestors, and that those are what the preconditioner works on. So the banner is
+the hierarchy itself: the gmsh base mesh, and the two refinements built from it,
+each one subdividing every cell of the one before and snapping the new boundary
+nodes back onto the bounding circles.
+
+Building three meshes at refinement 0, 1 and 2 gives exactly the three levels a
+single `refinement=2` mesh holds in `dm_hierarchy` -- same base mesh, same
+refinement callback -- and does it through the supported API rather than
+reaching into the DMPlex objects.
+
+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 underworld3 as uw
+import pyvista as pv
+
+RADIUS_INNER, RADIUS_OUTER = 0.5, 1.0
+CELL_SIZE = 0.25          # the BASE mesh; refinement takes it from there
+OUT = "articles/setting-up-full-multigrid/figures/banner.png"
+
+pv.global_theme.allow_empty_mesh = True
+pv.global_theme.background = "white"
+
+levels = []
+for refinement in (0, 1, 2):
+    mesh = uw.meshing.Annulus(radiusInner=RADIUS_INNER, radiusOuter=RADIUS_OUTER,
+                              cellSize=CELL_SIZE, qdegree=2,
+                              refinement=refinement)
+    pvm = uw.visualisation.mesh_to_pv_mesh(mesh)
+    print("refinement %d: %6d cells, %d hierarchy level(s)"
+          % (refinement, pvm.n_cells, len(mesh.dm_hierarchy)))
+    levels.append(pvm)
+
+# Wide and short: a banner is cropped hard on a narrow screen, so the three
+# panels have to read at a glance and nothing may depend on fine detail.
+pl = pv.Plotter(shape=(1, 3), window_size=(2100, 700), off_screen=True,
+                border=False)
+for col, pvm in enumerate(levels):
+    pl.subplot(0, col)
+    pl.set_background("white")
+    # Line width drops as the cells get smaller, so the finest panel reads as a
+    # texture rather than as a block of ink.
+    pl.add_mesh(pvm, color="white", show_edges=True, edge_color="#1a1a1a",
+                line_width=(1.6, 1.1, 0.7)[col], lighting=False)
+    pl.enable_parallel_projection()
+    pl.view_xy()
+    pl.reset_camera(bounds=(-1.02, 1.02, -1.02, 1.02, -0.05, 0.05))
+    # reset_camera leaves a wide margin; a banner is short and gets cropped, so
+    # the meshes have to carry the strip rather than float in it.
+    pl.camera.zoom(1.35)
+pl.screenshot(OUT)
+print("wrote", OUT)

articles/setting-up-full-multigrid/examples/layer_figure.py

@@ -0,0 +1,104 @@
+"""Figure: what happens to the comparison when the viscosity structure is
+concentrated rather than spread through the box.
+
+Run from the repository root:
+
+    python3 articles/setting-up-full-multigrid/examples/layer_figure.py
+
+The numbers are those printed by `fmg-vs-gamg.py solkz layer` and are repeated
+here so the figure has a single, checkable source. Both panels are relative to
+the same baseline -- FMG on the constant-viscosity problem -- so the two
+profiles can be read against each other. At contrast 1 the two profiles ARE the
+same problem (eta = 1 everywhere) and the runs agree to the last digit, which is
+why the curves start together.
+
+Colour carries the preconditioner and line style carries the viscosity profile,
+so neither is distinguished by colour alone. The palette is checked for
+colour-vision separation (worst pair dE 24.7, against a target of 8).
+
+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 matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+OUT = "articles/setting-up-full-multigrid/figures/layer-vs-smooth.png"
+
+CONTRAST = [1e0, 1e2, 1e4, 1e6]
+
+# Relative to FMG at contrast 1. `None` = did not converge: GAMG on the band at
+# 1e6 ran past 20 000 multigrid cycles per velocity solve without reaching the
+# tolerance, so there is no cost to plot -- reporting the work it managed to do
+# before being stopped would be reporting a cap, not a method.
+WORK = {
+    ("fmg", "smooth"):  [1.0, 1.3, 1.6, 1.8],
+    ("gamg", "smooth"): [2.2, 3.3, 4.0, 4.0],
+    ("fmg", "band"):    [1.0, 1.9, 2.1, 2.5],
+    ("gamg", "band"):   [2.2, 8.1, 28.1, None],
+}
+TIME = {
+    ("fmg", "smooth"):  [1.0, 1.1, 1.2, 1.3],
+    ("gamg", "smooth"): [1.2, 1.5, 1.7, 1.6],
+    ("fmg", "band"):    [1.0, 1.3, 1.4, 1.5],
+    ("gamg", "band"):   [1.2, 2.7, 7.7, None],
+}
+
+COLOUR = {"fmg": "#2a78d6", "gamg": "#eb6834"}      # validated categorical 1, 2
+STYLE = {"smooth": "-", "band": "--"}
+INK = "#0b0b0b"
+INK_MUTED = "#52514e"
+GRID = "#e4e3df"
+
+
+def panel(ax, data, title):
+    for (pc, profile), y in data.items():
+        xs = [x for x, v in zip(CONTRAST, y) if v is not None]
+        ys = [v for v in y if v is not None]
+        ax.plot(xs, ys, STYLE[profile], color=COLOUR[pc], linewidth=2.0,
+                marker="o", markersize=5, markerfacecolor="white",
+                markeredgewidth=1.6, zorder=3,
+                label="%s, %s" % (pc.upper(), profile))
+        # An open marker at the last point that converged, with the failure
+        # named rather than left as a gap the reader has to interpret.
+        if y[-1] is None:
+            # Say what the gap means. An unexplained stop reads as missing data.
+            ax.annotate("does not converge\nbeyond this point",
+                        xy=(xs[-1], ys[-1]), xytext=(10, -2),
+                        textcoords="offset points", ha="left", va="top",
+                        fontsize=8.5, color=INK_MUTED, linespacing=1.35)
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlabel("viscosity contrast", fontsize=9.5, color=INK_MUTED)
+    ax.set_title(title, fontsize=10.5, color=INK, pad=10)
+    ax.grid(True, which="major", color=GRID, linewidth=0.8, zorder=0)
+    ax.set_axisbelow(True)
+    for side in ("top", "right"):
+        ax.spines[side].set_visible(False)
+    for side in ("left", "bottom"):
+        ax.spines[side].set_color(GRID)
+    ax.tick_params(colors=INK_MUTED, labelsize=9)
+    # Same tick VALUES in both panels so the two can be read against each
+    # other, but each panel framed on its own data rather than padded to a
+    # shared top -- the time differences are genuinely smaller and should look
+    # it without half the panel being empty.
+    ax.set_yticks([1, 1.5, 2, 3, 5, 10, 20, 30])
+    ax.set_yticklabels(["1", "1.5", "2", "3", "5", "10", "20", "30"])
+    finite = [v for y in data.values() for v in y if v is not None]
+    ax.set_ylim(min(finite) * 0.88, max(finite) * 1.45)
+    ax.minorticks_off()
+
+
+fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.1), sharex=True)
+fig.patch.set_facecolor("white")
+panel(axes[0], WORK, "work per unknown, relative to FMG at contrast 1")
+panel(axes[1], TIME, "time per unknown, relative to FMG at contrast 1")
+# handlelength long enough that the dashed entries are visibly dashed -- the
+# line style is half the encoding, so a legend that hides it breaks identity.
+axes[0].legend(frameon=False, fontsize=9, labelcolor=INK_MUTED,
+               loc="upper left", handlelength=3.4, borderaxespad=0.2)
+fig.tight_layout()
+fig.savefig(OUT, dpi=200, facecolor="white")
+print("wrote", OUT)

articles/setting-up-full-multigrid/examples/relax_annulus_figure.py

@@ -0,0 +1,72 @@
+"""Figure: what the coarse mesh leaves in the fine one, and what relax() does.
+
+Run from the repository root:
+
+    python3 articles/setting-up-full-multigrid/examples/relax_annulus_figure.py
+
+The annulus is the clearer subject than a box: its base mesh is coarse relative
+to the curvature, so the coarse cells' edges survive as visible seams across the
+refined mesh. It also shows that the boundary snapping holds -- every node that
+starts on a bounding circle is still exactly on it afterwards, because relax()
+moves interior coordinates only.
+
+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 underworld3 as uw
+import pyvista as pv
+
+RADIUS_INNER, RADIUS_OUTER = 0.5, 1.0
+CELL_SIZE = 0.25          # the BASE mesh; refinement takes it from there
+REFINEMENT = 2
+OUT = "articles/setting-up-full-multigrid/figures/relax-annulus.png"
+
+pv.global_theme.allow_empty_mesh = True
+pv.global_theme.background = "white"
+
+mesh = uw.meshing.Annulus(radiusInner=RADIUS_INNER, radiusOuter=RADIUS_OUTER,
+                          cellSize=CELL_SIZE, qdegree=2, refinement=REFINEMENT)
+before = uw.visualisation.mesh_to_pv_mesh(mesh)
+print("cells %d, hierarchy levels %d" % (before.n_cells, len(mesh.dm_hierarchy)))
+
+
+def report(pvm, tag):
+    """Element quality, and whether the bounding circles are still exact."""
+    cells = pvm.cells.reshape(-1, 4)[:, 1:]
+    p = pvm.points[:, :2]
+    a = np.linalg.norm(p[cells[:, 1]] - p[cells[:, 2]], axis=1)
+    b = np.linalg.norm(p[cells[:, 2]] - p[cells[:, 0]], axis=1)
+    c = np.linalg.norm(p[cells[:, 0]] - p[cells[:, 1]], axis=1)
+    s = 0.5 * (a + b + c)
+    area = np.sqrt(np.maximum(s * (s - a) * (s - b) * (s - c), 0.0))
+    q = 2.0 * (area / np.maximum(s, 1e-30)) / np.maximum(
+        a * b * c / np.maximum(4.0 * area, 1e-30), 1e-30)
+    r = np.linalg.norm(p, axis=1)
+    near = (np.abs(r - RADIUS_OUTER) < 2e-3) | (np.abs(r - RADIUS_INNER) < 2e-3)
+    exact = (np.abs(r - RADIUS_OUTER) < 1e-6) | (np.abs(r - RADIUS_INNER) < 1e-6)
+    print("%-7s q median %.3f  10th pct %.3f  min %.3f | on the circle: %d/%d"
+          % (tag, np.median(q), np.percentile(q, 10), q.min(),
+             exact.sum(), near.sum()))
+
+
+report(before, "before")
+mesh.relax()
+after = uw.visualisation.mesh_to_pv_mesh(mesh)
+report(after, "after")
+
+pl = pv.Plotter(shape=(1, 2), window_size=(1900, 950), off_screen=True, border=False)
+for col, (m, title) in enumerate(((before, "Annulus, refinement=2"),
+                                  (after, "after relax()"))):
+    pl.subplot(0, col)
+    pl.set_background("white")
+    pl.add_mesh(m, color="white", show_edges=True, edge_color="#1a1a1a",
+                line_width=1.1, lighting=False)
+    pl.add_text(title, position="upper_edge", font_size=14, color="black")
+    pl.enable_parallel_projection()
+    pl.view_xy()
+    pl.reset_camera(bounds=(-1.05, 1.05, -1.05, 1.12, -0.05, 0.05))
+pl.screenshot(OUT)
+print("wrote", OUT)

articles/setting-up-full-multigrid/metadata.yml

@@ -0,0 +1,38 @@
+# 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 <slug>.md, that canonical_path
+# matches the slug, and that no legacy DOI is ever paired with a new registrant.
+id: UWTN 2026-014
+slug: setting-up-full-multigrid
+title: Setting Up Full Multigrid
+article_type: technical-note
+status: published
+authors:
+  - name: Louis Moresi
+    orcid: 0000-0003-3685-174X
+    affiliation: Australian National 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: /setting-up-full-multigrid/
+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:
+  - mantle-convection
+methods:
+  - solvers
+  - meshing
+  - parallel-hpc
+ghost_tags:
+  - Underworld Code
+figures: 2
+# Generated from the model rather than a stock photograph, so there is nobody to
+# credit. `figures` counts the figures in the body; the banner is not one.
+banner: figures/banner.png
+banner_credit: null
+source: native