Solving PDE inside a subdomain defining new measures

Hi!

I am trying to simulate the incompressible Navier-Stokes equation throughout a pipe using Chorin's projection method as the solver, so my mesh (generated with Gmsh) consists of two cylinders, one inside the other one, and with a longitude of 1 in the z-axis. I just want to study the velocity and the pressure in the inner cylinder (with radius r2), so I created a subdomain and then defined new measures so later in the solver I could integrate just in that inner cylinder:

r1=1.0 #radius of big cylinder
r2 = 0.5715 #radius of small cylinder

# Set parameter values
dt = 0.001 #time-step
T =  0.1 #length of time interval
nu = 6.602E-7 #kinematic viscosity

# Define pressure boundary conditions
p_in = Constant(10.0) 

# Define velocity boundary conditions
noslip = Constant((0.0, 0.0, 0.0))  #With 0 < z < 1

# bottom z=0
bottom = AutoSubDomain(lambda x, on_boundary: on_boundary and near(x[2], 0.) and (x[0]**2 + x[1]**2 < r2**2))

# boundary of inner cylinder
insidewall = AutoSubDomain(lambda x, on_boundary: near(x[0]**2 + x[1]**2, r2**2))

# Define boundary conditions
bc3 = DirichletBC(Q, p_in, bottom)
bc1 = DirichletBC(V, noslip, insidewall, "pointwise")

bcp = [bc3]
bcu = [bc1]

# inner pipe
class Pipe(SubDomain):
    def inside(self, x, on_boundary):
        return ((x[0]**2 + x[1]**2) < r2**2)
pipe = Pipe()

# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
pipe.mark(domains, 1)

# Define new measures associated with the interior domains
dx = Measure("dx")[domains]

# Create functions
u0 = Function(V)
u0 = interpolate(Constant((0., 0., 0.)), V)
u1 = Function(V)
p1 = Function(Q)

# Define coefficients
k = Constant(dt)
f = Constant((0., 0., 1.))

# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx(1) + inner(grad(u0)*u0, v)*dx(1) 
+ nu*inner(grad(u), grad(v))*dx(1) - inner(f, v)*dx(1)
a1 = lhs(F1)
L1 = rhs(F1)

# Pressure update
a2 = inner(grad(p), grad(q))*dx(1)
L2 = -(1/k)*div(u1)*q*dx(1)

# Velocity update
a3 =  inner(u, v)*dx(1)
L3 =  inner(u1, v)*dx(1) - k*inner(grad(p1), v)*dx(1)

but I am getting the following error

*** Error:   Unable to set given rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  unknown

What am I doing wrong?

Thank you in advance!

I have changed my code adding keep_diagonal=True when assembling matrices:

# Assemble matrices keeping diagonal entries in sparsity pattern
A1 = assemble(a1, keep_diagonal=True)
A2 = assemble(a2, keep_diagonal=True)
A3 = assemble(a3, keep_diagonal=True)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

t = dt
while t < T + DOLFIN_EPS:
    # Compute tentative velocity step
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "gmres", prec)
    end()

    # Pressure correction
    begin("Computing pressure correction")
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    solve(A2, p1.vector(), b2, "gmres", prec)
    end()

    # Velocity correction
    begin("Computing velocity correction")
    b3 = assemble(L3)
    [bc.apply(A3, b3) for bc in bcu]
    solve(A3, u1.vector(), b3, "gmres", prec)
    end()

and now it allows me to integrate just with the index 1 so the simulation will take place in the inner cylinder , but now I am getting:

*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76.
*** Where:   This error was encountered inside /build/buildd/dolfin-1.4.0+dfsg/dolfin /la/PETScKrylovSolver.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  unknown

Did I do a good use of keep_diagonal?

Thanks!