Adaptive nonlinear solver in 2d doesn't work if only components of field are specified with bcs

Hi,

I am trying to use the adaptive non linear solver in 2d to compute the flow field of two counter-rotating vortices inside an area and I have the following issue:

I want to use a 'free' or 'slip' or invicid boundary condition for the centerline of the chamber between the vortices since the flow field will be symmetric at that point.

I can do this in the normal nonlinear solver by initializing one of the components to 0. i.e:

bcslip = DirichletBC(W.sub(0).sub(0), Constant(0.0), Lid())

However when I try using this with the adaptive solver I get this:

Error: Unable to create Dirichlet boundary condition.
*** Reason: Expecting a vector-valued boundary value but given function is scalar

Another way of doing this is to enforce the boundary condition in the variational form with something like this (added last term):

n = FacetNormal(mesh)
F = (nuinner(grad(u), grad(v)) - div(v)p + qdiv(u) + inner(grad(u)u, v))dx() + nuinner(grad(u) * n, v)*ds(2)

Where ds(2) corresponds to the Lid() boundary.

This allows me to not have the component of W.sub(0) set to 0, which solves the adaptive solver error. However the solution has velocity vectors that are not tangent to the Lid() boundary. Is the last term correct?

Do you know if the error I got is due to a boundary condition the adaptive solver can't use?

What do you suggest I use for a inviscid boundary condition that will work with the adaptive nonlinear solver? I would prefer to have a term in F since it would allow me to have a free slip boundary in any geometry in the future

The code is here:

from dolfin import *

#Area measurements
rad = 0.01
vfr = 0.004

#Peak vortex speed
Upeak = 1.0

meshtol = 1e-5

class InR(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (sqrt(x[0]*x[0] + x[1]*x[1]) < (vfr + meshtol))

class PPoint(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < -rad + meshtol and x[1] < -rad + meshtol

class NoSlip(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (sqrt(x[0]*x[0] + x[1]*x[1]) > vfr) and x[0] > -rad

class Lid(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < -rad + meshtol


# Use compiler optimizations
parameters["form_compiler"]["cpp_optimize"] = True

# Allow approximating values for points that may be generated outside
# of domain (because of numerical inaccuracies)
parameters["allow_extrapolation"] = True

# Choose refinement algorithm (needed for solver to choose right refinement algorithm)
parameters["refinement_algorithm"] = "plaza_with_parent_facets"


# Material parameters
nu = Constant(0.0002)

# Mesh
mesh = Mesh('twovortexmesh.xml')

# Create boundary subdomains
bc_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bc_markers.set_all(0)
InR().mark(bc_markers, 1)
Lid().mark(bc_markers, 2)

# Define new measure with associated subdomains
ds = Measure("ds")[bc_markers]

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

# Define unknown and test function(s)
(v, q) = TestFunctions(W)
w = Function(W)
(u, p) = (as_vector((w[0], w[1])), w[2])

# Define variational forms for Navier-Stokes
n = FacetNormal(mesh)
F = (nu*inner(grad(u), grad(v)) - div(v)*p + q*div(u) + inner(grad(u)*u, v))*dx()# + nu*inner(grad(u) * n, v)*ds(2)

# Define boundary conditions
coreu  = Expression(("Upeak*x[1]/vfr","-Upeak*x[0]/vfr"), Upeak=Upeak, vfr=vfr)
bcu    = DirichletBC(W.sub(0), coreu, InR())
bcslip = DirichletBC(W.sub(0).sub(0), Constant(0.0), Lid())
bcp    = DirichletBC(W.sub(1), Constant(0.0), PPoint(), "pointwise")
noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), NoSlip())
#bc = [bcu, bcp, noslip]
bc = [bcu, bcp, bcslip, noslip]

# Define goal
M = p*ds(1)

# Define error tolerance (with respect to goal)
tol = 1e-07

# Compute Jacobian form
J = derivative(F, w)

# Define variational problem
pde = NonlinearVariationalProblem(F, w, bc, J)

# Define solver
#solver = AdaptiveNonlinearVariationalSolver(pde, M)
solver = NonlinearVariationalSolver(pde)

# Solve to given tolerance
#solver.solve(tol)
solver.solve()

# Show all timings
list_timings()

# Solutions on coarsest and finest mesh:
mesh0 = mesh.root_node()
mesh1 = mesh.leaf_node()
plot(mesh0, title="Initial mesh")
plot(mesh1, title="Final mesh")

(u0, p0) = w.root_node().split()
(u1, p1) = w.leaf_node().split()
plot(u0, title="Velocity on initial mesh")
plot(u1, title="Velocity on final mesh")
plot(p0, title="Pressure on initial mesh")
plot(p1, title="Pressure on final mesh")
#(u, p) = w.split()
#plot(u, title="Velocity field")
#plot(p, title="Pressure")
interactive()

Thank you

Alex