BUG: Changing boundary condition order results in divergence

Describe the bug When solving using a list of boundary conditions, changing their order in the list results in divergence.

Steps to Reproduce Minimal example:

import firedrake as fd

meshDim = 50
mesh = fd.UnitSquareMesh(meshDim,meshDim)

P2 = fd.VectorElement("P", mesh.ufl_cell(), 2)
P1 = fd.FiniteElement("P", mesh.ufl_cell(), 1)
TH = fd.MixedElement([P2, P1])

W = fd.FunctionSpace(mesh, TH)

# Define boundary conditions
bc_left   = fd.DirichletBC(W.sub(0), (0, 1), 1)
bc_right  = fd.DirichletBC(W.sub(0), (0, -1), 2)
bc_bottom = fd.DirichletBC(W.sub(0), (0, 0), 3)
bc_top    = fd.DirichletBC(W.sub(0), (0, 0), 4)

# bcW = [bc_left, bc_right, bc_top, bc_bottom]  # Working order of boundary conditions
bcW = [bc_bottom, bc_right, bc_top, bc_left]  # Diverging order of boundary conditions

v, q = fd.TestFunctions(W)
w    = fd.Function(W)
u, p = fd.split(w)

# Define the variational forms (Navier Stokes)
F = fd.Constant(1 / 50) * fd.inner(fd.grad(u), fd.grad(v)) * fd.dx \
    + fd.dot(fd.dot(fd.grad(u), u), v) * fd.dx \
    - p * fd.div(v) * fd.dx - q * fd.div(u) * fd.dx

# Solve the problem
fd.solve(F == 0, w, bcs=bcW)

Expected behavior I expect the order of the boundary conditions to not matter.

Error message

---------------------------------------------------------------------------
ConvergenceError                          Traceback (most recent call last)
Cell In[86], line 28
     23 F = fd.Constant(1 / 50) * fd.inner(fd.grad(u), fd.grad(v)) * fd.dx \
     24     + fd.dot(fd.dot(fd.grad(u), u), v) * fd.dx \
     25     - p * fd.div(v) * fd.dx - q * fd.div(u) * fd.dx
     27 # Solve the problem
---> 28 fd.solve(F == 0, w, bcs=bcW)

File petsc4py/PETSc/Log.pyx:115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/adjoint_utils/solving.py:57, in annotate_solve.<locals>.wrapper(*args, **kwargs)
     54     tape.add_block(block)
     56 with stop_annotating():
---> 57     output = solve(*args, **kwargs)
     59 if annotate:
     60     if hasattr(args[1], "create_block_variable"):

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/solving.py:128, in solve(*args, **kwargs)
    126 # Call variational problem solver if we get an equation
    127 if isinstance(args[0], ufl.classes.Equation):
--> 128     _solve_varproblem(*args, **kwargs)
    129 else:
    130     # Solve pre-assembled system
    131     return _la_solve(*args, **kwargs)

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/solving.py:180, in _solve_varproblem(*args, **kwargs)
    173 # Create solver and call solve
    174 solver = vs.NonlinearVariationalSolver(problem, solver_parameters=solver_parameters,
    175                                        nullspace=nullspace,
    176                                        transpose_nullspace=nullspace_T,
    177                                        near_nullspace=near_nullspace,
    178                                        options_prefix=options_prefix,
    179                                        appctx=appctx)
--> 180 solver.solve()

File petsc4py/PETSc/Log.pyx:115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/adjoint_utils/variational_solver.py:89, in NonlinearVariationalSolverMixin._ad_annotate_solve.<locals>.wrapper(self, **kwargs)
     86     tape.add_block(block)
     88 with stop_annotating():
---> 89     out = solve(self, **kwargs)
     91 if annotate:
     92     block.add_output(self._ad_problem._ad_u.create_block_variable())

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/variational_solver.py:282, in NonlinearVariationalSolver.solve(self, bounds)
    280     work.copy(u)
    281 self._setup = True
--> 282 solving_utils.check_snes_convergence(self.snes)
    284 # Grab the comm associated with the `_problem` and call PETSc's garbage cleanup routine
    285 comm = self._problem.u.function_space().mesh()._comm

File ~/Code/firedrake/firedrake/src/firedrake/firedrake/solving_utils.py:138, in check_snes_convergence(snes)
    136         else:
    137             msg = reason
--> 138         raise ConvergenceError(r"""Nonlinear solve failed to converge after %d nonlinear iterations.
    139 Reason:
    140    %s""" % (snes.getIterationNumber(), msg))

ConvergenceError: Nonlinear solve failed to converge after 32 nonlinear iterations.
Reason:
   DIVERGED_DTOL

Environment:

OS: MacOS M1 pro
Python version: 3.11.5
Output of firedrake-status:

Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: False
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    tinyasm: False
    torch: False
    petsc_int_type: int32
    cache_dir: /Users/tlynn/Code/firedrake/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: /opt/homebrew/opt/openblas
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |47f6c37   |False     |
|PyOP2               |master                        |8e4ad930  |False     |
|fiat                |master                        |8c66270   |False     |
|firedrake           |master                        |35f2b3cce |False     |
|h5py                |firedrake                     |6cc4c912  |False     |
|libspatialindex     |master                        |4768bf3   |True      |
|libsupermesh        |master                        |b145b65   |False     |
|loopy               |main                          |8158afdb  |False     |
|petsc               |firedrake                     |9364cb008b|False     |
|pyadjoint           |master                        |42959ef   |False     |
|pytest-mpi          |main                          |a478bc8   |False     |
|tsfc                |master                        |6f72c9c   |False     |
|ufl                 |master                        |55c281d7  |False     |
---------------------------------------------------------------------------

Any relevant environment variables or modifications: N/A

Additional Info N/A

Handling boundary conditions in a lid-driven cavity (or I guess a shear-driven cavity like you have here) is actually quite delicate. As specified, the boundary data are discontinuous at the corners.

Let's just consider the top-left corner:

x--x--x--
|    
x  
|
x
|

If you provide a list of boundary conditions, Firedrake applies them in turn. If the "top" (zero) bc is last, then the top-left corner dof will be (0, 0). On the other hand, if the "left" (non-zero) bc is last, then the top-left corner dof will be (0, 1).

This introduces a singularity in the numerical scheme that isn't there in the physical system. To handle this, one needs to treat the singularity somehow (Section 3.3 of this paper has a brief overview of techniques). Probably the simplest numerical approach is to pick a regularisation scheme that replaces the discontinuous data with $C^\infty$ data that is $1$ over most of the left edge and then falls to zero quickly at the corners. For a reasonable choice here (and how to change it as $\text{Re}$ is varied) see eq (6) of Lopez et al. (2017) (their problem is in 3D and drives the $z = 1$ lid, but it should be straightforward to modify for 2D).

Placing the zero Dirichlet boundary conditions last has, from the point of view of implementation, the same effect as this kind of mathematical regularisation. The reason being if we have a function that is one everywhere on the left boundary and decays to zero quickly near the corners, then at the dof that is on the midpoint of the edge near the top-left corner, the function will still be one. If we pick regularised boundary data, interpolating it to the finite element space will give use dof values:

0--x---x
|
1
|
1

Which is exactly the effect obtained by applying the top (zero) data last.

firedrakeproject / firedrake

BUG: Changing boundary condition order results in divergence #3093