Repository navigation
SUPG and Levelset #657
Description
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
Activity
lmoresi commentedon Sep 3, 2026
Thanks for this, NengLu. We measured the
levelsetbranch 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 Rform works exactly as advertised on PetscDS with P2 elements, and the class runs unchanged on currentdevelopment. 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.
level_set_SUPG.LevelSetSolver.solvecallsself._adv_solver.solve(dt)positionally after the solver'ssolvebecame keyword-only, so the LeVeque script cannot run without a one-line fix.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.dtis 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.phi_old. PR Eulerian SUPG advection-diffusion with symbolic BDF and Adams-Moulton orders (#657 follow-on) #673 builds the scheme from theEulerianDDt 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.shapelyis 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 undersrc/and use argparse where the repo usesuw.Params.The level-set wrappers are not in #673. Once it merges,
LevelSetSolvercan take the new class withintegrator="am", order=1, theta=0.5to 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 atdocs/developer/design/eulerian-supg-transport.md.lmoresi commentedon Sep 3, 2026
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 positionalsolve(dt)call, makesshapelyoptional, 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.Eulerian advection-diffusion with SUPG: BDF and Adams-Moulton orders …
Navier-Stokes with Eulerian SUPG momentum transport, and a partition-…