Skip to content

Add CitcomS predictor-corrector DDt manager for composed transport - #689

Open
gthyagi wants to merge 36 commits into
underworldcode:developmentfrom
gthyagi:feature/zhong2008-supg-integration
Open

Add CitcomS predictor-corrector DDt manager for composed transport#689
gthyagi wants to merge 36 commits into
underworldcode:developmentfrom
gthyagi:feature/zhong2008-supg-integration

Conversation

@gthyagi

@gthyagi gthyagi commented Sep 5, 2026

Copy link
Copy Markdown
Contributor

Summary

Add the Zhong/CitcomS predictor-corrector as a transport manager for the
composing uw.systems.AdvDiffusion solver introduced in #688.

transport = uw.systems.ddt.EulerianSUPGPC(
    mesh, T, U.sym, method="citcoms", temperature_rate_field=Tdot,
)
thermal = uw.systems.AdvDiffusion(mesh, T, U.sym, DuDt=transport)
thermal.constitutive_model.Parameters.diffusivity = 1.0
thermal.solve(timestep=thermal.estimate_dt())

The solver assembles diffusion, source, and the manager's transport
contributions. The manager owns the persistent rate, steady directional
SUPG parameter, correction loop, stability estimate, and restart metadata.
The scalar solver supplies residual assembly through a small execution hook.
Implicit CN/BDF continue to use the existing EulerianSUPG manager.

Numerical Behavior

Manager selection Purpose Accuracy contract
method="citcoms" Fixed two-correction Zhong/CitcomS update Matches the independently assembled recurrence; approximately first order for nonuniform diffusion
method="pc_converged" Consistent-residual correction to tolerance Second-order reference with adv_gamma=0.5; raises on non-convergence

Both PC modes require continuous scalar P1 fields on triangles or tetrahedra.
Their timestep is 0.9 * min(dt_adv, dt_diff). They reuse PETSc vectors and
geometry arrays and allocate no unused implicit history.

Restart and Integration

  • Manager snapshots preserve rate initialization, numerical configuration,
    correction diagnostics, and timestep; mesh snapshots carry the rate field.
  • Solver snapshots preserve the timestep and field-change estimator.
  • Implicit cell-size fields are registered before fresh-process restoration.
  • A PC manager must bind to its exact temperature field and mesh.
  • Checkpoints must be restored into a matching object graph. Snapshots from
    the old standalone SUPG solver are not a migration format.

Validation

Validated locally on 8 September 2026 against the merged development
architecture. Times are elapsed suite times, not CPU-hours. Independent
suites sometimes ran concurrently; rows are not an additive wall-time total.

Check Result Time
Focused serial API, manager mutations and snapshots 33 passed 32.58 s
Frozen f41bcd2f source versus manager, both PC modes in 2-D/3-D 4 passed 29.32 s
Serial time-order, fresh-process restart and workspace reuse 20 passed 185.39 s
Additional transport/transient tests 14 passed 301.80 s
Extended analytical pulse/rotation/spherical diffusion suite 26 passed 1226.94 s
Eight-rank API/equivalence/time-order/workspace/manager suite 49 passed per rank 284.75 s
Eight-rank disk replay, all five methods 5 passed per rank 24.60 s
Eight-rank separate-process restart, four methods 4 passed 158.37 s
Full repository Level 1, serial 1797 passed; 41 skipped; 2 expected failures 1633.89 s

The initial MPI group had five obsolete checkpoint-capability skips, caused
by inspecting the enhanced-variable wrapper signature. That guard is now
removed and all five disk cases pass in the separate row above. Six
200-update RSS soak cases remain opt-in Level 3 tests; the default workspace
checks run eight updates. Optional frozen-baseline cases need the old source
file and were run separately in serial and included in the MPI group.

Fixed PC2 retains order approximately 1.01; composed CN and converged PC
measure order 2.00 on triangles and tetrahedra. Converged-PC agreement with
the independent trapezoidal map is below 5.2e-14. MPI fresh-process PC
restarts are exact; CN/BDF2 maximum field differences are below 9e-16.

The full Level 1 invocation (pytest tests -m level_1) also deselected 1362
other-level cases. Its 27 min 13.89 s duration is the repository-wide suite,
not a focused SUPG subtotal. GitHub CI on this new revision is pending;
previous-revision CI is not counted as validation.

Scope

Targets development following #688. Generic mesh fixes #691/#692 are already
upstream. No Stokes, geoid, Navier-Stokes, or dependency changes are included.
The implementation and small tests do not require a long coupled A1 run.

Checklist

  • Predictor-corrector moved into the DDt transport interface.
  • Fixed and converged methods explicitly distinguished.
  • Independent analytical, restart, and workspace tests migrated.
  • Serial and eight-rank migration validation completed.
  • CI passed on this revision.

