Unable to access vector of degrees of freedom fro function.

I am iterating a non-linear problem using simple 'Picard' iteration (as per the example
at http://fenicsproject.org/documentation/tutorial/nonlinear.html)
Incidentally, this example does not work with the current version of FEniCS - and needs updating to reflect the
change to the VariationalProblem() interface.

Anyway, my problem is as follows. I have a MixedFunctionSpace, and I wish to access one of the vectors, to calculate
the norm.

#!/usr/bin/python

from dolfin import *
import numpy

# Parameters
Ra=100.0
XX = 2.0

# uniform mesh
mesh=Rectangle(0,0,XX,1.0,int(20*XX),20)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
X = MixedFunctionSpace([V, Q, Q])

u,p,T = TrialFunctions(X)
w,q,S = TestFunctions(X)

u0=Function(V)

# bilinear form - try to put trial function first...
a = ( inner(grad(w),sym(grad(u))) \
    - div(w)*p \
    - w[1]*T*Ra \
    - q*div(u) \
    + S*dot(u0,grad(T)) \
    - dot(grad(S),grad(T)) \
    )*dx

f = Constant(0.0)
L = q*f*dx

# Define boundaries
def bound0(x,on_boundary): # x=0
    return x[0] < DOLFIN_EPS and on_boundary
def bound1(x,on_boundary): # x=XX
    return x[0] > XX*(1.0-DOLFIN_EPS) and on_boundary
def bound2(x,on_boundary): # z=0
    return x[1] < DOLFIN_EPS and on_boundary
def bound3(x,on_boundary): # z=1
    return x[1] > (1.0-DOLFIN_EPS) and on_boundary

## Define BCs for velocity
vzero = Constant((0.0,0.0))
zero = Constant(0.0)
Hsource=Expression("exp(-pow((x[0]-1.0)/a,2))",a=1.0)
Vbc0 = DirichletBC(X.sub(0),vzero,bound0) # v=0 on x=0
Vbc1 = DirichletBC(X.sub(0),vzero,bound1) # v=0 on x=XX
Vbc2 = DirichletBC(X.sub(0),vzero,bound2) # v=0 on z=0
Vbc3 = DirichletBC(X.sub(0).sub(1),zero,bound3) # vz=0 on z=1
Vbc4 = DirichletBC(X.sub(1),zero,bound3) # p=0 on z=1
Vbc5 = DirichletBC(X.sub(2),zero,bound3) # T=0 on z=1
Vbc6 = DirichletBC(X.sub(2),Hsource,bound2) # T=e^-(r/a)^2 on z=0
Vbcs = [ Vbc0, Vbc1, Vbc2, Vbc3, Vbc4, Vbc5, Vbc6 ]

U=Function(X)
u,p,T=U.split()

PP = LinearVariationalProblem(a,L,U,bcs=Vbcs)
solver= LinearVariationalSolver(PP)

#iterate, changing the velocity field each step
for i in range(10):
    solver.solve()
    diff = u.vector().array() - u0.vector().array()
    eps = numpy.linalg.norm(diff, ord=numpy.Inf)
    print eps
    u0.assign(u)

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://answers.launchpad.net/dolfin
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom fro function.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** -------------------------------------------------------------------------