Fix SL trace-back FE overshoot + add monotone limiter to DDt#186
Merged
Conversation
…an DDt
The SL trace-back in SemiLagrangian.update_pre_solve was calling
uw.function.global_evaluate(...) at two sites (velocity midpoint and
psi_star endpoint) without propagating its enclosing evalf kwarg, so
the trace-back always used FE Lagrange shape-function interpolation
regardless of any user-facing RBF flag. P3 Lagrange overshoots at
non-nodal upstream points in cells with sharp gradients (Runge-like
local non-monotonicity), producing psi_star values outside the data
range that feed the implicit SNES and ignite catastrophic ringing on
free-surface convection at high Ra.
Two changes in ddt.py SemiLagrangian:
1. Plumb evalf through both global_evaluate calls (velocity midpoint
and psi_star endpoint). When evalf=True the trace-back now routes
through Shepard 3-NN inverse-distance interpolation (bounded by
neighbour values), as the calling solver could already request via
AdvDiffusionSLCN.solve(_evalf=True).
2. Add monotone_mode instance attribute + per-call kwarg with two
limiter options applied after the trace-back evaluation:
"clamp" (B.2): clip the trace-back result to [nbr_min, nbr_max]
of the k=dim+1 nearest psi_star DOFs. Bit-identical
to FE in smooth regions; saves the run at the
catastrophe step. Cost: one extra kdtree query per
update_pre_solve.
"pick" (B.1): keep the FE result if it lies within nbr bounds,
else re-evaluate via RBF at out-of-bounds DOFs.
More conservative than clamp; doubles trace-back
cost when triggered.
Default (None) is unchanged behaviour (pure FE, can overshoot).
Settable via ``adv_diff.DuDt.monotone_mode = "clamp"`` (same
pattern as ``theta``).
Tested against the RK4-full free-surface convection canonical case
(Ra=1e5, ρg=1e5, P3 T, CN diffusion, 35 steps; default FE rings
catastrophically at step 35 with post_raw=[-1.12, +2.56], 432 bad
DOFs). With ``monotone_mode = "clamp"`` the run is bit-identical to
FE through step 30 (no overshoots, no DOFs out of nbr bounds; the
clamp is a no-op), then at step 35 it produces a clean physical
solution: post_raw=[0, 1.000], Nu=143 (vs FE artefact 245), vrms=81,
h_max=0.151. Wall time per step unchanged (~19s).
Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR updates the semi-Lagrangian DDt trace-back in SemiLagrangian_DDt to (a) correctly propagate the caller’s evalf choice into global_evaluate calls, and (b) add an optional monotonicity limiter (monotone_mode) to prevent FE overshoot at sharp gradients.
Changes:
- Plumb
evalf=evalfthrough the velocity midpoint and upstreampsi_starglobal_evaluatecalls inupdate_pre_solve. - Add
monotone_modeas an instance attribute and per-call kwarg, with"clamp"and"pick"limiter implementations applied after trace-back evaluation. - Implement neighbour-range detection via KDTree over
psi_starDOF coordinates to compute local[nbr_min, nbr_max]bounds.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+2020
to
+2033
| # Plain-numpy coords for kdtree (handle pint/unit-aware) | ||
| if hasattr(end_pt_coords, "magnitude"): | ||
| epc_nd = np.asarray(end_pt_coords.magnitude) | ||
| else: | ||
| epc_nd = np.asarray(end_pt_coords) | ||
| psi_coords_nd = np.asarray(self.psi_star[i].coords_nd) | ||
| if hasattr(psi_coords_nd, "magnitude"): | ||
| psi_coords_nd = np.asarray(psi_coords_nd.magnitude) | ||
| nnn = self.mesh.dim + 1 | ||
| kdt = uw.kdtree.KDTree( | ||
| np.ascontiguousarray(psi_coords_nd)) | ||
| _, idxs = kdt.query( | ||
| np.ascontiguousarray(epc_nd), k=nnn, | ||
| sqr_dists=False) |
Comment on lines
+2047
to
+2053
| # B.1 "pick": re-evaluate via RBF where FE was out of bounds. | ||
| value_rbf = uw.function.global_evaluate( | ||
| expr_to_evaluate, end_pt_coords, evalf=True) | ||
| vrbf_flat = np.asarray(value_rbf).reshape(nbr_min.shape) | ||
| out_of_bounds = ((veep_flat < nbr_min) | ||
| | (veep_flat > nbr_max)) | ||
| veep_lim = np.where(out_of_bounds, vrbf_flat, veep_flat) |
1. Non-dimensionalise end_pt_coords for the kdtree query. In unit-aware runs end_pt_coords is dimensional (e.g. metres) while psi_star[i].coords_nd is [0,1] non-dimensional, so without explicit conversion the kdtree picks wrong neighbours and the limiter would clamp/pick against unrelated DOFs. Match the same conversion pattern used earlier in the function for psi_star_0_coords → psi_star_0_coords_nd. 2. In pick mode, only evaluate the RBF for the subset of upstream coords whose FE result is out-of-bounds. The previous code unconditionally ran a second global_evaluate(..., evalf=True) over ALL coords, paying ~2x trace-back cost on every call even when zero DOFs needed re-evaluation. Now the cost is dominated by the cheap FE pass when most DOFs are in-bounds. Underworld development team with AI support from Claude Code
lmoresi
added a commit
that referenced
this pull request
May 14, 2026
2 tasks
jcgraciosa
pushed a commit
to jcgraciosa/underworld3
that referenced
this pull request
May 20, 2026
PR underworldcode#186 landed the monotone limiter on SemiLagrangian_DDt. To use it via the high-level solver, users had to construct the solver and then mutate the DDt instances after the fact: adv_diff = uw.systems.AdvDiffusionSLCN(mesh, u_Field=T, V_fn=v.sym) adv_diff.DuDt.monotone_mode = "clamp" adv_diff.DFDt.monotone_mode = "clamp" This change adds a one-line constructor idiom: adv_diff = uw.systems.AdvDiffusionSLCN( mesh, u_Field=T, V_fn=v.sym, monotone_mode="clamp") The new kwarg forwards to both internally-constructed SemiLagrangian DDts (DuDt and DFDt). When the caller supplies a pre-built DuDt explicitly, that DDt's monotone_mode is preserved as the source of truth — the kwarg only takes effect for the internally-constructed DFDt in that case (documented in the class docstring). Default unchanged (None = pure FE trace-back). No behaviour change for code that doesn't pass the kwarg. Regression test in tests/test_1054_advdiff_monotone_mode_kwarg.py covers: default None, clamp forwarded to both DuDt and DFDt, pick forwarded to both, and the explicit-DuDt-preserves-its-own-mode case. Underworld development team with AI support from Claude Code
Merged
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.
Summary
evalfkwarg through bothglobal_evaluatecalls insideSemiLagrangian.update_pre_solve(velocity midpoint and psi_star endpoint). Previously the calling solver's_evalf=Truehad no effect on the trace-back — it always used FE Lagrange shape functions regardless.monotone_modeinstance attribute + per-call kwarg onSemiLagrangian_DDtwith two limiter options applied after the trace-back evaluation. Default unchanged (None= pure FE).Background
P3 (and higher) Lagrange FE shape functions are non-monotone: they can produce values outside
[data_min, data_max]at non-nodal points inside a cell with sharp gradient. The SL trace-back samplespsi_star.symat upstreamend_pt_coords— exactly the non-nodal points where the overshoot fires. On free-surface convection at Ra ≥ 1e5 the resulting psi_star values outside[0, 1]feed the implicit SNES and ignite catastrophic ringing (T → [-1.12, +2.56] in one step at step 35 of the canonical RK4-full case, with 432 DOFs out of bounds and the mode-5 plume structure destroyed by global noise pepper).A
--rbf-advection-style flag had been wired into the calling solver but never reached the actual trace-back, so the failure was opaque until we instrumented psi_star pre/post the trace-back call and saw it producing values like[0, +1.1473]from input data strictly in[0, 1].The two limiter options
None(default)\"clamp\"(B.2)[nbr_min, nbr_max]of thek = dim + 1nearest psi_star DOFs. Bit-identical to FE in smooth regions (clamp is a no-op); saves the run at the catastrophe step.update_pre_solve\"pick\"(B.1)Set via the same pattern as
theta:Test (manual, free-surface convection)
Canonical RK4-full free-surface convection: Ra=1e5, ρg=1e5, P3 T, CN diffusion, 35 steps. Default FE behaviour rings catastrophically at step 35.
monotone_mode=None)[-1.12, +2.56]monotone_mode=\"clamp\"[0, +1.000]monotone_mode=\"pick\"[0, +1.000]For steps 5–30 the
\"clamp\"mode is bit-identical to pure FE (no DOFs out of nbr bounds, clamp is a no-op), so accuracy in the smooth regime is preserved.Test plan
monotone_mode=\"clamp\"and confirms the canonical case completes without ringingFiles changed
src/underworld3/systems/ddt.py(+81 lines, no deletions)Underworld development team with AI support from Claude Code