diff --git a/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh b/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh new file mode 100644 index 0000000..6ba684c --- /dev/null +++ b/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh @@ -0,0 +1,433 @@ +$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 98 1 98 +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 7 +5 +6 +7 +8 +9 +10 +11 +0.1249999999997738 0 0 +0.2499999999994124 0 0 +0.3749999999990479 0 0 +0.4999999999986942 0 0 +0.6249999999990126 0 0 +0.7499999999993417 0 0 +0.8749999999996708 0 0 +1 12 0 7 +12 +13 +14 +15 +16 +17 +18 +0.8749999999995011 1 0 +0.7500000000003465 1 0 +0.6250000000012137 1 0 +0.500000000002059 1 0 +0.3750000000015605 1 0 +0.2500000000010404 1 0 +0.1250000000005201 1 0 +1 13 0 7 +19 +20 +21 +22 +23 +24 +25 +1 0.1249999999997738 0 +1 0.2499999999994124 0 +1 0.3749999999990479 0 +1 0.4999999999986942 0 +1 0.6249999999990126 0 +1 0.7499999999993417 0 +1 0.8749999999996708 0 +1 14 0 7 +26 +27 +28 +29 +30 +31 +32 +0 0.8749999999995011 0 +0 0.7500000000003465 0 +0 0.6250000000012137 0 +0 0.500000000002059 0 +0 0.3750000000015605 0 +0 0.2500000000010404 0 +0 0.1250000000005201 0 +2 1 0 66 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +0.4375000000018098 0.8917468245270954 0 +0.1212591409528901 0.4311211703685352 0 +0.5624999999988534 0.1082531754730395 0 +0.8892750608248819 0.5640948631063587 0 +0.1082531754723038 0.6875000000007802 0 +0.690108146438289 0.8929849573910127 0 +0.312810865984173 0.1044931860863129 0 +0.8817373148985181 0.3307740260777834 0 +0.890514807813104 0.7966579344024687 0 +0.1880296711609206 0.8941534943677034 0 +0.812499999999506 0.1082531754733398 0 +0.1207695367638072 0.1869093576240868 0 +0.3125882785282735 0.892147936167085 0 +0.3777617321236122 0.7851464931730501 0 +0.5004602886882805 0.7837691230743562 0 +0.4407840225037774 0.6771478521171512 0 +0.564960049482756 0.6767503730786034 0 +0.5004884414966647 0.567277882976253 0 +0.3869428580759623 0.5672991787102901 0 +0.4430103817492976 0.4568493155218519 0 +0.5634998038750219 0.4584684189290819 0 +0.5010850309385081 0.3501225286940344 0 +0.3795377339355259 0.3481218582070821 0 +0.6251266299094012 0.3504430388690966 0 +0.6876877389643827 0.4586835206367175 0 +0.4375000000020022 0.2422277716919826 0 +0.75000000000057 0.3504809471648015 0 +0.6866823184207802 0.2476289266480872 0 +0.3141186395982236 0.2496554432804577 0 +0.2441861548444152 0.3577361238201676 0 +0.8124289467221131 0.4550024831196128 0 +0.7500000000003217 0.5669872981100748 0 +0.8221267806304088 0.6767859519483694 0 +0.6985159857017829 0.6821170180061178 0 +0.7624318058360666 0.7868767533111305 0 +0.2353327497992248 0.7709339938935139 0 +0.1082531754734316 0.312500000001173 0 +0.2225773588041955 0.4734472164134149 0 +0.1371230538040118 0.5617982471836261 0 +0.2180925345221853 0.6548089282261794 0 +0.6284495566844617 0.7850749770380402 0 +0.6874999999992675 0.1082531754732872 0 +0.4374999999991104 0.1082531754729983 0 +0.5632759514069947 0.8921152691542631 0 +0.6275253365868215 0.5683807519561415 0 +0.5625000000014364 0.2422277716920833 0 +0.1040995142213909 0.8120678220969118 0 +0.191853721181113 0.09292564867358062 0 +0.8965842169429771 0.1933620068391708 0 +0.8088666438741593 0.8971730864423774 0 +0.2076726160145944 0.2562891844220704 0 +0.7989110605604646 0.2500733495358726 0 +0.3260733379406406 0.6776011426060606 0 +0.9154642987621083 0.4374999999988711 0 +0.3281595767787979 0.4577243358118324 0 +0.9226624454758982 0.6874999999991773 0 +0.09150635094649401 0.9084936490535113 0 +0.9084936490536325 0.9084936490536784 0 +0.09150635094641726 0.09150635094661723 0 +0.9084936490535547 0.09150635094629651 0 +0.2698281199876322 0.5654465081585671 0 +0.6249999999989277 0.1752404735813565 0 +0.4999999999993456 0.1752404735815081 0 +0.3750000000023208 0.1752404735838395 0 +0.7398391675636558 0.1753067707867434 0 +0.2510355733930947 0.1746449973844669 0 +$EndNodes +$Elements +5 194 1 194 +1 11 1 8 +1 1 5 +2 5 6 +3 6 7 +4 7 8 +5 8 9 +6 9 10 +7 10 11 +8 11 2 +1 12 1 8 +9 4 12 +10 12 13 +11 13 14 +12 14 15 +13 15 16 +14 16 17 +15 17 18 +16 18 3 +1 13 1 8 +17 2 19 +18 19 20 +19 20 21 +20 21 22 +21 22 23 +22 23 24 +23 24 25 +24 25 4 +1 14 1 8 +25 3 26 +26 26 27 +27 27 28 +28 28 29 +29 29 30 +30 30 31 +31 31 32 +32 32 1 +2 1 2 162 +33 37 68 79 +34 68 37 72 +35 62 34 69 +36 59 40 63 +37 34 62 70 +38 63 40 86 +39 28 29 71 +40 65 41 67 +41 67 38 73 +42 41 65 88 +43 29 34 71 +44 37 28 71 +45 64 36 65 +46 63 36 64 +47 66 49 77 +48 66 67 73 +49 68 42 79 +50 34 70 71 +51 64 66 77 +52 49 66 73 +53 71 70 93 +54 72 71 93 +55 62 69 83 +56 40 59 84 +57 64 65 66 +58 43 81 84 +59 43 84 97 +60 68 72 85 +61 37 71 72 +62 66 65 67 +63 38 67 82 +64 33 45 46 +65 33 16 45 +66 34 30 69 +67 15 16 33 +68 45 42 68 +69 29 30 34 +70 73 38 76 +71 38 14 76 +72 17 42 45 +73 15 33 76 +74 35 9 74 +75 22 23 36 +76 31 44 69 +77 8 9 35 +78 14 15 76 +79 16 17 45 +80 39 7 75 +81 33 46 47 +82 30 31 69 +83 8 35 75 +84 33 47 76 +85 10 43 74 +86 7 8 75 +87 47 73 76 +88 17 18 42 +89 46 45 68 +90 9 10 74 +91 24 25 41 +92 31 32 44 +93 10 11 43 +94 6 7 39 +95 47 46 48 +96 57 64 77 +97 61 55 62 +98 57 63 64 +99 20 21 40 +100 27 28 37 +101 50 48 51 +102 54 52 55 +103 47 49 73 +104 57 59 63 +105 53 56 57 +106 50 51 52 +107 47 48 49 +108 54 55 58 +109 13 14 38 +110 57 56 59 +111 49 48 50 +112 53 52 54 +113 50 52 53 +114 53 54 56 +115 58 55 61 +116 53 57 77 +117 49 50 77 +118 50 53 77 +119 59 56 60 +120 54 58 78 +121 56 54 78 +122 60 56 78 +123 70 62 87 +124 44 80 98 +125 83 44 98 +126 26 27 79 +127 5 6 80 +128 19 20 81 +129 36 63 86 +130 12 13 82 +131 13 38 82 +132 65 36 88 +133 27 37 79 +134 6 39 80 +135 20 40 81 +136 87 51 93 +137 67 41 82 +138 51 48 85 +139 46 68 85 +140 22 36 86 +141 40 21 86 +142 62 55 87 +143 52 51 87 +144 36 23 88 +145 24 41 88 +146 3 26 89 +147 4 12 90 +148 1 5 91 +149 32 1 91 +150 2 19 92 +151 11 2 92 +152 18 3 89 +153 25 4 90 +154 51 85 93 +155 61 62 83 +156 69 44 83 +157 59 60 84 +158 70 87 93 +159 48 46 85 +160 81 40 84 +161 21 22 86 +162 55 52 87 +163 85 72 93 +164 19 81 92 +165 26 79 89 +166 5 80 91 +167 23 24 88 +168 42 18 89 +169 41 25 90 +170 44 32 91 +171 43 11 92 +172 12 82 90 +173 60 94 97 +174 39 75 96 +175 94 78 95 +176 35 74 94 +177 35 94 95 +178 78 58 95 +179 75 35 95 +180 60 78 94 +181 58 61 96 +182 95 58 96 +183 75 95 96 +184 94 74 97 +185 74 43 97 +186 79 42 89 +187 81 43 92 +188 80 44 91 +189 96 61 98 +190 39 96 98 +191 82 41 90 +192 80 39 98 +193 84 60 97 +194 61 83 98 +$EndElements diff --git a/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh.h5 b/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh.h5 new file mode 100644 index 0000000..2d4aeee Binary files /dev/null and b/articles/setting-up-full-multigrid/examples/.meshes/uw_simplexbox_minC(0.0, 0.0)_maxC(1.0, 1.0)_csize0.125_regFalse.msh.h5 differ diff --git a/articles/setting-up-full-multigrid/examples/banner.py b/articles/setting-up-full-multigrid/examples/banner.py new file mode 100644 index 0000000..b4defc7 --- /dev/null +++ b/articles/setting-up-full-multigrid/examples/banner.py @@ -0,0 +1,61 @@ +"""Banner: the three levels a `refinement=2` mesh carries. + +Run from the repository root: + + python3 articles/setting-up-full-multigrid/examples/banner.py + +The subject of the note is that a mesh built with `refinement` keeps its coarser +ancestors, and that those are what the preconditioner works on. So the banner is +the hierarchy itself: the gmsh base mesh, and the two refinements built from it, +each one subdividing every cell of the one before and snapping the new boundary +nodes back onto the bounding circles. + +Building three meshes at refinement 0, 1 and 2 gives exactly the three levels a +single `refinement=2` mesh holds in `dm_hierarchy` -- same base mesh, same +refinement callback -- and does it through the supported API rather than +reaching into the DMPlex objects. + +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that +identifies what these numbers came from. +""" +import underworld3 as uw +import pyvista as pv + +RADIUS_INNER, RADIUS_OUTER = 0.5, 1.0 +CELL_SIZE = 0.25 # the BASE mesh; refinement takes it from there +OUT = "articles/setting-up-full-multigrid/figures/banner.png" + +pv.global_theme.allow_empty_mesh = True +pv.global_theme.background = "white" + +levels = [] +for refinement in (0, 1, 2): + mesh = uw.meshing.Annulus(radiusInner=RADIUS_INNER, radiusOuter=RADIUS_OUTER, + cellSize=CELL_SIZE, qdegree=2, + refinement=refinement) + pvm = uw.visualisation.mesh_to_pv_mesh(mesh) + print("refinement %d: %6d cells, %d hierarchy level(s)" + % (refinement, pvm.n_cells, len(mesh.dm_hierarchy))) + levels.append(pvm) + +# Wide and short: a banner is cropped hard on a narrow screen, so the three +# panels have to read at a glance and nothing may depend on fine detail. +pl = pv.Plotter(shape=(1, 3), window_size=(2100, 700), off_screen=True, + border=False) +for col, pvm in enumerate(levels): + pl.subplot(0, col) + pl.set_background("white") + # Line width drops as the cells get smaller, so the finest panel reads as a + # texture rather than as a block of ink. + pl.add_mesh(pvm, color="white", show_edges=True, edge_color="#1a1a1a", + line_width=(1.6, 1.1, 0.7)[col], lighting=False) + pl.enable_parallel_projection() + pl.view_xy() + pl.reset_camera(bounds=(-1.02, 1.02, -1.02, 1.02, -0.05, 0.05)) + # reset_camera leaves a wide margin; a banner is short and gets cropped, so + # the meshes have to carry the strip rather than float in it. + pl.camera.zoom(1.35) +pl.screenshot(OUT) +print("wrote", OUT) diff --git a/articles/setting-up-full-multigrid/examples/fmg-vs-gamg.py b/articles/setting-up-full-multigrid/examples/fmg-vs-gamg.py new file mode 100644 index 0000000..4d7b925 --- /dev/null +++ b/articles/setting-up-full-multigrid/examples/fmg-vs-gamg.py @@ -0,0 +1,323 @@ +"""FMG against GAMG on the SolKz problem: timing only. + +SolKz is the standard smooth-viscosity Stokes benchmark: on the unit box, + + eta(z) = exp(2 B z) viscosity varying exponentially with depth + f = (0, sin(m pi z) cos(n pi x)) + free slip on all four walls + +and the viscosity contrast across the box is exp(2B). It has a closed-form +solution, but nothing here uses it: every number below is a TIME, and the point +is how the cost of the solve behaves, not how accurate the answer is. + +A second viscosity profile is measured alongside it. SolKz spreads its +viscosity variation over the whole box, which is the structure algebraic +coarsening reads best, so on its own it cannot say whether the two methods +separate when the structure is CONCENTRATED. The `layer` profile puts the same +total contrast into a band a twentieth of the box thick, + + eta(z) = 1 + (contrast - 1) exp(-((z - 1/2) / w)^2), w = 0.05 + +with the forcing and the boundary conditions unchanged, so the viscosity is the +only thing that differs. It is smooth, so nothing here is testing how either +method copes with an under-resolved jump. + +Two things are deliberately fixed across every run: + + * the solver tolerance. Holding it constant while the mesh refines is the + usual way to show a multigrid scaling, and it is a choice rather than a + neutral one -- a finer mesh has a smaller discretisation error, so an + argument could be made for tightening the solver tolerance alongside it, + which would make the work per unknown grow. Fixed tolerance answers "what + does it cost to solve this system", not "what does it cost to reach the + accuracy the mesh can support". + + * the velocity iteration cap, raised well above the default 200 so that + neither preconditioner is truncated rather than allowed to converge. + +The velocity block is solved with fgmres preconditioned by multigrid, so one +Krylov iteration is one multigrid cycle. Counting them needs care: the Schur +factorisation invokes the velocity solve SIXTEEN times per Stokes solve, so a +naive total folds in a factor that is identical for both preconditioners. What +is reported here is therefore cycles PER VELOCITY SOLVE, which is the +like-for-like number. Even so it counts cycles, not work -- a geometric cycle +and an algebraic one are different amounts of it, which is what the timings +are for. + +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that +identifies what these numbers came from. +""" +import math +import os +import subprocess +import sys +import time + +from petsc4py import PETSc + +PETSc.Log.begin() # before any solving, or the counts start from nowhere + +import sympy +import underworld3 as uw + +BASE_CELL = 1 / 8 +QDEG = 3 +TOL = 1.0e-6 +VEL_CAP = 20000 # high enough that NOTHING is truncated at any size + # measured here. At 2000, GAMG at the largest problem is + # cut off mid-solve, the outer solve compensates with twice + # as many velocity solves, and the reported work jumps 35% + # for a reason that has nothing to do with the method. +# Each configuration is measured in a FRESH PROCESS. +# +# Measured alone, every configuration here is exactly reproducible: the same +# flop count to five figures, and sixteen velocity solves per Stokes solve. +# Run back to back in one process they are not -- one sweep reported a doubled +# invocation count and a 35% larger GAMG figure, which no configuration +# reproduces on its own. Something survives between solves. Rather than trust a +# number that depends on what ran before it, the parent process below invokes +# this file once per configuration and parses one line back. +N_X, M_Z = 3, 2 # the SolKz wavenumbers +LAYER_W = 0.05 # half-width of the concentrated viscosity band + + +def viscosity(profile, contrast, z): + """The viscosity field, as a sympy expression in the depth coordinate.""" + if profile == "solkz": + B = 0.5 * math.log(contrast) if contrast > 1.0 else 0.0 + return sympy.exp(2 * B * z) + if profile == "layer": + return 1 + (contrast - 1) * sympy.exp(-(((z - 0.5) / LAYER_W) ** 2)) + raise ValueError("unknown viscosity profile %r" % profile) + + +def solve_once(pc, refinement, contrast, profile="solkz"): + """Time one solve. Returns wall seconds, flops, cycles and unknowns.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), + cellSize=BASE_CELL, qdegree=QDEG, refinement=refinement) + x, z = mesh.X.coords_sym if hasattr(mesh.X, "coords_sym") else mesh.X + + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = viscosity( + profile, contrast, z) + stokes.bodyforce = sympy.Matrix( + [0, sympy.sin(M_Z * sympy.pi * z) * sympy.cos(N_X * sympy.pi * x)]) + + # Free slip: the wall-normal component is held, the tangential one is free. + for wall, comp in (("Left", (0,)), ("Right", (0,)), + ("Bottom", (1,)), ("Top", (1,))): + stokes.add_dirichlet_bc((0.0,) * mesh.dim, wall, comp) + + stokes.preconditioner = pc + stokes.tolerance = TOL + stokes.petsc_options.setValue("fieldsplit_velocity_ksp_max_it", VEL_CAP) + + # The OUTER Krylov, set explicitly for both preconditioners so this script + # measures the same thing whatever Underworld's default happens to be. + # + # It has to be flexible. Both sub-blocks are Krylov solves run to a + # tolerance, so inside the Schur factorisation the operator the outer method + # applies differs from one iteration to the next, and plain gmres assumes it + # does not. On a smooth problem the outer converges in about two iterations + # and the choice is invisible; on a concentrated viscosity structure the + # inner solves are genuinely inexact and it is not. Leaving this to the + # default would measure that choice rather than the preconditioner. + stokes.petsc_options.setValue("ksp_type", "fgmres") + if pc == "gamg": + # PETSc's own default cycle for GAMG, set explicitly so that this + # script measures the same thing whatever Underworld ships. + stokes.petsc_options.setValue( + "fieldsplit_velocity_pc_mg_type", "multiplicative") + + # Build everything; not timed, and deliberately not converged. This solve + # exists to make PETSc assemble the operators and set up the fieldsplit, so + # that the block size can be corrected below -- the ANSWER is thrown away. + # + # It has to be capped. For GAMG this runs at the operator's default block + # size of 1, which is the configuration the correction below exists to undo, + # and on a concentrated viscosity structure at high contrast that costs + # HOURS while the measured solve costs seconds. Capped at one iteration it + # costs nothing and builds exactly the same objects. + stokes.petsc_options.setValue("fieldsplit_velocity_ksp_max_it", 1) + stokes.petsc_options.setValue("ksp_max_it", 1) + stokes.petsc_options.setValue("snes_max_it", 1) + try: + stokes.solve() + except (RuntimeError, PETSc.Error): + pass # not converging is the point + stokes.petsc_options.setValue("fieldsplit_velocity_ksp_max_it", VEL_CAP) + stokes.petsc_options.delValue("ksp_max_it") + stokes.petsc_options.delValue("snes_max_it") + + if pc == "gamg": + # Tell GAMG the field has two components per node, so it aggregates + # nodes rather than scalars. The operator is built by PETSc from the + # DM, so this has to be set after it exists and the PC rebuilt. + sub_v = stokes.snes.getKSP().getPC().getFieldSplitSubKSP()[0] + A, P = sub_v.getOperators() + A.setBlockSize(mesh.dim) + if P.handle != A.handle: + P.setBlockSize(mesh.dim) + sub_v.getPC().reset() + + # Count iterations across the WHOLE solve, not the last linear solve. A + # linear Stokes SNES does its work in the first solve and the second + # converges immediately, so getIterationNumber() afterwards reports the + # trivial one and makes every resolution look identical. + ksp = stokes.snes.getKSP() + sub = ksp.getPC().getFieldSplitSubKSP()[0] + # Record the ITERATION INDEX, not just a count: every time it restarts at + # zero the Schur factorisation has invoked the velocity solve again, which + # is what separates "cycles per solve" from "cycles in total". + seen = [] + sub.setMonitor(lambda k, i, rn: seen.append(i)) + + v.array[...] = 0.0 + p.array[...] = 0.0 + + # Flops as well as time. The flop count is deterministic and independent of + # the machine, so it measures the METHOD; the time measures this machine. + # Their ratio is the rate, which is where the two disagree. + f0 = PETSc.Log.getFlops() + t0 = time.time() + stokes.solve(zero_init_guess=True) + wall = time.time() - t0 + flops = PETSc.Log.getFlops() - f0 + + ok = stokes.snes.getConvergedReason() > 0 and sub.getConvergedReason() > 0 + starts = [j for j, i in enumerate(seen) if i == 0] + per_solve = [len(seen[a:b]) for a, b in zip(starts, starts[1:] + [len(seen)])] + return dict(wall=wall, flops=flops, ok=ok, + ndof=stokes.snes.getSolution().getSize(), + invocations=len(starts), velocity=len(seen), + cycles=(min(per_solve), max(per_solve)) if per_solve else (0, 0)) + + +def _measure_and_print(pc, refinement, contrast, profile): + """Child-process entry point: measure one configuration, print one line.""" + r = solve_once(pc, refinement, contrast, profile) + print("RESULT %s %d %g %s %d %r %.6e %.6f %d %d %d" + % (pc, refinement, contrast, profile, r["ndof"], r["ok"], r["flops"], + r["wall"], r["invocations"], r["cycles"][0], r["cycles"][1])) + + +def measure(pc, refinement, contrast, profile="solkz"): + """Run one configuration in a fresh interpreter and read the result back.""" + out = subprocess.run( + [sys.executable, os.path.abspath(__file__), + "--one", pc, str(refinement), repr(contrast), profile], + capture_output=True, text=True) + for line in out.stdout.splitlines(): + if line.startswith("RESULT "): + f = line.split() + # Echo each configuration as it lands. The tables are assembled at + # the end, and a sweep that only prints then loses everything if it + # is interrupted -- which matters when one configuration can run for + # an hour. + print(" " + line, flush=True) + return dict(pc=f[1], ndof=int(f[5]), ok=f[6] == "True", + flops=float(f[7]), wall=float(f[8]), + invocations=int(f[9]), + cycles=(int(f[10]), int(f[11]))) + sys.stderr.write(out.stdout + out.stderr) + raise RuntimeError("no result from %s at refinement %d, contrast %g (%s)" + % (pc, refinement, contrast, profile)) + + +def contrast_table(profile="solkz"): + """Cost against viscosity contrast, at one resolution.""" + rows = [(c, pc, measure(pc, 2, c, profile)) + for c in (1.0, 1.0e2, 1.0e4, 1.0e6) for pc in ("fmg", "gamg")] + # Work per unknown, relative to FMG on the easiest problem. A ratio does + # not depend on the machine it was measured on. + bw = next(r["flops"] / r["ndof"] for c, pc, r in rows if c == 1.0 and pc == "fmg") + bt = next(r["wall"] / r["ndof"] for c, pc, r in rows if c == 1.0 and pc == "fmg") + print("\n%s: cost against viscosity contrast" % profile) + print("\n| preconditioner | viscosity contrast | relative work per unknown " + "| relative time per unknown | cycles per velocity solve |") + print("|---|---|---|---|---|") + for c, pc, r in sorted(rows, key=lambda k: (k[1] != "fmg", k[0])): + lo, hi = r["cycles"] + cycles = str(lo) if lo == hi else "%d–%d" % (lo, hi) + if r["ok"]: + print("| %s | 10^%d | %.1f | %.1f | %s |" + % (pc.upper(), round(math.log10(c)), + (r["flops"] / r["ndof"]) / bw, (r["wall"] / r["ndof"]) / bt, + cycles)) + else: + # Not converged. Say so rather than reporting the cost of a solve + # that did not finish -- and report the cycle count anyway, because + # a count sitting exactly on VEL_CAP means truncation, not stall. + print("| %s | 10^%d | did not converge | did not converge | %s |" + % (pc.upper(), round(math.log10(c)), cycles)) + + +def scaling_table(): + """Cost against problem size, at constant viscosity.""" + # BOTH run to refinement 5. GAMG at that size is expensive, but stopping it + # a level short leaves the largest FMG row with nothing to compare against, + # and the wall-clock comparison is the one that turns on the largest size. + out = {pc: [measure(pc, r, 1.0) for r in refs] + for pc, refs in (("fmg", (1, 2, 3, 4, 5)), ("gamg", (1, 2, 3, 4, 5)))} + + bw = out["fmg"][0]["flops"] / out["fmg"][0]["ndof"] + bt = out["fmg"][0]["wall"] / out["fmg"][0]["ndof"] + print("\n| preconditioner | unknowns | relative work per unknown " + "| relative time per unknown | Gflop/s |") + print("|---|---|---|---|---|") + for pc in ("fmg", "gamg"): + for r in out[pc]: + print("| %s | %d | %.1f | %.1f | %.1f |" + % (pc.upper(), r["ndof"], (r["flops"] / r["ndof"]) / bw, + (r["wall"] / r["ndof"]) / bt, r["flops"] / r["wall"] / 1e9)) + + # Cycles at EVERY size, not the first few: the claim being made is that the + # count stops growing, and a table that stops short of the largest problem + # cannot support it. + print("\n| preconditioner | %s |" + % " | ".join("%d unknowns" % r["ndof"] for r in out["fmg"])) + print("|---|" + "---|" * len(out["fmg"])) + for pc in ("fmg", "gamg"): + cells = [] + for r in out[pc]: + lo, hi = r["cycles"] + cells.append(str(lo) if lo == hi else "%d\u2013%d" % (lo, hi)) + print("| %s | %s |" % (pc.upper(), " | ".join(cells))) + invocations = sorted({r["invocations"] for rows in out.values() for r in rows}) + print("\n(velocity solves per Stokes solve: %s)" % invocations) + + print() + for pc, rows in out.items(): + steps = [math.log(b["flops"] / a["flops"]) / math.log(b["ndof"] / a["ndof"]) + for a, b in zip(rows, rows[1:])] + xs = [math.log(r["ndof"]) for r in rows] + ys = [math.log(r["wall"]) for r in rows] + mx, my = sum(xs) / len(xs), sum(ys) / len(ys) + power = (sum((x - mx) * (y - my) for x, y in zip(xs, ys)) + / sum((x - mx) ** 2 for x in xs)) + print("%s: time ~ N^%.2f, flops ~ N^%s per step" + % (pc.upper(), power, [round(s, 2) for s in steps])) + + +if __name__ == "__main__": + if len(sys.argv) > 1 and sys.argv[1] == "--one": + _measure_and_print(sys.argv[2], int(sys.argv[3]), float(sys.argv[4]), + sys.argv[5]) + elif len(sys.argv) > 1: + # One table at a time, so a re-measurement does not mean re-running + # everything: `fmg-vs-gamg.py solkz`, `layer`, or `scaling`. + for name in sys.argv[1:]: + if name == "scaling": + scaling_table() + else: + contrast_table(name) + else: + contrast_table("solkz") + contrast_table("layer") + scaling_table() diff --git a/articles/setting-up-full-multigrid/examples/layer_figure.py b/articles/setting-up-full-multigrid/examples/layer_figure.py new file mode 100644 index 0000000..5180ca0 --- /dev/null +++ b/articles/setting-up-full-multigrid/examples/layer_figure.py @@ -0,0 +1,104 @@ +"""Figure: what happens to the comparison when the viscosity structure is +concentrated rather than spread through the box. + +Run from the repository root: + + python3 articles/setting-up-full-multigrid/examples/layer_figure.py + +The numbers are those printed by `fmg-vs-gamg.py solkz layer` and are repeated +here so the figure has a single, checkable source. Both panels are relative to +the same baseline -- FMG on the constant-viscosity problem -- so the two +profiles can be read against each other. At contrast 1 the two profiles ARE the +same problem (eta = 1 everywhere) and the runs agree to the last digit, which is +why the curves start together. + +Colour carries the preconditioner and line style carries the viscosity profile, +so neither is distinguished by colour alone. The palette is checked for +colour-vision separation (worst pair dE 24.7, against a target of 8). + +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that +identifies what these numbers came from. +""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +OUT = "articles/setting-up-full-multigrid/figures/layer-vs-smooth.png" + +CONTRAST = [1e0, 1e2, 1e4, 1e6] + +# Relative to FMG at contrast 1. `None` = did not converge: GAMG on the band at +# 1e6 ran past 20 000 multigrid cycles per velocity solve without reaching the +# tolerance, so there is no cost to plot -- reporting the work it managed to do +# before being stopped would be reporting a cap, not a method. +WORK = { + ("fmg", "smooth"): [1.0, 1.3, 1.6, 1.8], + ("gamg", "smooth"): [2.2, 3.3, 4.0, 4.0], + ("fmg", "band"): [1.0, 1.9, 2.1, 2.5], + ("gamg", "band"): [2.2, 8.1, 28.1, None], +} +TIME = { + ("fmg", "smooth"): [1.0, 1.1, 1.2, 1.3], + ("gamg", "smooth"): [1.2, 1.5, 1.7, 1.6], + ("fmg", "band"): [1.0, 1.3, 1.4, 1.5], + ("gamg", "band"): [1.2, 2.7, 7.7, None], +} + +COLOUR = {"fmg": "#2a78d6", "gamg": "#eb6834"} # validated categorical 1, 2 +STYLE = {"smooth": "-", "band": "--"} +INK = "#0b0b0b" +INK_MUTED = "#52514e" +GRID = "#e4e3df" + + +def panel(ax, data, title): + for (pc, profile), y in data.items(): + xs = [x for x, v in zip(CONTRAST, y) if v is not None] + ys = [v for v in y if v is not None] + ax.plot(xs, ys, STYLE[profile], color=COLOUR[pc], linewidth=2.0, + marker="o", markersize=5, markerfacecolor="white", + markeredgewidth=1.6, zorder=3, + label="%s, %s" % (pc.upper(), profile)) + # An open marker at the last point that converged, with the failure + # named rather than left as a gap the reader has to interpret. + if y[-1] is None: + # Say what the gap means. An unexplained stop reads as missing data. + ax.annotate("does not converge\nbeyond this point", + xy=(xs[-1], ys[-1]), xytext=(10, -2), + textcoords="offset points", ha="left", va="top", + fontsize=8.5, color=INK_MUTED, linespacing=1.35) + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlabel("viscosity contrast", fontsize=9.5, color=INK_MUTED) + ax.set_title(title, fontsize=10.5, color=INK, pad=10) + ax.grid(True, which="major", color=GRID, linewidth=0.8, zorder=0) + ax.set_axisbelow(True) + for side in ("top", "right"): + ax.spines[side].set_visible(False) + for side in ("left", "bottom"): + ax.spines[side].set_color(GRID) + ax.tick_params(colors=INK_MUTED, labelsize=9) + # Same tick VALUES in both panels so the two can be read against each + # other, but each panel framed on its own data rather than padded to a + # shared top -- the time differences are genuinely smaller and should look + # it without half the panel being empty. + ax.set_yticks([1, 1.5, 2, 3, 5, 10, 20, 30]) + ax.set_yticklabels(["1", "1.5", "2", "3", "5", "10", "20", "30"]) + finite = [v for y in data.values() for v in y if v is not None] + ax.set_ylim(min(finite) * 0.88, max(finite) * 1.45) + ax.minorticks_off() + + +fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.1), sharex=True) +fig.patch.set_facecolor("white") +panel(axes[0], WORK, "work per unknown, relative to FMG at contrast 1") +panel(axes[1], TIME, "time per unknown, relative to FMG at contrast 1") +# handlelength long enough that the dashed entries are visibly dashed -- the +# line style is half the encoding, so a legend that hides it breaks identity. +axes[0].legend(frameon=False, fontsize=9, labelcolor=INK_MUTED, + loc="upper left", handlelength=3.4, borderaxespad=0.2) +fig.tight_layout() +fig.savefig(OUT, dpi=200, facecolor="white") +print("wrote", OUT) diff --git a/articles/setting-up-full-multigrid/examples/relax_annulus_figure.py b/articles/setting-up-full-multigrid/examples/relax_annulus_figure.py new file mode 100644 index 0000000..4f07e13 --- /dev/null +++ b/articles/setting-up-full-multigrid/examples/relax_annulus_figure.py @@ -0,0 +1,72 @@ +"""Figure: what the coarse mesh leaves in the fine one, and what relax() does. + +Run from the repository root: + + python3 articles/setting-up-full-multigrid/examples/relax_annulus_figure.py + +The annulus is the clearer subject than a box: its base mesh is coarse relative +to the curvature, so the coarse cells' edges survive as visible seams across the +refined mesh. It also shows that the boundary snapping holds -- every node that +starts on a bounding circle is still exactly on it afterwards, because relax() +moves interior coordinates only. + +Run against underworld3 `development` at commit `0addec15` +(0addec1595f8d7a59b99e15b42455267a73dab86, 2026-08-15). `uw.__version__` +reports 0.0.0 for every build, so the commit is the only thing that +identifies what these numbers came from. +""" +import numpy as np +import underworld3 as uw +import pyvista as pv + +RADIUS_INNER, RADIUS_OUTER = 0.5, 1.0 +CELL_SIZE = 0.25 # the BASE mesh; refinement takes it from there +REFINEMENT = 2 +OUT = "articles/setting-up-full-multigrid/figures/relax-annulus.png" + +pv.global_theme.allow_empty_mesh = True +pv.global_theme.background = "white" + +mesh = uw.meshing.Annulus(radiusInner=RADIUS_INNER, radiusOuter=RADIUS_OUTER, + cellSize=CELL_SIZE, qdegree=2, refinement=REFINEMENT) +before = uw.visualisation.mesh_to_pv_mesh(mesh) +print("cells %d, hierarchy levels %d" % (before.n_cells, len(mesh.dm_hierarchy))) + + +def report(pvm, tag): + """Element quality, and whether the bounding circles are still exact.""" + cells = pvm.cells.reshape(-1, 4)[:, 1:] + p = pvm.points[:, :2] + a = np.linalg.norm(p[cells[:, 1]] - p[cells[:, 2]], axis=1) + b = np.linalg.norm(p[cells[:, 2]] - p[cells[:, 0]], axis=1) + c = np.linalg.norm(p[cells[:, 0]] - p[cells[:, 1]], axis=1) + s = 0.5 * (a + b + c) + area = np.sqrt(np.maximum(s * (s - a) * (s - b) * (s - c), 0.0)) + q = 2.0 * (area / np.maximum(s, 1e-30)) / np.maximum( + a * b * c / np.maximum(4.0 * area, 1e-30), 1e-30) + r = np.linalg.norm(p, axis=1) + near = (np.abs(r - RADIUS_OUTER) < 2e-3) | (np.abs(r - RADIUS_INNER) < 2e-3) + exact = (np.abs(r - RADIUS_OUTER) < 1e-6) | (np.abs(r - RADIUS_INNER) < 1e-6) + print("%-7s q median %.3f 10th pct %.3f min %.3f | on the circle: %d/%d" + % (tag, np.median(q), np.percentile(q, 10), q.min(), + exact.sum(), near.sum())) + + +report(before, "before") +mesh.relax() +after = uw.visualisation.mesh_to_pv_mesh(mesh) +report(after, "after") + +pl = pv.Plotter(shape=(1, 2), window_size=(1900, 950), off_screen=True, border=False) +for col, (m, title) in enumerate(((before, "Annulus, refinement=2"), + (after, "after relax()"))): + pl.subplot(0, col) + pl.set_background("white") + pl.add_mesh(m, color="white", show_edges=True, edge_color="#1a1a1a", + line_width=1.1, lighting=False) + pl.add_text(title, position="upper_edge", font_size=14, color="black") + pl.enable_parallel_projection() + pl.view_xy() + pl.reset_camera(bounds=(-1.05, 1.05, -1.05, 1.12, -0.05, 0.05)) +pl.screenshot(OUT) +print("wrote", OUT) diff --git a/articles/setting-up-full-multigrid/figures/banner.png b/articles/setting-up-full-multigrid/figures/banner.png new file mode 100644 index 0000000..7905614 Binary files /dev/null and b/articles/setting-up-full-multigrid/figures/banner.png differ diff --git a/articles/setting-up-full-multigrid/figures/layer-vs-smooth.png b/articles/setting-up-full-multigrid/figures/layer-vs-smooth.png new file mode 100644 index 0000000..b8417af Binary files /dev/null and b/articles/setting-up-full-multigrid/figures/layer-vs-smooth.png differ diff --git a/articles/setting-up-full-multigrid/figures/relax-annulus.png b/articles/setting-up-full-multigrid/figures/relax-annulus.png new file mode 100644 index 0000000..06b9d58 Binary files /dev/null and b/articles/setting-up-full-multigrid/figures/relax-annulus.png differ diff --git a/articles/setting-up-full-multigrid/metadata.yml b/articles/setting-up-full-multigrid/metadata.yml new file mode 100644 index 0000000..bc250dd --- /dev/null +++ b/articles/setting-up-full-multigrid/metadata.yml @@ -0,0 +1,38 @@ +# Validated in CI against schemas/article-metadata.schema.json. +# `pixi run validate` checks this and the cross-file invariants a schema cannot +# express -- that the article file is named .md, that canonical_path +# matches the slug, and that no legacy DOI is ever paired with a new registrant. +id: UWTN 2026-014 +slug: setting-up-full-multigrid +title: Setting Up Full Multigrid +article_type: technical-note +status: published +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliation: Australian National University +publication_date: 2026-08-17 +version: 1.0.0 +# The deposit writes archive_doi and repository_record_id when the note is +# published; leave them out until then. `doi` and `doi_registrant` were here +# once and are not fields the schema knows -- every note made from this +# template failed `pixi run validate` on all three of them. +license: CC-BY-4.0 +canonical_path: /setting-up-full-multigrid/ +legacy_paths: [] +# Facets, from vocabulary.yml. Both keys must be present even when empty: a +# note with no subject is normal -- many are purely about method. +subjects: + - mantle-convection +methods: + - solvers + - meshing + - parallel-hpc +ghost_tags: + - Underworld Code +figures: 2 +# Generated from the model rather than a stock photograph, so there is nobody to +# credit. `figures` counts the figures in the body; the banner is not one. +banner: figures/banner.png +banner_credit: null +source: native diff --git a/articles/setting-up-full-multigrid/setting-up-full-multigrid.md b/articles/setting-up-full-multigrid/setting-up-full-multigrid.md new file mode 100644 index 0000000..fb87464 --- /dev/null +++ b/articles/setting-up-full-multigrid/setting-up-full-multigrid.md @@ -0,0 +1,301 @@ +--- +title: Setting Up Full Multigrid +description: >- + Underworld can precondition a Stokes solve with geometric multigrid built + from a real hierarchy of meshes, rather than with the algebraic multigrid it + falls back to. What the difference is, how to build a mesh that has a + hierarchy, and when the choice is worth making. +date: 2026-08-14 +authors: + - name: Louis Moresi + orcid: 0000-0003-3685-174X + affiliations: + - Australian National University +license: CC-BY-4.0 +banner: figures/banner.png +keywords: + - Underworld Code + - Tricks of the Trade + - development +exports: + - format: typst + logo: ../../static/uwtn-logo.png + series: "Underworld Technical Notes" + origin_url: https://www.underworldcode.org/setting-up-full-multigrid/ + template: ../../templates/pdf + output: setting-up-full-multigrid.pdf + article_id: UWTN 2026-014 + article_version: 1.0.0 + software_version: underworld3 development @ 0addec15 +--- +
+ +Most of the effort in a geodynamics model goes into solving the Stokes +equations, and most of that is spent on the velocity solver. How that solver is +preconditioned sets how long a model takes to run. Multigrid methods accelerate elliptic +solvers using a hierarchy of mesh resolutions. Underworld can build the +preconditioner from a real hierarchy of meshes or directly from the matrix +alone, and this note is about how to build that hierarchy and what it is worth. + +## What multigrid does + +An iterative solver reduces error unevenly. Simple relaxation methods +(smoothers) are good at removing error that varies rapidly from cell to cell, +and poor at removing error that varies smoothly across the whole domain. The +smooth part is then what costs: it takes many sweeps to carry information from +one side of a fine mesh to the other. + +Multigrid turns that around. Error that is smooth on a fine mesh is *not* smooth +on a mesh twice as coarse, where it spans half as many cells, so a few sweeps +there can quickly reduce error that the fine grid struggles with. + +The method is a recursion on that idea: smooth on the fine mesh, transfer what is left to a +coarser one, solve there with another recursive step, transfer the correction +back, smooth again. Each *level* handles the band of error it can see, +and the work per unknown flattens out as the model gets bigger. + +The two things we need: the coarse levels, and the operators that move +between them. + +## Two ways to build the coarse levels + +**Algebraic multigrid** (AMG) builds the coarse levels from the matrix +representation of the problem directly. It inspects the operator's connectivity, +groups unknowns that are strongly +coupled, and constructs coarse representations and transfer operators from that grouping +alone. Because it asks nothing of the mesh, it works everywhere. +PETSc's implementation is GAMG, and it is what Underworld +uses by default. + +**Geometric multigrid** (GMG) uses actual coarser meshes. If the fine mesh was built +by ***refining*** a coarse one, those coarser meshes already exist, and the transfers +follow from the refinement relation: each fine node either coincides with a +coarse node or lies inside a known coarse cell. In PETSc this is `pc_type=mg`, +and Underworld runs it as a *full* multigrid cycle: solve on the coarsest mesh, +interpolate that solution up to the next mesh as a starting point, and repeat, +cycling up and down through the levels at each stage to clear the error the +interpolation leaves behind. We'll refer to this as FMG. + +The difference between them shows up when the operator's connectivity is a poor +guide to the geometry, which can happen when there are jumps or steep gradients in material +properties. The algebraic method sees a matrix whose couplings are +dominated by the stiff region and groups unknowns accordingly; the geometric +method uses the grids it was given and is indifferent to the coefficients. Which +of those is the better bet is not obvious in advance, so we'll measure it and see what works. +Switching the geometric solver on comes first, and we start with +the hierarchy of meshes. + +## Making a mesh with a hierarchy + +A mesh only carries a hierarchy if it was built with one. That is the +`refinement` argument: + +```python +mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.05, refinement=2) +len(mesh.dm_hierarchy) # 3: the base and two refinements +``` + +The mesh you get back is the finest level, and the coarser ones are stored +within it as PETSc `DMPlex` objects in `mesh.dm_hierarchy`. Those are what the +preconditioner uses. They are not Underworld meshes, there is no `Mesh` object +for a coarse level, and making one takes work. The hierarchy is something +the solver knows about rather than something you handle. + +Build the same mesh at `cellSize=0.0125` with no refinement and you get about +the same resolution with no hierarchy at all, and so no geometric multigrid. +Algebraic multigrid runs perfectly well on a mesh that has a hierarchy; it just +does not use it. + +The solver picks it up on its own: + +```python +stokes = uw.systems.Stokes(mesh) +stokes.solve() # preconditioner defaults to "auto" + +stokes.preconditioner = "fmg" # or "gamg", or back to "auto" +stokes.solve() +``` + +`"auto"` uses geometric multigrid when `len(mesh.dm_hierarchy) > 1` and GAMG +when it does not. You can also ask directly, with `"fmg"` or `"gamg"`. Two +behaviours are deliberate: `"auto"` only ever adds geometric multigrid to an +untouched configuration, so it will not overwrite preconditioner options you +have set yourself; and where a request cannot be honoured, you can find the +reasoning in `stokes.pc_fallbacks`. `stokes.preconditioner_settings` reports +what was actually applied. + +## Curved boundaries need help + +Refining a mesh puts new nodes at the midpoints of existing edges. On a +straight boundary that works well. On a curved one it is not enough: the coarse mesh +approximates a circle by a polygon, and the midpoint of a chord lies inside the +circle rather than on it. Those midpoints have to be moved back onto the +boundary each time the mesh is refined. + +Underworld's curved meshes carry a refinement *callback* that snaps +boundary-labelled nodes back onto the true surface after each refinement. For +the annulus it is a few lines: take the nodes labelled `Upper` and `Lower` and +rescale each to the correct radius: + +```python +coords[upper] *= radiusOuter / R[upper] +coords[lower] *= radiusInner / R[lower] +``` + +This happens automatically for the built-in meshes, so in normal use there is +nothing you need to do. It does matter if you build a mesh of your own with curved +boundaries: without a callback of this kind, the finest mesh still has the +coarsest mesh's faceting, and the boundary it represents is the polygon rather +than the circle. + +## What the coarse mesh leaves behind + +Because the fine mesh is made by subdividing the coarse one, the coarse mesh's +own triangulation is still visible in it. Its edges survive as continuous lines +across the fine mesh, and its vertices remain places where an unusual number of +elements meet and the element density is inherited from the coarse triangulation. + +`mesh.relax()` moves the nodes to improve element shapes while keeping the +resolution and the topology, and it helps to loosen that imprint. + +```{figure} figures/relax-annulus.png +:alt: Two wireframe views of a twice-refined annulus mesh. On the left, long straight seams from the coarse base mesh run across the refined mesh and meet at vertices where an unusual number of elements converge. On the right the same mesh after relaxation, with the seams much less apparent and the element sizes more even. + +A twice-refined annulus, before and after `relax()`. The seams on the left are +the coarse mesh's own edges, still visible after two refinements. Median +element quality goes from 0.940 to 0.972 and the tenth percentile from 0.839 to +0.927; the worst single cell goes the other way, 0.827 to 0.763, because the +mover minimises a global energy and will trade a few cells for many. Every node +that began on a bounding circle is still exactly on it: relaxation moves +interior coordinates only. +``` + +Relaxation moves coordinates but not topology, so the refinement relation between +the levels is untouched and the hierarchy survives it. Moving a mesh without +rebuilding it is the subject of a [related technical note](/moving-the-mesh-without-remaking-it/). + +## When the choice matters + +The test is SolKz, a standard Stokes benchmark on the unit box: viscosity +varying exponentially with depth, $\eta = e^{2Bz}$, forced by +$\mathbf{f} = (0,\; \sin(m\pi z)\cos(n\pi x))$ with $n = 3$, $m = 2$, free slip +on all four walls, Taylor–Hood $P_2$–$P_1$ elements. The viscosity contrast +across the box is $e^{2B}$. SolKz has a closed-form solution, and Underworld +carries it in `uw.analytic` — putting it to its usual use, checking that a +solver converges to the right answer, is the subject of a +[companion note](/testing-a-solver-against-exact-solutions/). Nothing here uses +it: these are cost measurements, not accuracy measurements. + +*Computational work* is reported as floating-point operations (counted by PETSc) per unknown and +relative to the constant-viscosity FMG run with a fixed target iteration tolerance across the board. Flops are a useful measure: they do not depend on the machine, and they are generally deterministic. Times are given alongside, because they are what you wait for and they reveal that different algorithms may be more or less efficient on particular hardware. + +First we examine the cost against viscosity contrast, at 11 727 unknowns: + +| preconditioner | viscosity contrast | relative work per unknown | relative time per unknown | +|---|---|---|---| +| FMG | 10⁰ | 1.0 | 1.0 | +| FMG | 10² | 1.3 | 1.1 | +| FMG | 10⁴ | 1.6 | 1.2 | +| FMG | 10⁶ | 1.8 | 1.3 | +| GAMG | 10⁰ | 2.2 | 1.2 | +| GAMG | 10² | 3.3 | 1.5 | +| GAMG | 10⁴ | 4.0 | 1.7 | +| GAMG | 10⁶ | 4.0 | 1.6 | + +Both solver configurations handle the whole range in the viscosity gradient sweep and there is almost +no change in the amount of work required or time taken to solve. GAMG does slightly more work but does so a little more efficiently. + +Next, let us look at the cost (per unknown) as we change the problem size, at constant viscosity: + +| preconditioner | unknowns | relative work per unknown | relative time per unknown | Gflop/s | +|---|---|---|---|---| +| FMG | 2 947 | 1.0 | 1.0 | 1.6 | +| FMG | 11 727 | 1.7 | 1.3 | 2.1 | +| FMG | 46 783 | 2.2 | 1.4 | 2.4 | +| FMG | 186 879 | 2.5 | 1.6 | 2.4 | +| FMG | 747 007 | 2.6 | 1.8 | 2.2 | +| GAMG | 2 947 | 2.7 | 1.3 | 3.4 | +| GAMG | 11 727 | 3.8 | 1.5 | 3.8 | +| GAMG | 46 783 | 4.0 | 1.6 | 3.9 | +| GAMG | 186 879 | 4.6 | 1.7 | 4.2 | +| GAMG | 747 007 | 5.5 | 1.9 | 4.6 | + +Both flatten (which is good). Step by step, FMG's flop count goes as $N^{1.38}$, $N^{1.20}$, +$N^{1.07}$, $N^{1.03}$; GAMG's as $N^{1.24}$, $N^{1.04}$, $N^{1.11}$, $N^{1.12}$. Both are +converging on linear, which is what multigrid of either kind is supposed to do: +the algebraic method reaches it as surely as the geometric one. + +So the geometric hierarchy is worth roughly a factor of two in throughput: FMG does about +half GAMG's arithmetic on this problem and holds on to that advantage through six orders of viscosity contrast and a 250-fold growth in problem size. + +In wall clock the advantage is smaller — 1.8 against 1.9 at the largest +size, which is no practical advantage at all. This is because GAMG runs at +3.4–4.6 Gflop/s where FMG runs at 1.6–2.4, and FMG's rate *falls* at the largest +size as the problem outgrows cache. + +The shape of the arithmetic matters as much as the amount. FMG's operators are +sparse and indirect and its coarse levels are too small to keep a processor +busy; GAMG's are denser and more regular and run closer to peak. How those +trade off is a property of the machine, so time both on yours and use whichever +wins. + +So far the viscosity has varied smoothly across the whole box, which is the +distribution algebraic coarsening reads best. The second test puts the same +total contrast into a band a twentieth of the box thick, + +$$\eta(z) = 1 + (\eta_0 - 1)\,e^{-\left((z - 1/2)/w\right)^2}, \qquad w = 0.05$$ + +and changes nothing else: same forcing, same boundary conditions, same mesh. At +contrast 1 the two problems are the same problem, and the runs agree to the last +digit. + +```{figure} figures/layer-vs-smooth.png +:alt: Two side-by-side log-log panels. Both have viscosity contrast on the horizontal axis, marked at 1, 10^2, 10^4 and 10^6. The left panel plots relative work per unknown, the right relative time per unknown. Each panel carries four curves with open circular markers at those four contrasts: blue for FMG and orange for GAMG, solid for the smooth SolKz viscosity and dashed for the concentrated band. Reading the left panel, FMG smooth runs 1.0, 1.3, 1.6, 1.8; FMG band runs 1.0, 1.9, 2.1, 2.5; GAMG smooth runs 2.2, 3.3, 4.0, 4.0. The dashed orange GAMG band curve leaves the others behind — 2.2, 8.1, 28.1 — and then stops after 10^4, with the text "does not converge beyond this point" beside its last marker. In the right panel the first three curves are packed into a narrow span between 1.0 and 1.7, while GAMG band again separates, running 1.2, 2.7, 7.7 before stopping at the same place. All four curves meet at the left-hand edge, because at contrast 1 the two viscosity distributions are the same problem. + +Cost against viscosity contrast for the two viscosity distributions, at 11 727 +unknowns. FMG is almost indifferent to the change. GAMG is not: on the band its +cost climbs with the contrast, and at $10^6$ it does not converge at all — 20 000 +multigrid cycles per velocity solve without reaching the tolerance. This is not +an artefact of an under-resolved band. At this resolution the band spans about +six cells, and resolving it better does not rescue GAMG: at four and sixteen +times the cell count, GAMG's disadvantage against FMG on the band grows from +5.5 to 7.9 to 8.0 times its disadvantage on the smooth problem at the same size. +``` + +The geometric hierarchy is indifferent to the coefficients because it never +consults them. That is the case for keeping a hierarchy when you can. + +:::{note} What this comparison does not address +The solver tolerance is held fixed as the mesh refines. That is the usual way +to show a multigrid scaling, but a finer mesh has a smaller +discretisation error, so an argument can be made for tightening the solver +tolerance alongside it, which would make the work per unknown grow. The table +should be read as "what does it cost to solve this system", not "what does it cost to +reach the accuracy the mesh can support". + +Both preconditioners are also run in serial, on one machine. GAMG's coarsening +is the part of it that changes most under parallel decomposition and is not measured here. +::: + +## Using it + +```python +# Build the hierarchy at mesh construction: base + two refinements. +mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.05, refinement=2) + +stokes = uw.systems.Stokes(mesh) +stokes.preconditioner = "auto" # FMG when a hierarchy exists +stokes.solve() + +print(len(mesh.dm_hierarchy)) # how many levels there are +print(stokes.preconditioner_settings) # what was applied +print(stokes.pc_fallbacks) # anything declined, and why +``` + +If a solve is slow and the mesh has no hierarchy, this is the first thing to +change: rebuild it with `refinement` rather than at a fine `cellSize`, and the +same resolution arrives with the coarse levels attached. Note that `cellSize` +then sets the *coarsest* mesh, and each refinement halves it — so `cellSize=0.05` +with `refinement=2` resolves like `cellSize=0.0125`. + +
Comments
Discussion of these notes happens in GitHub Discussions, so it stays with the source and is searchable alongside it.
Read the discussionStart one
diff --git a/classification.yml b/classification.yml index 68d9c8e..65201be 100644 --- a/classification.yml +++ b/classification.yml @@ -317,3 +317,8 @@ introducing-the-technical-notes: article_type: news subjects: [] methods: [] + +setting-up-full-multigrid: + article_type: technical-note + subjects: [mantle-convection] + methods: [solvers, meshing, parallel-hpc]