SUPG and Levelset #657

Activity

1. NengLu

   added the

   Feature

   issue type on Aug 27, 2026
2. NengLu

   added

   enhancementNew feature or requestNew feature or request

   on Aug 27, 2026
3. lmoresi commented on Sep 3, 2026

   

   lmoresi

   on Sep 3, 2026

   Member

   More actions

   Thanks for this, NengLu. We measured the levelset branch against SLCN and then built the production version of the solver on top of your weak-form idea. PR #673 has the result; the numbers and the code findings are here so the discussion stays with the issue.

   What the SUPG prototype gets right. The F1 += tau u R form works exactly as advertised on PetscDS with P2 elements, and the class runs unchanged on current development. On a rotating Gaussian at Courant 0.5 it is far more accurate than SLCN on the same mesh (0.6% against 21% round-trip error after one revolution) and 4 to 6 times cheaper per step. On your LeVeque level-set comparison (T = 2, 64 by 64) it gives a shape error of 0.028 against 0.051 for SLCN at half the wall time. The implicit scheme is stable at any cell Courant number: a band refined to h/9 across the orbit, at local Courant 13, changes the error in the third digit only. That is the property the transport design note asks for.

   What limits it. Accuracy is set by u dt against the feature width, as dt^2 for Crank-Nicolson: 9.8% at Courant 2, 66% with undershoots to -0.35 at Courant 8, while SLCN stays at 8% until its trace-back arc exceeds about 10 degrees per step. So the two schemes cross over near Courant 2 on the feature's own scale.

   Findings in the branch as pushed.

   1. level_set_SUPG.LevelSetSolver.solve calls self._adv_solver.solve(dt) positionally after the solver's solve became keyword-only, so the LeVeque script cannot run without a one-line fix.
   2. The discontinuity-capturing term's "residual" is u . grad(phi_old) alone, without the time-derivative part, so it is O(1) on any cleanly moving feature and adds first-order diffusion of size h|u| everywhere. With the default coefficient it erodes the Gaussian peak to 0.34 and loses 28% of the mass at Courant 2. Reinitialisation hides this on a tanh level set; on temperature it does not.
   3. dt is captured as a float literal in the kernels, so every change of timestep recompiles (1.2 s against 0.03 s per step). The fix is a UW expression for the timestep, which is what the SLCN solver does; PR Eulerian SUPG advection-diffusion with symbolic BDF and Adams-Moulton orders (#657 follow-on) #673 does that.
   4. Crank-Nicolson only, with a hand-rolled phi_old. PR Eulerian SUPG advection-diffusion with symbolic BDF and Adams-Moulton orders (#657 follow-on) #673 builds the scheme from the Eulerian DDt history instead, which gives BDF 1-3 and Adams-Moulton 1-3 from one code path; the study there shows Adams-Moulton above order 1 blows up from Courant 1 on advection and BDF2 is the robust default, with Crank-Nicolson 3 to 4 times more accurate below Courant 2.
   5. shapely is a hard import of the level-set modules and is not in the pixi environment; the two level-set files duplicate about 700 lines; the scripts live under src/ and use argparse where the repo uses uw.Params.

   The level-set wrappers are not in #673. Once it merges, LevelSetSolver can take the new class with integrator="am", order=1, theta=0.5 to keep Crank-Nicolson behaviour, and the discontinuity-capturing term can be redone with the full lagged residual (it needs phi at n-1, which the order-2 history provides).

   Measurement scripts and CSVs: ~/+Simulations/supg_vs_slcn_657/ on the ANU machine; design note in the PR at docs/developer/design/eulerian-supg-transport.md.
4. lmoresi commented on Sep 3, 2026

   

   lmoresi

   on Sep 3, 2026

   Member

   More actions

   Follow-on for the level-set part: PR #677 brings your level-set pipeline onto the shared solver interface as uw.systems.LevelSetSolver(psi, velocity=..., epsilon=..., advection='supg'|'slcn'), with your commits kept at the base of the branch and credit in the description, module and docs. It fixes the positional solve(dt) call, makes shapely optional, retires the prototype SUPG solver in favour of the one in #673 (which reproduces it to four digits), moves the LeVeque comparison to the examples, and adds a rotation test for both transport choices. It is opened against #673 and will be retargeted at development once that merges; comments on the API naming welcome.

   👍React with 👍1Reacted by Neng Lu
5. gthyagi

   added a commit that references this issue on Sep 5, 2026

   Eulerian advection-diffusion with SUPG: BDF and Adams-Moulton orders …

   ...

   b686530
6. lmoresi

   added a commit that references this issue on Sep 8, 2026

   Navier-Stokes with Eulerian SUPG momentum transport, and a partition-…

   ...

   1d5e8c1
7. lmoresi

   added a commit that references this issue on Oct 23, 2026

   Merge pull request #677 from underworldcode/feature/levelset-supg

   ...

   7222caf
8. lmoresi

   added

   fixed-in-developmentFix merged to development; the issue closes when it is released to mainFix merged to development; the issue closes when it is released to main

   on Oct 28, 2026
9. lmoresi commented on Oct 28, 2026

   

   lmoresi

   on Oct 28, 2026

   Member

   More actions

   Delivered by PR #677, merged to development as 7222caf3 on 2026-10-23.
   Labelling fixed-in-development.

   Probed rather than taken on the PR's word. src/underworld3/level_set.py and
   tests/test_1105_levelset_rotation.py were verified byte-identical to
   development's copies first (blob 1f9cedcc for the module), so this run
   exercises development's level set and not a variant of it:

   test_circle_returns_after_one_revolution[supg]                PASSED
   test_circle_returns_after_one_revolution[slcn]                PASSED
   test_the_geometry_helpers_refuse_clearly_without_shapely      PASSED
   test_initialise_from_polygon                                  SKIPPED  (shapely absent)
   test_material_property_field_blends                           PASSED
   test_sync_velocity_from_reaches_the_solve_time_vector         PASSED
   test_sync_velocity_from_refuses_a_bare_expression             PASSED
   
   6 passed, 1 skipped in 260.84s
   

   uw.systems.LevelSetSolver and uw.level_set both import and the rotating
   circle returns after one revolution under each of the two transports, which is
   what this issue asked for.

   Two things we are recording rather than glossing. The single skip is
   interface_geometry="polygon", which needs shapely — a package we never
   declare, so that path is uncovered on every run; the behaviour when it is
   missing is tested. Filed as #864, which also covers a worse instance of the
   same pattern. And level_set.py:372 builds its advection solver through
   systems.AdvDiffusionSLCN, so under the transport renaming in #800 it will
   emit a FutureWarning until it is migrated.

   Left open under the standing policy: issues close when development is
   released to main, not on merge.