Parallel bug - Calculating a boundary integral causes the next Stokes solve to hang

[jcgraciosa]

Describe the bug
If we calculate a boundary integral then do a Stokes solve afterwards, the Stokes solve "hangs."

• What version code?
  development branch

• Steps to reproduce?
  Simple code below; bug can be reproduced with nproc >= 2:

import sys
import petsc4py
petsc4py.init(sys.argv)

import underworld3 as uw
from underworld3.systems import Stokes
import numpy as np
import sympy

Ra       = 1e4
k        = 1.0
viscosity = 1.0
tol      = 1e-6

# --- Mesh ---
meshbox = uw.meshing.UnstructuredSimplexBox(
    minCoords=(0.0, 0.0),
    maxCoords=(1.0, 1.0),
    cellSize=1.0 / 8,
    qdegree=3,
)


x, z = meshbox.X

# --- Variables ---
v_soln = uw.discretisation.MeshVariable("U", meshbox, meshbox.dim, degree=2)
p_soln = uw.discretisation.MeshVariable("P", meshbox, 1, degree=1, continuous=True)
t_soln = uw.discretisation.MeshVariable("T", meshbox, 1, degree=3)

meshbox.dm.view()
#meshbox.view()
#meshbox.view_parallel()

# --- Stokes solver ---
stokes = Stokes(meshbox, velocityField=v_soln, pressureField=p_soln)
stokes.tolerance = tol
stokes.petsc_options["snes_converged_reason"] = None
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.constitutive_model.Parameters.shear_viscosity_0 = viscosity
stokes.saddle_preconditioner = 1.0 / viscosity
stokes.add_dirichlet_bc((sympy.oo, 0.0), "Bottom")
stokes.add_dirichlet_bc((sympy.oo, 0.0), "Top")
stokes.add_dirichlet_bc((0.0, sympy.oo), "Left")
stokes.add_dirichlet_bc((0.0, sympy.oo), "Right")
stokes.bodyforce = sympy.Matrix([0, Ra * t_soln.sym[0]])

# --- Initial temperature: linear + sinusoidal perturbation ---
pertStrength = 0.1
T_init = (1.0 - z) + pertStrength * sympy.cos(sympy.pi * x) * sympy.sin(sympy.pi * z)
with meshbox.access(t_soln):
    t_soln.data[:, 0] = np.clip(
        uw.function.evaluate(T_init, t_soln.coords).squeeze(), 0.0, 1.0
    )

# this will reproduce the error with the Stokes solve not finishing
uw.pprint("Initial boundary integral ... ")
lw_int = uw.maths.BdIntegral(meshbox, t_soln.sym[0], "Bottom").evaluate()
uw.pprint("Solving stokes ... ")
stokes.solve(zero_init_guess=True)

[lmoresi]

@jcgraciosa — PR #99 (merged to development) fixes a related BdIntegral MPI issue: ranks with empty local boundary strata now correctly participate in the Allreduce instead of skipping it, which was causing hangs.

Could you test your reproduction script against current development to see if #99 resolves this? The fix is in UW_DMPlexComputeBdIntegral in petsc_compat.h — it now checks DMLabelGetStratumSize before entering the assembly path.

If the hang persists, the root cause may be different (e.g. the Stokes DM/DS state being corrupted by the BdIntegral evaluation). Let us know either way.

[lmoresi]

@jcgraciosa — PR #99 (merged to development) fixes a related BdIntegral MPI issue: ranks with empty local boundary strata now correctly participate in the Allreduce instead of skipping it, which was causing hangs.

Could you test your reproduction script against current development to see if #99 resolves this? The fix is in UW_DMPlexComputeBdIntegral in petsc_compat.h — it now checks DMLabelGetStratumSize before entering the assembly path.

If the hang persists, the root cause may be different (e.g. the Stokes DM/DS state being corrupted by the BdIntegral evaluation). Let us know either way.

[jcgraciosa]

Thanks for letting me know of the potential fix.
I re-checked this with the latest development branch, but the hang still persists with nproc >= 2.

[lmoresi]

@jcgraciosa - See PR #115 for a potential fix for the problem.