The analytic solution suite: uw.analytic with the full Velic family, an oracle-free validation contract, and four fixed defects

docs/api/analytic.md

@@ -0,0 +1,25 @@
+# Analytic Solutions
+
+```{eval-rst}
+.. automodule:: underworld3.analytic
+   :members:
+   :show-inheritance:
+```
+
+## The contract
+
+Every solution satisfies the same contract, so a validation run reads the same way
+whichever one you use.
+
+```{eval-rst}
+.. automodule:: underworld3.analytic._base
+   :members:
+   :show-inheritance:
+```
+
+## See also
+
+- {doc}`solvers` — the solvers these solutions validate.
+- `docs/developer/subsystems/analytic-solutions.md` — the implementation form,
+  the validation protocol every transcription must pass, and the provenance of
+  each vendored reference kernel.

docs/api/function.md

@@ -78,8 +78,5 @@ Factory function for creating UWQuantity objects with units.
 
 ## Analytic Functions
 
-```{eval-rst}
-.. automodule:: underworld3.function.analytic
-   :members:
-   :show-inheritance:
-```
+The analytic solutions have moved to {doc}`analytic` — `underworld3.function.analytic`
+is a deprecation shim. Use `uw.analytic.SolCx(mesh, ...)`.

docs/api/index.md

@@ -29,6 +29,7 @@ model
 utilities
 visualisation
 adaptivity
