[BUG] - Navier-stokes solution unstable after first timestep

[bknight1]

Describe the bug
The numerical solution for the Navier-Stokes solver is unstable after the first solve.

To debug, I created a Stokes solver which has been extended to include the NS terms. After debugging, I believe the error is due to the SemiLagrangian history manager for DuDt, as when switching to the Eulerian history manager the solution is stable. I believe the flux term may also need to be switched to Eulerian, although this hasn't been tested as in my current implementation I re-construct the flux term from the velocity history. Interestingly, the bug only appears at lower cell sizes (<0.1).

There is also an error from the Eulerian history manager during the pre-solve update, but doesn't seem to be altering the numerical solution:
Callback error in <function _BaseMeshVariable._create_canonical_data_array.<locals>.canonical_data_callback at 0x32a94eac0>: cannot reshape array of size 1969 into shape (2)
Screenshots
SL after first timestep

Eulerian after first timestep

[lmoresi]

@bknight1 — Thanks for the detailed diagnosis. The SemiLagrangian vs Eulerian comparison is useful, though we suspect the Eulerian DDt may be masking the problem through numerical diffusion rather than actually fixing the underlying issue (the Eulerian advection term is only a cursory implementation and effectively lowers the Reynolds number).

To dig into this properly, could you share:

1. A reproducer script (or point us to the test/notebook you're running) — the solver setup, mesh, BCs, and timestepping parameters
2. Which branch/commit you're testing against — there are DDt changes in progress from the viscoelastic work that may touch the same code paths
3. The timestep size and Reynolds number — to check whether the CFL condition is being respected for the SemiLagrangian transport

The resolution dependence (only cellSize < 0.1) is interesting — finer meshes mean more swarm particles and more migration events, which could point to a swarm data synchronisation issue in the SemiLagrangian history manager. The reshape error in the Eulerian path (cannot reshape array of size 1969 into shape (2)) also looks like a data cache bug worth tracking separately.

Underworld development team with AI support from Claude Code

[bknight1]

Sure, I've attached a script that allows changing the DuDt and DFDt that is passed into the NS solver to test the Eulerian V SemiLagrangian history managers.

NS_issue_example.py

Tested against the latest development branch using PETSc 3.24.4

[lmoresi]

Fix: SemiLagrangian departure point clamping (commit 4d522d8)

Root cause: The SemiLagrangian DDt computed upstream departure points (mid_pt_coords and end_pt_coords) without clamping them to the domain boundary. Near corners — especially the velocity discontinuities in the lid-driven cavity — departure points landed outside the mesh. global_evaluate then extrapolated using the P2 basis, producing values of order 1000x the actual field maximum. This injected a massive spurious acceleration that blew the solver up on step 2.

Fix: Clamp both mid_pt_coords and end_pt_coords using mesh.return_coords_to_bounds() in SemiLagrangian.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 test: 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, cellSize=0.04, t=5.0:

Re=400, cellSize=0.04, t=10.0:

(Images to be uploaded)

Next steps

@bknight1 — The fix is on development. Could you:

1. Validate with your original test case
2. Test with order=2 (second-order DuDt) — this is an important capability for NS
3. The Eulerian DDt callback reshape error (cannot reshape array of size 1969 into shape (2)) is a separate bug worth tracking

PR #83 has been closed since the Eulerian workaround is no longer needed.

Underworld development team with AI support from Claude Code