SUPG and Levelset

[NengLu]

I've created a branch: https://github.com/underworldcode/underworld3/tree/levelset

which includes several tests in the test_tem directory. The SUPG code is saved in the solver_supg.py file.
And the levelset implementations are saved as level_set_SLCN and level_set_SUPG.

Please take a look and feel free to commit any changes. @lmoresi @gthyagi

[lmoresi]

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.

[lmoresi]

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.