diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index b3c0064..01e8040 100644 --- a/CONTRIBUTING.md +++ b/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, diff --git a/acknowledgements.yml b/acknowledgements.yml index 7cfba8c..46f300f 100644 --- a/acknowledgements.yml +++ b/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 diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.png b/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.png new file mode 100644 index 0000000..01850d5 Binary files /dev/null and b/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.png differ diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.svg b/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.svg new file mode 100644 index 0000000..655c5a5 --- /dev/null +++ b/articles/symbolic-fem-and-adjoint-forms/figures/adjoint-loop.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/banner.jpg b/articles/symbolic-fem-and-adjoint-forms/figures/banner.jpg new file mode 100644 index 0000000..f57d4f8 Binary files /dev/null and b/articles/symbolic-fem-and-adjoint-forms/figures/banner.jpg differ diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.png b/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.png new file mode 100644 index 0000000..34d1a75 Binary files /dev/null and b/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.png differ diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.svg b/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.svg new file mode 100644 index 0000000..79323d4 --- /dev/null +++ b/articles/symbolic-fem-and-adjoint-forms/figures/misfit-surface.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/src/adjoint-loop.typ b/articles/symbolic-fem-and-adjoint-forms/figures/src/adjoint-loop.typ new file mode 100644 index 0000000..d309cde --- /dev/null +++ b/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") +}) diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/src/generate-misfit-surface-data.py b/articles/symbolic-fem-and-adjoint-forms/figures/src/generate-misfit-surface-data.py new file mode 100644 index 0000000..77eb7b4 --- /dev/null +++ b/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)) diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface-data.json b/articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface-data.json new file mode 100644 index 0000000..d8a9cd4 --- /dev/null +++ b/articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface-data.json @@ -0,0 +1 @@ +{"lines": [[[0.0, 0.6899220000000001], [0.030250000000000006, 0.6932024200000001], [0.06050000000000001, 0.6973104800000001], [0.09075000000000001, 0.70224618], [0.12100000000000002, 0.7080095200000001], [0.15125000000000002, 0.7146005000000002], [0.18150000000000002, 0.7220191200000001], [0.21175000000000005, 0.7302653800000001], [0.24200000000000005, 0.7393392800000002], [0.27225, 0.7492408200000001], [0.30250000000000005, 0.7599700000000003], [0.3327500000000001, 0.7715268200000001], [0.36300000000000004, 0.7839112800000001], [0.3932500000000001, 0.79712338], [0.4235000000000001, 0.8111631200000001], [0.45375, 0.8260305000000001], [0.4840000000000001, 0.84172552], [0.5142500000000001, 0.8582481800000001], [0.5445, 0.8755984800000001], [0.5747500000000001, 0.8937764200000001], [0.6050000000000001, 0.9127820000000002], [0.6352500000000001, 0.93261522], [0.6655000000000002, 0.9532760800000002], [0.6957500000000001, 0.9747645800000001], [0.7260000000000001, 0.9970807200000001], [0.7562500000000002, 1.0202245000000003], [0.7865000000000002, 1.04419592], [0.8167500000000001, 1.0689949800000003], [0.8470000000000002, 1.0946216800000002], [0.8772500000000002, 1.12107602], [0.9075, 1.148358], [0.9377500000000001, 1.1764676200000002], [0.9680000000000002, 1.2054048800000001], [0.9982500000000002, 1.23516978], [1.0285000000000002, 1.2657623200000003], [1.05875, 1.2971825000000001], [1.089, 1.3294303200000002], [1.11925, 1.3625057800000002], [1.1495000000000002, 1.39640888], [1.1797500000000003, 1.4311396200000004], [1.2100000000000002, 1.4666980000000005]], [[0.26666666666666666, 0.521858], [0.29691666666666666, 0.5263000200000001], [0.32716666666666666, 0.5315696800000002], [0.35741666666666666, 0.53766698], [0.3876666666666667, 0.5445919200000001], [0.4179166666666667, 0.5523445000000001], [0.4481666666666667, 0.5609247200000002], [0.4784166666666667, 0.5703325800000001], [0.5086666666666667, 0.5805680800000002], [0.5389166666666667, 0.59163122], [0.5691666666666667, 0.6035220000000001], [0.5994166666666667, 0.6162404200000001], [0.6296666666666667, 0.6297864800000001], [0.6599166666666667, 0.6441601800000001], [0.6901666666666668, 0.6593615200000001], [0.7204166666666667, 0.6753905], [0.7506666666666668, 0.6922471200000001], [0.7809166666666667, 0.7099313800000001], [0.8111666666666666, 0.72844328], [0.8414166666666667, 0.7477828200000001], [0.8716666666666668, 0.7679500000000001], [0.9019166666666667, 0.7889448200000001], [0.9321666666666668, 0.8107672800000001], [0.9624166666666667, 0.8334173800000001], [0.9926666666666668, 0.8568951200000001], [1.022916666666667, 0.8812005000000002], [1.0531666666666668, 0.9063335200000001], [1.0834166666666667, 0.9322941800000001], [1.1136666666666668, 0.9590824800000002], [1.143916666666667, 0.9866984200000002], [1.1741666666666666, 1.015142], [1.2044166666666667, 1.04441322], [1.2346666666666668, 1.07451208], [1.264916666666667, 1.1054385800000002], [1.2951666666666668, 1.1371927200000003], [1.3254166666666667, 1.1697745000000002], [1.3556666666666666, 1.20318392], [1.3859166666666667, 1.2374209800000002], [1.4161666666666668, 1.2724856800000002], [1.446416666666667, 1.3083780200000004], [1.4766666666666668, 1.3450980000000003]], [[0.5333333333333333, 0.38553800000000005], [0.5635833333333333, 0.3911416200000001], [0.5938333333333333, 0.39757288000000013], [0.6240833333333333, 0.4048317800000001], [0.6543333333333333, 0.4129183200000001], [0.6845833333333333, 0.42183250000000005], [0.7148333333333333, 0.43157432000000007], [0.7450833333333333, 0.4421437800000001], [0.7753333333333334, 0.4535408800000001], [0.8055833333333333, 0.4657656200000001], [0.8358333333333334, 0.4788180000000001], [0.8660833333333334, 0.49269802000000007], [0.8963333333333334, 0.5074056800000001], [0.9265833333333334, 0.52294098], [0.9568333333333334, 0.53930392], [0.9870833333333333, 0.5564945000000001], [1.0173333333333334, 0.5745127200000001], [1.0475833333333333, 0.5933585800000001], [1.0778333333333334, 0.6130320800000001], [1.1080833333333335, 0.6335332200000001], [1.1383333333333334, 0.6548620000000002], [1.1685833333333333, 0.67701842], [1.1988333333333334, 0.7000024800000001], [1.2290833333333335, 0.7238141800000001], [1.2593333333333334, 0.7484535200000001], [1.2895833333333335, 0.7739205000000002], [1.3198333333333334, 0.8002151200000001], [1.3500833333333335, 0.8273373800000001], [1.3803333333333336, 0.8552872800000002], [1.4105833333333335, 0.8840648200000002], [1.4408333333333334, 0.91367], [1.4710833333333335, 0.9441028200000001], [1.5013333333333336, 0.9753632800000002], [1.5315833333333335, 1.0074513800000002], [1.5618333333333334, 1.0403671200000002], [1.5920833333333335, 1.0741105000000002], [1.6223333333333332, 1.1086815200000002], [1.6525833333333333, 1.1440801800000002], [1.6828333333333334, 1.1803064800000003], [1.7130833333333335, 1.2173604200000003], [1.7433333333333336, 1.2552420000000004]], [[0.7999999999999999, 0.2809620000000001], [0.8302499999999999, 0.28772722000000006], [0.8604999999999999, 0.2953200800000001], [0.8907499999999999, 0.3037405800000001], [0.9209999999999999, 0.3129887200000001], [0.9512499999999999, 0.3230645000000001], [0.9814999999999999, 0.3339679200000001], [1.01175, 0.3456989800000001], [1.042, 0.3582576800000001], [1.07225, 0.3716440200000001], [1.1025, 0.3858580000000001], [1.1327500000000001, 0.4008996200000001], [1.163, 0.41676888000000006], [1.19325, 0.4334657800000001], [1.2235, 0.4509903200000001], [1.25375, 0.46934250000000005], [1.284, 0.48852232000000007], [1.31425, 0.5085297800000002], [1.3445, 0.52936488], [1.3747500000000001, 0.5510276200000002], [1.405, 0.5735180000000002], [1.43525, 0.59683602], [1.4655, 0.6209816800000001], [1.4957500000000001, 0.6459549800000002], [1.526, 0.6717559200000001], [1.5562500000000001, 0.6983845000000002], [1.5865, 0.7258407200000001], [1.6167500000000001, 0.75412458], [1.6470000000000002, 0.7832360800000002], [1.6772500000000001, 0.8131752200000002], [1.7075, 0.843942], [1.7377500000000001, 0.8755364200000001], [1.7680000000000002, 0.9079584800000001], [1.7982500000000001, 0.9412081800000002], [1.8285, 0.9752855200000002], [1.8587500000000001, 1.0101905], [1.8889999999999998, 1.04592312], [1.91925, 1.0824833800000002], [1.9495, 1.1198712800000001], [1.9797500000000001, 1.1580868200000003], [2.0100000000000002, 1.1971300000000002]], [[1.0666666666666667, 0.20813000000000004], [1.0969166666666668, 0.21605682000000004], [1.1271666666666667, 0.22481128000000006], [1.1574166666666668, 0.23439338000000004], [1.1876666666666666, 0.24480312000000007], [1.2179166666666668, 0.25604050000000006], [1.2481666666666666, 0.26810552000000004], [1.2784166666666668, 0.2809981800000001], [1.3086666666666666, 0.29471848000000006], [1.3389166666666665, 0.30926642000000004], [1.3691666666666666, 0.3246420000000001], [1.3994166666666668, 0.34084522000000006], [1.4296666666666666, 0.35787608000000004], [1.4599166666666668, 0.37573458000000004], [1.4901666666666666, 0.39442072000000006], [1.5204166666666667, 0.41393450000000004], [1.5506666666666669, 0.43427592000000004], [1.5809166666666667, 0.4554449800000001], [1.6111666666666666, 0.47744168000000003], [1.6414166666666667, 0.5002660200000001], [1.6716666666666669, 0.5239180000000001], [1.7019166666666667, 0.54839762], [1.7321666666666669, 0.5737048800000001], [1.7624166666666667, 0.59983978], [1.7926666666666669, 0.6268023200000001], [1.822916666666667, 0.6545925000000001], [1.8531666666666669, 0.68321032], [1.8834166666666667, 0.7126557800000001], [1.9136666666666668, 0.7429288800000001], [1.943916666666667, 0.7740296200000001], [1.9741666666666666, 0.805958], [2.0044166666666667, 0.8387140200000001], [2.034666666666667, 0.8722976800000002], [2.064916666666667, 0.9067089800000001], [2.0951666666666666, 0.9419479200000003], [2.1254166666666667, 0.9780145], [2.155666666666667, 1.01490872], [2.1859166666666665, 1.0526305800000002], [2.216166666666667, 1.0911800800000002], [2.2464166666666667, 1.1305572200000003], [2.276666666666667, 1.1707620000000003]], [[1.3333333333333333, 0.16704200000000005], [1.3635833333333334, 0.17613042000000007], [1.3938333333333333, 0.18604648000000004], [1.4240833333333334, 0.19679018000000004], [1.4543333333333333, 0.20836152000000008], [1.4845833333333334, 0.22076050000000008], [1.5148333333333333, 0.23398712000000005], [1.5450833333333334, 0.24804138000000006], [1.5753333333333333, 0.26292328000000004], [1.6055833333333331, 0.27863282000000006], [1.6358333333333333, 0.29517000000000004], [1.6660833333333334, 0.31253482000000005], [1.6963333333333332, 0.33072728000000007], [1.7265833333333334, 0.34974738000000005], [1.7568333333333332, 0.36959512000000005], [1.7870833333333334, 0.3902705], [1.8173333333333335, 0.41177352000000006], [1.8475833333333334, 0.4341041800000001], [1.8778333333333332, 0.45726248], [1.9080833333333334, 0.48124842000000007], [1.9383333333333335, 0.5060620000000001], [1.9685833333333334, 0.53170322], [1.9988333333333335, 0.5581720800000001], [2.0290833333333333, 0.5854685800000001], [2.0593333333333335, 0.6135927200000001], [2.0895833333333336, 0.6425445000000002], [2.1198333333333332, 0.6723239200000001], [2.1500833333333333, 0.7029309800000001], [2.1803333333333335, 0.7343656800000001], [2.2105833333333336, 0.7666280200000002], [2.2408333333333332, 0.7997179999999999], [2.2710833333333333, 0.8336356200000001], [2.3013333333333335, 0.8683808800000001], [2.3315833333333336, 0.9039537800000002], [2.3618333333333332, 0.9403543200000002], [2.3920833333333333, 0.9775825], [2.4223333333333334, 1.01563832], [2.452583333333333, 1.0545217800000002], [2.4828333333333337, 1.0942328800000003], [2.5130833333333333, 1.1347716200000004], [2.5433333333333334, 1.1761380000000003]], [[1.5999999999999999, 0.15769800000000003], [1.63025, 0.16794802000000003], [1.6604999999999999, 0.17902568000000002], [1.69075, 0.19093098000000003], [1.7209999999999999, 0.20366392000000005], [1.75125, 0.21722450000000004], [1.7814999999999999, 0.23161272], [1.81175, 0.24682858000000002], [1.8419999999999999, 0.26287208], [1.8722499999999997, 0.27974322], [1.9024999999999999, 0.29744200000000004], [1.93275, 0.31596842], [1.9629999999999999, 0.33532248000000003], [1.99325, 0.35550418], [2.0235, 0.37651352000000005], [2.05375, 0.3983505], [2.084, 0.42101512], [2.11425, 0.44450738000000006], [2.1445, 0.46882728], [2.17475, 0.49397482000000004], [2.205, 0.5199500000000001], [2.2352499999999997, 0.5467528199999999], [2.2655000000000003, 0.57438328], [2.29575, 0.6028413800000001], [2.326, 0.63212712], [2.35625, 0.6622405000000001], [2.3865, 0.69318152], [2.41675, 0.72495018], [2.447, 0.7575464800000001], [2.47725, 0.7909704200000001], [2.5075, 0.8252219999999999], [2.53775, 0.8603012200000001], [2.568, 0.8962080800000001], [2.59825, 0.9329425800000001], [2.6285, 0.9705047200000002], [2.65875, 1.0088945], [2.689, 1.04811192], [2.7192499999999997, 1.0881569800000002], [2.7495000000000003, 1.1290296800000001], [2.77975, 1.1707300200000004], [2.81, 1.2132580000000002]], [[1.8666666666666667, 0.180098], [1.8969166666666668, 0.19150962], [1.9271666666666667, 0.20374888000000002], [1.9574166666666668, 0.21681578], [1.9876666666666667, 0.23071032000000002], [2.017916666666667, 0.24543250000000005], [2.048166666666667, 0.26098232], [2.0784166666666666, 0.27735978000000006], [2.1086666666666667, 0.29456488000000003], [2.138916666666667, 0.31259762], [2.169166666666667, 0.33145800000000003], [2.199416666666667, 0.35114602000000006], [2.2296666666666667, 0.37166168000000005], [2.259916666666667, 0.39300498], [2.290166666666667, 0.41517592000000003], [2.3204166666666666, 0.4381745], [2.3506666666666667, 0.46200072000000003], [2.380916666666667, 0.48665458000000006], [2.4111666666666665, 0.51213608], [2.441416666666667, 0.5384452200000001], [2.4716666666666667, 0.5655820000000001], [2.501916666666667, 0.59354642], [2.532166666666667, 0.62233848], [2.5624166666666666, 0.65195818], [2.5926666666666667, 0.6824055200000001], [2.622916666666667, 0.7136805000000002], [2.653166666666667, 0.74578312], [2.683416666666667, 0.77871338], [2.7136666666666667, 0.8124712800000001], [2.743916666666667, 0.8470568200000002], [2.774166666666667, 0.88247], [2.8044166666666666, 0.91871082], [2.834666666666667, 0.9557792800000001], [2.864916666666667, 0.9936753800000001], [2.895166666666667, 1.0323991200000002], [2.925416666666667, 1.0719505], [2.9556666666666667, 1.11232952], [2.9859166666666668, 1.1535361800000001], [3.016166666666667, 1.1955704800000002], [3.046416666666667, 1.2384324200000003], [3.076666666666667, 1.2821220000000002]], [[2.1333333333333333, 0.23424199999999998], [2.1635833333333334, 0.24681522], [2.1938333333333335, 0.26021607999999996], [2.224083333333333, 0.27444458], [2.2543333333333333, 0.28950072000000004], [2.2845833333333334, 0.3053845], [2.3148333333333335, 0.32209592], [2.345083333333333, 0.33963498000000003], [2.3753333333333333, 0.35800168], [2.4055833333333334, 0.37719601999999997], [2.4358333333333335, 0.397218], [2.4660833333333336, 0.41806762], [2.4963333333333333, 0.43974488], [2.5265833333333334, 0.46224978], [2.5568333333333335, 0.48558232], [2.587083333333333, 0.5097425], [2.6173333333333333, 0.53473032], [2.6475833333333334, 0.56054578], [2.677833333333333, 0.58718888], [2.7080833333333336, 0.6146596200000001], [2.7383333333333333, 0.642958], [2.7685833333333334, 0.6720840199999999], [2.7988333333333335, 0.7020376800000001], [2.829083333333333, 0.73281898], [2.8593333333333333, 0.76442792], [2.8895833333333334, 0.7968645000000001], [2.9198333333333335, 0.83012872], [2.9500833333333336, 0.86422058], [2.9803333333333333, 0.8991400800000001], [3.0105833333333334, 0.9348872200000001], [3.0408333333333335, 0.9714619999999999], [3.071083333333333, 1.00886442], [3.1013333333333337, 1.0470944800000002], [3.1315833333333334, 1.08615218], [3.1618333333333335, 1.1260375200000001], [3.1920833333333336, 1.1667505], [3.2223333333333333, 1.20829112], [3.2525833333333334, 1.25065938], [3.2828333333333335, 1.2938552800000003], [3.3130833333333336, 1.3378788200000002], [3.3433333333333337, 1.3827300000000002]], [[2.4, 0.32012999999999997], [2.43025, 0.3338648199999999], [2.4605, 0.3484272799999999], [2.49075, 0.36381737999999997], [2.521, 0.38003511999999995], [2.55125, 0.39708049999999995], [2.5815, 0.41495351999999996], [2.61175, 0.43365418], [2.642, 0.45318247999999994], [2.67225, 0.47353841999999996], [2.7025, 0.494722], [2.7327500000000002, 0.5167332199999999], [2.763, 0.5395720799999999], [2.79325, 0.5632385799999999], [2.8235, 0.58773272], [2.85375, 0.6130544999999998], [2.884, 0.6392039199999999], [2.91425, 0.66618098], [2.9444999999999997, 0.6939856799999999], [2.9747500000000002, 0.72261802], [3.005, 0.752078], [3.03525, 0.7823656199999999], [3.0655, 0.81348088], [3.09575, 0.84542378], [3.126, 0.87819432], [3.15625, 0.9117925], [3.1865, 0.9462183200000001], [3.21675, 0.9814717799999999], [3.247, 1.0175528800000002], [3.27725, 1.05446162], [3.3075, 1.0921979999999998], [3.3377499999999998, 1.13076202], [3.3680000000000003, 1.17015368], [3.39825, 1.21037298], [3.4285, 1.2514199200000002], [3.45875, 1.2932945], [3.489, 1.3359967199999998], [3.51925, 1.37952658], [3.5495, 1.42388408], [3.57975, 1.4690692200000002], [3.6100000000000003, 1.515082]], [[0.0, 0.6899220000000001], [0.06, 0.6493399200000002], [0.12, 0.61036488], [0.18, 0.57299688], [0.24, 0.5372359200000001], [0.3, 0.5030820000000001], [0.36, 0.4705351200000001], [0.42000000000000004, 0.43959528], [0.48, 0.4102624800000001], [0.5399999999999999, 0.3825367200000001], [0.6, 0.3564180000000001], [0.6599999999999999, 0.33190632000000014], [0.72, 0.3090016800000001], [0.78, 0.2877040800000001], [0.8400000000000001, 0.26801352], [0.9, 0.24993000000000004], [0.96, 0.23345352000000003], [1.02, 0.21858408000000007], [1.0799999999999998, 0.2053216800000001], [1.1400000000000001, 0.19366632], [1.2, 0.18361800000000003], [1.26, 0.17517672000000004], [1.3199999999999998, 0.16834248000000004], [1.38, 0.16311528000000006], [1.44, 0.15949512000000002], [1.5, 0.157482], [1.56, 0.15707592000000004], [1.6199999999999999, 0.15827688000000004], [1.6800000000000002, 0.16108488000000004], [1.7399999999999998, 0.16549992000000002], [1.8, 0.171522], [1.8599999999999999, 0.17915112], [1.92, 0.18838728000000002], [1.98, 0.19923048], [2.04, 0.21168072], [2.1, 0.22573800000000002], [2.1599999999999997, 0.24140231999999995], [2.2199999999999998, 0.2586736799999999], [2.2800000000000002, 0.2775520800000001], [2.34, 0.29803752], [2.4, 0.32012999999999997]], [[0.13444444444444448, 0.7108366666666667], [0.19444444444444448, 0.6714161866666667], [0.2544444444444445, 0.6336027466666667], [0.31444444444444447, 0.5973963466666667], [0.37444444444444447, 0.5627969866666667], [0.43444444444444447, 0.5298046666666667], [0.49444444444444446, 0.4984193866666667], [0.5544444444444445, 0.4686411466666667], [0.6144444444444445, 0.4404699466666668], [0.6744444444444444, 0.4139057866666667], [0.7344444444444445, 0.3889486666666668], [0.7944444444444444, 0.3655985866666668], [0.8544444444444445, 0.34385554666666673], [0.9144444444444445, 0.3237195466666667], [0.9744444444444446, 0.3051905866666667], [1.0344444444444445, 0.28826866666666673], [1.0944444444444446, 0.2729537866666667], [1.1544444444444446, 0.25924594666666667], [1.2144444444444442, 0.2471451466666667], [1.2744444444444447, 0.23665138666666669], [1.3344444444444443, 0.22776466666666667], [1.3944444444444444, 0.2204849866666667], [1.4544444444444444, 0.2148123466666667], [1.5144444444444445, 0.2107467466666667], [1.5744444444444445, 0.20828818666666668], [1.6344444444444446, 0.20743666666666666], [1.6944444444444446, 0.20819218666666667], [1.7544444444444443, 0.21055474666666668], [1.8144444444444447, 0.2145243466666667], [1.8744444444444444, 0.22010098666666664], [1.9344444444444444, 0.22728466666666666], [1.9944444444444445, 0.23607538666666666], [2.0544444444444445, 0.24647314666666664], [2.1144444444444446, 0.2584779466666667], [2.1744444444444446, 0.27208978666666667], [2.2344444444444447, 0.28730866666666666], [2.2944444444444443, 0.3041345866666666], [2.3544444444444443, 0.32256754666666654], [2.414444444444445, 0.3426075466666667], [2.4744444444444444, 0.3642545866666666], [2.5344444444444445, 0.3875086666666666]], [[0.26888888888888896, 0.7480997777777779], [0.32888888888888895, 0.7098408977777779], [0.38888888888888895, 0.6731890577777777], [0.44888888888888895, 0.6381442577777778], [0.508888888888889, 0.6047064977777779], [0.568888888888889, 0.5728757777777778], [0.6288888888888889, 0.5426520977777778], [0.688888888888889, 0.5140354577777778], [0.7488888888888889, 0.48702585777777785], [0.8088888888888889, 0.4616232977777779], [0.8688888888888889, 0.4378277777777778], [0.9288888888888889, 0.41563929777777786], [0.9888888888888889, 0.39505785777777785], [1.048888888888889, 0.37608345777777785], [1.108888888888889, 0.3587160977777778], [1.168888888888889, 0.3429557777777778], [1.228888888888889, 0.3288024977777778], [1.288888888888889, 0.3162562577777778], [1.3488888888888888, 0.30531705777777784], [1.408888888888889, 0.2959848977777778], [1.468888888888889, 0.28825977777777784], [1.528888888888889, 0.2821416977777778], [1.5888888888888888, 0.27763065777777785], [1.6488888888888888, 0.2747266577777778], [1.708888888888889, 0.2734296977777778], [1.768888888888889, 0.2737397777777778], [1.828888888888889, 0.2756568977777778], [1.8888888888888888, 0.2791810577777778], [1.9488888888888891, 0.2843122577777778], [2.0088888888888885, 0.2910504977777778], [2.068888888888889, 0.2993957777777778], [2.1288888888888886, 0.3093480977777778], [2.188888888888889, 0.32090745777777774], [2.2488888888888887, 0.33407385777777776], [2.308888888888889, 0.3488472977777778], [2.368888888888889, 0.3652277777777778], [2.4288888888888884, 0.3832152977777777], [2.488888888888889, 0.40280985777777767], [2.5488888888888894, 0.4240114577777779], [2.608888888888889, 0.44682009777777776], [2.6688888888888886, 0.4712357777777777]], [[0.4033333333333334, 0.8017113333333334], [0.4633333333333334, 0.7646140533333334], [0.5233333333333334, 0.7291238133333333], [0.5833333333333334, 0.6952406133333334], [0.6433333333333333, 0.6629644533333334], [0.7033333333333334, 0.6322953333333334], [0.7633333333333334, 0.6032332533333333], [0.8233333333333335, 0.5757782133333333], [0.8833333333333333, 0.5499302133333335], [0.9433333333333334, 0.5256892533333335], [1.0033333333333334, 0.5030553333333334], [1.0633333333333332, 0.48202845333333344], [1.1233333333333333, 0.4626086133333334], [1.1833333333333333, 0.4447958133333334], [1.2433333333333334, 0.42859005333333333], [1.3033333333333335, 0.4139913333333334], [1.3633333333333333, 0.40099965333333343], [1.4233333333333333, 0.38961501333333337], [1.4833333333333332, 0.37983741333333343], [1.5433333333333334, 0.3716668533333334], [1.6033333333333333, 0.3651033333333334], [1.6633333333333333, 0.3601468533333334], [1.7233333333333332, 0.35679741333333337], [1.7833333333333332, 0.3550550133333334], [1.8433333333333333, 0.35491965333333336], [1.9033333333333333, 0.3563913333333334], [1.9633333333333334, 0.3594700533333334], [2.0233333333333334, 0.36415581333333336], [2.0833333333333335, 0.3704486133333334], [2.143333333333333, 0.37834845333333333], [2.2033333333333336, 0.38785533333333333], [2.263333333333333, 0.39896925333333333], [2.3233333333333333, 0.41169021333333333], [2.3833333333333333, 0.42601821333333334], [2.4433333333333334, 0.44195325333333335], [2.5033333333333334, 0.45949533333333337], [2.563333333333333, 0.4786444533333333], [2.623333333333333, 0.49940061333333324], [2.6833333333333336, 0.5217638133333334], [2.743333333333333, 0.5457340533333332], [2.8033333333333332, 0.5713113333333333]], [[0.5377777777777779, 0.8716713333333335], [0.597777777777778, 0.8357356533333336], [0.6577777777777779, 0.8014070133333333], [0.7177777777777778, 0.7686854133333334], [0.7777777777777779, 0.7375708533333334], [0.837777777777778, 0.7080633333333335], [0.8977777777777779, 0.6801628533333334], [0.957777777777778, 0.6538694133333334], [1.017777777777778, 0.6291830133333334], [1.077777777777778, 0.6061036533333335], [1.137777777777778, 0.5846313333333334], [1.1977777777777778, 0.5647660533333334], [1.2577777777777779, 0.5465078133333334], [1.317777777777778, 0.5298566133333334], [1.377777777777778, 0.5148124533333334], [1.4377777777777778, 0.5013753333333334], [1.4977777777777779, 0.48954525333333343], [1.557777777777778, 0.4793222133333334], [1.6177777777777778, 0.4707062133333334], [1.677777777777778, 0.4636972533333334], [1.7377777777777779, 0.4582953333333334], [1.797777777777778, 0.4545004533333334], [1.8577777777777778, 0.4523126133333334], [1.9177777777777778, 0.4517318133333334], [1.9777777777777779, 0.4527580533333334], [2.037777777777778, 0.45539133333333337], [2.097777777777778, 0.4596316533333334], [2.1577777777777776, 0.46547901333333336], [2.217777777777778, 0.4729334133333334], [2.2777777777777777, 0.48199485333333336], [2.337777777777778, 0.49266333333333334], [2.397777777777778, 0.5049388533333333], [2.457777777777778, 0.5188214133333333], [2.517777777777778, 0.5343110133333333], [2.577777777777778, 0.5514076533333334], [2.637777777777778, 0.5701113333333334], [2.6977777777777776, 0.5904220533333333], [2.7577777777777777, 0.6123398133333333], [2.817777777777778, 0.6358646133333334], [2.8777777777777778, 0.6609964533333333], [2.937777777777778, 0.6877353333333333]], [[0.6722222222222223, 0.9579797777777779], [0.7322222222222223, 0.9232056977777778], [0.7922222222222223, 0.8900386577777778], [0.8522222222222222, 0.8584786577777779], [0.9122222222222223, 0.8285256977777778], [0.9722222222222223, 0.8001797777777779], [1.0322222222222224, 0.7734408977777778], [1.0922222222222224, 0.7483090577777778], [1.1522222222222223, 0.7247842577777779], [1.212222222222222, 0.7028664977777779], [1.2722222222222221, 0.6825557777777779], [1.3322222222222222, 0.6638520977777779], [1.3922222222222222, 0.6467554577777779], [1.4522222222222223, 0.6312658577777779], [1.5122222222222224, 0.6173832977777778], [1.5722222222222224, 0.6051077777777778], [1.6322222222222222, 0.5944392977777778], [1.6922222222222223, 0.5853778577777778], [1.7522222222222221, 0.5779234577777779], [1.8122222222222224, 0.5720760977777778], [1.8722222222222222, 0.5678357777777778], [1.9322222222222223, 0.5652024977777779], [1.9922222222222221, 0.5641762577777778], [2.0522222222222224, 0.5647570577777778], [2.112222222222222, 0.5669448977777778], [2.1722222222222225, 0.5707397777777778], [2.232222222222222, 0.5761416977777778], [2.292222222222222, 0.5831506577777779], [2.352222222222222, 0.5917666577777778], [2.412222222222222, 0.6019896977777778], [2.4722222222222223, 0.6138197777777779], [2.532222222222222, 0.6272568977777778], [2.5922222222222224, 0.6423010577777778], [2.652222222222222, 0.6589522577777778], [2.7122222222222225, 0.6772104977777778], [2.772222222222222, 0.6970757777777778], [2.8322222222222218, 0.7185480977777776], [2.8922222222222222, 0.7416274577777777], [2.9522222222222227, 0.7663138577777779], [3.0122222222222224, 0.7926072977777777], [3.072222222222222, 0.8205077777777777]], [[0.8066666666666668, 1.060636666666667], [0.8666666666666667, 1.0270241866666667], [0.9266666666666667, 0.9950187466666667], [0.9866666666666668, 0.9646203466666667], [1.0466666666666669, 0.9358289866666667], [1.1066666666666667, 0.9086446666666668], [1.1666666666666667, 0.8830673866666667], [1.2266666666666668, 0.8590971466666667], [1.2866666666666666, 0.8367339466666668], [1.3466666666666667, 0.8159777866666668], [1.4066666666666667, 0.7968286666666669], [1.4666666666666668, 0.7792865866666668], [1.5266666666666668, 0.7633515466666667], [1.586666666666667, 0.7490235466666668], [1.646666666666667, 0.7363025866666667], [1.7066666666666668, 0.7251886666666667], [1.7666666666666666, 0.7156817866666667], [1.8266666666666667, 0.7077819466666667], [1.8866666666666667, 0.7014891466666667], [1.9466666666666668, 0.6968033866666667], [2.006666666666667, 0.6937246666666668], [2.066666666666667, 0.6922529866666668], [2.1266666666666665, 0.6923883466666667], [2.1866666666666665, 0.6941307466666667], [2.2466666666666666, 0.6974801866666668], [2.3066666666666666, 0.7024366666666667], [2.3666666666666667, 0.7090001866666668], [2.4266666666666667, 0.7171707466666667], [2.486666666666667, 0.7269483466666667], [2.5466666666666664, 0.7383329866666667], [2.606666666666667, 0.7513246666666668], [2.6666666666666665, 0.7659233866666667], [2.7266666666666666, 0.7821291466666667], [2.7866666666666666, 0.7999419466666667], [2.8466666666666667, 0.8193617866666667], [2.9066666666666667, 0.8403886666666668], [2.9666666666666663, 0.8630225866666665], [3.0266666666666664, 0.8872635466666666], [3.086666666666667, 0.9131115466666668], [3.1466666666666665, 0.9405665866666666], [3.2066666666666666, 0.9696286666666667]], [[0.9411111111111112, 1.1796420000000003], [1.0011111111111113, 1.14719112], [1.0611111111111113, 1.11634728], [1.1211111111111112, 1.0871104800000002], [1.1811111111111112, 1.05948072], [1.2411111111111113, 1.0334580000000002], [1.301111111111111, 1.00904232], [1.3611111111111112, 0.9862336800000001], [1.4211111111111112, 0.9650320800000002], [1.4811111111111113, 0.9454375200000001], [1.5411111111111113, 0.9274500000000001], [1.6011111111111112, 0.9110695200000001], [1.6611111111111112, 0.8962960800000002], [1.7211111111111113, 0.8831296800000001], [1.7811111111111113, 0.8715703200000001], [1.8411111111111111, 0.8616180000000001], [1.9011111111111112, 0.8532727200000001], [1.9611111111111112, 0.8465344800000001], [2.021111111111111, 0.8414032800000001], [2.0811111111111114, 0.8378791200000001], [2.141111111111111, 0.8359620000000001], [2.2011111111111115, 0.8356519200000001], [2.261111111111111, 0.8369488800000001], [2.321111111111111, 0.8398528800000001], [2.381111111111111, 0.84436392], [2.4411111111111112, 0.8504820000000001], [2.5011111111111113, 0.8582071200000001], [2.561111111111111, 0.86753928], [2.6211111111111114, 0.8784784800000001], [2.681111111111111, 0.89102472], [2.7411111111111115, 0.905178], [2.801111111111111, 0.9209383200000001], [2.861111111111111, 0.93830568], [2.921111111111111, 0.9572800800000001], [2.9811111111111113, 0.97786152], [3.0411111111111113, 1.00005], [3.101111111111111, 1.02384552], [3.161111111111111, 1.04924808], [3.2211111111111115, 1.07625768], [3.281111111111111, 1.10487432], [3.341111111111111, 1.135098]], [[1.0755555555555558, 1.3149957777777779], [1.1355555555555559, 1.283706497777778], [1.195555555555556, 1.2540242577777778], [1.2555555555555558, 1.2259490577777779], [1.3155555555555558, 1.1994808977777778], [1.3755555555555559, 1.174619777777778], [1.4355555555555557, 1.151365697777778], [1.4955555555555557, 1.1297186577777778], [1.5555555555555558, 1.109678657777778], [1.6155555555555559, 1.091245697777778], [1.675555555555556, 1.074419777777778], [1.7355555555555557, 1.059200897777778], [1.7955555555555558, 1.045589057777778], [1.8555555555555558, 1.0335842577777778], [1.915555555555556, 1.0231864977777778], [1.9755555555555557, 1.014395777777778], [2.035555555555556, 1.0072120977777779], [2.095555555555556, 1.001635457777778], [2.155555555555556, 0.9976658577777779], [2.215555555555556, 0.9953032977777779], [2.275555555555556, 0.9945477777777779], [2.3355555555555556, 0.9953992977777779], [2.3955555555555557, 0.9978578577777779], [2.4555555555555557, 1.0019234577777778], [2.5155555555555558, 1.0075960977777778], [2.575555555555556, 1.014875777777778], [2.635555555555556, 1.023762497777778], [2.695555555555556, 1.0342562577777779], [2.755555555555556, 1.0463570577777779], [2.8155555555555556, 1.0600648977777778], [2.8755555555555556, 1.0753797777777778], [2.9355555555555557, 1.0923016977777777], [2.9955555555555557, 1.1108306577777778], [3.055555555555556, 1.130966657777778], [3.115555555555556, 1.1527096977777778], [3.175555555555556, 1.176059777777778], [3.2355555555555555, 1.2010168977777778], [3.2955555555555556, 1.2275810577777777], [3.355555555555556, 1.255752257777778], [3.4155555555555557, 1.2855304977777777], [3.4755555555555557, 1.3169157777777778]], [[1.2100000000000002, 1.4666980000000005], [1.2700000000000002, 1.4365703200000004], [1.33, 1.4080496800000002], [1.3900000000000001, 1.3811360800000003], [1.4500000000000002, 1.3558295200000003], [1.5100000000000002, 1.3321300000000003], [1.5700000000000003, 1.3100375200000003], [1.6300000000000003, 1.2895520800000002], [1.6900000000000002, 1.2706736800000002], [1.75, 1.2534023200000002], [1.81, 1.2377380000000002], [1.87, 1.2236807200000004], [1.9300000000000002, 1.2112304800000002], [1.9900000000000002, 1.2003872800000002], [2.0500000000000003, 1.1911511200000002], [2.1100000000000003, 1.1835220000000002], [2.17, 1.1774999200000003], [2.2300000000000004, 1.1730848800000002], [2.29, 1.1702768800000003], [2.3500000000000005, 1.1690759200000003], [2.41, 1.1694820000000004], [2.47, 1.1714951200000003], [2.5300000000000002, 1.1751152800000002], [2.59, 1.1803424800000002], [2.6500000000000004, 1.1871767200000003], [2.71, 1.1956180000000003], [2.7700000000000005, 1.2056663200000002], [2.83, 1.2173216800000002], [2.8900000000000006, 1.2305840800000003], [2.95, 1.24545352], [3.0100000000000002, 1.2619300000000002], [3.0700000000000003, 1.2800135200000002], [3.13, 1.2997040800000002], [3.1900000000000004, 1.3210016800000002], [3.25, 1.3439063200000003], [3.3100000000000005, 1.3684180000000001], [3.37, 1.39453672], [3.4299999999999997, 1.42226248], [3.49, 1.4515952800000003], [3.55, 1.48253512], [3.6100000000000003, 1.515082]]], "point": [0.851, 0.4147144800000001], "foot": [0.851, 0.1932], "patch": [[0.5409999999999999, 0.44743288000000014], [0.9410000000000001, 0.2955416800000001], [1.161, 0.38199608000000007], [0.761, 0.5338872800000001]], "foot_x": [0.62, 0.0], "foot_y": [0.231, 0.1932], "arrow": [[0.851, 0.4147144800000001], [1.3292250075392191, 0.329572099329455]], "base_arrow": [[0.851, 0.1932], [1.3292250075392191, 0.2975995777193376]], "path": [[0.851, 0.4147144800000001], [1.3648532000000002, 0.370882397845248], [1.6366441856, 0.39377742678565375], [1.7811214861683202, 0.42251808729536683], [1.8584605990377883, 0.44536575445770304], [1.9002615325622112, 0.46191173814986136], [1.9231510806390935, 0.4736913037996944]], "minimum": [1.955, 0.5060000000000001], "minimum_foot": [1.955, 0.5060000000000001], "axes": {"o": [0.0, 0.0], "x": [2.55, 0.0], "y": [1.2925000000000002, 1.0810000000000002], "z": [0.0, 0.576]}, "base": [[0.0, 0.0], [2.4, 0.0], [3.6100000000000003, 1.0120000000000002], [1.2100000000000002, 1.0120000000000002]]} \ No newline at end of file diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface.typ b/articles/symbolic-fem-and-adjoint-forms/figures/src/misfit-surface.typ new file mode 100644 index 0000000..402de68 --- /dev/null +++ b/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") +}) diff --git a/articles/symbolic-fem-and-adjoint-forms/figures/src/uwfig.typ b/articles/symbolic-fem-and-adjoint-forms/figures/src/uwfig.typ new file mode 100644 index 0000000..d538951 --- /dev/null +++ b/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) +} diff --git a/articles/symbolic-fem-and-adjoint-forms/metadata.yml b/articles/symbolic-fem-and-adjoint-forms/metadata.yml new file mode 100644 index 0000000..99af402 --- /dev/null +++ b/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 .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 Logan Voss / Unsplash" diff --git a/articles/symbolic-fem-and-adjoint-forms/references.bib b/articles/symbolic-fem-and-adjoint-forms/references.bib new file mode 100644 index 0000000..f374028 --- /dev/null +++ b/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} +} diff --git a/articles/symbolic-fem-and-adjoint-forms/symbolic-fem-and-adjoint-forms.md b/articles/symbolic-fem-and-adjoint-forms/symbolic-fem-and-adjoint-forms.md new file mode 100644 index 0000000..5c7215f --- /dev/null +++ b/articles/symbolic-fem-and-adjoint-forms/symbolic-fem-and-adjoint-forms.md @@ -0,0 +1,410 @@ +--- +title: Symbolic FEM and Adjoint Forms +description: >- + Every implicit solve in Underworld3 is a residual written in SymPy, and its + Newton solver assembles the tangent stiffness from four symbolic derivative + blocks. We show that the same blocks give the adjoint: one transposed linear + solve, and an integral for each parameter, yield the derivative of a misfit + with respect to every parameter of the model at once. +date: 2026-10-25 +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliations: + - Australian National University +license: CC-BY-4.0 +banner: figures/banner.jpg +bibliography: + - references.bib +keywords: + - Underworld Code + - adjoint + - inversion + - data assimilation +exports: + - format: typst + logo: ../../static/uwtn-logo.png + series: "Underworld Technical Notes" + origin_url: https://www.underworldcode.org/symbolic-fem-and-adjoint-forms/ + template: ../../templates/pdf + output: symbolic-fem-and-adjoint-forms.pdf + article_id: UWTN 2026-020 + article_version: 1.0.0 + software_version: underworld3 development @ f414da63 +--- +
Photo by Logan Voss / Unsplash
+ +In this article, we will see how the underworld3 sympy formulation combined with the PETSc finite element, strong-form framework builds model derivatives which can solve non-linear problems using a Newton method. This requires us to calculate the local derivative of the model with respect to the unknowns so that we know the direction to move towards a better solution. + +If we differentiate the prediction a model makes with respect to one of its parameters, then we can determine how sensitive the prediction is to changes in that parameter and therefore whether the model can be used to constrain the value of that parameter. Similarly, model derivatives are used in data assimilation to determine how to adjust a model's parameters so it better fits the observations we have [@Bunge_2003; @Liu_2008]. + +## The forward problem, solved by Newton + +Underworld3 starts from the strong form of the equations which is usually how they are written in a textbook or a paper. These are symbolic, sympy objects that pass directly to compiled objects in PETSc [@Moresi_2025; @Knepley_2013]: everything +inside a divergence term goes into the flux $\mathbf F_1$, everything else into the +force $F_0$. The discrete solution is written in the finite-element basis, +$u = \sum_j u_j \phi_j$, and the weak form is one equation for each basis +function $\phi_i$: the integral of the two templates against it and its +gradient, + +```{math} +:label: eq-residual +R_i(u;\, m) \;=\; +\int_\Omega F_0(u, \nabla u;\, m)\,\phi_i +\;+\; \mathbf F_1(u, \nabla u;\, m) \cdot \nabla\phi_i \;=\; 0 . +``` +$R$, the residual {eq}`eq-residual`, is the **out-of-balance force**: what the +assembled system has left over when the current $u$ is put into the governing +equation. A solution is a $u$ that makes it vanish. Here $m$ is a set of parameters of the model (for example a Rayleigh number, or a constitutive modulus) that the user typically would adjust to run an example. + +A Newton non-linear solver needs the derivative of $R$ with respect to the unknowns +[@Moresi_2026_sympy]: + +```{math} +:label: eq-blocks +G_0 = \frac{\partial F_0}{\partial u}, \quad +G_1 = \frac{\partial F_0}{\partial (\nabla u)}, \quad +G_2 = \frac{\partial \mathbf F_1}{\partial u}, \quad +G_3 = \frac{\partial \mathbf F_1}{\partial (\nabla u)} , +``` + +which underworld computes directly from sympy, compiles and registers with PETSc. The Jacobian is the derivative of each equation $R_i$ with respect to each nodal unknown $u_j$. Since $\partial u/\partial u_j = \phi_j$ and $\partial (\nabla u)/\partial u_j = \nabla\phi_j$, the chain rule turns the four blocks into $K_{ij} = \partial R_i/\partial u_j$, assembled element by element, + +```{math} +:label: eq-jacobian +K_{ij} = \int_\Omega +\phi_i \, G_0 \, \phi_j +\;+\; \phi_i \, G_1 \cdot \nabla\phi_j +\;+\; \nabla\phi_i \cdot G_2 \, \phi_j +\;+\; \nabla\phi_i \cdot G_3 \cdot \nabla\phi_j . +``` + +$K = \partial R/\partial u$, assembled by {eq}`eq-jacobian`, is the **tangent +stiffness matrix**. For a linear problem it is just the stiffness matrix, +and + +```{math} +:label: eq-stiffness +R = Ku - f +``` + +For a nonlinear problem, $K(u)$ is the stiffness of the material at +the current solution state. Newton uses it to reduce the out-of-balance +force: expanding $R(u + \delta u) \approx R(u) + K\,\delta u$ and asking for +the left side to vanish gives the correction + +```{math} +:label: eq-newton +K\,\delta u = -R(u), \qquad u \leftarrow u + \delta u , +``` + +repeated until $R$ is zero to the tolerance asked for. For a linear problem +{eq}`eq-stiffness` makes the expansion exact, and one step gives the solution. + +## The adjoint problem + +Now suppose we want to choose the parameters $m$ so that the model reproduces +some observations as well as it can. The observations are values +$d_i^{\mathrm{obs}}$ — a surface velocity at a station, or a stress at a point — +and the model predicts the same quantities, $d_i(u)$, from its solution. The +quality of the fit is the misfit + +```{math} +:label: eq-misfit +J(u) = \frac12 \sum_i \big(d_i(u) - d_i^{\mathrm{obs}}\big)^2 , +``` + +and the task is to minimise this by adjusting $m$, under the +constraint that $u$ always solves the model: $R(u;\, m) = 0$ at every $m$ we +try. To do that we need the gradient $\mathrm dJ/\mathrm dm$ (which way, and +how far, to move each parameter) and the adjoint approach provides this for +all the parameters at once - at the expense of one additional *linear* solve. + +```{figure} figures/misfit-surface.svg +:name: fig-misfit-surface +:width: 70% +:alt: A bowl-shaped surface drawn as a grid above the parameter plane with axes m_1 to the right, m_2 towards the rear, and J upward. A red dot on the surface labelled J(m) marks the current parameters, with a dashed vertical line down to the plane beneath it. The tangent plane touches the surface at the dot, and a rust arrow labelled minus grad J points down the slope from it. A dashed line with small open circles runs from the dot down the bowl to a black dot at its lowest point, labelled m-star. + +The misfit $J$ as a surface over two parameters (in this case). A forward solve locates the model at one point on this surface. The adjoint solve gives the tangent plane at that location. Following the gradient downhill reaches $m^*$, where the model best fits the observations. +``` + +The constraint is added to the minimisation with a Lagrange multiplier +$\mu$ [@Plessix_2006]: + +```{math} +:label: eq-lagrangian +\mathcal L(u, \mu;\, m) = J(u;\, m) + \mu^{T} R(u;\, m) , +``` + +which equals $J$ for any $\mu$ whenever the model is solved, because $R = 0$ +there. We are free to choose $\mu$, and we choose it so that $\mathcal L$ does +not change, to first order, when the solution changes: +$\partial\mathcal L/\partial u = \partial J/\partial u + \mu^{T}K = 0$. +Transposed, that is the **adjoint problem**, + +```{math} +:label: eq-adjoint +K^{T}\mu = -\Big(\frac{\partial J}{\partial u}\Big)^{T}, +``` + +a linear solve for a field $\mu$ on the same mesh as $u$. It is linear because +$u$ is already known: $K$ and $\partial J/\partial u$ are evaluated at the +converged solution and held fixed, and $\mu$ appears only multiplied by them. +The operator is the tangent stiffness $K$ at that solution, transposed. The +force is the +misfit's derivative with respect to the solution, + +```{math} +:label: eq-djdu +\frac{\partial J}{\partial u_j} += \sum_i \big(d_i(u) - d_i^{\mathrm{obs}}\big)\,\frac{\partial d_i}{\partial u_j} , +``` + +one entry per nodal unknown: each misfit $d_i - d_i^{\mathrm{obs}}$, placed at +the nodes its prediction depends on. This is **the misfit written as a force**, +assembled like any finite-element force vector, and it is the only place the +observations enter. With $\mu$ chosen this way, the solution's response to $m$ drops out of +$\mathrm d\mathcal L/\mathrm dm$, and what is left is the gradient for each +parameter, + +```{math} +:label: eq-gradient +\frac{\mathrm d J}{\mathrm d m_k} = \frac{\partial J}{\partial m_k} + \mu^{T}\frac{\partial R}{\partial m_k} , +``` + +where $\partial R/\partial m_k$ is a **pseudo-force**: the out-of-balance +force that appears when $m_k$ is nudged and the solution held still. It is +an integral for each parameter, not a solve. + +[The appendix](#sec-derivation) reaches {eq}`eq-adjoint` and +{eq}`eq-gradient` a second way, by differentiating the constraint directly, +which shows what each term means. + +Reducing $J$ is then a loop, [](#fig-adjoint-loop): (1) solve the model, +(2) solve the adjoint problem, (3) form the gradient for every parameter, +(4) update the parameters and go round again. In practice a quasi-Newton method such as +L-BFGS chooses the step $\alpha$ and builds up curvature from successive +gradients. + +```{figure} figures/adjoint-loop.svg +:name: fig-adjoint-loop +:width: 78% +:alt: Four rounded boxes in a square, numbered 1 to 4 and joined by arrows running clockwise. Top left, 1. Solve the model, R(u; m) = 0 by Newton. An arrow labelled K leads right to 2. Solve the adjoint, K-transpose mu equals minus the transpose of dJ/du, into which a heavier rust-coloured arrow labelled dJ/du comes down from the words observations d-obs above. A teal arrow labelled mu leads down to 3. Form the gradient, dJ/dm_k as an integral. An arrow labelled dJ/dm leads left to 4. Update the parameters, m becomes m minus alpha grad J, and an arrow labelled m leads up to close the loop. A dashed rust-coloured arrow labelled dR/dm_k runs diagonally from the model box to the gradient box. + +One pass of an inversion. Teal carries a field, rust a force, black a quantity handed on. The tangent stiffness $K$ assembled by the forward solve is the operator of the adjoint solve, and the pseudo-force $\partial R/\partial m_k$ comes from the same residual; the observations enter only as the adjoint's load. +``` + +Each pass costs one nonlinear solve and one linear one, irrespective of the number of +parameters there are, and the fit improves until the gradient vanishes, at a +minimum of $J$ where no small change of $m$ brings the model closer to the +observations. The minimum may be a local one: the gradient describes the +surface only where the model is. In +Underworld3, with the misfit `J` written as an integrand over the mesh, +steps 2 and 3 are a single call, + +```python +dJdm = stokes.gradient(J, parameters=[m1, m2, m3])["parameters"] +dJdm[m1] # dJ/dm1 +``` + +and step 4 is whatever optimiser one prefers. + +## Re-using the forward solve + +The adjoint problem {eq}`eq-adjoint` asks for two things: the transposed stiffness $K^{T}$, and the pseudo-force $\partial R/\partial m$. Neither is a new object. Both are the blocks {eq}`eq-blocks` that Newton already uses, read differently. + +**The transposed stiffness is a relabelling.** $K^{T}$ is the bilinear form of +{eq}`eq-jacobian` with trial and test exchanged. Exchanging $\phi_i$ and +$\phi_j$ leaves the two diagonal terms in place, transposed on their paired +indices, and swaps the two mixed ones — the $G_1$ term becomes a $G_2$ term +and the other way about, with their components exchanged: + +```{math} +:label: eq-transpose +G_0 \;\rightarrow\; G_0^{T}, \qquad +G_3 \;\rightarrow\; G_3^{T}, \qquad +G_1 \;\rightarrow\; G_2^{T}, \qquad +G_2 \;\rightarrow\; G_1^{T} . +``` + +So the adjoint operator is assembled from the *same four compiled kernels*, +registered in a different arrangement. We do that wherever the kernels carry +the whole operator, rather than transposing an assembled matrix, which matters +for Stokes, where the $(u,p)$ and $(p,u)$ blocks exchange places by the same +rule and the composite operator stays a composite operator the Krylov solver +can be preconditioned against. Nothing is differentiated a second time, and the +assembled operator agrees with the explicit transpose to machine precision. + +**The parameter derivative is the same `diff`, on a different symbol.** +A named quantity $m$ — a friction coefficient, a reference viscosity, a +continuation parameter — survives into $F_0$ and $\mathbf F_1$ as a symbol, +because the solver carries the constitutive model's own expression rather than +the number it currently holds. Differentiating {eq}`eq-residual` with respect +to $m$ instead of $u$ is the same operation that produced {eq}`eq-blocks`: + +```{math} +:label: eq-drdm +\frac{\partial R}{\partial m} \;=\; +\int_\Omega \frac{\partial F_0}{\partial m}\,\phi +\;+\; \frac{\partial \mathbf F_1}{\partial m}\cdot\nabla\phi . +``` + +## Accuracy of the gradient + +Three things decide the accuracy of the gradient. The forward solve must be converged before +the adjoint is taken, because $\mu$ is the multiplier for the constraint +$R = 0$ and has meaning only where that holds. The adjoint solve itself can be +run to a loose tolerance: every pass re-solves the forward problem and +recomputes $J$, so a direction that is slightly off costs an extra iteration +and the next misfit shows whether it helped. And the matrix must be the true +derivative $\partial R/\partial u$. With a frozen-viscosity Picard matrix, $\mu$ +solves its own equation exactly and the gradient it gives belongs to no +function, so the optimiser can stop where the true gradient is still large. +The tangent is a property of the code rather than of the iterate, so it is +checked once, by finite differences, and assembling it symbolically from the +same residual is what makes it right by construction. + ++++ {"part": "appendix"} + +(sec-derivation)= +## Derivation of the adjoint + +The adjoint problem {eq}`eq-adjoint` and the gradient {eq}`eq-gradient` also +follow from differentiating the constraint directly, without the Lagrangian +{eq}`eq-lagrangian`, and this route shows what each term means. + +We explicitly wrote $R$ in terms of $m$ to demonstrate that there is an entire family of solutions to $R=0$ that we can access by changing the model's free parameters. If the model makes a prediction, $J(u)$, we can ask how changing $m$ affects that prediction (the sensitivity of the prediction to $m$). That is + +```{math} +:label: eq-question +\frac{\mathrm d J}{\mathrm d m} +``` + +If $J$ is an underworld function then we can form a symbolic expression. First we note that $u$ and $m$ are not independent: the solution is constrained to the family of solutions of {eq}`eq-residual`, $R = 0$ holds *for every* $m$ and, because of this, $dR/dm$ is zero too + +```{math} +:label: eq-total +\frac{\mathrm d R}{\mathrm d m} += \frac{\partial R}{\partial u}\frac{\mathrm d u}{\mathrm d m} ++ \frac{\partial R}{\partial m} = 0 . +``` + +The first factor is the tangent stiffness $K$ again, so + +```{math} +:label: eq-implicit +K\,\frac{\mathrm d u}{\mathrm d m} = -\frac{\partial R}{\partial m}, +\qquad \text{so} \qquad +\frac{\mathrm d u}{\mathrm d m} = -K^{-1}\frac{\partial R}{\partial m}, +``` + +$\partial R/\partial m$ is the pseudo-force of {eq}`eq-gradient`: the out-of-balance force that +appears if the parameter is nudged and the solution is held still. It is +assembled exactly as $R$ is, with $\partial F_0/\partial m$ and +$\partial \mathbf F_1/\partial m$ in place of $F_0$ and $\mathbf F_1$. + +We interpret {eq}`eq-implicit` as: **the sensitivity of the +solution to a parameter is the response of the structure to that parameter's +pseudo-force.** + +The prediction depends on $m$ in the same two ways: directly, and through the +solution. With $\partial J/\partial u$ the row of derivatives +$\partial J/\partial u_j$ with respect to the nodal unknowns — the prediction +written as a force — the chain rule and {eq}`eq-implicit` give + +```{math} +:label: eq-chain +\frac{\mathrm d J}{\mathrm d m} += \frac{\partial J}{\partial m} ++ \frac{\partial J}{\partial u}\frac{\mathrm d u}{\mathrm d m} += \frac{\partial J}{\partial m} +- \frac{\partial J}{\partial u}\,K^{-1}\frac{\partial R}{\partial m} . +``` + +### Many parameters at once + +A real model has many parameters, $m_1, \dots, m_N$, and we want the +sensitivity for all of them. We first write {eq}`eq-chain` like this: + +```{math} +:label: eq-chain-k +\frac{\mathrm d J}{\mathrm d m_k} += \frac{\partial J}{\partial m_k} +- \frac{\partial J}{\partial u}\,K^{-1}\,\frac{\partial R}{\partial m_k}, +\qquad k = 1, \dots, N . +``` + +These are $N$ equations that differ only in their last factor, the pseudo-force $\partial R/\partial m_k$. The row in front of it, $\big(\partial J/\partial u\big)K^{-1}$, contains the prediction and the stiffness but no parameter; it is the same for every $k$ and we only need to compute it once. + +Now define $\mu$ through + +```{math} +:label: eq-mu-row +\mu^{T} = -\frac{\partial J}{\partial u}\,K^{-1} . +``` + +We do not have access to $K^{-1}$ so we multiply both sides by $K$ on +the right, take the transpose and substitute into {eq}`eq-chain-k` +to recover {eq}`eq-adjoint`: + +```{math} +K^{T}\mu = -\Big(\frac{\partial J}{\partial u}\Big)^{T}, \qquad +\frac{\mathrm d J}{\mathrm d m_k} = \frac{\partial J}{\partial m_k} + \mu^{T}\frac{\partial R}{\partial m_k} . +``` + +One linear solve for $\mu$, whatever $N$ is, and then one product for each +parameter. Written as fields, with $\mu = \sum_i \mu_i \phi_i$, that product is + +```{math} +:label: eq-mu-test +\mu^{T}\frac{\partial R}{\partial m_k} += \int_\Omega \frac{\partial F_0}{\partial m_k}\,\mu +\;+\; \frac{\partial \mathbf F_1}{\partial m_k}\cdot\nabla\mu , +``` + +which is the pseudo-force's weak form with $\mu$ in the place of the test function. + +$\mu$, the **adjoint state**, is the response of the same structure to the +misfit's force: an influence field, in the sense of an influence line. It says +how strongly an out-of-balance force applied anywhere in the domain would tell +on the misfit. With it in hand, the derivative for a parameter is +$\mu^{T}\partial R/\partial m$ — the virtual work that parameter's pseudo-force +does against the influence field — and that is an integral, not a solve. + +Solving {eq}`eq-implicit` for each $\mathrm d u/\mathrm d m_k$ and contracting it +with $\partial J/\partial u$ gives the same numbers, at the cost of $N$ solves. +That the two agree is Maxwell–Betti reciprocity: the +work done by one force through the displacements of another equals the work +done by the second through the displacements of the first. The transpose is +what makes it true for a problem whose stiffness is not symmetric. + +So the adjoint is cheap in parameters and dear in misfits. One solve for ten +parameters or a thousand; another solve for every additional misfit. That +suits a model with many dials and one thing you care about, which is the usual +shape of the problem. + +### One stiffness, three forces + +The Newton update, the forward sensitivity and the adjoint are three linear +solves with the same matrix and different right-hand sides: + +```{math} +:label: eq-family +K\,\delta u = -R, \qquad\quad +K\,\frac{\mathrm d u}{\mathrm d m} = -\frac{\partial R}{\partial m}, +\qquad\quad +K^{T}\mu = -\Big(\frac{\partial J}{\partial u}\Big)^{T} . +``` + +Newton has already assembled and preconditioned that matrix, so the adjoint +costs one more linear solve. + ++++ + + + +*The Underworld project is supported by AuScope and the Australian Government through the National Collaborative Research Infrastructure Strategy (NCRIS). Source code:* [*github.com/underworldcode/underworld3*](https://github.com/underworldcode/underworld3) + + + +
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 60749b4..9bf4824 100644 --- a/classification.yml +++ b/classification.yml @@ -331,3 +331,7 @@ boundary-conditions-on-non-planar-boundaries: article_type: technical-note subjects: [surface-processes] methods: [solvers, finite-elements, meshing] +symbolic-fem-and-adjoint-forms: + article_type: technical-note + subjects: [] + methods: [solvers, symbolic-codegen] diff --git a/pages-src/submit.md b/pages-src/submit.md index f42c987..c329764 100644 --- a/pages-src/submit.md +++ b/pages-src/submit.md @@ -76,7 +76,9 @@ the working version. That is deliberate, and it is why the series is on a repository that versions. **Commit your figure sources, not just the pictures.** If a figure is drawn from -data, the script and the data belong beside it. Sixteen figures in the older +data, the script and the data belong beside it, in +`articles//figures/src/` — not in `examples/`, which is for what a reader +runs. Sixteen figures in the older material on this site went dark when the server they were linked from was retired — and not one of them was recoverable from the Internet Archive. They came back only because copies existed elsewhere, which was luck rather than diff --git a/templates/article-template/ARTICLE.md b/templates/article-template/ARTICLE.md index 3865815..ec3e559 100644 --- a/templates/article-template/ARTICLE.md +++ b/templates/article-template/ARTICLE.md @@ -107,8 +107,10 @@ PDF, and produce one references section. ## Examples -Notebooks and small data files go in `examples/` beside this file. They are -deposited with the note, so the archived record carries working code. +Notebooks a reader is meant to run — the ones a Binder link opens — and the +small data files they read go in `examples/` beside this file. The scripts that +draw the figures go in `figures/src/`. Both are deposited with the note, so the +archived record carries working code. When a release breaks one, update it and publish a **new version** of the same deposit rather than a new note — the DOI already in circulation then resolves to diff --git a/templates/pdf/lapreprint.typ b/templates/pdf/lapreprint.typ index eef1e34..bb789d7 100644 --- a/templates/pdf/lapreprint.typ +++ b/templates/pdf/lapreprint.typ @@ -65,6 +65,8 @@ // The path to a bibliography file if you want to cite some external works. bibliography-file: none, bibliography-style: "apa", + // Printed after the references, as a paper sets it. + appendix: none, // The paper's content. body ) = { @@ -209,7 +211,7 @@ #show: smallcaps #v(20pt, weak: true) #if it.numbering != none and not is-ack { - numbering(heading-numbering, ..levels) + numbering(it.numbering, ..levels) [.] h(7pt, weak: true) } @@ -221,7 +223,7 @@ #set text(style: "italic") #v(10pt, weak: true) #if it.numbering != none { - numbering(heading-numbering, ..levels) + numbering(it.numbering, ..levels) [.] h(7pt, weak: true) } @@ -230,7 +232,7 @@ ] else [ // Third level headings are run-ins too, but different. #if it.level == 3 { - numbering(heading-numbering, ..levels) + numbering(if it.numbering != none { it.numbering } else { heading-numbering }, ..levels) [. ] } _#(it.body):_ @@ -420,4 +422,16 @@ if (abstracts != none and abstracts.len() > 0) { show bibliography: set text(8pt) bibliography(bibliography-file, title: text(10pt, "References"), style: bibliography-style) } + + // Lettered, and counted afresh, so that A.1 cannot be read as a section of + // the body. The show rule above prints `it.numbering` for this reason. + if appendix != none { + counter(heading).update(0) + set heading(numbering: (..n) => if n.pos().len() == 1 { + "Appendix " + numbering("A", ..n) + } else { + numbering("A.1", ..n) + }) + appendix + } } diff --git a/templates/pdf/template.typ b/templates/pdf/template.typ index c382753..1a2bb67 100644 --- a/templates/pdf/template.typ +++ b/templates/pdf/template.typ @@ -56,6 +56,11 @@ [# endif #] [# if doc.bibtex #] bibliography-file: "[-doc.bibtex-]", +[# endif #] +[# if parts.appendix #] + appendix: [ +[-parts.appendix-] + ], [# endif #] authors: ( [# for author in doc.authors #] diff --git a/templates/pdf/template.yml b/templates/pdf/template.yml index 2980885..33ba500 100644 --- a/templates/pdf/template.yml +++ b/templates/pdf/template.yml @@ -93,6 +93,12 @@ parts: abstract will prepare readers for the detailed information to follow, communicate the essence of the article and help readers take away and remember key points. max_words: 500 + - id: appendix + description: > + Material the article points to but does not need inline, such as a + formal derivation. Printed after the references, with lettered + headings. Marked in the source with a `+++ {"part": "appendix"}` block, + not a heading, so no ordinary section is ever lifted into it. # summary / acknowledgements / availability are NOT declared. # # A declared part is filled from a body section with a matching heading, so diff --git a/tests/test_migration.py b/tests/test_migration.py index efea3c2..0d83387 100644 --- a/tests/test_migration.py +++ b/tests/test_migration.py @@ -211,12 +211,14 @@ def test_a_labelled_equation_is_referenceable(): import re for path in sorted((ROOT / "articles").glob("*/*.md")): body = path.read_text(encoding="utf-8") - # Two forms label an equation in this corpus, and both are valid MyST: - # `$$ ... $$ (label)`, and amsmath's \label{} inside \begin{equation}. - # A \label{} inside a plain $$ block is NEITHER -- that is LaTeX the + # Three forms label an equation in this corpus, and all are valid MyST: + # `$$ ... $$ (label)`, amsmath's \label{} inside \begin{equation}, and + # the `:label:` option of a {math} directive, which native notes use. + # A \label{} inside a plain $$ block is NONE of them -- that is LaTeX the # renderer never sees, and it is the mistake this test exists to catch. labels = set(re.findall(r"(?m)^\$\$\s*\(([^)\s]+)\)\s*$", body)) labels |= set(re.findall(r"\\begin\{equation\}\s*\\label\{([^}]+)\}", body)) + labels |= set(re.findall(r"(?m)^\s*:label:\s*(\S+)\s*$", body)) referenced = set(re.findall(r"\{eq\}`([^`]+)`", body)) missing = referenced - labels assert not missing, (