+analytic
 ```
 
 ## Quick Links
@@ -50,6 +51,9 @@ adaptivity
 - **{doc}`scaling`** - Units, quantities, and non-dimensionalisation
 - **{doc}`maths`** - Mathematical operations and integrals
 
+### Validation
+- **{doc}`analytic`** - Exact solutions for benchmarking and convergence testing
+
 ### Infrastructure
 - **{doc}`model`** - Model management and configuration
 - **{doc}`utilities`** - I/O, mesh import, and helper functions

docs/developer/index.md

@@ -37,6 +37,7 @@ same topic are reference or historical material subordinate to the governing doc
 | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) |
 | Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) |
 | Rotated free-slip & wall-normal datum | [subsystems/rotated-freeslip.md](subsystems/rotated-freeslip.md) |
+| Analytic & benchmark solutions | [subsystems/analytic-solutions.md](subsystems/analytic-solutions.md) |
 | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) |
 | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) |
 | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) |
@@ -197,6 +198,7 @@ subsystems/constitutive-models-anisotropy
 subsystems/swarm-system
 subsystems/data-access
 subsystems/interpolation
+subsystems/analytic-solutions
 subsystems/expressions-functions
 subsystems/containers
 subsystems/checkpointing-system

docs/examples/fluid_mechanics/advanced/Ex_Stokes_Annulus_Benchmark_Kramer.py

@@ -72,7 +72,12 @@
 import numpy as np
 import sympy
 import os
-import assess
+# `assess` (Kramer et al. 2021) is an optional dependency: pip install
+# "underworld3[benchmarks]". Imported through the wrapper so a missing install
+# reports what to do about it rather than a bare ModuleNotFoundError.
+from underworld3.analytic.kramer import require_assess
+
+assess = require_assess()
 import h5py
 from enum import Enum
 

docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_SolC.py

@@ -14,38 +14,41 @@
 
 # %% [markdown]
 """
-# Stokes Benchmark SolCx
+# Stokes Benchmark SolC
 
 **PHYSICS:** fluid_mechanics
 **DIFFICULTY:** advanced
 
 ## Description
 
-The SolCx benchmark tests the Stokes solver with a sharp viscosity contrast.
-A vertical step in viscosity at x = 0.5 creates a challenging test case for
-iterative solvers. Compares Dirichlet and natural boundary conditions.
+SolC is the *isoviscous* benchmark: a dense column occupying half the box drives
+flow through a fluid of uniform viscosity. The difficulty is the discontinuous
+forcing, not a material contrast — the pressure has a kink at the column edge
+that a discretisation has to resolve.
+
+(For the sharp **viscosity** contrast, see SolCx. This file was previously
+titled and documented as SolCx while solving SolC, which is worth knowing if
+you are comparing old output.)
 
 ## Key Concepts
 
-- **Viscosity jump**: Step function viscosity contrast (10^6)
-- **Benchmark validation**: Comparison with analytical solution (if UW2 available)
-- **Natural vs Dirichlet BCs**: Free-slip implemented different ways
-- **Piecewise functions**: Using sympy.Piecewise for sharp interfaces
-- **Multigrid preconditioning**: Essential for high viscosity contrast
+- **Discontinuous body force**: buoyancy steps at x = 0.5, viscosity is uniform
+- **Benchmark validation**: against `uw.analytic.SolC`, in-tree and validated
+- **Truncated series**: the exact solution is a Fourier sum, `modes` terms of it
+- **Free slip on all walls**, with a pressure null space
 
 ## Mathematical Formulation
 
-Viscosity step function:
-$$\\eta(x) = \\begin{cases} 10^6 & x > 0.5 \\\\ 1 & x \\le 0.5 \\end{cases}$$
-
-Buoyancy forcing:
-$$f_y = -\\cos(\\pi x) \\sin(2\\pi y)$$
+Uniform viscosity $\\eta = 1$, with buoyancy stepping at the column edge. Both
+the forcing and the exact velocity come from `uw.analytic.SolC`, so they cannot
+disagree about sign or side — see the note in the validation section below.
 
 ## Parameters
 
 - `uw_resolution`: Mesh resolution
 - `uw_refinement`: Mesh refinement level
-- `uw_viscosity_contrast`: log10 of viscosity contrast
+- `uw_modes`: Fourier modes in the SolC analytic solution
+- `uw_viscosity_contrast`: log10 contrast, used by the SolCx section
 """
 
 # %% [markdown]
@@ -87,14 +90,15 @@
 params = uw.Params(
     uw_resolution = 4,              # Base mesh resolution
     uw_refinement = 2,              # Mesh refinement levels
-    uw_viscosity_contrast = 6,      # log10 of viscosity contrast
     uw_use_simplex = 1,             # Use simplex mesh (1) or quad (0)
     uw_penalty = 100,               # Stokes penalty parameter
+    uw_modes = 40,                  # Fourier modes in the SolC analytic solution
+    uw_viscosity_contrast = 6,      # log10 contrast, for the SolCx section below
 )
 
 # Derived parameters
-eta_ratio = 10 ** params.uw_viscosity_contrast
 use_simplex = bool(params.uw_use_simplex)
+eta_ratio = 10 ** params.uw_viscosity_contrast
 
 # %% [markdown]
 """
@@ -159,20 +163,19 @@
 """
 
 # %%
-eta_0 = 1
 x_c = sympy.Rational(1, 2)
-f_0 = 1
+
+# The exact solution supplies the forcing as well as the answer. Writing the
+# body force out by hand here is how this file came to solve a mirrored,
+# sign-flipped problem from the one it compared against: SolC's buoyancy is
+# negative on x < x_c, and the Piecewise previously used was +1 on x > x_c.
+# Nobody noticed, because the comparison sat behind an `import underworld` that
+# always failed.
+solC = uw.analytic.SolC(mesh, x_c=x_c, modes=int(params.uw_modes))
 
 stokes.penalty = params.uw_penalty
-stokes.bodyforce = sympy.Matrix(
-    [
-        0,
-        Piecewise(
-            (f_0, x > x_c),
-            (0.0, True),
-        ),
-    ]
-)
+stokes.constitutive_model.Parameters.shear_viscosity_0 = solC.fn_viscosity
+stokes.bodyforce = solC.fn_bodyforce
 
 # Free-slip boundary conditions (Dirichlet form)
 stokes.add_dirichlet_bc((sympy.oo, 0.0), "Top")
@@ -208,6 +211,25 @@
 # %%
 stokes.solve()
 
+# %% [markdown]
+"""
+## Validation against the analytic solution
+
+This has to happen **here**, not at the end of the file. Everything below
+reconfigures the same solver for other experiments, and `v` then holds those
+answers rather than this one — which is exactly how the old check came to
+compare a SolC analytic solution against a SolCx solve with penalty boundary
+conditions.
+"""
+
+# %%
+# `error` is a global reduction, so this is the same number on any rank count.
+solC_velocity_error = solC.error("velocity", v)
+solC_pressure_error = solC.error("pressure", p)
+
+uw.pprint(f"SolC relative velocity error: {solC_velocity_error:.6e}")
+uw.pprint(f"SolC relative pressure error: {solC_pressure_error:.6e}")
+
 # %% [markdown]
 """
 ## SolCx Benchmark Configuration
@@ -231,7 +253,7 @@
 timing.reset()
 timing.start()
 stokes.solve(zero_init_guess=True)
-timing.print_table(display_fraction=0.999)
+timing.print_table()  # see #499: display_fraction was removed from the API
 
 # Save solution with Dirichlet BCs
 v0.data[...] = v.data[...]
@@ -256,7 +278,7 @@
 timing.reset()
 timing.start()
 stokes.solve()
-timing.print_table(display_fraction=0.999)
+timing.print_table()  # see #499: display_fraction was removed from the API
 
 v1.data[...] = v.data[...]
 
@@ -278,7 +300,7 @@
 timing.reset()
 timing.start()
 stokes.solve()
-timing.print_table(display_fraction=0.999)
+timing.print_table()  # see #499: display_fraction was removed from the API
 
 # %% [markdown]
 """
@@ -330,32 +352,19 @@
         show_scalar_bar=False,
     )
 
-    pl.show(cpos="xy")
-
-# %% [markdown]
-"""
-## Validation Against UW2 (if available)
-"""
+    # Only when there is somewhere to show it. Guarded on mpi.size alone, this
+    # blocks a script run forever waiting on a window that never opens — which
+    # is why running this file to completion was not something anyone had done.
+    if uw.is_notebook:
+        pl.show(cpos="xy")
 
 # %%
-try:
-    import underworld as uw2
-
-    solC = uw2.function.analytic.SolC()
-    vel_soln_analytic = solC.fn_velocity.evaluate(mesh.X.coords)
-    from mpi4py import MPI
-    from numpy import linalg as LA
-
-    comm = MPI.COMM_WORLD
-
-    num = function.evaluate(v.fn, mesh.X.coords)
-    if comm.rank == 0:
-        print(f"Velocity difference norm: {LA.norm(v.data - vel_soln_analytic):.6e}")
-    comm.barrier()
-except ImportError:
-    import warnings
-
-    warnings.warn("Unable to validate against UW2 analytical solution (UW2 not available).")
-
-# %%
-print(f"SolCx benchmark complete: resolution {n_els}, refinement {refinement}")
+uw.pprint(
+    f"Complete: resolution {n_els}, refinement {refinement}, "
+    f"modes {int(params.uw_modes)}"
+)
+uw.pprint(f"  SolC velocity error (validated above): {solC_velocity_error:.6e}")
+uw.pprint(
+    "  The SolCx and natural-BC solves that follow it are BC experiments, "
+    "not benchmarks — nothing here compares them against an exact solution."
+)

docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_SolNL.py

@@ -78,24 +78,19 @@
 # %%
 # NL problem
 # Create solution functions
-from underworld3.function.analytic import (
-    AnalyticSolNL_velocity,
-    AnalyticSolNL_bodyforce,
-    AnalyticSolNL_viscosity)
+from underworld3 import analytic as A
 
 x, y = mesh.X
 
-r = mesh.r
 eta0 = 1.0
 n = 1
 r0 = 1.5
-params = (eta0, n, r0)
-sol_bf_ijk = AnalyticSolNL_bodyforce(*params, *r)
-sol_vel_ijk = AnalyticSolNL_velocity(*params, *r)
 
-sol_bf = mesh.vector.to_matrix(sol_bf_ijk)
-sol_vel = mesh.vector.to_matrix(sol_vel_ijk)
-sol_visc = AnalyticSolNL_viscosity(*params, *r)
+sol = A.SolNL(mesh, eta_0=eta0, n=n, r=r0)
+
+sol_vel = sol.fn_velocity
+sol_bf = sol.fn_bodyforce
+sol_visc = sol.fn_viscosity
 
 # debug - are problems just because there is no analytic solution module on mac
 # The solNL case is a MMS force term (complicated) designed to produce a specific

docs/examples/fluid_mechanics/advanced/Ex_Stokes_Spherical_Benchmark_Kramer.py

@@ -79,7 +79,12 @@
 import numpy as np
 import sympy
 import os
-import assess
+# `assess` (Kramer et al. 2021) is an optional dependency: pip install
+# "underworld3[benchmarks]". Imported through the wrapper so a missing install
+# reports what to do about it rather than a bare ModuleNotFoundError.
+from underworld3.analytic.kramer import require_assess
+
+assess = require_assess()
 import h5py
 import sys
 

docs/examples/fluid_mechanics/advanced/Ex_Stokes_Spherical_Benchmark_Thieulot.py

@@ -81,7 +81,12 @@
 import numpy as np
 import sympy as sp
 import os
-import assess
+# `assess` (Kramer et al. 2021) is an optional dependency: pip install
+# "underworld3[benchmarks]". Imported through the wrapper so a missing install
+# reports what to do about it rather than a bare ModuleNotFoundError.
+from underworld3.analytic.kramer import require_assess
+
+assess = require_assess()
 import h5py
 import sys
 from petsc4py import PETSc

docs/examples/snesfas_investigation/benchmark_3way.py

@@ -19,7 +19,7 @@
 import numpy as np
 import sympy
 import underworld3 as uw
-from underworld3.function import analytic as A
+from underworld3 import analytic as A
 
 fv = importlib.util.module_from_spec(
     importlib.util.spec_from_file_location("fv", os.path.join(os.path.dirname(os.path.abspath(__file__)), "fas_vanka.py")))

docs/examples/utilities/advanced/Ex_Stokes_Annulus_Benchmark_Kramer_etal.py

@@ -44,7 +44,12 @@
 import os
 import matplotlib.pyplot as plt
 import cmcrameri.cm as cmc
-import assess
+# `assess` (Kramer et al. 2021) is an optional dependency: pip install
+# "underworld3[benchmarks]". Imported through the wrapper so a missing install
+# reports what to do about it rather than a bare ModuleNotFoundError.
+from underworld3.analytic.kramer import require_assess
+
+assess = require_assess()
 # -
 
 os.environ["SYMPY_USE_CACHE"] = "no"

docs/examples/utilities/advanced/Ex_Stokes_Annulus_Benchmark_Thieulot.py

@@ -44,7 +44,12 @@
 import os
 import matplotlib.pyplot as plt
 import cmcrameri.cm as cmc
-import assess
+# `assess` (Kramer et al. 2021) is an optional dependency: pip install
+# "underworld3[benchmarks]". Imported through the wrapper so a missing install
+# reports what to do about it rather than a bare ModuleNotFoundError.
+from underworld3.analytic.kramer import require_assess
+
+assess = require_assess()
 # -
 
 os.environ["SYMPY_USE_CACHE"] = "no"