Maintainer review of the manager API remains requested.

lmoresi and others added 21 commits September 2, 2026 17:29
A solver that assembles its own weighted sum of history terms (an Eulerian
scheme applying a multistep rule to a spatial operator) needs the
constants-routed coefficient expressions, not just their current values.
Read-only accessors; no behaviour change.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…vars

A MeshVariable that is dropped and garbage-collected (the default Model
holds the only strong reference; uw.reset_default_model() releases it, and
the statistics helpers delete temporaries deliberately) leaves its PETSc
field in the DM. Mesh.update_lvec zipped mesh.vars.values() against the
field decomposition by position, and the JIT's petsc_a[] offsets were a
running count over the live variables, so every later variable was packed
into, and read from, the wrong slots. Measured: a P0 cell-size field landing
in a P2 slot as garbage, NaN residuals in one run and a subtly wrong answer
in the next, depending on when the collector ran.

update_lvec now packs by field name and zeroes an orphaned field; the JIT
reads component offsets from the DM's own field list and patches each
variable from its field_id. Regression test: 2 of its 3 checks fail without
the fix.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…r for scalar variables

A Gaussian carried round the origin by rigid rotation while diffusing is
exact at every time (rotation commutes with the Laplacian), so a transport
scheme's error can be measured directly and the round trip after one
revolution is an absolute check. AnalyticSolution.error(norm='integral')
added a 1x1 Matrix symbol to a scalar expression and had never been
exercised on a scalar variable.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…fusion solvers

The cell-crossing / diffusion-time reduction (isotropic or direction-aware,
minimum or percentile) becomes a module-level helper so the Eulerian
solver can call it rather than carrying a copy. SLCN behaviour unchanged.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…nditioner

A solver with no managed option block (_pc_option_prefix is None) sets its
own PC; installing the adapt child's PCMG hierarchy on it segfaulted inside
PETSc (additive-Schwarz PC, PCMG calls). The gate now treats that state as
the explicit choice it is, alongside preconditioner='gamg' and the user
override latch.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…from the symbolic history

