UWTN 2026-015: Testing a solver against exact solutions

articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0)_maxC(1.0, 1.0).msh

@@ -0,0 +1,125 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$PhysicalNames
+5
+1 11 "Bottom"
+1 12 "Top"
+1 13 "Right"
+1 14 "Left"
+2 99999 "Elements"
+$EndPhysicalNames
+$Entities
+4 4 1 0
+1 0 0 0 0 
+2 1 0 0 0 
+3 0 1 0 0 
+4 1 1 0 0 
+11 0 0 0 1 0 0 1 11 2 1 -2 
+12 0 1 0 1 1 0 1 12 2 4 -3 
+13 1 0 0 1 1 0 1 13 2 2 -4 
+14 0 0 0 0 1 0 1 14 2 3 -1 
+1 0 0 0 1 1 0 1 99999 4 11 13 12 14 
+$EndEntities
+$Nodes
+9 25 1 25
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+1 0 0
+0 3 0 1
+3
+0 1 0
+0 4 0 1
+4
+1 1 0
+1 11 0 3
+5
+6
+7
+0.2499999999994109 0 0
+0.4999999999986921 0 0
+0.7499999999993406 0 0
+1 12 0 3
+8
+9
+10
+0.7500000000003471 1 0
+0.5000000000020595 1 0
+0.2500000000010405 1 0
+1 13 0 3
+11
+12
+13
+1 0.2499999999994109 0
+1 0.4999999999986921 0
+1 0.7499999999993406 0
+1 14 0 3
+14
+15
+16
+0 0.7500000000003471 0
+0 0.5000000000020595 0
+0 0.2500000000010405 0
+2 1 0 9
+17
+18
+19
+20
+21
+22
+23
+24
+25
+0.2499999999998184 0.2500000000006332 0
+0.2500000000002257 0.5000000000012177 0
+0.2500000000006332 0.7500000000000954 0
+0.4999999999995339 0.2500000000002257 0
+0.5000000000003759 0.5000000000003757 0
+0.5000000000012175 0.7499999999998438 0
+0.7499999999995922 0.2499999999998184 0
+0.749999999999844 0.4999999999995339 0
+0.7500000000000953 0.7499999999995921 0
+$EndNodes
+$Elements
+5 32 1 32
+1 11 1 4
+1 1 5 
+2 5 6 
+3 6 7 
+4 7 2 
+1 12 1 4
+5 4 8 
+6 8 9 
+7 9 10 
+8 10 3 
+1 13 1 4
+9 2 11 
+10 11 12 
+11 12 13 
+12 13 4 
+1 14 1 4
+13 3 14 
+14 14 15 
+15 15 16 
+16 16 1 
+2 1 3 16
+17 1 5 17 16 
+18 16 17 18 15 
+19 15 18 19 14 
+20 14 19 10 3 
+21 5 6 20 17 
+22 17 20 21 18 
+23 18 21 22 19 
+24 19 22 9 10 
+25 6 7 23 20 
+26 20 23 24 21 
+27 21 24 25 22 
+28 22 25 8 9 
+29 7 2 11 23 
+30 23 11 12 24 
+31 24 12 13 25 
+32 25 13 4 8 
+$EndElements

articles/testing-a-solver-against-exact-solutions/examples/.meshes/uw_structuredQuadBox_minC(0.0, 0.0, 0.0)_maxC(1.0, 1.0, 1.0).msh

