updated checkpoint_xdmf method save discontinuous fields properly by gthyagi · Pull Request #1 · underworldcode/underworld3

gthyagi · 2025-06-30T04:35:28Z

I have updated checkpoint_xdmf method to save discontinuous fields properly.

Copilot

Pull Request Overview

This PR updates the checkpoint_xdmf method to properly save discontinuous fields. Key changes include introducing a helper function to retrieve cell field sizes from HDF5 files, adjusting the file naming scheme for mesh variables, and dynamically determining the data centering and dimensions based on the field’s continuity.

Comments suppressed due to low confidence (3)

src/underworld3/discretisation.py:2897

[nitpick] Consider renaming the helper function to more clearly indicate that it returns the size of cell fields. For example, 'fetch_cell_field_size' could better express that a file read operation is performed.

 def get_cell_field_size(h5_filename, mesh_var):

src/underworld3/discretisation.py:2904

Verify that the updated file naming pattern (including the leading dot) is consistent with other parts of the codebase and intended file structure.

 var_filename = filename + f".mesh.{var.clean_name}.{index:05}.h5"

src/underworld3/discretisation.py:2927

[nitpick] The variable name 'numItems' could be more descriptive; consider renaming it to something like 'field_size' to clarify that it represents the number of data elements being used.

 Dimensions="1 {numItems} {var.num_components}"

Updated headline status from "Phase D delivers BDF-class fault mechanics" to honest scope: Phase B ships for VE / mild VEP; Phase D and Phase E investigations land on the branch as documented experiments but are NOT shipping features for deep-yield VEP. TL;DR rewritten to call out the three regime distinctions and the structural physics mismatch that ETD-as-designed cannot close without rebuilding from the ground up. API section split into "production" (Phase B integrator='etd' on existing classes) and "experimental — DO NOT USE for production" (Phase D split class, Phase E hybrid mode). Added Phase E section under Implementation phasing (DONE), with killer-test numbers and the slip-drift diagnosis. Renamed "generic TimeIntegrator refactor" to Phase F (deferred). New lessons: 8. The diagnostic that mattered: |σ_∥|, not |σ_xy| 9. The structural BDF-vs-ETD slip-rate difference is physics 10. Phase D recommendations checklist (rec #1, #4, #5 implemented; #2, #3 not; raw-E rejected; back-derived η rejected; Min cap rejected) 11. Phase E hybrid drifts because of shared σ* history 12. Investigation closed — ETD for VE / mild VEP, BDF for deep yield Underworld development team with AI support from Claude Code (https://claude.com/claude-code)

- _rss_mb(): switch from resource.ru_maxrss (peak/high-water RSS) to psutil.Process().memory_info().rss (current RSS), with /proc/self/statm Linux fallback and ru_maxrss kept only as a last resort. Current RSS is preferred so freed memory shows as a negative delta. (Copilot #1, #6) - diff(): apply 0.01 MiB threshold before including rss_mb in the delta. Smaller noise prints as "+0.00 MiB" and defeats the "no change" fast path. (Copilot #2) - dump_petsc_leaks_at_finalize(): drop leading "-" from PETSc option keys (codebase convention is no dash; "-malloc_dump" was a no-op), and align docstring with what the function actually sets. (Copilot #3) - Remove unused `os` import from test module. (Copilot #4) - Soften ckdtree.pyx comment claiming Cython "guarantees deterministic destruction" — that's CPython refcounting, which can lag for objects trapped in reference cycles until cyclic GC runs. (Copilot #5) - Test: assert KDTree count deltas relative to immediate before/after, never to absolute baselines. Earlier tests in the same pytest session can leave KDTree refs alive that the cyclic GC may collect at any time, shifting the absolute count. CI surfaced this with "assert 332 == 338" — six trees from earlier tests collected during our gc.collect(). Underworld development team with AI support from Claude Code

…an oracle-free validation contract, and four fixed defects (#571) * Add uw.analytic: one namespace and one contract for exact solutions 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 Code * Move the analytic extension out of uw.function, behind a shim underworld3.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 * Transcriber for the Maple-generated kernels, and SolCx under validation 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 * SolCx as validated SymPy: the transcription replaces the compiled kernel 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 * Validation harness: the checks a transcription must pass, as reusable 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 transcribed, and the reader bug a second kernel exposed 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 * Schmid & Podladchikov inclusion: potentials derived and verified, not 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 * Schmid & Podladchikov elliptical inclusion, validated 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 * SolKx: exponentially varying viscosity, validated by the equations 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 * SolDB2d and SolDB3d: polynomial manufactured solutions, and the first 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 * SolKz: depth-dependent viscosity, and the stress convention that is not 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 * Unify the solution interface, and a conformance suite that covers all 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 * SolA and SolB: the two isoviscous solutions 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 * SolM, and a published stress that is wrong 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 * SolC: the first truncated-series solution, and the sign that only momentum 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 * docs: record what remains untranscribed, and the specific obstacles 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 * docs: record what a SolDA transcription runs into 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 * SolDA: the hardest of the family, and the one that combines all the others 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 * SolH completes the Velic family: 3D, a double series, and a corrected 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 * Transport, Richards and the optional Kramer wrapper complete the suite 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 * Switch the four scalar tests onto uw.analytic, removing the inline copies 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 * Link the analytic subsystem doc into the developer index 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 * Install assess and actually validate CylindricalStokes 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 * Make the SolC example a real benchmark; fix MeshVariable.clone (#498) 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 * Make the t=0 singularity refuse rather than warn, and give it a resolution 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 * Analytic suite: unify the convention audit, fix four defects, complete 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 * Record what set_fields was given, rather than scraping the source for it 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 Code * Split the analytic suite into a per-PR tier and a full-family tier 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

gthyagi added 2 commits June 27, 2025 11:18

updated checkpoint_xdmf method save discontinuous fields properly
17a12bf

added var as argument to def
aa65291

CopilotAI review requested due to automatic review settings June 30, 2025 04:35

CopilotAI reviewed Jun 30, 2025
View reviewed changes

julesghub approved these changes Jun 30, 2025
View reviewed changes

julesghub merged commit faed46c into underworldcode:developmentJun 30, 2025

gthyagi deleted the fix_checkpoint_xdmf branch April 8, 2026 01:47

lmoresi mentioned this pull request Jul 5, 2026
chore(wave-a): quality-campaign deletions & dead code (D1, D7, D13, D16 + loose-ends verdicts) #325
Merged

Navigation Menu

Uh oh!

updated checkpoint_xdmf method save discontinuous fields properly - #1

updated checkpoint_xdmf method save discontinuous fields properly#1
julesghub merged 2 commits into
underworldcode:developmentfrom
gthyagi:fix_checkpoint_xdmf

gthyagi commented Jun 30, 2025

Uh oh!

CopilotAI left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Uh oh!

Conversation

gthyagi commented Jun 30, 2025

Uh oh!

CopilotAI left a comment

Choose a reason for hiding this comment

Pull Request Overview

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants