Free-surface hardening: material-boundary gauge, moving-mesh transport fix, full-density (Boussinesq) formulation#435
Merged
Conversation
…ace integrator
Package the validated free-slip + prognostic dynamic-topography recipe as a
single user-facing object, uw.systems.FreeSurface. The user configures one
stress-free-top Stokes solve (constitutive model, driving body force, wall BCs,
options); FreeSurface reads it and owns the rest:
- derives the held rotated free-slip lid (constraint reaction sigma_nn -> h_inf
via dynamic_topography) and the consistent solve (prescribed surface-normal
velocity, so the surface stays a material boundary);
- runs the three-number L-stable exponential surface step;
- transports a composition (mesh level set / temperature) with enclosed-area
conservation;
- carries the surface increment inward with a Laplacian diffuser and deforms.
Surface nodes come from the DMPlex boundary stratum (label-based, parallel-correct)
rather than a |y-H| radial heuristic; estimate_dt adds a surface-motion cap on top
of the advective CFL; continuous pressure is checked (CBF recovery requirement).
Reproduces the reference driver (crameri_study Case-2, res48) to < 0.2% on surface
amplitude over two steps. Serial and np2 regression tests included.
v1 scope (documented in code): Cartesian vertical topography, and linear /
field-dependent viscosity (a strain-rate-dependent law is not yet supported on the
derived solves). Tangent-advection is deferred. These are the next work items.
Underworld development team with AI support from Claude Code
The derived held/consistent solves inherited two linear-only assumptions from the
Crameri driver; both are FreeSurface-local (the nonlinear solver subsystem is
unchanged and correct):
1. snes_type was hardcoded to "ksponly" on the derived solves — a single linear
solve. Now mirror the free solve's nonlinear configuration (snes_type +
consistent_jacobian): a nonlinear free solve gives nonlinear derived solves, a
linear one keeps its ksponly one-shot. (The rotated held solve additionally
auto-detects nonlinearity and dispatches its own linear/Newton path.)
2. The constitutive-model copy carried the free solve's velocity and velocity-
gradient symbols into the derived solves, so a strain-rate-dependent viscosity
computed itself from the free solution. This is a solver-copy / model-retarget
problem: a constitutive expression is tagged with its source solver's variable
identity. _copy_constitutive_model now unwraps the copied viscosity and rebinds
the source u and Unknowns.L atoms onto each derived solve's own (extracted as
_velocity_rebind_map; a no-op for linear / field-dependent viscosity). The derived
held solve is now bit-identical to an independently built self-consistent solve.
The retarget pitfall (and its adjoint-operator relevance) is documented in
docs/developer/subsystems/constitutive-models.md ("Retargeting a model to another
solver") so it is not re-hit.
Regression test (test_1070): a power-law free surface whose FreeSurface held solve
must equal a hand-built held solve using its own velocity — fails without the rebind,
passes with it. Linear Crameri result and np2 parallel test unchanged.
Underworld development team with AI support from Claude Code
…lter
Additions toward moving-plate / curved-geometry free surfaces (all default-off, so
the validated Cartesian path is unchanged):
- tangent_advect ("shape" | "fields"): first-pass along-surface semi-Lagrangian
transport of the surface fields — barycentric (1-D, bounded) interpolation at the
tangential launch point s - v_t*dt, applied before the exponential relaxation. The
SL-transport + exp-relaxation composition is the exact integrating factor along the
characteristic; a standalone check reproduces the exact translating-decaying cosine
solution to ~4% (SL diffusion). On the flat Cartesian benchmark the correction is
~0.1%/step (negligible); it becomes first-order on the annulus (rotation), which is
where it will be validated and where the parametrisation swaps x -> theta + periodic.
- surface_mask: h_dot -> h_dot * mask, pins the surface rate near a driven wall so a
stress-free top does not fight an imposed wall velocity (the outflow-corner spike).
- surface_filter: Taubin low-pass of the recovered equilibrium topography over the
extracted surface graph (parallel-safe), attenuating the facet-scale CBF boundary
node-noise that appears under strong forcing. Closes a v1 gap for actively-forced
free surfaces.
- _apply_walls now replays a NON-zero essential wall velocity (a plate push) verbatim
(mutable-Matrix wrap for add_essential_bc), not just no-slip.
Underworld development team with AI support from Claude Code
…gent transport FreeSurface gains annulus/cylindrical support and reads the surface rate ḣ (and the consistent-solve normal) along the DEFORM direction — radial on an annulus, vertical in Cartesian — rather than the FE facet normal of the deformed surface. On a deforming/rotating surface the tilted facet normal projects the large tangential velocity into ḣ, an unstable explicit advection that grows the surface out of round-off (measured: a no-load rotating surface blew up 0 -> 7 km over 16 steps; the deform- direction normal holds it flat at ~7e-3 km). The along-surface (tangential) motion is carried once, by the semi-Lagrangian transport (shape/fields modes, optional diffusion- free spectral launch-point interpolation), not double-counted through the normal read. Adds the surface-geometry helpers (_surface_height, _normal_direction, _along_surface, _launch_interp, _fourier_interp) and a Taubin surface-filter hook. Small-slope / near- radial-deformation assumption is documented on _surface_normal_velocity. Validated through the annulus free-surface + rigid-rotation (TPW) study. Underworld development team with AI support from Claude Code
… free-slip
The rotated strong free-slip constraint imposed u.n = 0 only; a non-zero normal datum
raised NotImplementedError, so the free-surface "consistent" solve fell back to an ill-
conditioned penalty natural BC (the term that tripped newtonls line search on an
otherwise linear solve).
Extend the LINEAR rotated solve to a prescribed normal datum:
- build_rotation takes datum_specs {boundary: value | sympy expr | field.sym} and
returns a signed datum_map, accounting for the SVD sign of Vt[0] = ±n̂;
- solve_rotated_freeslip reads solver._rotated_freeslip_datum, imposes it via
zeroRowsColumns(x, b), and restores it on the solution after the iterative path's
v_n=0 cleanup;
- a zero datum is dropped from the map so it takes the untouched free-slip RHS path
(matches free-slip to iterative round-off); the nonlinear path refuses a datum loudly.
One rotated matrix, a switchable datum: 0 -> free-slip / held, ũ_n -> prescribed radial
(the material-surface "consistent" velocity), differing only in the constraint RHS.
Tests (tests/test_1018): datum-zero-matches-freeslip (round-off) and prescribed u.n =
cos(theta) imposed to machine precision at the surface nodes.
Remaining (tracked separately): wire the FreeSurface consistent solve to use this and
drop the penalty; add a datum argument to add_rotated_freeslip_bc (pyx); np2/np4;
nonlinear path.
Underworld development team with AI support from Claude Code
… the penalty (#403) The consistent (material-surface) solve now prescribes u.n = ũ_n through the STRONG rotated free-slip constraint (datum = the un_target field) instead of the penalty natural BC. Two wins: - machine-precision u.n = ũ_n (no penalty leak), along the same rotated per-node normal the rate was measured along; - the solve routes through the rotated LINEAR solver (direct KSP), so a linear (isoviscous) flow no longer thrashes newtonls the way the penalty did — the annulus FS-convection driver runs clean under the DEFAULT newtonls (0 DIVERGED_LINE_SEARCH, T in [0,1], vrms/h_max match the ksponly baseline). No manual snes_type=ksponly needed. Closes #403 item 1. Existing free-surface + rotated level_1 tests green (16). Underworld development team with AI support from Claude Code
…#403 item 3) The datum surgery (_set_rows_local, the xhat ownership loop, per-node datum eval in build_rotation) is the np>1 crash class if any of it indexed the local slice with global rows. New tests/parallel/test_1064: prescribe u.n = cos(theta) on the annulus surface and check, with parallel-safe reductions, that the datum is imposed (integral relL2 ~ 4e-5) and ∫v·v matches the serial reference to iterative round-off. Verified: np1 relL2 3.978674e-05 / ∫v·v 6.6559607579; np2 IDENTICAL relL2, ∫v·v to 9 digits — fully partition-independent. (Skips serial per the mpi(min_size=2) convention.) Underworld development team with AI support from Claude Code
…4 to PR #404) PR #404 (bugfix/3d-boundary-flux-desmear -> development) adds a NEW tests/parallel/test_1064_rotated_freeslip_parallel.py. Rename ours to 1066 to keep test numbers unique once both land. No content change; the datum coverage is unchanged. Underworld development team with AI support from Claude Code
…nsport (annulus seam) On an annulus the surface filter and the tangential transport ran on the LOCAL x-sorted ring. Two failures fell out of that: 1. SEAM at 3 and 9 o'clock (θ=0, θ=π). The filter mapped the ring to an extract_surface submesh by KDTree, but the ring (106 nodes) and the submesh (104 DOFs) disagreed, so two pairs of distinct ring nodes — at θ≈0 and θ≈π, the x-extremes where dx/dθ→0 — collapsed onto one submesh DOF. array[map]=values is last-write-wins, so both nodes read back the same smoothed value → a 5 km radial jump exactly at those two angles (measured; median node jump ~1.2 km). Transport was never involved. 2. NOT PARALLEL. Both the filter (submesh-mapped) and the transport (launch-point interp) were rank-local; a launch point or smoothing stencil crossing a partition boundary was mishandled. Replace both with a globally-gathered, along-surface-ordered ring (θ on an annulus, x on a box): _build_ring_gather caches the concat→s-sorted permutation once (the surface deforms only along the normal, so the order is invariant); _ring_gather/_ring_scatter move a per-local-node array to/from the full sorted ring; the filter is a periodic 1-D Taubin on that ring; the transport interpolates launch points against the full ring. No submesh, no KDTree collision, parallel-correct by construction. Validated: serial max node-to-node jump 5.3→1.79 km (uniform, seam gone); np1 vs np2 give the same ring (max_jump 1.79/1.80, ∮dr² match to solver round-off) with no partition- boundary seam; test_1070 free-surface + test_1018 rotated tests green. Underworld development team with AI support from Claude Code
…lus seam) test_freesurface_annulus_ring_filter_no_seam (level_1): filter a smooth ring field and assert the max node-to-node jump is uniform around the ring (no outlier at theta=0/pi, where the old extract_surface KDTree mapping collided 2 ring nodes onto 1 submesh DOF and seamed by ~5 km), and that _ring_gather -> _ring_scatter is the identity. Underworld development team with AI support from Claude Code
…tep, keep the surface cap The semi-Lagrangian transport is unconditionally stable, so estimate_dt gains an advect_scale that multiplies the ADVECTIVE limit only; the surface-motion cap (max_surface_cfl) still binds. Motivation: with a tiny advection-CFL dt the surface barely moves per step and the flow crawls through a near-flat spin-up transient where the free/held/consistent solves are ill-conditioned — measured as 12-to-73-minute non- converging solves (fgmres hitting its 300 max_it, each iteration an expensive coarse solve). A larger advective step blows through that transient in a few steps and the solves stay well-conditioned (~2s each) — no solver-code change needed. Crucially advect_scale does NOT scale the surface cap: multiplying the FINAL estimate_dt instead (an earlier driver bug) let a soft surface over-deform and the surface-flow feedback run away (h_max 27% of radius). Scaling only the advective part keeps the mesh safe: stiff surface takes big steps (advection-bound), soft surface stays capped (surface-bound), both stable. Underworld development team with AI support from Claude Code
…n the free surface Two joined fixes for hot material appearing in the COLD thermal boundary layer of a free-surface convection run, at the inflow/outflow (surface convergence) points. 1. return_coords_to_bounds is now deform-aware (Mesh). The mesh-specific ANALYTIC restore (annulus radial, box faces) is captured at construction with the ORIGINAL radii/extent. Once mesh.deform() moves the geometry that closure tests the WRONG surface: where a free surface has dipped inward, a point beyond the true boundary is not detected as outside, is never restored, and falls through to evaluate()'s RBF/Shepard fallback, which averages distant DOFs — injecting hot interior material into the cold TBL. deform() now sets _geometry_deformed (and release()s any rigid radial/plane BoundingSurface), and the property switches to a general facet restore: nearest point on the CURRENT boundary facets (_nearest_on_facets_2d/3d) with an outward-normal side test built from _boundary_facets' opposite-vertex orientation, so interior points are returned untouched. Cached against the coordinate array so it refreshes after each deform. Non-simplicial meshes keep the analytic closure. Verified: on an annulus deformed 1.0 -> 0.92, a point at r=0.95 (outside the TRUE surface but inside the ORIGINAL radius, which the analytic closure misses entirely) is restored to 0.92, while an interior point at 0.80 is untouched. Normal orientation is 100% correct on both arcs. 2. FreeSurface's composition transport now uses monotone_mode="clamp" (was None, a level-set default). With old_frame_traceback=True the departure point is deliberately NOT clamped, and SemiLagrangian._trace_back documents the limiter as the ONLY thing bounding a foot that lands outside the old mesh — so None was a contract violation. Measured (isoviscous annulus, gap-Ra 3e4, rho_g 5e5, step 60): analytic restore + no limiter : T [-0.361, 1.000] Nu 5.48 h 164.1 km + clamp only : T [-0.313, 1.000] Nu 5.34 h 162.7 km + clamp + deform-aware restore: T [-0.034, 1.000] Nu 2.97 h 165.1 km Undershoot down ~10x; the inflated Nu was the injected hot material, not physics. The topography is unchanged, i.e. the free-surface machinery was never the problem. Underworld development team with AI support from Claude Code
…ibed flux must be weak) Reverts the consistent-solve half of 7134d69. A prescribed-flux datum on a CLOSED boundary must be imposed WEAKLY: the penalty tolerates the small discrete incompatibility between the demeaned target and incompressibility (∮u.n dA = 0), whereas a strong rotated constraint over-determines it. The rotated prescribed-normal datum primitive itself (rotated_bc.py, 1a08d39) is retained and still tested — it is correct where discrete compatibility is guaranteed; it is the wrong tool for this particular BC. MEASUREMENT CAVEAT (honest): this revert is NOT yet demonstrated to fix the material-boundary error. At matched steps the manager's consistent solve shows |u.n - u~n|/|u~n| ~ 13-21% with EITHER the datum or the reverted penalty, while a stripped-down standalone penalty solve achieves 1.8-2.8% on the same problem. So the third solver is under-enforcing its material-boundary condition in the manager, and the cause is something OTHER than (or in addition to) the datum swap — candidates are the wall replay (_apply_walls) and the solver options inherited via _new_stokes. Tracked as the next diagnostic; the revert stands because it restores the previously validated design and the datum demonstrably did not help here. Underworld development team with AI support from Claude Code
…aint is the right tool The consistent (third) solve exists to hold the surface as a MATERIAL boundary by imposing u.n = ũ_n. It was measured under-enforcing at 13-21%, and 6067b26 reverted it from the strong rotated datum to a penalty on the reasoning that a prescribed flux on a closed boundary must be imposed WEAKLY. Both the measurement and the reasoning were wrong; this commit corrects them. The measurement was wrong. It sampled the advection velocity at fs._surf_coords -- the surface node coordinates captured once in __init__ and never refreshed after mesh.deform(). By step 12 that point sits 0.55 cells inside the moving surface. Sampled at the LIVE surface the same error is 3.8e-2, not 1.4e-1. The residue is below the 7.7e-2 gap between the raw P2 velocity and the smoothed P1 projection the target is built from, so the "5-10x worse than a bare solve" gap was mostly the metric. The reasoning was wrong in its premise. The prescribed datum is NOT meaningfully incompatible: measured with exact FE boundary quadrature the net flux was 8e-3 of the gross, and all of that came from _demean removing the NODAL mean of ũ_n. Surface nodes are unevenly spaced, so the node-count mean is the wrong quadrature; an incompressible interior over a closed no-slip base requires the ARC-LENGTH mean. Weighting the demean by ring arc length drops the net flux to 8e-5, a 100x reduction. With a genuinely flux-free datum the strong constraint is decisively better, not worse (annulus, gap-Ra 3e4, rho_g 5e5, dtscale 10, at matched steps): constraint |u.n - ũ_n|/|ũ_n| net volume flux through the surface penalty 3.8e-2 4.4e-2 strong 5.9e-3 7.0e-5 So the penalty was the leak: it misses the prescribed rate by 6.4x more AND passes 4% net volume through the surface, which is precisely the material-boundary violation that puts semi-Lagrangian departure feet outside the domain. Changes: - _surface_mean weights by ring arc length (live radii), not node count. Applies to the h / h_inf datum too, making it volume-preserving rather than node-count-preserving. - _ring_coords becomes a property reading LIVE mesh coordinates; _surface_normal_velocity and _tangential_velocity sample there. _current_shape already read geometry live, so velocity and geometry sampling were inconsistent. Low impact on the trajectory here (the P1 projection is smoothed over ~5 cells) but it was wrong. - consistent_constraint={"strong","penalty"} keyword, default "strong"; the penalty path is kept as a switchable fallback with consistent_penalty. - _build_ring_gather moves up to the geometry block, since _surface_mean now needs it. Tests: test_freesurface_prescribed_rate_is_flux_free (validated to FAIL at 9.0e-3 with the old nodal demean, threshold 5e-3) and test_freesurface_strong_constraint_beats_penalty. 19 pre-existing rotated-BC and free-surface tests still pass. Underworld development team with AI support from Claude Code
…a nonlinear rheology The rotated prescribed-normal datum is implemented on the LINEAR rotated path only, and the solve dispatches by a measured residual probe — so with consistent_constraint="strong" (now the default) a genuinely nonlinear rheology plus a composition field raises NotImplementedError from the primitive, where the penalty previously worked. Catch it at the manager boundary and name the knob (consistent_constraint="penalty") instead of leaving the caller with the primitive's message. Not covered by the existing nonlinear test, which calls solve() but never advance() and passes no composition, so the consistent solve never runs there. Also documents consistent_constraint / consistent_penalty in the class docstring. Underworld development team with AI support from Claude Code
… tests on #421 The two new datum tests fail on np2 for a reason they do not diagnose, so they are skipped there with the measurement in the reason string, and the half of the path that IS parallel-correct gets its own test that runs on every rank count. The datum field carries ~6-7% net flux on np2 versus ~0.06% in serial, INDEPENDENT of how the surface mean is removed -- so this predates the arc-length gauge and affects the penalty path too (measured: nodal demean 8.96e-3 serial / 7.06e-2 np2; arc-length 5.63e-4 serial / 5.75e-2 np2). It was never caught because the datum's net flux had not been measured in parallel, and a weak penalty absorbs the incompatibility quietly. The ring machinery is NOT the cause and is now pinned by test_freesurface_ring_quadrature_is_exact_in_parallel: on np2 the gathered ring holds 66 entries for 64 distinct nodes (the 2 seam vertices appear on both ranks) and the trapezoidal rule splits their weight, giving sum(w) = 2*pi*r_out to machine precision and weighted means of cos(k*theta) at ~1e-16. The loss is in the surface-array <-> field round trip across the seam. Filed as #421 with the suspects ranked. Underworld development team with AI support from Claude Code
…ing-mesh amplifier (#423) The serial T-unboundedness in annulus free-surface convection (T reaching [-1.6, 3.2] by ~8% deformation) was characterised stage-by-stage and localised by ablation. It is NOT a landing-point problem, NOT the monotone limiter, NOT the flux-history frame, and NOT the timestep: it is an exponential instability of the SL record->trace->solve loop whose necessary ingredients are per-step mesh motion, a sharp front in squeezed cells, and old_frame_traceback=True. Evidence chain (probes in ~/+Simulations/.../teaching/, full tables in #423): - Stage attribution: the clamp invariant range(psi*) within range(T_pre) held at every step (traceback CONTRACTS; feet exonerated) while the FE solve grew both extremes ~10%/step; the worst violator sat mesh-locked at r=0.952 for 22 straight steps. - Levers: BE (theta_flux=1) delays then saturates a period-2 limit cycle -- accelerant, not engine (per Louis: a step towards the diagnosis, not the cure). Flux-history old-frame: no effect. TRUE dt halving: WORSE at matched deformation -- the gain is per loop CYCLE, not per unit time. - Minimal reproducer (no Stokes, no free surface; prescribed 4-cell velocity, sharp cold layer, 9% squeezed annulus): static mesh bounded; +/-0.1%R/step wobble -> T=160 in 60 steps. The free surface is exonerated as physics; movers/adapt are equally exposed. - 2x2 localisation under motion: oldframe ON -> runaway at theta 0.5 AND 1.0; oldframe OFF (standard ALE) -> bounded at both. - Verification on the real problem: rho_g 2e5 with old_frame_traceback=False holds T in [-0.001, 1.000] for 60 steps through 16.9% deformation (previously unusable beyond ~5.4%); 5e5 exactly [0,1]. The irony: old-frame was introduced as the FIX for the earlier high-Ra blow-up, validated on a Cartesian box at small deformation. The hazards that motivated it are now covered by fixes landed since (deform-aware return_coords_to_bounds for the feet, the monotone clamp on the sample), and its old-geometry sampling leg is the amplifier. WHY it amplifies remains open at source level -- tracked in #423; the facility stays in ddt.py with its pathology documented there. Adds test_1071_slcn_moving_mesh.py: the wobbling-squeezed-annulus reproducer in the manager's configuration (standard ALE + clamp) must stay bounded -- a 12-step guard that reaches T>1.7 under the old-frame path. Underworld development team with AI support from Claude Code
…aluate The no-composition timestep path evaluated the free velocity at every mesh vertex through point location, which legitimately fails (loud NaN policy) on a deformed mesh and poisoned dt -> advance(NaN). A max-speed bound needs no interpolation: read the nodal data directly. Found by the topographic-relaxation A/B (teaching/relaxation_test.py), which also demonstrates that the reduced (driving-only) body force carries NO restoring force for surface deformation (h frozen at machine-noise velocities), while the full Boussinesq body force rho = rho_0(1 - alpha*dT) relaxes an imposed 4% topography exponentially at the Cathles scale (measured 2.1e4 vs rho_0*g/(2*eta*k) = 2.5e4) through the manager UNCHANGED. Formulation adoption for convection runs is a separate design step. Underworld development team with AI support from Claude Code
…REQUIRED on a deforming mesh Louis's ruling, demonstrated twice: the reduced (driving-only) body force has no restoring force for surface deformation anywhere in the momentum system (relaxation test: imposed 5% topography FROZEN), while rho = rho_0(1-alpha*dT) relaxes it at the Cathles rate (2.1e4 vs 2.5e4) and in convection holds the topography at its supported amplitude (restoring_force_demo.png) — through the manager unchanged. Records the alpha*dT <= ~0.3 coupling constraint (softness is not a free knob), the tolerance guidance for the hydrostatic-dominated RHS, and the #423 old-frame failure-mode row. Underworld development team with AI support from Claude Code
Louis's questions drove this to ground: at the full-density "equilibrium" the surface parked at 24% of the isostatic m=4 amplitude while the free solve pumped rms(u.n)=143 (~vrms) through a stationary surface. His framing — "if the surface is near equilibrium, release it and it does not move; why does the held lid feel a need to move?" — is exactly right in continuum, and the defect is NOT the held solve: it is the CBF reaction RECOVERY, which on a deformed boundary returns only ~0.16 of the current-shape load rho_0*g*h (probe A: pure statics, u=0, recovery 0.158 of the exact reaction) while recovering the driving component exactly (probe B: full vs reduced on the undeformed mesh, ratio 1.0000). The conversion then builds h_inf from a reaction that has forgotten 84% of "where we are". Fix: background_buoyancy=<background body force>. The manager runs a TWIN held solve with the background force alone on the same geometry; both recoveries carry the identical defective load leg, so the DIFFERENCE isolates the driving support sigma' exactly, and h_inf = sigma'/(rho_0*g) is the true equilibrium topography. Default None preserves the reduced-formulation path bit-for-bit (test_1070 green). Validated: - Relaxation test (T=0, 5% topography): corrected h_inf ~ 0; decay 264 km -> EXACTLY 0 (the old 5 km floor was the corrupted target's residue), early rate 1.7e4 vs Cathles 2.5e4. RETRACTION: the earlier "2.1e4" agreement of the uncorrected run was partly the scheme contracting geometrically to its wrong target. - Equilibrium restart (equil_full_1e6@300): through-flow ratio 2.15 -> 0.075 in ONE step; h_m4 climbs 0.00233 -> 0.00953 vs isostatic referee 0.00975 within 25 steps; T exact [0,1]. Cost note (Louis): the twin solve is a stopgap. The load is ANALYTIC (rho_0*g*h from the geometry); once the CBF load-leg recovery defect is fixed (kappa=0.16 -> 1, to be filed), h_inf = h + (sigma_rec - rho_0*g*h)/(rho_0*g) needs NO extra solve. The twin also depends on geometry only, so it can be cached/refreshed every N steps meanwhile. Underworld development team with AI support from Claude Code
…-density h_inf; CBF recovery EXONERATED Correction of the mechanism behind 997a6b0, with the label-based probes redone: the CBF reaction recovery is essentially EXACT in every leg (load on a 5%-deformed ring: 0.976 of the analytic reaction; driving deformed/undeformed: 0.972; full = exact superposition). The 0.158 "recovery defect" in issue #431 was an artifact of the probe's radius mask, which selected only crest arcs on a deformed ring — the label-based redo overturns it. The REAL defect is a sign/semantics error in the manager: under full density the recovered field contains the self-load component +h_current, and the reduced-form negation flips it, giving h_inf = -h_current + driving/rho_g. That map's fixed point sits at HALF the equilibrium (measured parking at a fraction of isostasy) and its -1 eigenvalue is the period-2 ringing observed in every full-density run. Since the recovery is faithful, the fix needs NO extra solve (Louis's requirement): background_buoyancy="analytic" subtracts the GEOMETRIC current height from the recovered field before the standard conversion, cancelling the self-load exactly (to the recovery's own 2.4% accuracy). Validated: - relaxation: 264 km -> 0.12 km, trajectory matching the two-reaction reference to <0.5% at every step (floor 0.12 km = the 2.4% residue vs exactly 0 for two-reaction); - equilibrium restart: bit-identical to the two-reaction mode at reported precision (step 15: h_m4 0.00747, through-flow ratio 0.245 in both). The two-reaction mode (background_buoyancy=<expr>) is retained as the exact reference; "analytic" is the recommended production setting for full-density runs. Underworld development team with AI support from Claude Code
…896114) Records the final ruling: the parked-surface/through-flow defect was the manager's reduced-form negation flipping the self-load in the recovered reaction; the CBF recovery is exact. 'analytic' subtracts the geometric height - no extra solve. Probe lesson recorded: select boundary DOFs by label, never a radius mask. Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
Member
Author
Adversarial self-review (house rule: posted before marking ready)Attacking the weakest points of this PR honestly: 1. The parallel exposure is real and is shipped as a known hole. The strong 2. The commit history contains a corrected mis-attribution. 997a6b0's message 3. The 'analytic' subtraction leans on the recovery being faithful (0.976). The 4. Old-frame retirement is manager-scope only. old_frame_traceback=True remains 5. Test-time cost. test_1070 grew two level_2 tests (~90s) and test_1071 adds 6. What is NOT validated here: 3D (all measurements are 2D annulus/box); None of these block the merge in my judgement; 1 and 2 deserve reviewer attention. Underworld development team with AI support from Claude Code |
This was referenced Jul 27, 2026
…datum fixed (#421) Closes the #421 defect chain, localised by three successive np2 probes: 1. u_n was multi-valued at seam nodes: _surface_normal_velocity point-evaluated the P1 velocity at its OWN nodes, and on-vertex point location mis-locates at process seams (the documented ddt band-aid class). Fixed by direct nodal data reads (the projection has already ghost-synced) — also faster. 2. With that fixed, h_inf remained ~50% multi-valued at seams: the gathered ring holds a seam vertex ONCE PER ADJACENT RANK, and the Taubin filter's roll-stencil treats the copies as two distinct neighbouring nodes (a phantom zero-length segment), smoothing each against a DIFFERENT stencil; each rank then scattered back its own diverged copy. Fix: deduplicate the ring at construction — gather averages seam copies (identical to round-off), all ring operations (filter, tangential transport, arc-length weights, datum gauge) act on unique physical nodes, scatter expands the unique value back to every copy. Ring fields are single-valued across ranks by construction. Measured: prescribed-datum net flux at np2 1.142e-3 vs serial 1.140e-3 (was 5.75e-2, 50x worse than serial). Seam probe shows zero multi-valued ring quantities. The two datum tests are un-skipped and pass at np2 (8 passed with --with-mpi). Also: NotImplementedError guard for dim != 2 (previously a 3D run would proceed with meaningless ring ordering and a stub de-smear — silently wrong; 3D is follow-on work via PR #404's boundary mass), and the background_buoyancy docstring entry. Underworld development team with AI support from Claude Code
Member
Author
|
Update to the adversarial review, point 1: the parallel datum defect (#421) is now FIXED in this PR (67d811d — ring seam deduplication + direct nodal velocity reads). The datum tests run and pass at np2; no parallel gating of the strong-constraint default is needed. Point 1 is withdrawn; points 2-6 stand. Also added in the same commit: an explicit NotImplementedError for 3D (previously would run silently wrong; follow-on work via #404's boundary mass), per discussion. Underworld development team with AI support from Claude Code |
Plain-language plan for the post-#435 work: A (land 2D, teaching rebuild, stationarity run), B (nonlinear datum + rotated-path unification, joint session, #438/#403), C (3D, structural, builds on #404 — suggested for Thyagi), plus the deferred/tracked list. 3D is assembly work, not new mathematics: the integrator, datum, carrier and full-density h_inf are dimension-free; the 2D-only pieces and their 3D replacements are tabulated with the #421 seam-dedup lesson carried over. Underworld development team with AI support from Claude Code
Member
Author
|
Roadmap for the follow-on work is now in the PR: docs/developer/design/FREE_SURFACE_ROADMAP_2026-07.md — three workstreams (A: land this PR + teaching rebuild + stationarity run; B: nonlinear datum + rotated-path unification, joint session, tracking #438; C: 3D, structural, builds on #404, suggested for @gthyagi), with sequencing and the deferred list. C is assembly work, not new mathematics — the integrator, datum constraint, deform carrier and full-density h_inf are all dimension-free; the 2D-only pieces and their 3D replacements are tabulated. Underworld development team with AI support from Claude Code |
…sidual report) Underworld development team with AI support from Claude Code
…rkstream C Underworld 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 free
to 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.
What this delivers
The free-surface toolkit (uw.systems.FreeSurface) taken from "stable only for stiff
surfaces at small deformation" to a validated capability: bounded temperature transport
to at least 17% surface deformation, a correct material-boundary constraint, and the
full-density (Boussinesq) formulation with its surface restoring force — the setting in
which topography relaxes at the physical rate and equilibrium states have closed
streamlines with no flow through the surface.
The four main fixes, each with its measurement
1. Material-boundary gauge (2c522c5, 8f49006). The rate prescribed to the
consistent solve must carry zero net flux on a closed boundary; the demean removed the
node-count mean instead of the arc-length mean, leaving ~1% net flux that a strong
constraint cannot absorb. With the arc-length gauge the strong rotated datum beats the
penalty on every measure: rate error 5.9e-3 vs 3.8e-2, net volume leak 7.0e-5 vs
4.4e-2, and the hot-injection Nu inflation disappears (3.06 vs 5.69 at matched steps).
2. Old-frame reach-back retired from the manager (b507aca, issue #423). The serial
temperature blow-up beyond ~5% deformation was localised by ablation to
old_frame_traceback=True — an exponential instability of the record/trace/solve loop
whose gain is per loop cycle (halving dt makes it worse). A minimal reproducer with no
free surface at all (prescribed velocity + a 0.1%/step mesh wobble) reproduces it;
with the standard ALE path the same configuration is bounded at any theta. The manager
now uses the standard ALE path; envelope measured from 2.6% to 17% deformation with
T in [0,1] to 1e-3. Guard test: test_1071_slcn_moving_mesh.py.
3. Timestep estimate on deformed meshes (84f8679). The no-composition CFL path
evaluated velocity by point location at mesh vertices, which legitimately NaNs on a
deformed mesh; it now reads nodal data directly.
4. Full-density formulation support (997a6b0, 9896114). With
rho = rho_0(1 - alpha dT) the free solve carries the surface restoring force
(topographic relaxation at the Cathles rate; without the background there is provably
no restoring force and topography grows without limit). The recovered held-lid
reaction then contains the self-load of the current shape, and the pre-existing sign
convention flips it — the surface parks at half its equilibrium with steady flow
through it. Fix: background_buoyancy="analytic" subtracts the geometric current
height (no extra solve); a twin-solve expression mode is retained as the exact
reference; both validated identical on the relaxation benchmark and an equilibrium
restart (surface climbs from 24% to the ~80%-of-isostasy dynamic equilibrium, with
the through-flow ratio collapsing 2.15 -> 0.075 in one step). The CBF recovery itself
was verified essentially exact (0.976/0.972 on a 5%-deformed ring, exact
superposition) — issue #431 records the retraction of an interim mis-attribution and
is rescoped to the benign 2.4% residual.
Default behaviour (no background_buoyancy, reduced formulation) is unchanged and
covered by the existing tests.
Tests
23 serial (test_1018, test_1070, test_1071) + 6 at np2 green. Two new datum tests are
skipped at np>1 against pre-existing issue #421 (parallel datum round trip), with the
measurement in the skip reason.
Issues touched (none auto-close; to be closed by hand where appropriate)
remains open.
An adversarial self-review follows as the first comment.
Underworld development team with AI support from Claude Code