uw.systems.AdvDiffusionSUPG(mesh, T, V_fn, order=N, integrator='bdf'|'am')
assembles the implicit weak form from the Eulerian DDt history: the BDF
stencil or the Adams-Moulton weights on the advective and diffusive terms
at every stored time level, plus the SUPG flux tau R u with the strong
residual of the same scheme. Timestep, multistep coefficients and the tau
weights are runtime constants of the compiled kernels, so a change of dt
costs nothing (the issue underworldcode#657 prototype recompiled on every change).
Diffusivity comes from the constitutive model like every scalar solver.

Measured on the rotating Gaussian: stable at any cell Courant number, error
set by u dt against the feature width (dt^2 for the second-order schemes),
unchanged to three digits by a band refined to h/9 at local Courant 13;
Crank-Nicolson reproduces the prototype's numbers to four digits.
Tests: API and no-recompile contract, temporal convergence (slopes 0.8/0.9
for BDF1, 1.9 for BDF2), band invariance, round trip, np=2 = serial.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…om the integrator study

Rotating-Gaussian study at res 32, Courant 0.25 to 8, pure advection and
kappa 1e-3: Adams-Moulton above order 1 blows up from Courant 1 (bounded
stability region), BDF3 fails from Courant 4, Crank-Nicolson is three to
four times more accurate than BDF2 at the same timestep but rings once the
feature is under-resolved in time, backward Euler carries 20-40% error at
any practical timestep. Cost per step is the same for every scheme. BDF2 is
the robust default; the note records the alternatives and when to pick them.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…Courant number

BDF2 and Crank-Nicolson track their res-32 errors at the same u dt. BDF3's
stability region misses the imaginary axis near the origin, so the
low-frequency modes of a finer mesh grow: 31x the exact field after 590 steps
at Courant 1. Safe only with diffusion, below Courant 2. Note and docstring
updated; the BDF2 default stands.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…-in replacement

The constructor, order, theta, f, V_fn, constitutive_model, delta_t,
estimate_dt and solve keep the meaning they have for AdvDiffusionSLCN, so a
script changes the class name and nothing else. order=1 with theta=0.5 is
Crank-Nicolson and the default, as for SLCN; order=2 takes theta=1 (BDF2)
unless 0.5 is asked for explicitly, which is refused for the reason the SLCN
documentation gives. The trace-back-only arguments (restore_points_func,
monotone_mode, old_frame_traceback, DFDt) are accepted and ignored with a
warning. integrator is inferred and only needs setting to reach the higher
Adams-Moulton rules. delta_t is settable and solve() reuses it; the notebook
viewer reports the scheme. User page docs/advanced/eulerian-advection-diffusion.md
with the swap table and the when-to-use-which guidance.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…e scheme

The only schemes the argument added were Adams-Moulton at orders 2 and 3,
which the integrator study shows blowing up on advection from Courant 1.
The multistep family now follows the order (the theta rule at order 1, BDF
above); the higher Adams-Moulton assembly stays in the code, reachable only
by switching the family on the instance, which is how the study measured it.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…the module and note

estimate_dt now returns the step at which the field changes by a fraction
(0.02) of its range: from the advective rate |u . grad phi| at the vertices
before the first solve, and from the rate the last step actually produced
after it. The cell-crossing time the semi-Lagrangian solver reports is not
a stability limit for this scheme and says nothing about its accuracy; it
stays available as basis='resolution'. The estimate is mesh-independent,
which the band test now checks (the resolution estimate collapses 3x on
the refined child, the accuracy estimate moves under 25%), and at the
default fraction Crank-Nicolson completes the rotating-Gaussian round trip
under one per cent. The advective rate uses the vertex Clement gradient
rather than a point evaluation of a derivative expression, which fails on a
mesh carrying many variables.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…oner="fmg" a real switch on the SUPG solver

The Eulerian SUPG step took two Newton iterations on a linear operator:
the Krylov default (rtol 1e-5) does not reach the SNES tolerance (1e-8),
and the second Jacobian assembly cost more than every linear solve of
the step. The Krylov tolerance is now 1e-9 and a step is one Newton
iteration: 1.54 s to 0.91 s per step at 256^2 in serial.

Measured against geometric multigrid at matched tolerances (design note,
"Preconditioner"), GMRES with additive-Schwarz ILU is the cheaper linear
solve at every Courant number from 1/2 to 32 and its iteration count is
the same on one and eight ranks; the multigrid's cycle count grows with
the Courant number nearly as fast, and a cycle costs about three Schwarz
iterations. Schwarz stays the default on every mesh.

preconditioner = "fmg" now hands the block to the managed multigrid
route (custom-P transfers over the refinement hierarchy or an adapt
child's coarse tail, flexible GMRES outside) for the rank count where a
one-level method runs out of coarse space. The solver's solve() builds
through the base _build, where a preconditioner choice is resolved; the
pre-run of the three setup stages marked the solver set up first, so the
request was silently inert. The semi-Lagrangian solvers share that
pattern and the defect (underworldcode#683).

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
…lows

The shipped convection examples set adv_diff.theta = 0.5 after building
the solver; the Eulerian drop-in refused it. The blend is a runtime
constant refreshed from the history manager before every solve, so the
setter updates it without a recompile (order 1 only, the constructor's
rule). Vector and tensor unknowns join the design note's deferred list:
the solver is scalar, where the semi-Lagrangian trace-back carries them.

Underworld development team with AI support from Claude Code

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL
Use one public AdvDiffusionSUPG class and one F0/F1 assembly for implicit CN/BE/BDF and explicit CitcomS updates. Preserve the P1 positive lumped mass, gamma=0.5 two-correction update, directional tau and conservative timestep estimate, cached PETSc vectors and geometry. Replace the former implementation with compatibility imports; no duplicate solver remains.

Register timestep-estimator and integrator snapshot state, omit unused CitcomS DDt history, retain explicit BDF selection and custom tau, and validate incompatible settings. Add frozen-source triangle/tetrahedron equivalence and in-memory/disk restart tests. Replace the MPI Gaussian cross-host golden value with a same-host serial reference without loosening the threshold.

Validation: rebuilt Mac worktree; 25 unified/API tests passed, 13 residual tests passed, seven MPI tests passed per rank including the Gaussian regression, and nine expanded two-rank migration/snapshot tests passed per rank. Style gate and diff whitespace checks passed. Gadi coupled benchmarks and memory gate are still pending; this is not a production-acceptance claim.

Underworld development team with AI support from OpenAI Codex.
Select backward Euler explicitly for quasi-steady manufactured and high-Peclet tests. Check the public generic tau expression instead of the removed P0 implementation detail. Compare spherical serial/MPI values on the same host and cached triangulation; keep manufactured convergence as the absolute accuracy gate. Apply the same near-machine-precision implicit replay tolerance as the unified snapshot tests. Syntax and whitespace checks pass; Gadi numerical validation remains pending.
The Gadi eight-rank migration gate failed before any timestep because the tiny triangle fixture leaves unsupported local simplex layouts on some ranks. Propagate layout errors to every rank before entering timestep reductions. Increase the equivalence fixture resolution and add a separate empty-partition rejection regression. Preserve numerical algorithms and tolerances. Rebuilt locally; Gadi validation pending.
Gadi preserved snapshot fields exactly but a rebuilt CN solve differed by 3.45e-10 with the default 1e-9 Krylov tolerance. Tighten only the restart-test KSP/SNES tolerances and retain the near-machine-precision replay assertion. Production solver defaults are unchanged.
Eight-rank PC2 migration and field replay now pass. The CN rate differed by 2.04e-12 because taking a timestep derivative divides field roundoff by dt=0.003. Propagate the unchanged 2e-14 field comparison bound through max(abs(delta T))/dt rather than using an unrelated fixed rate tolerance. Snapshot restoration remains exact and production settings are unchanged.
Exercise the explicit P1 gamma=0.5 two-correction path against smooth published pulse solutions on triangles/tetrahedra, the existing rotating-Gaussian oracle, and exact radial diffusion in a spherical shell. Keep SUPG active for advection and prescribe velocity without Stokes.

Bound finite-domain tails, use common fixed timesteps for spatial refinement, check absolute errors and convergence, and optionally retain compact HDF5 diagnostics for same-host serial/MPI comparison. Verify the spherical equation symbolically. Solver implementation unchanged. Syntax checks pass; numerical Gadi validation is pending.
Refine the spherical diffusion pair from 1/4-1/8 to 1/8-1/16 without relaxing the 8 percent absolute L2 criterion. Add fixed-1/8 dt, dt/2 and dt/4 trajectories, restoring identical T/Tdot/startup state and comparing full FE fields.

Record input Gmsh SHA256 and short mesh filenames for strict same-mesh serial/MPI comparison. Reuse a focused advance/measurement helper instead of duplicating the diagnostic code. Solver implementation unchanged. Syntax checks pass; bounded Gadi numerical follow-up pending.
@gthyagi
gthyagi marked this pull request as ready for review September 5, 2026 21:22
@gthyagi
gthyagi requested a review from lmoresi as a code owner September 5, 2026 21:22
@gthyagi

gthyagi commented Sep 5, 2026

Copy link
Copy Markdown
Contributor Author

@lmoresi This PR is ready for a focused review once the current CI run completes. Two decisions remain intentionally open:

  1. Generic mesh fixes: should e4d0eaab (collective cell-family inference on empty MPI ranks) and fa75557f (cell-local stabilization length, fixing mesh.cell_size() is partition-dependent: the per-cell radius comes from a kd-tree over the rank's centroids, so cells at partition boundaries get a different size on each rank count #687) remain in this integration PR, or should I move either fix to an independent PR before merge? Both have focused serial/MPI regressions and are required by the validation cases in their current form.
  2. CitcomS PC2 contract: the fixed two-correction algorithm is preserved as a compatibility integrator. Isolated diffusion tests measure temporal order near 1.01, while UW3 Crank-Nicolson measures 2.00 against the same exact discrete mode. The documentation therefore does not claim second-order time accuracy for PC2. Is retaining that exact compatibility behavior as an optional public integrator the intended contract, or should it remain benchmark-internal pending a later converged-correction variant?

The former 200-step RSS soak is now opt-in (UW_RUN_SUPG_MEMORY_SOAK=1). Routine Level 2 CI uses six short lifecycle cases (three transport methods in 2-D and 3-D), while preserving the full soak for release/HPC validation.

Keep test 1077's default cellsize 1/8, numerical setup, 5% analytical bound and 1e-8 serial/MPI comparison unchanged. Expose only the internal helper's mesh size for the requested eight-rank 1/4, 1/8 and 1/16 Mac investigation after an empty-partition mesh.

All three resolutions completed after mesh fix 7e17bec. Same-host serial/MPI errors agreed within 8.47e-9 at 1/8 and 2.46e-9 at 1/16. Preserve and document the independent 1/4 discrepancy of 7.00e-7; it persists with tighter solves on identical mesh hashes. Do not relax the original gate.
Add tiny 2D triangle and 3D tetrahedron tests with independent analytical P1 element matrices and exact semidiscrete eigenmode references. Compare every UW3 T/Tdot update against a closed-form two-correction map without using a fine numerical solution, Stokes, or A1.

Demonstrate that consistent residual mass with lumped correction mass and two fixed corrections approaches (2I-D^-1 M)D^-1 K, with first-order timestep differences. Separate the startup-rate error using exactly solved CN and genuinely lumped-residual controls, which establish why uniform scalar decay is insufficient.

Document the limitation without changing CitcomS semantics, solver defaults, or prior tolerances. Four tests pass on Mac serial (17.46 s) and eight MPI ranks (32.91 s); matching mesh hashes and update discrepancies below 4e-14. Style gate and whitespace checks pass.
Exercise the actual shared Eulerian CN solver on the same tiny triangle/tetrahedron meshes and timestep sequence as the PC2 diagnosis. Check the exact CN amplification map and second-order error against independently assembled generalized eigenmodes.

Two tests pass in serial (26.61 s) and on eight MPI ranks (13.69 s), with order 2.00 in both dimensions. Finest relative errors are approximately 5.4e-9 and 5.7e-9. No Stokes/A1 run, numerical-method changes, or relaxed error tolerances.
A tiny three-process restart regression exposed a missing _h_cell field: implicit SUPG previously created its generic-tau geometry field only during residual construction, while snapshots require the same fields registered before loading. Materialize this existing dependency at solver construction, only for automatic generic tau; do not relax snapshot schema checks or change the time integration.

Add independent full/write/resume worker tests for PC2, CN and BDF2 on a 370-tetrahedron mesh with varying timesteps and velocity. Require exact restored fields/state, then compare step-12 continuation with uninterrupted runs. The CN case demonstrably failed before the four-line constructor fix.

After rebuild: serial 3 passed in 48.95 s; eight-rank workers 3 passed in 58.98 s. PC2 continuation is exact; CN/BDF2 fields agree within 8e-16. Update the user guide with the fresh-interpreter workflow. No Stokes/A1 run or checkpoint tolerance relaxation.
Reduce timestep-control changes, FMG comparison norms, and adapted-mesh pulse maxima globally instead of assuming every partition contains the Gaussian peak. Make equality and finite-value failures collective before subsequent solver calls.

Keep the existing 1e-3 change, 1e-6 relative error, and 0.9-1.01 peak thresholds. No equations or solver tolerances change. The prior rank-local failures stranded peers in later PETSc collectives.

Validation: all 48 API/residual/migration/diffusion checks pass on eight Mac ranks in 152.40 s. The three modified API tests also pass in the final 12-test serial lifecycle sequence (98.14 s). Global FMG discrepancy is 3.04e-10 with amplitude 0.96026 on eight ranks.
…checks

Exercise PC2, CN and BDF2 for 200 changing-velocity updates on tiny triangles and tetrahedra, without Stokes, checkpoint output, reaction diagnostics or forced GC. Record per-rank current RSS after warm-up and require stable solver/vector handles and PC2 workspaces. Keep preset per-rank bounds of 16 MiB growth and 0.05 MiB/step late slope, with collective temperature-finiteness checks.

Delegate fresh-process restart execution caps and descendant cleanup to the existing MPI supervisor instead of duplicating process-group termination. Document actual UW3 CN second-order discrete diffusion validation and the distinction between small lifecycle regression and production-scale acceptance.

Validation: final serial API/restart/memory sequence 12 passed in 98.14 s; final eight-rank restart 3 passed in 60.59 s; final eight-rank memory 6 passed per rank in 60.56 s. Maximum resumed field discrepancy 8.89e-16, PC2 exact. Repeated memory checks pass unchanged limits; no universal zero-leak or second-order PC2 claim. Style gate and git diff --check pass.
Run singleton fresh-process phases directly when the target branch predates the MPI supervisor merged in development by underworldcode#678. Keep MPI restart validation conditional on that supervisor so parallel descendants remain bounded and diagnosable.

This lets the SUPG feature branch validate PC2, CN, and BDF2 restart state without importing the unrelated 691-line supervisor change into this review.

Underworld development team with AI support from Claude Code.
Detect the same-layout checkpoint API added by upstream underworldcode#674 before asserting MPI disk replay. Older transport-branch checkouts continue to test serial disk restore and MPI in-memory restore; rebased development checkouts automatically exercise the full distributed disk path.

This keeps the SUPG PR focused while making its dependency on the already-merged checkpoint fix explicit.

Underworld development team with AI support from Claude Code.
Keep six fast Level 2 workspace-reuse checks for PC2, CN, and BDF2 on triangles and tetrahedra. They run eight updates and retain the deterministic object-identity, finite-field, and boundedness assertions.

Reclassify the 200-update RSS and late-slope regression as an opt-in Level 3 slow test selected with UW_RUN_SUPG_MEMORY_SOAK=1. Preserve its warm-up, sampling, thresholds, and six-case matrix, and document both execution paths.

Validation: default serial 6 passed/6 skipped in 16.47 s wall; default eight-rank 6 passed/6 skipped per rank in 35.03 s wall; representative opt-in PC2 triangle soak passed in 10.83 s wall. Deprecated-pattern and whitespace checks pass.

Underworld development team with AI support from Claude Code.
@gthyagi
gthyagi force-pushed the feature/zhong2008-supg-integration branch from b739d68 to 048ee09 Compare September 5, 2026 22:03
@gthyagi

gthyagi commented Sep 5, 2026

Copy link
Copy Markdown
Contributor Author

Follow-up: the two generic mesh fixes have now been removed from this PR and submitted independently:

The #689 history was rewritten and force-pushed without either commit. A clean rebuild succeeded. Its broad serial SUPG suite then produced 73 passed, 9 expected skips, and one expected prerequisite failure: without #692, the existing high-Peclet SUPG error is 0.13349, above the unchanged 0.10290 threshold. I have documented #691/#692 as prerequisites rather than weakening that assertion.

The earlier request to decide whether the generic fixes should stay in #689 is therefore resolved. The remaining review question is only the intended public contract for fixed two-correction PC2 and its documented temporal-order limitation.

@gthyagi
gthyagi marked this pull request as draft September 5, 2026 22:09
Preserve the fixed two-correction CitcomS compatibility mode while adding an explicitly selected pc_converged integrator. Iterate the full Petrov-Galerkin rate residual at startup and each timestep, using the lumped mass only as a reusable correction preconditioner and failing clearly when configured tolerances are not reached.

Persist correction controls in solver snapshots, expose convergence diagnostics, and extend workspace and fresh-process restart coverage. Add independent finite-element diffusion tests proving second-order convergence and exact trapezoidal-map agreement in 2D/3D serial and MPI without relying on the Zhong A1 model.

Document the mathematical distinction from fixed PC2 and the current repeated-residual performance cost.
@gthyagi
gthyagi marked this pull request as ready for review September 6, 2026 11:49
@gthyagi

gthyagi commented Sep 6, 2026

Copy link
Copy Markdown
Contributor Author

Marked ready for review without changing the requested target: feature/eulerian-supg-transport.

The head is 17 commits ahead of that target and contains no unrelated target divergence. Generic mesh fixes remain in independent PRs #691 and #692 rather than being duplicated in this SUPG integration diff. The current red CI result is the known partition-dependent Mesh.cell_size() check fixed by #692; the pc_converged tests did not fail.

DMPlexIsSimplex returns false on an empty partition. Using that local value for coordinate FE construction mixed simplex and tensor bases across the same communicator, consuming different PETSc message tags. A subsequent mesh HDF5 labelsLoad then blocked in PetscSFSetUp_Basic/MPI_Waitall.

Gather cell-family decisions from populated ranks before FE setup and use the same classification for element metadata. Preserve the constructor hint for an entirely empty mesh; do not change SUPG algorithms, solver tolerances or MPI providers.

Add single-cell triangle/tet/quad/hex regressions that require empty ranks and subsequently load another mesh, checking P2 volume and boundary integrals. The triangle regression failed before the fix. The original SUPG migration/partition sequence plus all four regressions passed on eight Mac ranks: 15 tests in 26.05 s, 32.30 s including launcher, 2.74 GiB peak process-tree RSS. The pre-fix sequence hung and the mesh-only reproducer also hung. Gadi rerun remains pending.
…#687)

Adapt only the mesh-size correction from lmoresi's 68e545f on feature/navier-stokes-supg; do not import Navier-Stokes or other branch changes. Cache _radii_own from current DM vertex coordinates and use it for mesh.cell_size(). Preserve the legacy kd-tree radius arrays and global timestep/mesh-motion consumers.

Use coordinate-section offsets and the full vertex stratum so the own-cell RMS definition also handles hexahedra, which have eight vertices but six faces. Correct the field documentation and Nitsche mechanism tests for the new definition; retain physical solve tolerances and use the exact nearest-centroid <= own-centroid ordering instead of an arbitrary approximate-equality tolerance.

Add a first-failing independent geometry/deformation regression for triangles, tetrahedra, quadrilaterals and hexahedra plus a regular-square analytical control. Before: four failures in serial and on eight ranks. After rebuild: 21 passed/one expected skip serial (22.90 s), 22 passed on eight ranks (40.45 s), covering Nitsche solves, radius accessors, frozen PC2 migration and memory/disk snapshots. Own-cell geometry error is zero in these tests; style and whitespace gates pass.
Rename the new per-cell geometric radius cache from _radii_own to _cell_radii so the name describes cell geometry rather than rank ownership. Update the focused Nitsche and deformation checks accordingly.\n\nAdd an enumerated parallel regression that gathers owned-cell centroid/radius pairs and compares the complete sorted table with a fresh single-rank run on the same cached Gmsh mesh. This directly guards the rank-count-independence claim at np=2, np=4 and np=8 instead of relying only on within-rank geometric identities.\n\nValidated locally with 9 focused serial tests and the new MPI test at 2, 4 and 8 ranks.
Remove the local_h=False workaround from the boundary-normal MPI regression now that Mesh.cell_size() is partition independent. The test again exercises the public local_h=True default and compares its Nitsche solve with a fresh serial process.\n\nRecord the user-visible consequence in the development changelog: the rank-local centroid kd-tree moved the default Nitsche velocity answer by 6.6e-3, while the cell-geometry replacement is identical cell by cell from one through eight ranks.\n\nValidated the focused Nitsche regression at 2, 4 and 8 Open MPI ranks (10.99 s, 7.31 s and 9.60 s respectively).
@gthyagi

gthyagi commented Sep 6, 2026

Copy link
Copy Markdown
Contributor Author

PR #689 now includes the complete commit stacks from #691 and #692 as explicit cherry-picks, without retargeting the PR or importing unrelated development commits.

Source Original commit Commit in #689
#691 8e93fc17 0e90dcb1
#692 443ffb58 0f6b78fa
#692 6fde1ac9 accf259b
#692 e6eaac28 f41bcd2f

Validation after rebuilding the combined branch:

Run Result Time
Focused serial cell-size, Nitsche and SUPG residual suite 24 passed 23.39 s
Combined 8-rank empty-partition, cell-size, SUPG partition and default-Nitsche batch 7 passed per rank 41.50 s

The 8-rank batch includes test_1077_advdiff_supg_parallel.py::test_error_is_partition_independent, which was the sole failure in the previous #689 CI run.

@gthyagi

gthyagi commented Sep 6, 2026

Copy link
Copy Markdown
Contributor Author

@lmoresi The remaining current-development integration gate is now complete.

I created a temporary merge of current upstream/development and the exact #689 head (f41bcd2f). The source merge was clean apart from retaining the current-development changelog side, the merged tree built successfully, and the supervised fresh-process disk-restart test passed all four public transport selections on eight MPI ranks:

Integration check Result Time
PC2, converged PC, CN, and BDF2 fresh-process restart, 8 ranks 4 passed 94.59 s

This exercises the exact same-layout field reload from #674 and MPI descendant supervision from #678 without duplicating either implementation in #689. The PR description and restart table now record this result. Upstream CI for the exact PR head is still running.

The remaining maintainer decision is the public API contract: preserve fixed two-correction time_integrator="citcoms" for Zhong/CitcomS compatibility, with its documented first-order nonuniform-diffusion limitation, and expose time_integrator="pc_converged" as the separately selected second-order reference. Please confirm whether that split is the intended public contract.

@gthyagi

gthyagi commented Sep 6, 2026

Copy link
Copy Markdown
Contributor Author

Final author-side gate update: upstream CI passed for the exact PR head. The PR is now clean and mergeable, and every author-action checklist item is complete. The requested maintainer decision on the public PC2 API remains pending review.

@lmoresi
lmoresi changed the base branch from feature/eulerian-supg-transport to development September 8, 2026 06:08
@lmoresi

lmoresi commented Sep 8, 2026

Copy link
Copy Markdown
Member

#688 is merged to development (1d5e8c1), so this PR is retargeted there and now shows conflicts: it was built on the closed #673 branch and its diff carries 65 uses of names that no longer exist (AdvDiffusionSUPG, SNES_AdvectionDiffusion_SUPG). On development the scalar solver is uw.systems.AdvDiffusion (the composing solver; the semi-Lagrangian class is AdvDiffusionSLCN), its transport is the uw.systems.ddt.EulerianSUPG manager (the time scheme, the SUPG knobs and the advecting velocity live there; the solver composes DuDt.time_derivative() + DuDt.advection() and DuDt.stabilisation_flux(R)), and docs/developer/design/eulerian-supg-transport.md, src/underworld3/analytic/* and the rotation-test example moved under it as well. A predictor-corrector time scheme most naturally lands as a manager (a subclass of EulerianSUPG overriding time_derivative/spatial_weights, or a pre-solve on it) rather than in the solver. Section "The DDt as the transport plugin" of the design note has the contract.

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

…manager

Resolve the PR underworldcode#688 architecture transition on existing PR underworldcode#689. Retain the upstream composing AdvDiffusion solver and EulerianSUPG implicit manager; move fixed and residual-converged predictor-corrector rate history, tau, stability estimate, corrections and reusable state into EulerianSUPGPC. Keep PDE assembly in the composing solver through small manager execution hooks.

Preserve matching-graph snapshots for manager history and solver timestep/change estimates. Reject cross-solver manager reuse; keep theta validation transactional and rewire changed velocity expressions. Retain explicit-tau anisotropic behavior. Exclude unrelated mantle Stokes, geoid and dependency changes.

Validate frozen-source equivalence, independent analytical recurrence and temporal order, live manager mutations, workspace reuse, and serial/eight-rank disk/fresh-process restarts. The focused MPI group passes 49 tests plus five formerly skipped disk replays; all four fresh-process MPI methods pass. CN and converged PC are order 2; fixed CitcomS remains approximately order 1. Full repository Level 1 run remains in progress and is reported separately.
@gthyagi gthyagi changed the title Integrate CitcomS predictor-corrector into unified Eulerian SUPG Add CitcomS predictor-corrector DDt manager for composed transport Sep 8, 2026
@gthyagi

gthyagi commented Sep 8, 2026

Copy link
Copy Markdown
Contributor Author

@lmoresi implemented your DDt-manager recommendation in cd155fc5 on this
existing PR branch, after merging development at 1d5e8c15. The conflicts
are resolved; this PR is mergeable and remains ready for review.

transport = uw.systems.ddt.EulerianSUPGPC(
    mesh, T, U.sym, method="citcoms", temperature_rate_field=Tdot,
)
thermal = uw.systems.AdvDiffusion(mesh, T, U.sym, DuDt=transport)
  • EulerianSUPGPC owns the rate, startup, correction loop, directional
    stabilization, stability estimate and restart state. It derives from
    _DDtBase and reuses the Eulerian SUPG advection/stabilization methods
    without allocating unused implicit history.
  • The composing solver still owns the PDE residual. Small optional manager
    execution/estimation hooks support rate-based corrections, since these
    require repeated residual assembly rather than one implicit SNES solve.
    CN/BDF continue through the existing EulerianSUPG path.
  • Removed obsolete public solver usage. The only old class reference is in
    an optional test loading a frozen copy of the previous PR source.
  • Added regressions for manager binding, failed theta mutation, changed
    velocity expressions and explicit-tau anisotropic diffusion. Snapshots
    capture both manager history and solver timestep/change estimates.
Validation Result Elapsed
Focused serial API/snapshots/mutations 33 passed 32.58 s
Serial time-order/restart/workspaces 20 passed 185.39 s
Eight-rank API/equivalence/time-order/workspaces/mutations 49 passed per rank 284.75 s
Eight-rank disk replay 5 passed per rank 24.60 s
Eight-rank fresh-process restart, four methods 4 passed 158.37 s
Full repository Level 1, serial 1797 passed; 41 skipped; 2 expected failures 1633.89 s

All disk-replay cases now run: an obsolete capability check was inspecting
the enhanced-variable wrapper signature and incorrectly skipping them.
PC fresh-process replays are exact; CN/BDF2 field differences are below
9e-16. Fixed CitcomS still measures approximately first order; converged
PC and actual CN measure second order in both 2-D and 3-D. We have not
relabelled fixed PC2 as second order or used a long A1 run as verification.

The PR description includes the extended analytical results, timing details
and checkpoint migration caveat. Please review the manager execution-hook
contract in particular: it keeps the numerical correction algorithm out of
the composing solver while allowing it to reuse solver-owned assembly.

@lmoresi

lmoresi commented Sep 8, 2026

Copy link
Copy Markdown
Member

Measured the manager on the rotating-Gaussian round trip we use for the transport schemes (one revolution, sigma 0.12, 32 cells across, serial; branch cd155fc on development 1d5e8c1). At its own stability step (dt 0.0130, Courant 0.60, 485 steps): citcoms (two passes) 51% round-trip error, peak 0.59 of 1, 0.030 s/step; pc_converged 18%, peak 0.77, 0.63 s/step. The implicit EulerianSUPG manager on the same P1 field at the same step: 6.7%, peak 0.94, 0.030 s/step. So the two-pass step costs what the implicit P1 solve costs here and is eight times less accurate; the converged mode reaches the same implicit system at twenty times the cost per step and three times the error, which I read as its steady tau (no transient cap, directional length) adding streamline diffusion. The rows are written up in a docs PR against development; driver and logs in ~/+Simulations/supg_vs_slcn_657/rotation_pc/ on this machine. One mesh, one flow; happy to be shown a case where the preset's CitcomS fidelity is the point.

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

@lmoresi

lmoresi commented Sep 8, 2026

Copy link
Copy Markdown
Member

Ruling from Louis on the comparison above: the predictor-corrector manager stays as a CitcomS-compatibility choice, not a general-purpose scheme, and we are not asking for P2 support (P2 halves its explicit step and the two-pass error comes from not converging the implicit system, not from the element order). For the PR that means: keep citcoms as the preset name for the CitcomS settings, say in the class docstring and the user page that it exists to reproduce CitcomS behaviour and that uw.systems.AdvDiffusion with its default manager is the recommended scheme, and credit the method to Brooks and Hughes (the predictor-multicorrector of Brooks' thesis and the 1982 SUPG paper, carried through ConMan and CitcomS) rather than to CitcomS alone.

🤖 Generated with Claude Code

https://claude.ai/code/session_018T2VHUGaZiQVJ95qQ4DiSL

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants