[BUG] - Navier Stokes solver hangs in parallel if the domain has negative coordinates

[jcgraciosa]

Describe the bug

• What is the bug?
  The Navier-Stokes solver hangs during solve when run in parallel if the domain has negative coordinates. This is the root of Salam's reported issue in Kaiju.
• What version of the code?
  Latest development branch
• What system are you running the code on?
  Can be replicated in a local Mac installation or in Kaiju
• Steps to reproduce?
  This stripped version of Couette flow can be used to replicate the issue:

import sympy
import underworld3 as uw

Re         = 100
resolution = 16
U_wall     = 1.0
rho        = 1.0
nu         = 1.0 / Re
mu         = rho * nu

width       = 4.0
height      = 1.0

# will cause solve() to hang
minX, maxX  = -width/2., width/2.
minY, maxY  = -height/2., height/2.

# will NOT cause solve() to hang
#minX, maxX  = 0., width
#minY, maxY  = 0., height
#minX, maxX  = 1., 1. + width
#minY, maxY  = 0.5, 0.5 + height

uw.pprint(minX, minY)
uw.pprint(maxX, maxY)

mesh = uw.meshing.UnstructuredSimplexBox(
    minCoords=(minX, minY),
    maxCoords=(maxX, maxY),
    cellSize=1. / resolution,
    qdegree=3,
    regular=False,
)
mesh.dm.view()

v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2)
p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=False)

ns = uw.systems.NavierStokesSLCN(
    mesh, velocityField=v, pressureField=p, rho=rho, order=1
)
ns.constitutive_model = uw.constitutive_models.ViscousFlowModel
ns.constitutive_model.Parameters.viscosity = mu
ns.bodyforce = sympy.Matrix([0, 0])

ns.add_dirichlet_bc((U_wall, 0.0), "Top")
ns.add_dirichlet_bc((0.0,   0.0), "Bottom")

ns.tolerance = 1.0e-6

ns.petsc_options["snes_monitor"] = None
ns.petsc_options["snes_converged_reason"] = None
ns.petsc_options["ksp_monitor"] = None

dx = mesh.get_min_radius()
dt_adv = 1.0 * dx / U_wall
dt_diff = 1.0 * dx**2 / nu
dt = min(dt_adv, dt_diff)

uw.pprint(f"dt: {dt}")

ns.solve(timestep=dt, zero_init_guess=True)

The issue appears with at least 3 procs:

pixi run mpirun -n 3 python3 couette-stripped.py

[lmoresi]

Update: original repro fixed, but a narrower case survives

Re-investigated on current development (np=3). The original repro in this issue — UnstructuredSimplexBox(regular=False) with negative coordinates — now completes cleanly (5 timesteps, no hang). It was fixed by #222 (global_evaluate: faithful parallel evaluate() — out-of-domain mislocation fix, merged 2026-06-12, after this issue was filed). The NS semi-Lagrangian trace-back evaluates upstream velocity via global_evaluate, and #222 corrected the parallel point-location of out-of-domain trace-back points.

However, a related case still hangs: StructuredQuadBox + negative coordinates. A single solve() never returns (confirmed with mpirun --timeout 240 → hard timeout; "entering solve" printed, "SOLVE RETURNED" never reached).

Reproduction matrix (np=3, current development)

mesh	coords	result
UnstructuredSimplexBox(regular=False)	negative	✅ completes (this issue's original repro)
UnstructuredSimplexBox(regular=True)	positive	✅ completes
StructuredQuadBox	positive	✅ completes
StructuredQuadBox	negative	❌ hangs (single solve never returns)

So the hang now requires StructuredQuadBox AND negative coordinates specifically (quad+positive works; simplex+negative works).

Minimal repro

import sympy, underworld3 as uw
mesh = uw.meshing.StructuredQuadBox(
    minCoords=(-2.0, -0.5), maxCoords=(2.0, 0.5), elementRes=(16, 16), qdegree=3)
v = uw.discretisation.MeshVariable('U', mesh, mesh.dim, degree=2)
p = uw.discretisation.MeshVariable('P', mesh, 1, degree=1, continuous=False)
ns = uw.systems.NavierStokesSLCN(mesh, velocityField=v, pressureField=p, rho=1.0, order=1)
ns.constitutive_model = uw.constitutive_models.ViscousFlowModel
ns.constitutive_model.Parameters.viscosity = 0.01
ns.bodyforce = sympy.Matrix([0, 0])
ns.add_dirichlet_bc((1.0, 0.0), 'Top'); ns.add_dirichlet_bc((0.0, 0.0), 'Bottom')
dx = mesh.get_min_radius()
ns.solve(timestep=min(dx, dx**2/0.01), zero_init_guess=True)  # hangs at np>=3

mpirun -np 3 python repro.py

Likely locus

The hang is inside the first solve() (the SL DuDt trace-back / global_evaluate point-location). Since the same global_evaluate path works for simplex+negative, the quad-specific difference (StructuredQuadBox DM/coordinate layout vs simplex) is where the remaining out-of-domain point-location deadlocks under MPI. Re-scoping this issue to the StructuredQuadBox + negative-coordinate case.

Underworld development team with AI support from Claude Code

[lmoresi]

#211 — status on current development

Environment

• Branch: development @ a7f0b11f (post-Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265, pre-[BUG] - Navier Stokes solver hangs in parallel if the domain has negative coordinates #211-fix)
• PETSc 3.25.2 (conda-forge, dev env with the bumped 3.25 pin)
• MPICH 5.0.1

Reproduction matrix

Setup: StructuredQuadBox(elementRes=(16,16)), NavierStokesSLCN, timestep
min(dx, dx²/ν), single solve(). Coords "pos" = [0,4]×[0,1]; "neg" =
[-2,2]×[-0.5,0.5].

np	pos coords	neg coords
2	✅ SOLVE RETURNED (7.1 s)	❌ MPI_Allreduce: Message truncated
3	✅ SOLVE RETURNED (4.5 s)	❌ *** stack smashing detected *** (SIGABRT)
4	✅ SOLVE RETURNED (4.5 s)	❌ MPI_Allreduce: Message truncated

So on current development the bug is still present at every parallel np; the
np=3 failure is the memory-corruption twin of the np=2 / np=4 MPI truncation
(same PETSc call, same root cause; np=3 loses the race in DMPlex before MPI
catches it).

Failure trace (np=4, negative coords)

DMInterpolationSetUp_UW   (src/underworld3/function/petsc_tools.c:108)
 → DMLocatePoints         (dm/interface/dmcoordinates.c:1273)
 → DMLocatePoints_Plex    (dm/impls/plex/plexgeometry.c:1300)
 → DMGetBoundingBox       (dm/interface/dmcoordinates.c:994)
 → MPIU_Allreduce(lmin, gmin, cdim=2, MPIU_REAL, MPIU_MIN, dm->comm)
       ↳ MPI error: "Message truncated"

DMInterpolationSetUp_UW is called from
petsc_interpolate.interpolate_vars_on_mesh (_function.pyx:1301) as part of
global_evaluate → the semi-Lagrangian trace-back
(ddt.SemiLagrangian.update_pre_solve → global_evaluate).

Why this is quad-specific (not simplex)

The evaluator switches between two point-location paths via
Mesh._eval_use_robust_location()
(src/underworld3/discretisation/discretisation_mesh.py:5482):

• parallel AND dm.isSimplex() OR manifold (dim != cdim) →
  barycentric-hint bypass in petsc_tools.c (line 71:
  if (owning_cell) …), skips DMLocatePoints entirely;
• otherwise → DMLocatePoints path (line 99: else …), which hits the
  DMGetBoundingBox Allreduce.

StructuredQuadBox is a non-simplex volume mesh, so it always falls through
to the second path. #222 fixed the simplex path — that is why simplex + neg
now works and quad + neg still fails.

Where the fix has to live

Two options; neither is a one-liner.

1. Widen _eval_use_robust_location() to accept StructuredQuadBox. Tried
   with a monkey-patch — the barycentric hint (_robust_owning_cells) is
   simplex-only and the run crashes with stack smashing detected on quad
   cells (heap overrun in a barycentric helper). Doing this properly means
   teaching the hint code to handle quad reference maps, or supplying a
   different cell-containment hint for tensor-product cells (parametric
   inversion + interval test in the reference cube). Substantial change.

2. Root-cause the PETSc DMGetBoundingBox Allreduce on a negative-coord
   DMPlex. The Allreduce itself is (cdim=2, MPIU_REAL, MPIU_MIN) on the
   DM comm — the count/type are consistent across ranks, so "Message
   truncated" is almost certainly a corrupted sendbuf/recvbuf from
   DMGetLocalBoundingBox (rank-local step, immediately before) rather than
   an MPI collective inconsistency. Likely candidates: a stale
   coordinateDM or Vec with the wrong local size after a coord shift
   into the negative octant during DMPlex construction. A minimal PETSc-
   only reproducer (build a 2D DMPlex on [-2,2]×[-0.5,0.5], distribute,
   call DMGetBoundingBox) would isolate this cleanly and is worth filing
   upstream.

Interim workarounds for users

• Shift the mesh into the positive octant (as the issue already notes) — the
  physics is translation-invariant here.
• Use UnstructuredSimplexBox instead of StructuredQuadBox — global_evaluate: faithful parallel evaluate() (fix out-of-domain mislocation) #222 covers
  the simplex path.
• Serial (np=1) is unaffected.

Suggested next step

I'd recommend attacking option 2 first: write a PETSc-only reproducer (no
UW3), confirm the bounding-box Allreduce misbehaves for negative-coord DMPlex
in parallel, and file it upstream. That unblocks both quad and any future
non-simplex volume mesh path. Option 1 (extend the bypass to quads) is worth
having anyway but is a larger change.

[lmoresi]

Narrowing the diagnosis: getBoundingBox alone is not the culprit

Attempted a petsc-only reproducer for the MPI_Allreduce truncation my earlier
comment traced to DMLocatePoints_Plex → DMGetBoundingBox → MPI_Allreduce.
Two negative results, both worth recording so the search does not repeat:

1. A plain petsc4py DMPlex.createBoxMesh(simplex=False) on [-2,2]×[-0.5,0.5] at np=3 does not crash

from petsc4py import PETSc
dm = PETSc.DMPlex().createBoxMesh(faces=(16,16),
                                  lower=(-2.0,-0.5), upper=(2.0,0.5),
                                  simplex=False, interpolate=True,
                                  comm=PETSc.COMM_WORLD)
dm.distribute()
dm.getBoundingBox()   # ((-2.0, 2.0), (-0.5, 0.5)) — no MPI error

2. UW3's StructuredQuadBox with the exact #211 configuration at np=3 also does not crash on getBoundingBox

mesh = uw.meshing.StructuredQuadBox(
    minCoords=(-2.0, -0.5), maxCoords=(2.0, 0.5),
    elementRes=(16, 16), qdegree=3)
mesh.dm.getBoundingBox()   # ((-2.0, 2.0), (-0.5, 0.5)) — no MPI error

Ran at np=3 (the reported failure size) — dm.distribute() returned the
migration SF, distribution went through cleanly, bounding box came back
with the expected negative-corner bounds.

What this tells us

The MPI truncation in the crash trace is not produced by getBoundingBox
in isolation — it is triggered by something upstream in DMLocatePoints_Plex
that runs before, alongside, or after the bounding-box collective. Candidates
(from reading src/dm/impls/plex/plexgeometry.c):

• the cellSF build immediately following the bounding-box query,
• a PetscSF reduction on the point-location results,
• or the cell-containment fallback path (BVH search / midpoint refinement)
  when a query point sits near a partition seam.

The fix on the UW side stays the same as my earlier comment recommended:
widen Mesh._eval_use_robust_location to route non-simplex volume meshes
through the barycentric-hint bypass (which skips DMLocatePoints entirely).
That is a real change — _robust_owning_cells currently stack-smashes on
quad reference maps — but it is the right long-term fix.

An upstream PETSc filing needs a C reproducer (petsc4py does not expose
DMLocatePoints on DMPlex in any dev-env build I have here), which is
a bigger task and probably belongs after we have a concrete failing pointset
that survives isolation from UW3.

Underworld development team with AI support from Claude Code

[lmoresi]

Not reproducible on current development

Re-ran the surviving case from the narrowing above — StructuredQuadBox,
NavierStokesSLCN, one solve(), sign of the domain the only variable — on
development @ 5998283. The positive domain is the control.

np	∫ v·v (pos, [0,4]×[0,1])	∫ v·v (neg, [-2,2]×[-0.5,0.5])
1	5.416053533624e-02	5.416053533629e-02
2	5.416053533321e-02	5.416053533326e-02
3	5.416053532628e-02	5.416053532634e-02
4	5.416053534602e-02	5.416053534607e-02

Every run returns. Positive and negative agree to ~1e-12 relative at every rank
count, and the spread across rank counts is 1.6e-10 — the solve is
translation-invariant, which is the property this issue reports as broken.
Against the July matrix: MPI_Allreduce: Message truncated at np=2 and np=4,
*** stack smashing detected *** at np=3, on every negative run.

One trap worth recording, since it nearly produced a second bug report: the
global vector norm v.vec.norm() reads 6.0039e+00 at np=2 and np=3 but
4.3093e+00 at np=4 on the answers tabulated above. That is the norm counting
shared DOFs once per rank that holds them. Partition-independence has to be
judged on integrals.

What we cannot say

We cannot name the fix. The diagnosis in the comment above reached
DMLocatePoints_Plex → DMGetBoundingBox → MPI_Allreduce, and #556 rewrote that
entry path (reject foreign points before the locator, locate once,
containment-check the cell hint) — but the environment has also moved
underneath: MPI is now Open MPI 5.0.10 where the failing matrix ran MPICH, and
the failure was an MPI collective truncation. "Does not reproduce here" is
therefore weaker than "the defect is understood and removed".

Closing on the measurement. @jcgraciosa — if this survives on Kaiju, please
reopen and name that environment; the script and the full matrix are at
~/+Simulations/issue211_negative_coords/.

Underworld development team with AI support from Claude Code