Symbolic FEM and Adjoint Forms (UWTN 2026-020)

CONTRIBUTING.md

@@ -55,7 +55,7 @@ what each is meant to explain, and which of them are grouped into a candidate
 paper. Notes are written as micro-preprints: two or three related ones become
 the basis of a peer-reviewed modelling paper, so a note that will feed a paper
 should carry its own evidence — figures generated by a script in the note's
-`examples/`, not screenshots of a run nobody can repeat.
+`figures/src/`, not screenshots of a run nobody can repeat.
 
 ```bash
 pixi run new --slug my-note-slug --title "My note" --author louis
@@ -187,10 +187,11 @@ the author. Second person it reports as a density rather than a limit, for the
 reason just given, and it is not part of `pixi run test`, because several of
 its checks are densities whose right value depends on what the document is.
 
-### The example code is part of the note
+### The code is part of the note
 
-A note's `examples/` scripts are read, not only run, so they are held to the same
-standard as the prose. Underworld3's own guides are the authority and should not
+A note's code — the examples in `examples/` and the figure sources in
+`figures/src/` — is read, not only run, so it is held to the same standard as
+the prose. Underworld3's own guides are the authority and should not
 be restated here — the
 [Style Charter](https://github.com/underworldcode/underworld3/blob/development/docs/developer/UW3_STYLE_CHARTER.md)
 and the
@@ -237,7 +238,15 @@ something works, not how hard it was to make work.
   general thing the reader can check, and leave the specific one to its own
   note.
 - **Every figure and table ships the script that made it**, in the note's
-  `examples/`. A number a reader cannot regenerate is an assertion.
+  `figures/src/`, beside the figures it draws. A number a reader cannot
+  regenerate is an assertion.
+- **`examples/` holds what a reader is meant to launch**: the notebooks a
+  Binder link opens, and the scripts and data they read. A script that only
+  draws a figure goes in `figures/src/`, so a reader who opens `examples/`
+  finds the example rather than the plotting code around it. Both folders are
+  deposited with the note. Notes published before this rule keep their figure
+  scripts in `examples/`; they move when the note next gains a real example,
+  in the same change.
 - **Alt text describes what the image shows, including the numbers.** It is read
   by people who never see the figure, and by anyone reading the markdown.
 - **Citations are pinned** in the note's `references.bib` and cited in one form,

acknowledgements.yml

@@ -19,3 +19,4 @@
 - finding-particles-in-a-distributed-unstructured-mesh
 - underworld3-published-in-journal-of-open-source-software
 - underworld3-come-and-get-it
+- symbolic-fem-and-adjoint-forms

articles/symbolic-fem-and-adjoint-forms/figures/src/adjoint-loop.typ

@@ -0,0 +1,51 @@
+// One pass of an inversion: the four steps as a cycle, with the quantity each
+// hands to the next. The forward solve's tangent stiffness K is what the
+// adjoint solve reuses; the observations enter only as the adjoint's load;
+// the gradient is an integral, not a solve. House style: uwfig.typ.
+#import "@preview/cetz:0.3.4"
+#import "uwfig.typ": *
+
+#set page(width: auto, height: auto, margin: 8pt)
+#set text(font: ("Noto Sans", "Helvetica", "Arial"), size: 8pt, fill: ink)
+
+#cetz.canvas(length: 1cm, {
+  import cetz.draw: *
+
+  let W = 4.0
+  let H = 1.15
+  let node(pos, title, sub, name, bg: tint-fluid) = {
+    let (x, y) = pos
+    rect((x - W / 2, y - H / 2), (x + W / 2, y + H / 2),
+         fill: bg, ..styles.interface, radius: 3pt, name: name)
+    content((x, y + 0.2), text(weight: "semibold", title))
+    content((x, y - 0.22), text(size: 7.5pt, sub))
+  }
+  // edge with a label beside it; which: "vvec" (field), "fvec" (load), "hand" (ink)
+  let hand = (stroke: (paint: ink, thickness: 0.6pt),
+              mark: (end: "stealth", fill: ink, scale: 0.6, stroke: (dash: none)))
+  let edge(a, b, lbl, which, off, dash: none) = {
+    let st = if which == "hand" { hand } else { styles.at(which) }
+    if dash != none { st = (stroke: (..st.stroke, dash: dash), mark: st.mark) }
+    line(a, b, ..st)
+    content((a, 50%, b), box(fill: white, inset: 1.5pt, lbl), anchor: off)
+  }
+
+  // ---- the four steps, clockwise from top left
+  node((0, 3.2), [1. Solve the model], [$R(u; m) = 0$ by Newton], "fwd")
+  node((8, 3.2), [2. Solve the adjoint], [$K^T mu = -(partial J \/ partial u)^T$], "adj")
+  node((8, 0), [3. Form the gradient], [$partial J \/ partial m_k + mu^T partial R \/ partial m_k$], "grad")
+  node((0, 0), [4. Update the parameters], [$m <- m - alpha nabla_m J$], "upd")
+
+  // ---- what each step hands on
+  edge("fwd.east", "adj.west", $K$, "hand", "south")
+  edge("adj.south", "grad.north", $mu$, "vvec", "west")
+  edge("grad.west", "upd.east", $dif J \/ dif m$, "hand", "north")
+  edge("upd.north", "fwd.south", $m$, "hand", "east")
+  // the pseudo-load also comes out of the forward problem
+  edge("fwd.south-east", "grad.north-west", $partial R \/ partial m_k$, "fvec", "south-west", dash: (3pt, 2.5pt))
+
+  // ---- the observations: the only external input, as the adjoint's load
+  content((8, 5.3), text(size: 8pt, [observations $d^"obs"$]), name: "obs")
+  line((8, 5.0), "adj.north", ..styles.fvec)
+  content((8.15, 4.45), text(size: 7.5pt, fill: force, $partial J \/ partial u$), anchor: "west")
+})

articles/symbolic-fem-and-adjoint-forms/figures/src/generate-misfit-surface-data.py

@@ -0,0 +1,61 @@
+"""Geometry for the misfit-surface figure (cetz draws, Python computes).
+
+J(m1, m2) is a bowl over the parameter plane, drawn as a grid of surface
+lines in an oblique projection. One point on it is where the model has been
+solved; the adjoint gives the gradient there, drawn as the tangent plane and
+the downhill direction in the parameter plane; a few gradient steps lead to
+the minimum, where the model fits the observations.
+"""
+import json, math, os
+
+A, B, C = 0.62, 0.38, 0.22          # bowl: J = A x^2 + B y^2 + C x y  (x, y from the minimum)
+MIN = (1.35, 1.1)                  # m* in the parameter plane
+X0, X1, Y0, Y1 = 0.0, 2.4, 0.0, 2.2
+NX, NY = 9, 9
+ZS = 0.36                          # vertical scale for J
+
+def J(m1, m2):
+    x, y = m1 - MIN[0], m2 - MIN[1]
+    return A * x * x + B * y * y + C * x * y
+
+def grad(m1, m2):
+    x, y = m1 - MIN[0], m2 - MIN[1]
+    return (2 * A * x + C * y, 2 * B * y + C * x)
+
+# oblique projection: m1 to the right, m2 back-and-up, J straight up
+def P(m1, m2, z):
+    return [m1 + 0.55 * m2, 0.46 * m2 + ZS * z]
+
+lines = []
+for i in range(NX + 1):
+    m1 = X0 + (X1 - X0) * i / NX
+    lines.append([P(m1, Y0 + (Y1 - Y0) * k / 40, J(m1, Y0 + (Y1 - Y0) * k / 40)) for k in range(41)])
+for j in range(NY + 1):
+    m2 = Y0 + (Y1 - Y0) * j / NY
+    lines.append([P(X0 + (X1 - X0) * k / 40, m2, J(X0 + (X1 - X0) * k / 40, m2)) for k in range(41)])
+
+m = (0.62, 0.42)                   # the current parameters
+z = J(*m); g = grad(*m)
+# tangent plane patch, +-h in each parameter
+h = 0.20
+patch = [P(m[0] + sx * h, m[1] + sy * h, z + g[0] * sx * h + g[1] * sy * h)
+         for sx, sy in ((-1, -1), (1, -1), (1, 1), (-1, 1))]
+# gradient-descent path in the parameter plane, drawn on the surface
+path = [m]; p = m
+for _ in range(6):
+    gx, gy = grad(*p)
+    p = (p[0] - 0.36 * gx, p[1] - 0.36 * gy)
+    path.append(p)
+path_surf = [P(a, b, J(a, b)) for a, b in path]
+# the descent direction at m, in the base plane (z = 0) and on the surface
+gn = math.hypot(*g); d = (-g[0] / gn * 0.42, -g[1] / gn * 0.42)
+out = dict(lines=lines, point=P(*m, z), foot=P(*m, 0), patch=patch,
+           foot_x=P(m[0], Y0, 0), foot_y=P(X0, m[1], 0),
+           arrow=[P(*m, z), P(m[0] + d[0], m[1] + d[1], z + g[0] * d[0] + g[1] * d[1])],
+           base_arrow=[P(*m, 0), P(m[0] + d[0], m[1] + d[1], 0)],
+           path=path_surf, minimum=P(*MIN, 0.0), minimum_foot=P(*MIN, 0),
+           axes=dict(o=P(X0, Y0, 0), x=P(X1 + 0.15, Y0, 0), y=P(X0, Y1 + 0.15, 0), z=P(X0, Y0, 1.6)),
+           base=[P(X0, Y0, 0), P(X1, Y0, 0), P(X1, Y1, 0), P(X0, Y1, 0)])
+here = os.path.dirname(os.path.abspath(__file__))
+json.dump(out, open(os.path.join(here, "misfit-surface-data.json"), "w"))
+print("J at m:", round(z, 3), "grad:", [round(v, 2) for v in g], "steps:", len(path))

articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface.typ

@@ -0,0 +1,43 @@
+// The misfit as a surface over the parameter plane. One forward solve puts
+// us at one point; the adjoint gives the tangent plane there, hence the
+// downhill direction; a few steps reach the minimum. House style: uwfig.typ.
+#import "@preview/cetz:0.3.4"
+#import "uwfig.typ": *
+
+#set page(width: auto, height: auto, margin: 8pt)
+#set text(font: ("Noto Sans", "Helvetica", "Arial"), size: 8pt, fill: ink)
+
+#let d = json("misfit-surface-data.json")
+#let pt(p) = (p.at(0), p.at(1))
+
+#cetz.canvas(length: 1.6cm, {
+  import cetz.draw: *
+  // the parameter plane and its axes
+  line(..d.base.map(pt), close: true, stroke: (paint: mesh-grey, thickness: 0.4pt), fill: tint-solid.transparentize(60%))
+  line(pt(d.axes.o), pt(d.axes.x), ..styles.leader, mark: (end: "stealth", fill: ink, scale: 0.5))
+  line(pt(d.axes.o), pt(d.axes.y), ..styles.leader, mark: (end: "stealth", fill: ink, scale: 0.5))
+  line(pt(d.axes.o), pt(d.axes.z), ..styles.leader, mark: (end: "stealth", fill: ink, scale: 0.5))
+  content(pt(d.axes.x), $m_1$, anchor: "north-west", padding: 2pt)
+  label(pt(d.axes.y), $m_2$, align-to: "south")
+  content(pt(d.axes.z), $J$, anchor: "south", padding: 2pt)
+  // the foot of the current point: coordinates read off the two axes
+  line(pt(d.foot_x), pt(d.foot), stroke: (paint: ink, thickness: 0.5pt, dash: (1.5pt, 2pt)))
+  line(pt(d.foot_y), pt(d.foot), stroke: (paint: ink, thickness: 0.5pt, dash: (1.5pt, 2pt)))
+  circle(pt(d.foot), radius: 0.03, fill: ink, stroke: none)
+  // the surface: a grid of lines, thin
+  for l in d.lines { line(..l.map(pt), stroke: (paint: velocity.transparentize(35%), thickness: 0.4pt)) }
+  // the minimum, and the foot of the current point
+  circle(pt(d.minimum), radius: 0.045, fill: ink, stroke: none)
+  let (qx, qy) = pt(d.minimum)
+  label((qx, qy + 0.10), $m^*$, align-to: "south")
+  // the tangent plane, and the descent path
+  line(..d.patch.map(pt), close: true, fill: tint-fluid.transparentize(25%), ..styles.interface)
+  line(..d.path.map(pt), stroke: (paint: ink, thickness: 0.6pt, dash: (2pt, 2pt)))
+  for p in d.path.slice(1) { circle(pt(p), radius: 0.03, fill: white, stroke: ink + 0.5pt) }
+  line(pt(d.arrow.at(0)), pt(d.arrow.at(1)), ..styles.fvec)
+  circle(pt(d.point), radius: 0.055, fill: force, stroke: white + 0.6pt)
+  let (px, py) = pt(d.point)
+  label((px, py + 0.26), $J(m)$, align-to: "south")
+  let (ax, ay) = pt(d.arrow.at(1))
+  label((ax + 0.06, ay - 0.06), text(fill: force, $-nabla_m J$), align-to: "north-west")
+})

articles/symbolic-fem-and-adjoint-forms/figures/src/uwfig.typ

@@ -0,0 +1,70 @@
+#import "@preview/cetz:0.3.4"
+// House style for underworld3 line diagrams (Typst + CeTZ).
+// The same roles, weights, dash patterns and colours as uwfig.py, so a figure
+// can be prototyped in matplotlib and redrawn here (or the reverse) without
+// changing its look. Import from a figure with
+//     #import "uwfig.typ": *
+// and spread a role into a draw call:  line(a, b, ..styles.boundary)
+//
+// Versions: written against cetz 0.3.4 / Typst 0.13 (the versions the
+// cetz-figures skill compiles with); check the marks syntax on a newer cetz.
+
+// ---- colour tokens (secondary channel: every role also has a weight and dash)
+#let ink = rgb("#1f2a33")        // walls, bodies, text, dimensions
+#let mesh-grey = rgb("#9aa3aa")  // mesh edges
+#let velocity = rgb("#1b6f6a")   // motion: velocity vectors, inflow profiles, streamlines
+#let force = rgb("#a1522a")      // forces, pressure, tractions, applied stress
+#let energy = rgb("#6b3d8f")     // heat, work, energy flux
+#let tint-fluid = rgb("#eef3f6")
+#let tint-solid = rgb("#e4e0da")
+
+// ---- line roles: three weights (1.2 / 0.8 / 0.5 pt at final column width)
+#let styles = (
+  boundary:   (stroke: (paint: ink, thickness: 1.2pt)),                          // walls, bodies
+  interface:  (stroke: (paint: ink, thickness: 0.8pt)),                          // free surface
+  streamline: (stroke: (paint: velocity, thickness: 0.5pt)),                     // flow lines
+  hidden:     (stroke: (paint: ink, thickness: 0.8pt, dash: "dashed")),
+  cvline:     (stroke: (paint: ink, thickness: 0.8pt, dash: (5pt, 2pt, 1pt, 2pt))), // control volume
+  axis:       (stroke: (paint: ink, thickness: 0.5pt, dash: (8pt, 2pt, 1pt, 2pt))), // centre line
+  mesh:       (stroke: (paint: mesh-grey, thickness: 0.25pt)),
+  leader:     (stroke: (paint: ink, thickness: 0.5pt)),
+  // ---- arrow classes: each physical quantity has its own head and weight
+  vvec:       (stroke: (paint: velocity, thickness: 0.9pt, cap: "round"),
+               mark: (end: "stealth", fill: velocity, scale: 0.7, stroke: (dash: none))),
+  fvec:       (stroke: (paint: force, thickness: 1.3pt, cap: "round"),
+               mark: (end: "stealth", fill: force, scale: 0.9, stroke: (dash: none))),
+  evec:       (stroke: (paint: energy, thickness: 1.1pt, cap: "round"),
+               mark: (end: "stealth", fill: energy, scale: 0.8, stroke: (dash: none))),
+  dim:        (stroke: (paint: ink, thickness: 0.5pt),
+               mark: (symbol: ">", fill: ink, scale: 0.5)),                         // both ends
+)
+
+// ---- 45 degree hatching for solids / cross-sections (needs modpattern)
+// #import "@preview/modpattern:0.1.0": modpattern
+// #let hatch = modpattern((3pt, 3pt),
+//   std.line(start: (0%, 0%), end: (100%, 100%), stroke: (paint: ink, thickness: 0.4pt)))
+
+// ---- a parabolic inflow profile outside a vertical edge, arrows following the flow
+// side: "left" = inlet (profile to the left of x0, arrows into the domain)
+//       "right" = outlet (profile to the right of x0, arrows out of the domain)
+#let inflow-profile(x0, y0, y1, umax, side: "left", n: 9, scale: 1.0) = {
+  import cetz.draw: *
+  let h = y1 - y0
+  let u(y) = 4 * umax * (y - y0) * (y1 - y) / (h * h) * scale
+  let sgn = if side == "left" { -1 } else { 1 }
+  let pts = range(0, 41).map(i => { let y = y0 + h * i / 40; (x0 + sgn * u(y), y) })
+  line(..pts, stroke: (paint: velocity, thickness: 0.9pt))
+  for i in range(1, n + 1) {
+    let y = y0 + h * i / (n + 1)
+    if side == "left" { line((x0 - u(y), y), (x0, y), ..styles.vvec) }
+    else { line((x0, y), (x0 + u(y), y), ..styles.vvec) }
+  }
+}
+
+// ---- a direct label with a white backing so it reads over mesh lines
+// `anchor:` would be shadowed by `import cetz.draw: *` inside (cetz-figures gotcha 1)
+#let label(pos, body, align-to: "center") = {
+  import cetz.draw: *
+  content(pos, box(fill: white.transparentize(15%), inset: 2pt, radius: 2pt, text(8pt, fill: ink, body)),
+          anchor: align-to)
+}

articles/symbolic-fem-and-adjoint-forms/metadata.yml

@@ -0,0 +1,31 @@
+# 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-020
+slug: symbolic-fem-and-adjoint-forms
+title: Symbolic FEM and Adjoint Forms
+article_type: technical-note
+status: published
+authors:
+  - name: Louis Moresi
+    orcid: 0000-0003-3685-174X
+    affiliation: Australian National University
+publication_date: 2026-10-25
+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: /symbolic-fem-and-adjoint-forms/
+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:
+ghost_tags:
+  - Underworld Code
+figures: 2
+source: native
+banner_credit: "Photo by <a href=\"https://unsplash.com/@loganvoss?utm_source=underworld-technical-notes&utm_medium=referral&utm_campaign=api-credit\">Logan Voss</a> / <a href=\"https://unsplash.com/?utm_source=underworld-technical-notes&utm_medium=referral&utm_campaign=api-credit\">Unsplash</a>"

articles/symbolic-fem-and-adjoint-forms/references.bib

@@ -0,0 +1,90 @@
+% Resolved from doi.org and pinned here on purpose. Letting MyST fetch these
+% at build time makes the build depend on doi.org answering, which fails on CI
+% for one publisher or another and produces no PDF. A deposited PDF cannot be
+% repaired after the fact, so the bibliography is part of the source.
+
+@article{Moresi_2025,
+  title   = {Underworld3: Mathematically Self-Describing Modelling in Python
+             for Desktop, HPC and Cloud},
+  volume  = {10},
+  ISSN    = {2475-9066},
+  url     = {http://dx.doi.org/10.21105/joss.07831},
+  DOI     = {10.21105/joss.07831},
+  number  = {112},
+  journal = {Journal of Open Source Software},
+  publisher = {The Open Journal},
+  author  = {Moresi, Louis and Mansour, John and Giordani, Julian and
+             Knepley, Matt and Knight, Ben and Graciosa, Juan Carlos and
+             Gollapalli, Thyagarajulu and Lu, Neng and Beucher, Romain},
+  year    = {2025},
+  month   = aug,
+  pages   = {7831}
+}
+
+@misc{Moresi_2026_sympy,
+  title   = {How Underworld3 Turns SymPy into C},
+  url     = {https://doi.org/10.6084/m9.figshare.33193572},
+  DOI     = {10.6084/m9.figshare.33193572},
+  author  = {Moresi, Louis and Mansour, John and Giordani, Julian},
+  series  = {Underworld Technical Notes},
+  number  = {UWTN 2026-008},
+  publisher = {figshare},
+  year    = {2026}
+}
+
+@inproceedings{Knepley_2013,
+  title     = {Achieving High Performance with Unified Residual Evaluation},
+  url       = {https://arxiv.org/abs/1309.1204},
+  DOI       = {10.48550/arXiv.1309.1204},
+  author    = {Knepley, Matthew G. and Brown, Jed and Rupp, Karl and Smith, Barry F.},
+  year      = {2013},
+  publisher = {arXiv}
+}
+
+@article{Plessix_2006,
+  title   = {A review of the adjoint-state method for computing the gradient
+             of a functional with geophysical applications},
+  volume  = {167},
+  ISSN    = {0956-540X},
+  url     = {http://dx.doi.org/10.1111/j.1365-246X.2006.02978.x},
+  DOI     = {10.1111/j.1365-246X.2006.02978.x},
+  number  = {2},
+  journal = {Geophysical Journal International},
+  publisher = {Oxford University Press (OUP)},
+  author  = {Plessix, R.-E.},
+  year    = {2006},
+  pages   = {495--503}
+}
+
+@article{Bunge_2003,
+  title   = {Mantle circulation models with variational data assimilation:
+             inferring past mantle flow and structure from plate motion
+             histories and seismic tomography},
+  volume  = {152},
+  ISSN    = {1365-246X},
+  url     = {http://dx.doi.org/10.1046/j.1365-246X.2003.01823.x},
+  DOI     = {10.1046/j.1365-246x.2003.01823.x},
+  number  = {2},
+  journal = {Geophysical Journal International},
+  publisher = {Oxford University Press (OUP)},
+  author  = {Bunge, Hans-Peter and Hagelberg, C. R. and Travis, B. J.},
+  year    = {2003},
+  month   = feb,
+  pages   = {280--301}
+}
+
+@article{Liu_2008,
+  title   = {Simultaneous inversion of mantle properties and initial
+             conditions using an adjoint of mantle convection},
+  volume  = {113},
+  ISSN    = {0148-0227},
+  url     = {http://dx.doi.org/10.1029/2008JB005594},
+  DOI     = {10.1029/2008jb005594},
+  number  = {B8},
+  journal = {Journal of Geophysical Research: Solid Earth},
+  publisher = {American Geophysical Union (AGU)},
+  author  = {Liu, Lijun and Gurnis, Michael},
+  year    = {2008},
+  month   = aug,
+  pages   = {B08405}
+}