@@ -0,0 +1,169 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$PhysicalNames
+7
+2 11 "Bottom"
+2 12 "Top"
+2 13 "Right"
+2 14 "Left"
+2 15 "Front"
+2 16 "Back"
+3 99999 "Elements"
+$EndPhysicalNames
+$Entities
+8 12 6 1
+1 0 0 0 0 
+2 1 0 0 0 
+3 0 1 0 0 
+4 1 1 0 0 
+5 0 0 1 0 
+6 1 0 1 0 
+7 0 1 1 0 
+8 1 1 1 0 
+1 0 0 0 1 0 0 0 2 1 -2 
+2 1 0 0 1 1 0 0 2 2 -4 
+3 0 1 0 1 1 0 0 2 4 -3 
+4 0 0 0 0 1 0 0 2 3 -1 
+5 0 0 1 1 0 1 0 2 5 -6 
+6 1 0 1 1 1 1 0 2 6 -8 
+7 0 1 1 1 1 1 0 2 8 -7 
+8 0 0 1 0 1 1 0 2 7 -5 
+9 0 0 0 0 0 1 0 2 5 -1 
+10 1 0 0 1 0 1 0 2 2 -6 
+11 0 1 0 0 1 1 0 2 7 -3 
+12 1 1 0 1 1 1 0 2 4 -8 
+11 0 0 0 1 1 0 1 11 4 1 2 3 4 
+12 0 0 1 1 1 1 1 12 4 5 6 7 8 
+13 1 0 0 1 1 1 1 13 4 10 6 -12 -2 
+14 0 0 0 0 1 1 1 14 4 9 -4 -11 8 
+15 0 0 0 1 0 1 1 15 4 1 10 -5 9 
+16 0 1 0 1 1 1 1 16 4 -3 12 7 11 
+1 0 0 0 1 1 1 1 99999 6 15 13 16 12 14 11 
+$EndEntities
+$Nodes
+27 27 1 27
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+1 0 0
+0 3 0 1
+3
+0 1 0
+0 4 0 1
+4
+1 1 0
+0 5 0 1
+5
+0 0 1
+0 6 0 1
+6
+1 0 1
+0 7 0 1
+7
+0 1 1
+0 8 0 1
+8
+1 1 1
+1 1 0 1
+9
+0.4999999999986921 0 0
+1 2 0 1
+10
+1 0.4999999999986921 0
+1 3 0 1
+11
+0.5000000000020595 1 0
+1 4 0 1
+12
+0 0.5000000000020595 0
+1 5 0 1
+13
+0.4999999999986921 0 1
+1 6 0 1
+14
+1 0.4999999999986921 1
+1 7 0 1
+15
+0.5000000000020595 1 1
+1 8 0 1
+16
+0 0.5000000000020595 1
+1 9 0 1
+17
+0 0 0.5000000000020595
+1 10 0 1
+18
+1 0 0.4999999999986921
+1 11 0 1
+19
+0 1 0.5000000000020595
+1 12 0 1
+20
+1 1 0.4999999999986921
+2 11 0 1
+21
+0.5000000000003758 0.5000000000003758 0
+2 12 0 1
+22
+0.5000000000003758 0.5000000000003758 1
+2 13 0 1
+23
+1 0.499999999998692 0.4999999999986922
+2 14 0 1
+24
+0 0.5000000000020595 0.5000000000020595
+2 15 0 1
+25
+0.4999999999986921 0 0.5000000000003758
+2 16 0 1
+26
+0.5000000000020595 1 0.5000000000003756
+3 1 0 1
+27
+0.5000000000003757 0.5000000000003756 0.5000000000003759
+$EndNodes
+$Elements
+7 32 1 32
+2 11 3 4
+1 1 9 21 12 
+2 12 21 11 3 
+3 9 2 10 21 
+4 21 10 4 11 
+2 12 3 4
+5 5 13 22 16 
+6 16 22 15 7 
+7 13 6 14 22 
+8 22 14 8 15 
+2 13 3 4
+9 2 18 23 10 
+10 10 23 20 4 
+11 18 6 14 23 
+12 23 14 8 20 
+2 14 3 4
+13 5 17 24 16 
+14 16 24 19 7 
+15 17 1 12 24 
+16 24 12 3 19 
+2 15 3 4
+17 1 9 25 17 
+18 17 25 13 5 
+19 9 2 18 25 
+20 25 18 6 13 
+2 16 3 4
+21 3 11 26 19 
+22 19 26 15 7 
+23 11 4 20 26 
+24 26 20 8 15 
+3 1 5 8
+25 1 9 21 12 17 25 27 24 
+26 17 25 27 24 5 13 22 16 
+27 12 21 11 3 24 27 26 19 
+28 24 27 26 19 16 22 15 7 
+29 9 2 10 21 25 18 23 27 
+30 25 18 23 27 13 6 14 22 
+31 21 10 4 11 27 23 20 26 
+32 27 23 20 26 22 14 8 15 
+$EndElements

articles/testing-a-solver-against-exact-solutions/examples/banner.py

@@ -0,0 +1,61 @@
+"""Banner: three of the viscosity structures the exact solutions cover.
+
+Run from the repository root:
+
+    python3 articles/testing-a-solver-against-exact-solutions/examples/banner.py
+
+The note says the family covers "a viscosity jump, an exponentially varying
+viscosity, a laterally oscillating one". This is that sentence: SolCx's step,
+SolKz's exponential gradient, SolM's lateral oscillation, each evaluated from
+the solution's own symbolic `fn_viscosity` onto a mesh variable and rendered.
+
+Run against underworld3 `development` at commit `0addec15`
+(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__`
+reports 0.0.0 for every build, so the commit is the only thing that identifies
+what this came from.
+"""
+import numpy as np
+import underworld3 as uw
+import pyvista as pv
+
+OUT = "articles/testing-a-solver-against-exact-solutions/figures/banner.png"
+RES = 64
+
+pv.global_theme.allow_empty_mesh = True
+pv.global_theme.background = "white"
+
+panels = []
+for name, build in (
+        ("SolCx",  lambda m: uw.analytic.SolCx(m, eta_B=1.0e4)),
+        ("SolKz",  lambda m: uw.analytic.SolKz(m)),
+        ("SolM",   lambda m: uw.analytic.SolM(m)),
+):
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(RES, RES))
+    solution = build(mesh)
+    pvm = uw.visualisation.mesh_to_pv_mesh(mesh)
+    # Evaluate the solution's own symbolic viscosity straight onto the render
+    # points. No MeshVariable in between: nothing here is being solved, so
+    # there is nothing to project.
+    values = np.asarray(
+        uw.function.evaluate(solution.fn_viscosity, pvm.points[:, :mesh.dim]),
+        dtype=float).reshape(-1)
+    # log10, because these span four orders and a linear map shows one band.
+    pvm.point_data["log10_eta"] = np.log10(np.maximum(values, 1e-30))
+    print("%-6s log10 eta in [%.2f, %.2f]"
+          % (name, pvm.point_data["log10_eta"].min(),
+             pvm.point_data["log10_eta"].max()))
+    panels.append((name, pvm))
+
+pl = pv.Plotter(shape=(1, 3), window_size=(2100, 700), off_screen=True,
+                border=False)
+for col, (name, pvm) in enumerate(panels):
+    pl.subplot(0, col)
+    pl.set_background("white")
+    pl.add_mesh(pvm, scalars="log10_eta", cmap="RdBu_r", show_edges=False,
+                show_scalar_bar=False, lighting=False)
+    pl.enable_parallel_projection()
+    pl.view_xy()
+    pl.reset_camera(bounds=(0.0, 1.0, 0.0, 1.0, -0.05, 0.05))
+    pl.camera.zoom(1.30)
+pl.screenshot(OUT)
+print("wrote", OUT)