Uh oh!
There was an error while loading. Please reload this page.
The analytic solution suite: uw.analytic with the full Velic family, an oracle-free validation contract, and four fixed defects - #571
Merged
Conversation
Underworld2 shipped twelve exact ("Velic") Stokes solutions as the code's
source of truth for benchmarking. Two reached UW3, and the module they
landed in cannot be extended: underworld3.function.analytic is a compiled
Cython extension, so the name is owned by a .so and nothing else can live
under it. Exact solutions have leaked elsewhere as a result -- Gardner in
utilities/retention_curves.py, erfc diffusion and Ogata-Banks inline in
tests, and an undeclared 'assess' dependency for the Kramer benchmarks.
This is the interface commit: a new underworld3.analytic package holding
the contract every solution satisfies, so later work adds solutions rather
than restructuring around them. Nothing moves yet -- SolCx is re-exported
from its current home (the same class object, not a copy), so both import
paths stay valid.
AnalyticSolution gives each solution the exact fields (fn_velocity,
fn_pressure, fn_stress, fn_strainrate, fn_viscosity, fn_bodyforce), the
LaTeX that documents the problem it poses, evaluate() at arbitrary points,
and error() against a computed field. The nodal error norm is the global
MPI reduction from velocity_error generalised over field name, and it now
has the test that pins it: a perturbation confined to x < 0.5, which a
rank-local norm would report differently on every rank (the #370 failure).
A uniform perturbation cannot detect that.
The two boundary-condition mixins configure the solver in place rather
than returning something the caller applies. FreeSlipWalls uses the strong
rotated constraint, not component masking -- the two agree on an
axis-aligned box, but only the rotated form still holds when a solution is
used to validate a curved or adapted mesh.
docs/developer/subsystems/analytic-solutions.md records the decisions this
suite is built on: solutions are pure SymPy (so they JIT, carry a Jacobian,
and work as Dirichlet values), reference C kernels are kept as independent
oracles rather than deleted, and every transcription must clear six gates
before it lands -- including a negative control, because a gate that passes
a deliberately broken input is measuring nothing.
Verified: 17 contract tests pass at np=1, 2 and 3; test_1015_analytic_solcx
and test_1062_constrained_solcx pass unchanged; the style gate is clean.
Underworld development team with AI support from Claude Codeunderworld3.function.analytic was a compiled extension, so the name was owned by a .so and nothing could live under it -- no submodules, no registry, no pure-sympy solutions. That is what has kept the exact-solution suite from growing. The extension is now underworld3.analytic._reference._velic, and the old name is a deprecation shim. No behaviour changes and no test-file edits: the ten existing consumers pass unmodified, which is the acceptance criterion for this commit. The shim is a package DIRECTORY, not a module file, and must stay one. An orphaned analytic.cpython-*.so from an earlier install cannot be removed by pip uninstall once it has fallen out of the wheel RECORD, and it would be imported in place of the shim -- silently restoring the old module and bypassing every redirect. Python's path finder checks for a package directory before an extension of the same name, which is the only thing that makes such an orphan harmless. Verified by planting one: the directory wins, and SolCx still resolves. test_legacy_namespace_is_a_package_not_an_extension guards it. Names resolve through the shim's __getattr__ rather than being imported at module level, so the deprecation warning fires when a name is USED. Ten test files import this module at collection time, where a warning is noise the reader cannot act on. Object identity is preserved -- the shim returns the same classes, so isinstance and pickle work across both paths. uw.function.__getattr__ uses importlib.import_module rather than `from . import analytic`: the latter resolves the submodule by calling getattr() on the parent, which lands straight back in __getattr__ and recurses until the stack runs out. Found by running it. package_data gains its own "underworld3.analytic._reference" key. The existing "underworld3" globs do not reach a directory that is itself a package, and the JIT adds only the module's own directory to its include path -- so a missing header fails at solve time, when the kernel compiles, not at import. Confirmed the .so and all three headers appear in the installed RECORD. Verified: clean rebuild (build/lib.*, build/temp.* removed first, since a stale tree repackages the old .so and Cython's cached .c embeds the old module name); 20 contract tests at np=1 and np=2; test_1015_analytic_solcx and test_1062_constrained_solcx pass unchanged; the other five consumers and both parallel test files collect (68 tests); both the attribute and from-import paths work; style gate clean. Underworld development team with AI support from Claude Code
Groundwork for delivering the analytic solutions as SymPy. Not wired in: uw.analytic.SolCx is still the reference kernel, and velic.py is not exported, because two validation gates are open. _transcribe.py reads the published straight-line C -- the t125 = 0.4e1*t81*t83 + ... form -- and rebuilds the same expression tree in SymPy. It preserves the generator's grouping term for term and never simplifies: these expansions carry products that are stable only in the arrangement Maple produced, and a re-derivation is a different arrangement that can lose eight digits in exactly the regime the benchmarks probe. Transcription happens at run time rather than generating a checked-in module. Measured on solCx.c -- 1500 lines, the largest in the family -- reading both arrangements and both spatial branches takes 0.25 s, and the resulting expressions are 2000-7000 operations. That is cheap enough to do on demand, and doing it on demand means the SymPy form cannot drift from the C it came from: they are one artefact, not two copies. The measurement is also what shows the all-SymPy target is viable for the whole family, which was not obvious before. Numeric literals become exact Rationals. Maple writes them as 0.4e1, i.e. exact small values, so nothing is lost and a solution can be evaluated at arbitrary precision -- which is how a transcription error is told apart from the kernel's own double-precision cancellation. Gates 1 and 2 on the _solCx_A arrangement pass at 1e-14 to 1e-16 across every regime tested: contrast 1e-6 to 1e8, both orderings, several x_c and n. The table is in docs/developer/subsystems/analytic-solutions.md. Two gates are open and are why nothing is exported: - The _solCx_B arrangement transcribes to a different answer. The kernel dispatches on viscosity ordering purely for conditioning, so the two arrangements should agree -- and the compiled ones do, since the _solCx_A transcription reproduces the dispatcher's output in the regime where the dispatcher runs _solCx_B. Transcribing _solCx_B directly is wrong by a factor of tens everywhere. Both parse to the same structure, so the reader treats them identically and one is still being read wrongly. Unexplained. - The isoviscous case returns zero where the kernel does not, pointing at a 0/0 at unit viscosity ratio; the closed form carries ZR - 1 in denominators. Method note recorded in the subsystem doc: the first gate run reported 1e12 errors, which was the metric, not the transcription. A pointwise relative error divides by the true value, and these fields pass through zero. Normalising by the field magnitude over the sample is the honest measure -- suspect the metric before the result when every case fails alike. Underworld development team with AI support from Claude Code
uw.analytic.SolCx is now SymPy rebuilt from the published Maple kernel rather than a call into it. Same class, same name, drop-in: velocity_error, evaluate_stress and topography_top are kept, so the ten existing consumers pass unmodified and object identity across both namespaces still holds. What this buys is what a compiled kernel cannot do. The fields carry an analytic Jacobian, compile into a residual through the normal JIT path, and can be used as a Dirichlet boundary value -- test_transcription_is_usable_by_the_solver pins that last one, since it is the capability the whole change exists for. The old per-point evalf loop is gone with it. Validated against the kernel it came from, over ratios 1e-6 to 1e8 in both directions, wavenumbers 1 to 3, and an off-centre interface: worst-case normalised error 1e-14 to 1e-16, sampled with 40 stratified points plus the viscosity interface from both sides, the walls and the corners. Agreement at a handful of points was 7e-18. Three things looked like transcription failures and were not. All three are now tests, because each cost real time to find: _solCx_B is not a second conditioning of the same formula -- it is the mirror, A(x,z) = -B(1-x,z). The source dispatches on eta_A > eta_B, which reads as a conditioning choice, so transcribing both and picking looked obviously right; doing that gave answers wrong by a factor of tens. Evaluated exactly at 50 digits the relationship is clean, and the sign is the forcing cos(pi x) being odd about x = 1/2. Only _solCx_A is transcribed and there is no dispatch, which is safe because the stated reason for the dispatch was measured rather than assumed -- the error never leaves 1e-14 anywhere in the range. test_arrangements_are_mirror_images records the evidence, so if it ever fails we know the reasoning needs revisiting. Equal viscosities are a removable singularity -- the closed form carries (ZR - 1) in denominators. SymPy cancels it evaluating symbolically at a point, verified to 40 digits, but not in the compiled form, where it survives as 0/0. SolCx now raises rather than returning nonsense; uniform viscosity is a different benchmark. Parameters are substituted as exact Rationals regardless, since that is what allows the cancellation at all and it costs nothing. The first validation run reported errors of 1e12 and looked catastrophic. It was the metric: a pointwise relative error divides by the true value, and these fields pass through zero. The values were close all along. The test normalises by field magnitude and says so. Verified: 37 tests pass serially (the four analytic files, including the unchanged Stokes convergence test now driven by the transcription), 32 at np=2, style gate clean. Underworld development team with AI support from Claude Code
… code The gates were demonstrated ad hoc while transcribing SolCx. This makes them underworld3.analytic._validation, so the remaining eleven solutions get them by calling four functions rather than reinventing them, and so a regression in any of them is a test failure rather than something nobody re-runs. Two of the checks need no reference at all, and those are the strongest here. incompressibility_residual and momentum_residual put the fields back into the equations they claim to solve: div(u) and div(sigma) + f, using the solution's own stress and body force. They catch the failure a convergence test structurally cannot -- if a transcription and the solver share a mistaken convention the solve converges neatly to the wrong answer -- and they are what settles the body-force sign, where UW2's documentation and UW3's convention disagree. Measured on SolCx at contrast 1e3: 3.6e-17 and 2.3e-16. strainrate_consistency sits between comparison and physics. The kernels derive velocity and stress separately, so differentiating one and checking it against the other is a real cross-check, and it exercises the derivatives -- which is what a solver consumes and where a transcription can be wrong while still matching pointwise. 8.5e-16. Gate 5, the negative control, stays in each solution's test rather than the harness, because what counts as a plausible slip is solution-specific. test_the_checks_reject_a_broken_transcription perturbs one velocity coefficient by a part in a thousand and requires both the comparison and the oracle-free residual to report it. Without it the other checks are unfalsified: a check that passes a deliberately broken input is measuring nothing. The harness evaluates through lambdify rather than uw.function.evaluate, and that is not an optimisation. These checks differentiate the fields, and a viscosity-jump solution puts a large Piecewise inside a stress derivative; the JIT path took so long to generate and compile that a three-regime run did not finish in 45 minutes -- the same blow-up already recorded for add_nitsche_bc on SolCx. Lambdified, each check is under half a second. The expressions are pure SymPy in the mesh coordinates so this is exact, not an approximation; the one subtlety is that mesh coordinates cannot be bound as lambdify arguments and must be swapped for plain symbols first. All nine regimes now run every check: 31 tests in 2m22s, worst case 1e-10. Underworld development team with AI support from Claude Code
SolNL had no convenience class, and three of its six published entry points -- pressure, stress, strain rate -- were compiled but never reachable from Python. It is now a transcribed SymPy solution on the contract, with all six fields, and it is the first nonlinear one: the viscosity depends on the second invariant of the strain rate the solution itself produces, so it exercises a nonlinear solver rather than a linear one. Putting a second kernel through the transcriber was the point, and it found two defects SolCx could not have. The reader took the last identifier before `=` as the assignment target. SolNL writes its results through a struct, `out.x = ...`, which that rule reads as an assignment to `x` -- silently rebinding the coordinate. Every later statement using x then got the velocity component instead. The result was not a crash or an obvious mess: fn_velocity came out as exp(velocity_x)*sin(pi n z), a perfectly plausible-looking expression that happens to be wrong. Targets now keep any struct prefix. SolCx was unaffected because it writes through arrays that the tail truncation already excluded, which is exactly why one validated transcription is not evidence the reader is correct. C statements also wrap freely across lines, and a wrapped Python expression with indented continuations is a syntax error. SolCx's statements happened to be single-line. Expressions are now folded before evaluation. Two smaller additions the kernel needed: functions whose result is returned rather than assigned (evaluate_expression, CSource.returned), and non-void signatures (CSource.function(..., returns=)). Validated the same way as SolCx -- agreement with the published kernel at 1e-12 or better across three parameter sets, and divergence-free with no oracle. Plus one check worth having because it is cheap: the published velocity is short enough to assert outright, which catches a mangled read instantly. Verified: 37 transcription tests in 2m25s; the analytic consumers pass unmodified (25 tests); style gate clean; uw.analytic.available() reports both. Underworld development team with AI support from Claude Code
… yet exported Groundwork for the elliptical-inclusion benchmark (GJI 155, 269-288). The physics is settled; the representation is not, so nothing is exported and uw.analytic is unchanged. The authors' reference MATLAB publishes pressure, deviatoric stress and the rotation rate but not the velocity, so the Muskhelishvili potentials have to be recovered from what is there. phi comes from the matrix pressure via p = -2 Re[phi'(z)], giving phi'(z) = A/(zeta^2 - 1). That reading is then checked against something independent: the stress expression contains a term that must equal phi''(z) derived from the same phi'. It does, identically -- sympy.simplify of the difference is exactly zero. psi' is the remaining bracket, and psi'(inf) = -BC, the constant far field it should be. The reconstructed velocity is divergence-free to 2e-14 in the matrix, which is the first real evidence the reconstruction is right rather than merely plausible. Two representation problems remain, both recorded in the module docstring. Inverting z = zeta + 1/zeta as sqrt(z**2 - 4) cuts along a ray, so left of the origin it selects the root inside the unit circle -- the wrong sheet -- and the far field comes out asymmetric, about three times the imposed shear at (-50, 20). Writing it sqrt(z-2)*sqrt(z+2) cuts along the segment [-2, 2], which is the slit the map already has, and is correct. But SymPy will not then push re/im through it, and differentiating gives an unevaluated Derivative(re(...)) no code printer can emit. Building the components as (w + conj w)/2 and (w - conj w)/2i sidesteps re/im entirely; untried. The interior velocity is also still missing -- pressure and viscosity are Piecewise across the boundary but the velocity is not, so it is wrong inside. The interior is a uniform velocity gradient, fixed by the interior deviatoric stress and the rotation rate, both already computed here. Both are finishable, and the validation for them is already in place: the momentum and incompressibility residuals need no oracle, so they will confirm or refute the result directly, with the published pressure, interface pressure and rotation rate as three further independent checks. Two SymPy traps found on the way, noted in the code because they cost time. Integrating psi with the numeric constants already substituted puts SymPy in a floating complex polynomial ring where the division algorithm cannot detect zero and integration fails outright; integrating the zeta-shape once with a bare symbol keeps it exact. And the complex expression must be built on real-declared symbols with the mesh coordinates substituted at the end -- mesh coordinates carry no reality assumption, so re/im cannot be distributed through them. Verified: 57 analytic tests still pass, style gate clean, uw.analytic.available() unchanged at SolCx and SolNL. Underworld development team with AI support from Claude Code
uw.analytic.EllipticalInclusion: a viscous ellipse in a matrix under far-field general shear, with closed-form velocity and pressure inside and outside and no restriction on the viscosity ratio. No body force -- the flow is driven entirely by the far field, so it tests how a solver handles a strong contrast on a curved interface rather than how it handles forcing. This one is derived, not transcribed. The authors' MATLAB publishes pressure, stress and the rotation rate but not velocity, so the Muskhelishvili potentials had to be recovered from the fields and the velocity built from those. With no kernel to compare velocity against, the validation is physics and internal consistency: Stokes residual eta lap(v) - grad(p) 1.4e-17 (v reconstructed, p published) incompressibility 1.7e-16 velocity continuity across the interface 1e-5 at a 1e-7 step far field vs the imposed shear few parts in 1e6 (the 1/r^2 tail) interior strain rate uniform 1e-12 The first three are cross-checks between things derived separately, not restatements: the pressure is the published closed form, and the interior field comes from the published interior stress and rotation rate while the exterior comes from the potentials. Two traps, both of which produced a plausible wrong answer rather than an obvious one. A purely imaginary constant in phi' is invisible to the published data -- pressure is -2 Re[phi'] and stress involves phi'' -- but it is a far-field rigid rotation. Omitting it gives a flow with exactly the right strain and no spin, so an imposed simple shear comes back as pure shear at the correct magnitude. Its value came from a different published expression: taken to a circle the rotation rate collapses to -gr/2 for every viscosity ratio. When reading potentials back out of fields, ask what the fields are blind to. Inverting z = zeta + 1/zeta as sqrt(z**2 - 4) cuts along a ray and picks the root inside the unit circle for x < 0 -- the wrong sheet -- making the far field asymmetric, about three times too fast on one side. sqrt(z-2)*sqrt(z+2) cuts on [-2, 2], the slit the map already has. The test samples negative x deliberately; positive-only sampling would have missed it. That correct branch then defeats SymPy's re()/im(), which survive into derivatives as an unprintable Derivative(re(...)). The components are built as (w + conj w)/2 and (w - conj w)/2i instead, with conjugation done by flipping the sign of I -- for an expression in real symbols that is exactly conjugation, and unlike sympy.conjugate it distributes through a square root. Verified against numpy.conj on both sides of the cut. Verified: 16 inclusion tests in 16s; 78 analytic tests overall; style gate clean; uw.analytic.available() now lists EllipticalInclusion, SolCx, SolNL. Underworld development team with AI support from Claude Code
uw.analytic.SolKx — Stokes flow with eta = exp(2Bx) on the unit box, free slip everywhere, forced by (0, sin(m pi z) cos(n pi x)). The companion to SolCx: same geometry and forcing shape, but the viscosity varies smoothly instead of jumping, and the two fail differently. A jump tests how a discretisation copes with a discontinuity inside an element; a gradient tests whether the operator stays conditioned while the contrast builds across every element. Over the unit box the total contrast is exp(2B), so B = 5 already spans four orders. Transcribed from PETSc's copy of the kernel rather than Underworld2's: it is self-contained, returns every field in one call, and is maintained upstream. The source text is vendored at analytic/_reference/solKx.c with its BSD-2 notice, as transcription input rather than built code, and package_data now ships the .c alongside the headers. Validated without an oracle, and that is a deliberate choice rather than a shortcut. The forcing and the boundary conditions are both known, so by uniqueness a field set satisfying Stokes with them IS the solution: |div(sigma) + f| / |f| 2.5e-16 |div(v)| 4.3e-19 |v.n| on all four walls 1.4e-19 The tests found a real footgun. PETSc notes that the kernel admits non-integral m, and the first draft passed that through. But the vertical velocity carries sin(m pi z), which vanishes at z = 1 only for integer m -- so a fractional value still solves the equations while silently ceasing to satisfy free slip on the top wall, and the benchmark quietly becomes a different problem. Every residual check would still pass. m is now required to be a positive integer, with the reason in the error, and the test asserts the refusal. Three small transcriber additions, all mechanical: the PETSc spellings of the maths functions (PetscExpReal and friends) alongside the plain-C ones, C cast stripping since (PetscReal)n is juxtaposition in Python, and array-valued inputs so a kernel that reads its coordinates from pos[] can be bound. One note for the next transcription: these expressions run to tens of thousands of operations, so lambdify once per expression over the whole point set, not once per point. Doing it per point turned a two-minute suite into one that did not finish. Verified: 11 SolKx tests in 2m09s; 77 analytic tests elsewhere still pass; style gate clean; uw.analytic.available() now lists EllipticalInclusion, SolCx, SolKx, SolNL. Underworld development team with AI support from Claude Code
… 3D one Two Dohrmann-Bochev solutions transcribed from Underworld2's headers. SolDB2d is isoviscous; SolDB3d (Burstedde et al. 2013) carries a smooth viscosity peaked in the interior, exp(1 - beta[x(1-x)+y(1-y)+z(1-z)]), and is the suite's first 3D solution. That 3D gap mattered. Several parts of a Stokes discretisation genuinely differ between two and three dimensions -- the pressure space, the null space, the tensor assembly -- and no 2D benchmark can see a term that is wrong only in the third. SolDB3d also varies its viscosity in every direction at once, which none of the others do. These are the easiest solutions here to be sure of. The fields are short enough that div(v) and div(sigma) + f reduce SYMBOLICALLY to zero rather than to something small, so the tests assert exact equality: no sampling, no tolerance, no conditioning question. One convention trap, now pinned by a test. Unlike SolCx and SolKx, these kernels publish the DEVIATORIC stress rather than the total, so the pressure has to go back in as sigma = tau - p I. Reading the deviator as the total would leave the momentum residual wrong by exactly grad(p) -- large, but structured, and easy to misread as a transcription error rather than a convention one. Two transcriber additions, both from these files being C++ headers rather than C: identifiers Python reserves are renamed (these kernels take coordinates as `const double* in`, and `in[0]` does not parse), and a declaration packing several declarators into one statement is split, since `double x=in[0],y=in[1];` would otherwise be read as a single assignment whose value runs past the comma. The tests also caught a packaging gap of the kind PR 1 warned about: the new .hpp files were not in package_data, so they built fine and then failed at run time in the installed tree. package_data now covers .h, .hpp and .c. Two test-side notes worth keeping. `simplify` will not combine exponentials written in mesh coordinates -- the beta = 0 cases reduced and the others did not, purely because of the symbol type -- so the residual is rewritten over plain symbols first. And an exact Rational 4 and a float 4.0 in an exponent are equal but SymPy will not cancel them, so the expected form has to be built the same way the solution substitutes. Verified: 13 SolDB tests in 7.6s; 101 analytic tests overall; style gate clean. uw.analytic.available() now lists EllipticalInclusion, SolCx, SolDB2d, SolDB3d, SolKx, SolNL. Underworld development team with AI support from Claude Code
…ot uniform uw.analytic.SolKz — Stokes flow with eta = exp(2Bz), free slip on the unit box. The vertical twin of SolKx and not a redundant one: a viscosity varying with depth stratifies the flow along the direction buoyancy acts, coupling pressure and vertical velocity through the varying coefficient in a way a horizontal gradient never does, and it is the closer analogue of a real mantle profile. Validated by the equations, as SolKx was: |div(sigma)+f|/|f| is 2.2e-16 to 3.8e-16 across four regimes, div(v) ~1e-18, free slip ~1e-19 on all four walls. Two traps here, and the second is the one worth carrying forward. SolKz transposes SolKx. Its modes run in x rather than z, and its u1 is the VERTICAL velocity where SolCx and SolKx use u1 for the horizontal. The mapping is taken from the kernel's own output section rather than assumed, because reading it with the SolCx convention would silently transpose the entire solution -- div(v) would still vanish, free slip would still hold, and only the momentum residual would notice. The stress convention is not uniform across this family. SolCx and SolKx publish the total Cauchy stress; SolKz publishes the DEVIATOR -- into an array it calls `total_stress`. Following the name leaves the momentum residual at order |f| and invents a horizontal body force in a benchmark that has none: large, structured, and easy to misread as a bad transcription. Two cheap signatures separate them, and both are now standing tests. A deviator is traceless, so its xx and zz entries are exact negatives -- which the kernel's output visibly was. And tau = 2 eta edot, where the strain rate comes from the velocity, a different output of the same kernel. On SolKz the shear component agreed with 2 eta edot to machine precision while the normal components agreed with nothing, which located it at once. Recorded in the subsystem doc as a table of which solution publishes which, with the instruction not to trust the name. Verified: 12 SolKz tests; 49 neighbouring analytic tests; style gate clean. uw.analytic.available() now lists EllipticalInclusion, SolCx, SolDB2d, SolDB3d, SolKx, SolKz, SolNL. Underworld development team with AI support from Claude Code
… of them The solutions had drifted apart. Each assembled its own fn_* attributes and each had its own test file, and that combination let a real error through: SolNL's kernel publishes the deviatoric stress, it was stored as the total, and its momentum residual was 1.06 rather than zero. Its test file checked agreement with the kernel and incompressibility -- both passed -- and nothing checked the momentum balance. Fixed here, and made structurally hard to repeat. Assembly happens once. A solution hands its components to AnalyticSolution.set_fields, which applies the conventions; whether the source publishes sigma or tau is a class-level declaration, stress_is_deviatoric, honoured in exactly one place. Four solutions previously did this by hand, two of them differently. Conformance is checked for every registered solution. tests/test_1024_analytic_conformance.py iterates over uw.analytic.available() and applies the same six checks to all seven: contract populated, metadata declared, incompressible, momentum balance, stress and strain rate consistent. A solution added later is covered the moment it is registered. 35 checks, 101 s. Where a solution differs it says so through the contract rather than being exempted -- sample_points is new for exactly this. The elliptical inclusion is not box-filling and its conformal map is singular at the foci, so the generic unit-box sampler lands on both; it now supplies rings in the matrix instead. Three harness assumptions surfaced only once every solution went through the same path, which is the point of doing it: - adversarial_points only ever made 2D points, so the 3D solution could not be sampled by it at all; - momentum_residual normalised by the body force, and the inclusion is driven entirely by its boundary and has none -- dividing by zero gave 4e+285. It now scales by the largest term being cancelled, which is the right yardstick for a cancellation anyway; - sample silently cast complex results to real. It now checks the imaginary part is round-off first, because a genuinely complex result would mean the construction is wrong and discarding it would hide that. The first version of that check compared imaginary against total magnitude and fired on residuals, where both parts are round-off and the ratio is meaningless -- an absolute floor comes first. Also found while chasing what looked like a slow test: an orphaned full-suite pytest from earlier in the session had been competing for CPU, which is what made several unrelated runs look pathological. The conformance file itself was slow for a real reason too -- building a Stokes solver per solution dominated it, while checking nothing the contract tests do not already cover -- so that check now lives only in test_1016. Verified: 35 conformance checks; 37 transcription tests; 20 contract tests; style gate clean. Underworld development team with AI support from Claude Code
uw.analytic.SolA and SolB — constant viscosity on the unit box, free slip, forced by (0, sigma * sin(m pi z) cos(n pi x)) and its sinh counterpart. Worth having precisely because they are the simplest. They remove the viscosity structure entirely, so a discrepancy is in the discretisation or the solve rather than in how a hard coefficient is handled: run SolA before concluding anything from SolCx or SolKx. SolB then concentrates the response near one boundary instead of filling the box, which probes resolution where the solution is steep rather than accuracy where it is smooth. Both passed the conformance checks on the first attempt, including the momentum balance — so the forcing conventions inferred from the kernels were right, which after SolKz was not a safe assumption. Their stress is the TOTAL, not the deviator, and SolA is the clearest case in the family to read: its source writes `u3 = 2*kn*ss_z - pp`, with the pressure subtracted in plain sight, where SolKz's writes the same quantity without it. The provenance table now records which solution publishes which, and that SolNL belongs on the deviatoric side. One transcriber addition: sinh, cosh and tanh, which SolB needs and no earlier kernel used. The failure was clean and immediate — a NameError from the generated expression, not a wrong answer — which is the right way for an unsupported function to fail. Note the conformance fixture builds every registered solution together, so one solution failing to construct errors all of them. That is the cost of the shared fixture and it is worth it: the alternative is each solution's checks living somewhere they can be forgotten. Verified: 45 conformance checks over nine solutions plus 20 contract tests, 65 passing after a clean rebuild; style gate clean. uw.analytic.available() now lists EllipticalInclusion, SolA, SolB, SolCx, SolDB2d, SolDB3d, SolKx, SolKz, SolNL. Underworld development team with AI support from Claude Code
uw.analytic.SolM — Stokes flow with a laterally oscillating viscosity, 1 + eta_0(1 + cos(r pi x)), free slip on the unit box. Worth having because its viscosity oscillates rather than jumping (SolCx) or varying monotonically (SolKx, SolKz), and its wavelength is independent of the flow's. It is the one solution here where the coefficient structure and the solution structure can be deliberately mismatched: choose r incommensurate with n and every element sees a different viscosity profile, which tests quadrature more sharply than a smooth gradient does. The kernel's published stress is wrong. It declares its viscosity as (1 + cos(kr x)) eta0 + 1 and then computes stress as 2 (eta - 1) edot -- the constant part is missing. That is a defect in the source, not a transcription slip: the difference from 2 (eta - 1) edot is EXACTLY zero, and using the published stress leaves the momentum residual at 0.21 where deriving it from the kernel's own strain rate and viscosity gives 1.7e-16. Everything else SolM publishes is mutually consistent, so the transcription supplies the strain rate and lets set_fields derive the stress. This is the case for a check that consults no reference. Comparing SolM against its own kernel would have reproduced the error faithfully and reported agreement; only the momentum residual could see it. One transcriber addition: assignment targets may now carry an [index] as well as a struct prefix, because these kernels return results through out.xx = ... or out[0] = ... depending on vintage. Without it `out[0] = ...` matched nothing -- a loud failure rather than a quiet one, but a failure. Also fixed: an over-broad edit had replaced the same block in SolNL, which shares its shape, leaving SolNL deriving its stress under a comment about SolM's viscosity. SolNL publishes a correct stress and a correct strain rate, so it supplies both and the conformance check compares them. Verified: 50 conformance checks over ten solutions plus 20 contract tests, 70 passing after a clean rebuild; style gate clean. Underworld development team with AI support from Claude Code
…entum sees uw.analytic.SolC — isoviscous flow on the unit box driven by a dense column, sigma for x < x_c and zero beyond, free slip everywhere. It pairs with SolCx: SolCx puts a jump in the operator, SolC puts one in the right-hand side. The response is smooth in both, so trouble in either is in how the discontinuity is integrated rather than in the flow itself. This is the first solution here that is a truncated Fourier series rather than a closed form, which the transcriber now supports: CSource.loop_body extracts the mode loop and the caller evaluates it once per mode with the index bound, summing in SymPy. The accumulation itself cannot be read, since it uses += and the sum has to happen symbolically anyway. The body force is the RESOLVED step rather than a sharp one, and that is deliberate. The fields solve the problem with the density the kernel actually summed, so the pair is exact and the residual checks mean what they say. Comparing against a sharp step would report the truncation error as a defect. Raising `modes` sharpens the step and slows evaluation, since the expression carries one term per mode; the residuals are unchanged at 20 and 40 modes, which confirms they are measuring the transcription rather than the truncation. The body force is also MINUS the density. Most kernels in this family negate internally -- they write rho = -sigma*sin*cos and force with +sigma*sin*cos -- but SolC accumulates the density itself. As summed the momentum residual is 1.8; negated it is 1.6e-16. Worth stating as a rule in the subsystem doc because that sign is invisible to everything else: incompressibility was 1.4e-17 and free slip 1.8e-17 either way. Only the momentum balance could see it, and only because it does not consult the solution's own derivation. Verified: 55 conformance checks over eleven solutions; div 1.4e-17, momentum 1.6e-16, free slip 1.8e-17 for SolC itself; style gate clean. Underworld development team with AI support from Claude Code
SolDA and SolH are the last two Velic solutions. Their sources are vendored and the mode-loop machinery SolC needed is proven, but each carries a complication worth knowing before starting: SolDA combines a viscosity jump with a rectangular forcing in the largest kernel of the family, and SolH is 3D with a double mode loop (900 terms at the published default), nested branches selecting the zero modes, six stress components, and a transposed output mapping. Underworld development team with AI support from Claude Code
Probed rather than attempted. Three obstacles, all specific: the loop opens with a chained assignment (del_rhoB = del_rhoA = del_rho) that the statement reader mis-parses and which needs splitting into individual targets; there are two sequential spatial if/else blocks inside the mode loop, each ~790 lines, so every mode contributes a Piecewise and the structure compounds with mode count; and the loop body is an order of magnitude larger than SolC's, so the per-mode expression size has to be measured before a default mode count can be chosen. Vendors solDA.c alongside the sources already staged for SolH. Underworld development team with AI support from Claude Code
…thers uw.analytic.SolDA — a rectangular density anomaly in a fluid whose viscosity jumps at z_c, free slip on the unit box. The most demanding solution here and the only one that combines what the others test separately: a discontinuous forcing (as SolC), a discontinuous viscosity (as SolCx), and a truncated series. The two discontinuities are perpendicular, so a scheme that handles either alone still has to get their interaction right. I had recorded this as too large to attempt, on an estimate. The estimate was wrong and measuring was cheap: one mode on one side is 0.07 s and about four thousand operations per field, which puts twenty modes in the same range as SolKz. Worth remembering — the obstacle I could actually name (chained assignment) turned out to be a ten-line fix, and the one I could only guess at (size) was not an obstacle at all. The transcriber gained chained assignment: `del_rhoB = del_rhoA = del_rho;` assigns to both, but read as one statement its value is `del_rhoA = del_rho`, which is not an expression. Chains are now split innermost-first so each target is bound before the next uses it; verified on double and triple chains. Both of SolDA's `if (z < zc)` blocks branch on the same condition, so each mode is evaluated along one side and then the other and combined into a Piecewise — the SolCx pattern applied per mode. Every convention had to be read from the source and all of them held first time: total stress, minus-the-density forcing as in SolC, and the transposed mapping where u1 is the vertical velocity. It is genuinely expensive: 20 s to build at 8 modes and 47 s at 16, against 2.4 s for SolC at 40, because every mode carries a Piecewise. The default is 8 for that reason, and the docstring says so. The residuals are unchanged between 8 and 16 modes, which confirms they measure the transcription and not the truncation. Verified: div 1.1e-17, momentum 1.5e-15, free slip 2.2e-18; 60 conformance checks over twelve solutions plus 20 contract tests, 80 passing after a clean rebuild; style gate clean. Only SolH now remains untranscribed. Underworld development team with AI support from Claude Code
… estimate uw.analytic.SolH — isoviscous flow in the unit cube driven by a rectangular density block, free slip everywhere. The 3D counterpart of SolC, and the only 3D solution here with a discontinuous forcing. Three-dimensional flow around a compact body is not the 2D problem with an axis added: the return flow can go around the anomaly rather than only over it. I had recorded SolH as expensive and hard, on the strength of the kernel's own warning that it "can become *very* expensive to compute" and a 900-term count. That warning is about a COMPILED kernel, which re-sums every mode at every evaluation point. For a transcription it is backwards: each mode is about ninety operations, the smallest in the family, and the sum is built once. It builds in 1-2 s and validated on the first attempt -- div 3.4e-17, momentum 1.7e-16. That is the second estimate this session that measuring overturned in minutes, after SolDA. The pattern in both: the obstacle I could NAME was cheap to fix, and the one I could only guess at was not an obstacle. Recorded in the subsystem doc. Two transcriber additions, both mechanical once looked at. The C ternary, since SolH guards its zero modes with `(n!=0 || m!=0) ? ... : ...`; parenthesised groups are rewritten first so nested conditionals resolve, verified on both. And resolve_branches, which is the one that matters. These kernels guard their zero modes with tests on the loop indices, and those are bound to integers before anything is evaluated, so the construct collapses to whichever branch the C would take. Left unresolved it is not a crash: evaluate_block reads every assignment in order, so each guarded variable keeps the LAST branch's value, which in SolH silently zeroes two velocity components and leaves a plausible-looking solution. Verified: 65 conformance checks over thirteen solutions plus 20 contract tests, 85 passing after a clean rebuild; style gate clean. Underworld development team with AI support from Claude Code
The Stokes family was already in uw.analytic. This brings across the scalar
solutions that were scattered elsewhere, and closes the plan's remaining items.
Transport (transport.py) — Poisson1D, TwoLayerDarcy, ErfcDiffusion,
AdvectedFront. All four were written inline in the tests that used them, where
nothing checked them against the equations they solve. They declare
solves = "transport" and carry fn_solution / fn_coefficient / fn_source instead
of the velocity-and-pressure pair.
Richards (richards.py) — GardnerSteady and GardnerTransient, the one nonlinear
scalar family, previously NumPy functions in utilities/retention_curves.py.
Those functions keep their signatures and now evaluate the same SymPy
expression the classes build, so there is one formula rather than two copies
that can drift; both reproduce the previous arithmetic to reassociation
(4e-14 and 4e-16).
Kramer (kramer.py) — CylindricalStokes, wrapping the external `assess`, now
declared as the `benchmarks` extra. Four scripts under docs/examples/ imported
it while nothing declared it, so on a normal install they failed with a bare
ModuleNotFoundError.
Three things worth recording, all of them checks rather than solutions:
A residual only means something if it is the residual of the right equation.
AdvectedFront reported 1.44 next to a column of zeros, which reads
unambiguously as a broken solution — but diffusion_residual was testing pure
diffusion against an advecting front. The advection term is now included as the
general case; a purely diffusive solution declares no velocity and it drops
out. The conformance suite no longer skips transients, which is what let this
sit unnoticed.
A residual can be degenerate rather than wrong. richards_residual first
normalised by the flux divergence, which *is* the residual, and reported
exactly 1.00 for a solution that is exact to the last bit. It now normalises by
the terms that have to cancel, kept separately.
Not every perturbation is a negative control. Scaling K by a constant leaves
the steady Richards residual at zero — correctly, since that is a genuine
symmetry of the equation. Three controls that do discriminate are asserted
(wrong alpha in K: 0.64; head scaled 1%: 0.0099; K independent of head: 1.00),
and so is the symmetry, so neither is left as a claim in prose.
Also: lambdify("numpy") has no erfc and falls back to the scalar math.erfc
without complaint, failing much later from generated code. Unseen until now
because differentiating an erfc removes it, and only the Richards head keeps
one inside a logarithm. _validation.sample now asks for ["scipy", "numpy"].
CylindricalStokes is an oracle, not a member of the family: assess is numeric,
so none of the six gates can reach it. Declared as symbolic = False rather than
described, with the conformance sweep excluding on the declaration and then
asserting what it excluded — an accidental exclusion fails the suite instead of
quietly shrinking it.
uw.analytic.available() now lists 20 solutions, 19 of them symbolic and swept.
274 passed across the analytic and MG suites; 126 passed, 6 skipped (assess
absent) on the files touched here.
Underworld development team with AI support from Claude Code…pies test_1000 (Poisson sinusoid), test_1004 (two-layer Darcy), test_1005 (erfc diffusion) and test_1100 (advecting top hat) each carried their own copy of an exact solution that nothing checked against the equation it claimed to solve. They now use the registered solutions, which the conformance suite verifies. Assertions and tolerances are unchanged. Two things came out of the Darcy migration. TwoLayerDarcy needed generalising: test_1004 is posed on y in (-1, 0), not the unit column, and runs the case twice — with and without gravity. It now takes the column extent and a gravity term S, and the profile is *derived* from constant flux q = -k(dp/dz + S) rather than transcribed from the closed form the test carried. That makes the agreement a check on both rather than a copy of one: 1.1e-16 in both gravity cases, residual exactly zero. Its permeability arguments are now k_lower / k_upper rather than k1 / k2. test_1004's k1 is the *upper* layer, and getting that backwards produces a smooth, plausible, wrong answer with nothing to flag it — the names should not leave that available. test_1100's mesh0 case xpasses, as it did intermittently before; its xfail is strict=False and its note says either outcome is acceptable pending a rework. Nothing here was tuned to change that, and it is left alone. 29 passed, 1 xpassed across the four files. Underworld development team with AI support from Claude Code
It was written but never added to the authority map or the toctree, so Sphinx built it as an orphan and nothing pointed at it. docs-build succeeds. Underworld development team with AI support from Claude Code
assess is a 12 kB pure-Python wheel on PyPI with dependencies we already have, so the "working path untested" caveat was avoidable. It is now a dev dependency in pixi.toml — pyproject.toml keeps it as the `benchmarks` extra for users — and all four Kramer cases construct and evaluate. The wrapper's API guesses were right: the four CylindricalStokesSolution* classes, their constructor signatures, and .velocity_cartesian / .pressure_cartesian all match what the example scripts implied. "Returns finite values" is not validation, so the solution is now checked by finite differences — the same idea as Gate 4 with a weaker instrument, which is all a numeric oracle admits: div(u)/|u| ~ 1e-9 in all four cases (the difference floor at h=1e-6) free slip: u.n ~ 1e-17 on both arcs, with |u| ~ 1e-2 there zero slip: |u| ~ 1e-17 on the walls, 1e-5 inside Each carries a control. The free-slip wall is demonstrably slipping, so u.n = 0 is not passing because everything is zero; the divergence probe is checked against u = (x, y) to confirm it reports 2 rather than reporting 0 for everything. Installing assess also silently deleted the coverage that mattered most. The missing-dependency path is what a normal install takes and the thing the previous arrangement got wrong, and it was tested by skipping when the package was present — which, now that pixi supplies it, means never, least of all in CI. Absence is therefore simulated rather than waited for, and the fixture has its own negative control asserting the simulation actually blocks. Without that, a change to import resolution would let those tests pass by importing the real package while appearing to cover a path they never touch. The five example scripts now import through kramer.require_assess() instead of a bare `import assess`, so they report what to install rather than ModuleNotFoundError. That function is public because they call it: they use assess for cases this wrapper does not cover. 143 passed across the conformance, transport, Richards and optional suites; 24 of those are this file, with no skips. docs-build succeeds. Underworld development team with AI support from Claude Code
The example's validation was a shell-out to Underworld2 inside a try/except ImportError. UW2 is not installed, so the branch never ran — and behind it were three separate ways the comparison was wrong: * The body force was written out by hand as +1 on x > x_c. SolC's buoyancy is negative on x < x_c, so the file solved a mirrored, sign-flipped problem from the one it compared against. * The comparison sat at the *end* of the file, but the file runs four solves: SolC, then SolCx, then two boundary-condition experiments. By then `v` held a SolCx solve with penalty BCs, not the SolC answer. * It computed `num = function.evaluate(v.fn, ...)` and then never used it, differencing `v.data` instead. Forcing and viscosity now come from uw.analytic.SolC, so they cannot disagree with the solution about sign or side, and the validation happens immediately after the SolC solve while `v` still holds it. Measured, not asserted: res velocity pressure rate_v rate_p 8 1.323e-03 5.086e-03 16 1.659e-04 1.472e-03 3.00 1.79 32 1.890e-05 2.320e-04 3.13 2.66 64 9.520e-07 3.358e-05 4.31 2.79 Third order for P2 velocity, second for P1 pressure, as expected. Getting the file to run at all turned up three pre-existing defects, all filed: #498 MeshVariable.clone was broken at both levels — EnhancedMeshVariable forwarded no arguments to a base that requires two, and the base itself referenced a bare `MeshVariable`, which is not a name in its module. So it raised whichever way it was called, and no test covered it. Fixed here (it blocks the example at line 141) with tests, including that it still rejects the no-argument form. Six shipped examples were aborting on this line. #499 timing.print_table no longer accepts display_fraction/group_by/ output_file; 18 example files still pass them. Only this file's two call sites are fixed — the rest need a decision about whether to restore the keywords or update the call sites, which is not mine to make here. pl.show() was guarded on `uw.mpi.size == 1` alone, so a script run blocked forever on a window that never opens. Now also requires uw.is_notebook. This is endemic across the examples rather than specific to this file. The header also described SolCx — a 10^6 viscosity contrast with cos/sin buoyancy — while the code was isoviscous with a step force. Corrected, with a note for anyone comparing against old output. Underlying all of it: nothing runs the examples, so API drift lands in them unnoticed. Worth a smoke job at trivial resolution; noted on #499. Example runs to completion, exit 0. 91 passed across the clone and conformance suites. Underworld development team with AI support from Claude Code
…ution floor Every diffusive similarity solution here is a step with unbounded gradient at t = 0 — a state no finite element space can hold. Prose alone was not going to stop anyone benchmarking there, so this is enforced: sol.at(t) refuses t <= 0 with an explanation, rather than returning the singular profile sol.singular_at_origin declared per solution sol.earliest_resolvable_time(h) the floor the *mesh* imposes, not the solution: t0 >= (n_el h / 2)^2 / D, from requiring the front to span n_el elements The floor falls as h^2, so refining buys an earlier start quickly. It is a number rather than advice, which matters because the honest answer depends on resolution and nobody was going to work it out per run. Two things the docstrings now say outright, because neither is guessable: * A transient benchmark is a PAIR of times, never one. You initialise at t0 and compare at t1, and the error depends on both. An error quoted at a single time is uninterpretable — start too early and what you attribute to the timestepper is mostly initial projection error. * Below the floor you are measuring interpolation, not the solver. That immediately diagnoses test_1100_AdvDiffCartesian, which has carried an xfail calling itself "not a great test" and asking for "an error-function IC starting at t > 0 with a meaningful transport distance". At its own parameters (res 24, kappa 1, u 1/24, t0 1e-4, t1 2e-4): earliest resolvable t0 3.5e-3 -> it starts 35x too early front width at its t0 0.68 elements, narrower than one cell transport over the run 4.2e-6 = 0.0001 elements It initialises a profile the mesh cannot represent, then advects it by a ten-thousandth of a cell. It measures neither advection nor diffusion, which is why it has always been sensitive to which path uw.function.evaluate takes. The test is NOT reworked here. A resolution-consistent setup at res 24 still shows 11% error in five steps, dominated by time discretisation, so fixing it needs a timestep convergence study rather than new constants. Recorded in the subsystem doc with the numbers. Also documents the format to ask contributors for, since solutions derived from analytics can be supplied in whatever form we specify: SymPy on mesh.X rather than callables; no simplify(), preserve the derivation's grouping; declare the stress convention; give the equation and not only the answer; declare singularities in time and space; state the valid parameter ranges so they become constructor validation instead of folklore. 133 passed; docs-build succeeds. Underworld development team with AI support from Claude Code
…nto feature/analytic-suite # Conflicts: # src/underworld3/discretisation/enhanced_variables.py
…e the rename Landing work on top of the merged analytic suite. CONVENTION AUDIT. Measured rather than read: every registered Stokes solution already obeys one convention on its exposed fn_* fields -- total Cauchy stress, pressure positive in compression, div(sigma) + f = 0 -- and those are UW3's own, fixed by SNES_Stokes.stress and by F0 = -bodyforce against F1 = stress. The non-uniformity is in the published sources and is absorbed at one declared boundary (stress_is_deviatoric, honoured only inside set_fields). No convention needed changing; what was missing was enforcement. FOUR DEFECTS, all found by the oracle-free residual. 1. SolA is wrong for any viscosity but 1. solA.c:156 computes the zz stress without the factor of Z that its own xx stress carries and that solB.c:140 carries. The shortfall is tau_zz*(1-Z)/Z, identically zero at Z=1 -- the default eta, and the only value the kernel's own disabled driver exercised. At eta=3 the momentum residual is 2.8e-1, the deviator trace 6.7e-1 and the strain-rate consistency 6.7e-1, where |1-3|/3 = 0.667 exactly. Repaired by restoring the factor on the term that lost it, declared per solution; the vendored source stays verbatim. Deriving sigma_zz from sigma_xx via tracelessness was rejected because it would make the deviator traceless BY CONSTRUCTION and retire one of the three gates that caught the defect. 2. EllipticalInclusion ignored matrix_viscosity. The potentials are normalised to unit matrix viscosity; rescaling eta scales the stress AND the pressure, but only fn_viscosity was scaled, so the two parts of sigma were in different units. Momentum residual 6.3e-1 at matrix_viscosity=3, now 3.9e-15. 3. SolNL(r=2) raised KeyError: 'ComplexInfinity'. alpha = 1/r - 1 = -1/2 makes the published pressure's denominator vanish identically for every wavenumber -- a pole of the solution. Now refuses with a ValueError naming the cause. 4. SolKz's deviator/total boundary was described wrongly in our own docs. The C function does convert to the total and says so; what is deviatoric is what the transcription captures, because the transcriber stops at the first mode accumulation, which precedes the conversion. Corrected in the subsystem doc. GUARDS. The uniformity is now enforced rather than observed: - tests/test_1028_analytic_parameter_sweep.py re-applies all three residual gates AWAY from the defaults, with a table every registered Stokes solution must appear in. This is the class of bug that hid #1. - the body-force-sign negative control is asserted, not assumed: flipping the sign must move the momentum residual to order unity. - strainrate_consistency -- the genuinely independent velocity-vs-stress check -- promoted from one solution's file into the family-wide sweep. - test_stress_and_strain_rate_agree is labelled as structural for the ten solutions that publish only one of stress or strain rate, and the list of the three that publish both is asserted against the sources rather than commented. RENAME. Every in-repo caller moved to underworld3.analytic; the deprecating shim at underworld3.function.analytic stays, because the old path is public and external scripts cannot be audited. test_1016 still imports the old path deliberately -- it is the shim's contract test. test_1015's low-level AnalyticSolCx_* imports now go to analytic._reference._velic, where they live: those are the vendored kernel, not part of the uw.analytic surface. Underworld development team with AI support from Claude Code
test_the_list_of_solutions_publishing_both_is_accurate decided whether a
solution publishes both a stress and a strain rate by walking the MRO for a
class whose source contains "set_fields(" and then looking for the parameter
names in it. That is wrong for the one solution it matters for:
EllipticalInclusion never calls set_fields at all, so the walk fell through to
AnalyticSolution and matched `stress=` and `strainrate=` in the base class's own
signature. It reported True and the assertion failed.
set_fields now records what it was actually handed, as
publishes_both_stress_and_strainrate, defaulting to False for any solution that
bypasses it. The test reads that instead of the source text.
This matters beyond the false positive: the flag is what says whether the
conformance check sigma + p I == 2 eta edot is evidence or bookkeeping. Where
only one of the two was supplied, set_fields derived the other from exactly that
identity, so the check re-reads a derivation.
Underworld development team with AI support from Claude CodeContributor
There was a problem hiding this comment.
Pull request overview
This pull request lands the new underworld3.analytic (uw.analytic) namespace as a unified, contract-driven suite of analytic/benchmark solutions (Velic Stokes family plus scalar transport and Richards cases), including a deprecating compatibility shim for the legacy underworld3.function.analytic path. It also updates existing tests/examples to use the new API, vendors/relocates the reference kernels behind underworld3.analytic._reference, and includes a fix plus regression test for MeshVariable.clone.
Changes:
- Introduces
underworld3.analyticwith registry/contract, adds new analytic solution modules (transport, Richards, Kramer wrapper), and relocates the compiled reference-kernel extension tounderworld3.analytic._reference._velic. - Updates tests and example scripts to use
uw.analytic(and the legacy shim where explicitly required), adding new validation-focused test suites for transport and Richards. - Fixes
MeshVariable.clone/EnhancedMeshVariable.clonebehavior and adds a dedicated regression test.
Reviewed changes
Copilot reviewed 67 out of 75 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/test_1100_AdvDiffCartesian.py | Switches handwritten advected-front expression to uw.analytic.AdvectedFront. |
| tests/test_1062_constrained_solcx.py | Updates imports from legacy analytic namespace to underworld3.analytic. |
| tests/test_1026_analytic_richards.py | Adds Richards (Gardner) analytic conformance/negative-control tests. |
| tests/test_1025_analytic_transport.py | Adds transport analytic conformance/negative-control tests (Poisson/Darcy/diffusion/advection-diffusion). |
| tests/test_1023_analytic_solkz.py | Adds SolKz analytic validation tests (viscosity varies with depth). |
| tests/test_1022_analytic_soldb.py | Adds SolDB2d/SolDB3d manufactured-solution tests (symbolic simplifications). |
| tests/test_1021_analytic_solkx.py | Adds SolKx analytic validation tests (smooth viscosity contrast). |
| tests/test_1020_analytic_inclusion.py | Adds EllipticalInclusion validation tests (multiple independent checks). |
| tests/test_1018_rotated_freeslip.py | Updates analytic imports to new namespace. |
| tests/test_1017_custom_mg_stokes.py | Updates analytic imports to new namespace. |
| tests/test_1016_analytic_contract.py | Adds contract + legacy-shim behavior tests for uw.analytic. |
| tests/test_1015_analytic_solcx.py | Updates to uw.analytic and imports low-level reference bindings from _reference._velic. |
| tests/test_1005_TransientDarcyCartesian.py | Switches transient diffusion profile to uw.analytic.ErfcDiffusion. |
| tests/test_1004_DarcyCartesian.py | Switches two-layer Darcy analytic profile to uw.analytic.TwoLayerDarcy. |
| tests/test_1000_poissonCart.py | Switches sinusoidal Poisson source/solution to uw.analytic.Poisson1D. |
| tests/test_0839_nvb_parallel_adapt.py | Updates analytic imports to new namespace. |
| tests/test_0836_nvb_graded_adapt.py | Updates analytic imports to new namespace. |
| tests/test_0835_sbr_adapt_on_top.py | Updates analytic imports to new namespace. |
| tests/test_0301_meshvariable_clone.py | Adds regression tests for MeshVariable.clone behavior (issue #498). |
| tests/parallel/test_1064_rotated_freeslip_parallel.py | Updates analytic imports to new namespace (MPI test). |
| tests/parallel/test_1017_custom_mg_parallel_mpi.py | Updates analytic imports to new namespace (MPI test). |
| src/underworld3/utilities/retention_curves.py | Refactors Gardner retention helpers to delegate to underworld3.analytic.richards. |
| src/underworld3/function/analytic/init.py | Adds deprecation shim package for legacy underworld3.function.analytic. |
| src/underworld3/function/init.py | Adds lazy __getattr__ to keep uw.function.analytic reachable without eager import. |
| src/underworld3/discretisation/enhanced_variables.py | Updates EnhancedMeshVariable.clone to delegate to base clone implementation. |
| src/underworld3/discretisation/discretisation_mesh_variables.py | Fixes MeshVariable.clone NameError by calling the public factory. |
| src/underworld3/analytic/richards.py | Adds Gardner steady/transient Richards solutions and shared helpers. |
| src/underworld3/analytic/kramer.py | Adds Kramer/assess-backed curved-geometry Stokes wrapper solution(s). |
| src/underworld3/analytic/_reference/solH.c | Vendors reference kernel source into new _reference package location. |
| src/underworld3/analytic/_reference/solCx.h | Vendors reference kernel header into new _reference package location. |
| src/underworld3/analytic/_reference/solC.c | Vendors reference kernel source into new _reference package location. |
| src/underworld3/analytic/_reference/solB.c | Vendors reference kernel source into new _reference package location. |
| src/underworld3/analytic/_reference/solA.c | Vendors reference kernel source into new _reference package location. |
| src/underworld3/analytic/_reference/AnalyticSolNL.h | Vendors reference header into new _reference package location. |
| src/underworld3/analytic/_reference/AnalyticSolM.hpp | Vendors reference header into new _reference package location. |
| src/underworld3/analytic/_reference/AnalyticSolDB2d.hpp | Vendors reference header into new _reference package location. |
| src/underworld3/analytic/_reference/AnalyticSolCx.h | Adds UW3 adapter header for SolCx reference kernel in new location. |
| src/underworld3/analytic/_reference/AnalyticSolCx.c | Adds UW3 adapter C implementation for SolCx reference kernel in new location. |
| src/underworld3/analytic/_reference/init.py | Documents _reference purpose (vendored kernels; not public API). |
| src/underworld3/analytic/init.py | Adds public uw.analytic namespace exports plus registry (available, describe, is_available). |
| src/underworld3/init.py | Imports underworld3.analytic into the top-level package namespace. |
| setup.py | Renames compiled extension target to underworld3.analytic._reference._velic and adds package data for reference headers/sources. |
| pyproject.toml | Adds optional dependency group benchmarks = ["assess"]. |
| pixi.toml | Adds assess to dev dependencies for running benchmark-related tests/examples. |
| docs/examples/utilities/advanced/Ex_Stokes_Annulus_Benchmark_Thieulot.py | Imports assess via require_assess() to provide a clearer missing-dependency message. |
| docs/examples/utilities/advanced/Ex_Stokes_Annulus_Benchmark_Kramer_etal.py | Imports assess via require_assess() to provide a clearer missing-dependency message. |
| docs/examples/fluid_mechanics/advanced/Ex_Stokes_Spherical_Benchmark_Thieulot.py | Imports assess via require_assess() to provide a clearer missing-dependency message. |
| docs/examples/fluid_mechanics/advanced/Ex_Stokes_Spherical_Benchmark_Kramer.py | Imports assess via require_assess() to provide a clearer missing-dependency message. |
| docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_SolNL.py | Switches SolNL example to use uw.analytic.SolNL instead of low-level bindings. |
| docs/examples/fluid_mechanics/advanced/Ex_Stokes_Cartesian_SolC.py | Reworks SolC example to use uw.analytic.SolC, fixes validation placement, updates narrative. |
| docs/examples/fluid_mechanics/advanced/Ex_Stokes_Annulus_Benchmark_Kramer.py | Imports assess via require_assess() to provide a clearer missing-dependency message. |
| docs/examples/snesfas_investigation/benchmark_3way.py | Updates analytic imports to new namespace. |
| docs/developer/index.md | Adds analytic-solutions subsystem doc to the developer index/toctree. |
| docs/api/index.md | Adds analytic to API docs index. |
| docs/api/function.md | Updates function docs to point analytic solutions at the new analytic docs page. |
| docs/api/analytic.md | Adds API docs page for underworld3.analytic and its contract. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
+419
to
+424
| coordinate = sympy.Symbol("y") | ||
| saturation, _ = gardner_steady_saturation(coordinate, psi_0, psi_L, L, alpha) | ||
| # Normalised steady-state flux q* = q / Ks | ||
| q_star = (u_L - u_0 * np.exp(-alpha * L)) / (1.0 - np.exp(-alpha * L)) | ||
| u = sympy.lambdify(coordinate, saturation, ["scipy", "numpy"])( | ||
| np.asarray(y, dtype=float) | ||
| ) |
Comment on lines
+259
to
+262
| if self.boundary == "zero": | ||
| for wall in self.boundaries: | ||
| solver.add_dirichlet_bc((0.0, 0.0), wall) | ||
| return |
PR #571 did not fail, it TIMED OUT: CI cancelled at 1h00m22s against a 60-minute cap, with zero FAILED lines and the analytic files at 81% and passing. The suite costs 26.6 minutes in CI, and the run without it was already 55m07s. PROFILED FIRST. The cost is not solving -- there is exactly one Stokes solve in the whole suite. It is symbolic: every residual gate differentiates the solution's expressions and runs common-subexpression elimination over the result, and five solutions produce expressions with tens of thousands of operations. Measured over 1010s: test_1028 parameter sweep 416s SolDA 187s SolCx 14s test_1023 SolKz 189s SolKz 187s SolB 2s test_1024 conformance 169s SolKx 72s SolNL 0.7s test_1019 transcription 100s SolH 52s SolDB3d 0.6s test_1021 SolKx 88s SolC 35s SolA 0.5s everything else 48s Elliptical 34s SolM 0.2s The five expensive ones are 565s; the other eight together are about 17s. The solutions that historically caught defects -- SolA, SolM, SolNL, SolDB2d/3d -- are all in the cheap group, which is what makes the split affordable. CHEAPER WITHOUT LOSING ANYTHING. `momentum_residual` sampled the SYMBOLIC sum of its terms and then sampled each term again for the scale. The sum is the expensive one -- CSE over several nearly-cancelling series expressions -- and it was redundant: the terms are sampled anyway, so the residual can be summed NUMERICALLY. Same numbers to round-off, eight orders of margin against the 1e-8 gate, negative control still 2.000. Applied to the momentum, incompressibility, transport, diffusion and strain-rate gates. The two single-solution files also now cache their constructions instead of rebuilding per test (SolKz costs 6s to build, SolDA 12s). THE SPLIT, BY FILE PLACEMENT. scripts/test.sh batches by file GLOB, not by marker, so a `slow` marker alone would need every batch line to remember to deselect it. tests/analytic_full/ is a subdirectory, and the globs do not recurse, so it is excluded by construction; the level_2 mark additionally keeps it out of `pytest -m "level_1 and tier_a" tests/`, which does recurse. Verified all three collection paths. Solutions declare their own side of it (`expensive_to_validate`), so the tiers partition the family from one source of truth. Two guards: test_1024 asserts every solution it skips is NAMED in the full-family file, and that file asserts its hand-written list matches the declarations. NO GUARD WEAKENED. The momentum residual, incompressibility, tracelessness, strain-rate consistency and the body-force negative control all run in BOTH tiers on every solution that tier covers. The reduction is solutions per run, never checks per solution. per-PR analytic 16m50s -> 4m07s 307 passed CI batch 101*/102* 20m30s -> 7m30s 446 passed full family (opt-in) 8m34s 189 passed Underworld development team with AI support from Claude Code
lmoresi
commented
Aug 15, 2026
MemberAuthor
The CI result was a 60-minute job timeout, not a failure — zero FAILED lines, analytic files at 81% and passing. Measured from the log, the suite cost 26.6 minutes against a run (#568, same day, without it) that already took 55m07s. Fixed in The profile overturned the obvious hypothesis. There is exactly one Stokes solve in the entire suite; the cost is symbolic — every residual gate differentiates the solution's expressions and runs CSE over the result. Five solutions are 565s of 1010s (SolDA 187, SolKz 187, SolKx 72, SolH 52, SolC 35); the other eight total ~17s. SolC sums forty modes; SolKx/SolKz carry exponential viscosities. A real redundancy, removed: Split by file placement, not marker — Solutions declare their own tier ( No guard weakened. Momentum residual, incompressibility, tracelessness,
Scaling by the observed CI/local factor (26.6 min CI vs 16m50s local, ≈1.6×), the suite's CI cost should land near 6.5 minutes. Running both tiers now costs less than the old single tier did. Full family, one documented command (in the module docstring and pixi run -e amr-dev python -m pytest tests/analytic_full/ -vUnderworld development team with AI support from Claude Code |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for freeto join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Land the analytic solution suite: the Velic family as
uw.analytic, with an oracle-free validation contractBranched from
development, merging the long-strandedbugfix/analytic-module-rename.That branch name is misleading — it long outgrew a rename. This is the suite
itself, plus a convention audit and four defect fixes found while landing it.
What this adds
A single namespace,
uw.analytic, holding nineteen exact solutions behind onecontract:
SolA,SolB,SolC,SolCx,SolDA,SolDB2d,SolDB3d,SolH,SolKx,SolKz,SolM,SolNL.transcribed — the authors publish pressure, deviatoric stress and rotation
rate but no velocity, so the Muskhelishvili potentials were recovered from the
fields and the velocity built from those.
Poisson1D,ErfcDiffusion,AdvectedFront,TwoLayerDarcy,GardnerSteady,GardnerTransient.Every one reads the same way:
available()anddescribe()list the family without constructing it.Published kernels are vendored under
analytic/_reference/and transcribedinto SymPy rather than compiled, so the solutions are differentiable,
JIT-compilable, and usable as solver input rather than only as point oracles.
SolKzis the immediate motivation: a smooth viscosity variation with aknown solution, which is what makes a resolution-convergence claim checkable.
developmentcurrently has onlySolCx, whose viscosity is a step.Validation posture — read this part
Most of the 11k lines are mathematics transcribed from published kernels, where
a sign error is invisible to a build and to any test that compares against the
same transcription. So the question that matters is not "do the tests pass" but
"can the tests fail".
The oracle is a residual that consults nothing.
$\nabla\cdot\sigma + \mathbf f$ by symbolically differentiating the solution's
momentum_residualformsown Cauchy stress and adding its own body force. It never touches the vendored
kernel and never touches the solver. The scalar families have counterparts
(
transport_residual,diffusion_residual,richards_residual), andincompressibility_residualandstrainrate_consistencycomplete the set.This is not tautological, and the design forbids it becoming so:
AnalyticSolution.set_fieldstakes the body force as an independentlytranscribed input, and refuses a solution that supplies neither a stress nor
a strain rate — the one configuration in which the residual would be true by
construction.
The negative control fires. Flipping the sign of the body force moves the$10^{-16}$ to between 0.5 and 2.0 for every solution that
momentum residual from
has one. That is now asserted, not merely observed
(
test_flipping_the_body_force_breaks_the_momentum_balance).EllipticalInclusionis excluded by name because it is boundary-driven and has no body force to flip;
that exclusion is itself asserted.
Per-solution strength. All nineteen registered solutions are swept by
test_1024_analytic_conformance.py, which iteratesavailable()so a solutionadded later is covered the moment it is registered. Two exclusions exist and
both are declared by the solution, not listed by name:
symbolic = False(anexternal numeric oracle with nothing to differentiate —
CylindricalStokes) andrequires = "..."for an uninstalled optional dependency. The sweep assertsthat every exclusion has one of those grounds, so a typo in a class attribute
fails the test rather than quietly shrinking coverage.
Where the validation is weaker, stated plainly.
reference_agreementcompares a transcription against the kernel it came from. That tests the
transcription and not the mathematics: if the kernel is wrong, agreement is the
wrong answer. SolM is the worked example — see the errata.
One gate was partly true by construction, and is now labelled.
$\sigma + p\mathbf I = 2\eta\dot\varepsilon$ .$\tfrac12(\nabla\mathbf u + \nabla\mathbf u^T)$ against the strain rate
test_stress_and_strain_rate_agreeassertsFor a solution that supplied only one of stress or strain rate, that is exactly
the identity
set_fieldsused to derive the other — structural, not evidential.Only
SolNL,SolDB2dandSolDB3dsupply both. The gate is kept (it wouldcatch
set_fieldsregressing) but the docstring now says what it is worth, thelist of solutions publishing both is asserted against the sources rather than
commented, and the genuinely independent check —
strainrate_consistency, whichcompares
derived from the stress, two separate kernel outputs — has been promoted from a
single solution's file into the family-wide sweep.
The convention audit
Louis asked for the sign and stress conventions to be unified before this lands,
on the timing argument that nothing depends on them yet. The audit's finding is
that the exposed conventions are already uniform — measured, not asserted — and
that the non-uniformity is entirely in the published sources, where it is
absorbed at a single declared boundary. The full table, the per-kernel
provenance, and the errata are in
analytic_errata.md, written as prose for thetechnical note. Summary:
UW3's own conventions win, and the suite already matches them:
SNES_Stokes.stressset_fieldsbuilds exactly thisF0 = -bodyforce,F1 = stressmomentum_residualforms precisely thisboundary_normal_traction(post-#561)The sources are not uniform: four of thirteen publish the deviator, the rest
the total; component order varies (
solKx.cuses[xx, xz, zz]where everysibling uses
[xx, zz, xz]); and several label the vertical velocityu1. Allof that is absorbed by the
stress_is_deviatoricdeclaration and an explicitcomponent mapping, honoured in one place. Two things are uniform across every
source and neither is documented in most of them — both had to be measured by
finite-differencing the kernels' own outputs: the momentum sign, and pressure
positive in compression.
Errata (each recorded with its citation, not silently conformed)
solA.c:156) —new, and it was live in this suite. See below.
AnalyticSolM.hpp) — viscosity isPublished stress → residual 0.2117; derived from the kernel's own strain
rate → 1.7e-16. Re-verified independently while landing this: the
published deviator differs from
from
defect in the source rather than a slip in transcription.
solC.c) — as summed 1.848,negated 2.5e-16.
own. The subsystem doc said
solKz.cwrites its deviator into an array namedtotal_stress. As a claim about the C function that is wrong: it convertsexplicitly (
u6 -= u5; /* get total stress */) and says so. What is true is aclaim about the transcription's cut point — the transcriber stops at the
first mode accumulation, which precedes the conversion, so what we capture is
the deviator.
stress_is_deviatoric = Trueis correct but describes the cut,not the kernel. Corrected in the errata.
Four defects fixed, all found by the oracle
Three were found by a single question: do the gates still hold away from the
default parameters? The conformance sweep builds every solution from a mesh
alone, which is right — the defaults are part of the interface — but it means a
coefficient that is unity by default multiplies a term nothing ever looks at.
1. SolA is wrong for any viscosity other than 1 (the significant one)
solA.c:156computes theu3 = 2.0*kn*ss_z - pp. The matchingline in
solB.c:140reads2.0*Z*kn*ss_z - pp, and solA's ownthe next line but one carries the
Z. The shortfall isidentically zero at
exercised, and the default
etain our transcription.etaThree independent gates, all silent at the default.
Repaired by restoring the missing factor on the term that lost it, declared per$\sigma_{zz}$ from $\sigma_{xx}$ via tracelessness would also work and$Z=3$ , solB's is —
solution (
_zz_stress_lost_the_viscosity), leaving the vendored source verbatim.Deriving
was rejected: it would make the deviator traceless by construction and
silently retire one of the three gates that caught the defect. The test asserts
both halves — solA's published deviator is not traceless at
so the control says this is a correction to a defect, not a convention imposed
on the family.
2.
EllipticalInclusionignoresmatrix_viscosityAt
$\eta\to\lambda\eta$ at fixed boundary velocity the stress and the pressure$\lambda$ . $\sigma$ were in different units — invisible to any gate that
matrix_viscosity = 3.0the momentum residual is 0.632 while tracelessnessand strain-rate consistency stay at machine precision — the combination that
localises it. The potentials are normalised to unit matrix viscosity; under
both scale by
fn_viscositywas scaled andfn_pressurewas not, sothe two parts of
looks at only one of them. Now 3.9e-15.
Worth noting why it escaped: this is the only solution that assigns
fn_stress,fn_pressureandfn_strainratedirectly instead of going throughset_fields, which is the one place the stress/pressure relationship is applied.3.
SolNL(r=2.0)raisedKeyError: 'ComplexInfinity'denominator (
AnalyticSolNL.c:52) vanishes identically — for everywavenumber, so it is a pole of the solution rather than a bad
It now raises a
ValueErrornaming the cause; previously it producedComplexInfinity in three fields and surfaced as a
KeyErrorfrom SymPy's codeprinter, a long way from the cause.
4.
MeshVariable.cloneis still aNameErrorondevelopmentUnrelated to the analytic suite, but it is in the diff and it is the one merge
conflict.
developmentfixed#498 via the quickfix batch (#532) withnewMeshVariable = MeshVariable(...)insidediscretisation_mesh_variables.py— where the only class defined is_BaseMeshVariableandMeshVariableis not a name. The branch's versionuses
uw.discretisation.MeshVariable(...), which works; the merge took it. TheEnhancedMeshVariablehalf now delegates to the base rather than duplicatingthe constructor call, since duplicating it is what let the two halves of #498
drift apart.
The rename, and the shim
underworld3.function.analytic→underworld3.analytic. The suite outgrewfunction: it now carries boundary conditions, error norms and a registry, noneof which belong under "symbolic function evaluation".
Decision: move every caller AND keep a deprecating shim. The old path is a
public namespace and users' scripts live outside this repo, so a hard break is
not something we can audit. The shim:
__getattr__, so theDeprecationWarningfires when aname is used, not when the module is imported — ten test files import it at
collection time and a warning there is noise the reader cannot act on;
isinstanceand pickling workacross both paths;
function.analyticwas a compiled extension, and an orphaned.sofrom anearlier install would otherwise be imported in preference to the shim. Python
checks for a package directory before an extension of the same name, so the
directory always wins.
test_1016_analytic_contract.pytests the shim, including that orphaned-.socase, and is therefore the one file that still imports the old path deliberately.
Callers moved:
tests/test_1015,test_1017,test_1018,test_0835,test_0836,test_0839,test_1062,tests/parallel/test_1017,tests/parallel/test_1064,docs/examples/snesfas_investigation/benchmark_3way.py,and
docs/examples/.../Ex_Stokes_Cartesian_SolNL.py(rewritten to useuw.analytic.SolNLrather than three low-level kernel bindings). Docs atdocs/api/function.md,docs/api/analytic.mdanddocs/api/index.mdmatch.One nuance:
test_1015uses the low-levelAnalyticSolCx_*bindings, which arethe vendored kernel and deliberately not part of the
uw.analyticpublicsurface — that namespace exposes solutions, not the functions they were
transcribed from. Those imports now go to
underworld3.analytic._reference._velic, where they actually live, rather thanthrough the shim.
Drift: 152 commits
The merge itself was one conflict (above). Everything else was checked by
running, not by reading:
no impact, and checked specifically.
CylindricalStokesis the only solutionon a curved boundary and applies
add_rotated_freeslip_bc(0.0, wall); thedatum is zero and
boundary_normal_tractionor
dynamic_topographyback, so no error measure consumes the Rotated free-slip: weight the nodal normal by the facet measure the assembly integrates over (#560) #561 sign.The one place the suite states a traction sign is
SolCx.topography_top,returning
— the same sign as UW3's
Consistent. One difference recorded:
boundary_normal_tractionreturns itsresult mean-removed and
topography_topdoes not, so a caller comparing themmust remove the mean itself (
test_1018does).options), Reuse the rotated free-slip solver workspace across solves — reworked with a constants-aware key (#417, supersedes #418) #543 (rotated workspace cache) — no breakage.
Nothing in the suite needed fixing for the drift. The defects fixed here are
older than the drift and were found by the audit, not caused by it.
Test evidence
All runs in the
analytic-suiteworktree against its own pixi environment(
amr-dev), sequentially, on a machine also running other sessions.The suite's own tests —
test_101[5-9]_analytic*andtest_102[0-8]_analytic*:The one failure was a guard added in this PR, and it was a genuine false
positive in the guard rather than in the code — see "Record what set_fields was
given" below. After the fix, the two affected tests pass (
14 passed, 84 deselected in 65.07s). Every other test in the suite passed on the first runafter the merge, before any of the fixes here: 270 passed across the eleven
pre-existing analytic files, which is the answer to "did a fortnight of API
change break it" — it did not.
The tier-A gate —
pytest -m "level_1 and tier_a" tests/:All 26 skips are
need --with-mpi option to runplus two self-documentingpoint-locator skips; the xfail is the tracked 1-manifold evaluate limitation.
Nothing related to this change.
scripts/test.sh— the full serial suite, end to end, because the tier-Agate and what CI actually runs are not the same set:
The
test_101*/test_102*batch is the one this PR changes most, and it is alsothe batch that previously matched no CI glob (fixed under issue #504). It now
carries the whole analytic suite plus the new parameter sweep, at 20:30.
scripts/test.sh --parallel-only --p 2— run because this PR edits two filesunder
tests/parallel/:The xfail/xpass are the pre-existing #564 free-slip partition-dependence family
(
test_1063,test_1064,test_1066), each carrying its own long explanationin the test file recording that it is not caused by #560/#561 and that it passes
locally on macOS/arm64 (hence
strict=False). Unrelated to this PR; noted hereso the counts are not read as new.
The off-default parameter sweep, run directly as a probe before and after
the fixes — this is the measurement that found three of the four defects:
eta=3.0eta=0.25eta=3, n=2, m=1.5matrix_viscosity=3.0r=2.0KeyError: 'ComplexInfinity'ValueErrornaming the pole38 other configurations across 13 solutions passed both before and after, so the
fixes are localised and did not perturb anything else.
One thing worth flagging separately
pytest -m "level_1 and tier_a"with no path argument — the exact command inCLAUDE.mdunder "Quick validation" — dies silently: exit 1, zero bytes ofoutput, no traceback. Adding
tests/fixes it (pytest otherwise collects fromthe repo root, including
.pixi/andpetsc-custom/). This is pre-existing andunrelated to this PR, but it is the "silent pytest death" signature and it means
the documented quick-validation command has been reporting failure without
saying why. Worth its own issue.
Not done
CylindricalStokesremainssymbolic = False— it is an external numericoracle, so the residual gates cannot apply to it. That is declared by the
class and asserted by the sweep, not silently skipped.
resolutions. Convergence-rate testing is a separate concern.
docs/docstrings/review_queue.mdstill lists the removedsrc/underworld3/function/analytic.pyx. It is generated byscripts/docstring_sweep.pyand should be regenerated rather than hand-edited.Underworld development team with AI support from Claude Code