Applying DirichletBC on component of a VectorFunction that is in a mixed function space does not work

[florianwechsung]

I tried to adapt this mail in the mailinglist to the Stokes equations. I would like to apply Dirichlet BCs on only one component of the velocity field. However, in that case the solver fails (failed to converge after 10000 iterations).

I have slightly adapted my use case to illustrate the problem.

from firedrake import *
# Load mesh
mesh = UnitSquareMesh(10, 10)

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

applyBcsComponentWise = True
if applyBcsComponentWise:
    bc0 = DirichletBC(W.sub(0).sub(0), 0, [3, 4])
    bc1 = DirichletBC(W.sub(0).sub(1), 0, [3, 4])
    bc2 = DirichletBC(W.sub(0).sub(0), 1, 1)
    bc3 = DirichletBC(W.sub(0).sub(1), 0, 1)
    bcs = [bc0, bc1, bc2, bc3]
else:
    bc0 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), [3, 4])
    bc1 = DirichletBC(W.sub(0), Constant((1.0, 0.0)), 1)
    bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

U = Function(W)
solve(a == L, U, bcs=bcs)

for both values of the bool applyBcsComponentWise this snippet should produce the same result. However, for applyBcsComponentWise = True, the solver runs out of iterations. According to @wence- , to many entries in the matrix are zero'd.

[wence-]

To elaborate, for these VectorFunctionSpace component BCs, we encode the component in the high bits of the map value. This is fragile. In particular we pass information about boundary conditions in the matrix through as follows:

• For normal (apply to all dofs on a node) bcs we encode this information as -(map_value + 1).
• For componentwise bcs we encode as -(map_value + 1 + sum(2**(30 - cmpt), components))

In pyop2 we can then undo this and pull out the relevant components. For nodes that are on the corner of a domain (i.e. more then one subdomain applies) this can result in us producing "negative" values that overflow (e.g. 2**30 + 2**30 in 32bit signed integer arithmetic is -2**31, not 2**31).

A proper fix would be to instead decorate the map with the components that the bcs should apply to, and we can then pass more information into pyop2 and generate the correct code directly.

[wence-]

We think this is fixed by #951.

[wence-]

And merged.

[florianwechsung]

I'm afraid there is still an issue, this test fails:

from firedrake import *
import numpy as np
n = 2
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

bcs_1=[]
bcs_1.append(DirichletBC(V.sub(1), 0, 3))
bcs_1.append(DirichletBC(V, Constant((0, 0)), 4))
bcs_1.append(DirichletBC(V.sub(0), 1, 1))
bcs_1.append(DirichletBC(V.sub(1), 0, 1))

bcs_2=[]
bcs_2.append(DirichletBC(V.sub(1), 0, 3))
bcs_2.append(DirichletBC(V, Constant((0, 0)), 4))
bcs_2.append(DirichletBC(V, Constant((1, 0)), 1))

a = inner(grad(u), grad(v)) * dx

def asarray(A):
    return A.M.handle[:, :]
A_1 = asarray(assemble(a, bcs=bcs_1, mat_type="aij"))
A_2 = asarray(assemble(a, bcs=bcs_2, mat_type="aij"))
assert(np.allclose(A_1, A_2))

[wence-]

I think I got it now.

[florianwechsung]

Perfect, thank you!