Fix SemiLagrangian DDt blowup: clamp departure points to domain#85
Fix SemiLagrangian DDt blowup: clamp departure points to domain#85lmoresi merged 2 commits intodevelopmentfrom
Conversation
…#82) The SemiLagrangian update_pre_solve computed upstream departure points (mid_pt_coords, end_pt_coords) without clamping them to the domain. Near boundaries — especially corner singularities in lid-driven cavity flow — departure points landed outside the mesh, causing global_evaluate to extrapolate wildly (values of order 1000x the actual field maximum). This injected a massive spurious acceleration that blew up the solver on the second timestep. Fix: clamp both mid_pt_coords and end_pt_coords using the mesh's existing return_coords_to_bounds function. This mirrors the treatment already present in the Lagrangian DDt's swarm advection path. Validated against Ghia et al. (1982) lid-driven cavity benchmark at Re=100 (RMS error 0.004) and Re=400. Stable across all tested combinations of Re (10, 100, 400, 1000) and cellSize (0.1, 0.05, 0.025). Underworld development team with AI support from Claude Code
Clean rewrite of Ex_NS_PipeFlow.py: Poiseuille flow between parallel plates with analytic steady-state comparison. Replaces the unmigrated HPC batch script with a self-contained notebook following style guide. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
Fixes a SemiLagrangian DDt stability issue by preventing upstream departure points from leaving the mesh domain, and modernizes the Navier–Stokes pipe/Poiseuille flow example into a simpler standalone notebook-style script.
Changes:
- Clamp SemiLagrangian midpoint and end-point departure coordinates to mesh bounds before evaluation.
- Rewrite
Ex_NS_PipeFlowinto a streamlined transient Poiseuille-flow example with analytic comparison.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
src/underworld3/systems/ddt.py |
Adds domain clamping for semi-Lagrangian departure points to avoid out-of-mesh extrapolation blowups. |
docs/examples/fluid_mechanics/intermediate/Ex_NS_PipeFlow.py |
Replaces migrated legacy example with a simpler, self-contained Poiseuille-flow notebook example and plotting/comparison workflow. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # Clamp midpoint coordinates to the domain boundary | ||
| if self.mesh.return_coords_to_bounds is not None: | ||
| mid_pt_coords = self.mesh.return_coords_to_bounds(mid_pt_coords) |
| if step % 5 == 0 or step == NSTEPS - 1: | ||
| # Sample centreline velocity at the outlet (x = max, y = 0) | ||
| vx_max = np.abs(v_soln.data[:, 0]).max() | ||
| print(f"step {step:3d} | max |vx| = {vx_max:.4e}") |
| sample_coords = np.column_stack([np.full(n_sample, x_sample), y_sample]) | ||
|
|
||
| delta_t = Cmax*delta_x/max_vel | ||
| vx_numerical = uw.function.evaluate(v_soln.sym[0], sample_coords) |
| delta_x = meshbox.get_min_radius() | ||
| max_vel = vel | ||
| # Sample the velocity along a vertical line near the outlet | ||
| x_sample = 0.4 * ASPECT_RATIO * HEIGHT * 0.5 # near the outlet |
|
Here are two quick "benchmark" runs at Re=100, Re=400 for the Ghia et al benchmark. Qualitatively, this looks promising with the counterflow vortices developing at approximately the correct size ratios. @bknight1 - needs more testing, but merge if you think the main bugs are killed. 
|
|
Looks good to me! Thanks @lmoresi. A quick suggestion would be allowing the user to pass a DFDt history manager, as we currently enforce the SL version. We currently do this for all solvers that require the flux history. |
Order-1 and Order-2 validation: Lid-driven cavity Re=100Tested both Results
Order 1
Order 2
Note on order-2: This does not really look stable, so we will keep this open for now and keep digging Underworld development team with AI support from Claude Code |
Clamp mid_pt_coords and end_pt_coords to domain bounds in SemiLagrangian.update_pre_solve() to prevent extrapolation blowup near boundary discontinuities. Underworld development team with AI support from Claude Code
|
Update: Order-2 results with working BDF/AM coefficients The earlier order-2 results were invalid. We discovered that BDF/AM coefficients were frozen at order 1 from the first solve — This has been fixed in PR #88 (BDF/AM coefficients routed through 
Both orders are stable with the departure-point clamping fix from this PR. The order-2 instabilities seen previously were entirely the frozen-coefficients bug, not a problem with the clamping.
|
|
I think we can merge this fix. Test passing, comments from @bknight1 approving 1st order variant of the PR |
Summary
Root cause
The SemiLagrangian
update_pre_solvecomputed upstream departure points (mid_pt_coordsandend_pt_coords) without clamping them to the domain boundary. Near corners -- especially velocity discontinuities in the lid-driven cavity -- departure points landed outside the mesh.global_evaluatethen extrapolated using the P2 basis, producing values ~1000x the actual field maximum. This injected a spurious acceleration that blew up the solver on the second timestep.Fix
Clamp both
mid_pt_coordsandend_pt_coordsusingmesh.return_coords_to_bounds()inSemiLagrangian.update_pre_solve(). The Lagrangian DDt already did this for swarm advection; the SemiLagrangian (nodal) path was missing it. Two lines of code.Validation
Stability: All combinations stable (Re = 10, 100, 400, 1000; cellSize = 0.1, 0.05, 0.025).
Ghia et al. (1982) benchmark at Re=100 (cellSize=0.04, 200 steps, t=10): RMS error = 0.004 against published centreline velocity data.
Benchmark images (Re=100 and Re=400 streamlines) to follow in comments.
Test plan
Underworld development team with AI support from Claude Code