Stalling least squares functional convergence

Question

I am seeing stalling convergence of the least squares functional when solving a least squares discretization of Poisson. This is part of a larger algorithm in which I solve the problem on a coarse mesh, then solve a residual equation (where the residual is calculated using the coarse solution) on a refined mesh. Furthermore, the right hand side for the residual problem is "chopped" such that it is non-zero only over a subset of the domain and set to zero elsewhere. As I refine the mesh for the residual equation and measure the least squares functional, I see convergence stall. Theory says that this should not happen. Have I implemented something incorrectly in my calculation of the least squares functional, or is Fenics doing some incorrect calculation somewhere here? A minimal example code which displays this stalling convergence behavior is below.

from dolfin import *

# Global options and files
#set_log_level(PROGRESS)
outfile = open('lsf.txt', 'w')
n_coarse = 2
n_refine = 7
order = 1

# Define expressions and subdomains used later
##############################################
##############################################

class HomeTurf(SubDomain):
    def inside(self, x, on_boundary):
        return True if (x[0] <= 0.5 and x[1] <= 0.5) else False

class Chopper(Expression):  
    def __init__(self, homeTurf):
        self.homeTurf = homeTurf
    def eval(self, values, x):
        if self.homeTurf.inside(x, False):
            values[0] = 1.0
        else:
            values[0] = 0.0

class allBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

class westeastBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

class southnorthBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Create boundary objects
boundarySN  = southnorthBoundary()
boundaryWE  = westeastBoundary()
boundaryALL = allBoundary()
home_turf   = HomeTurf()


# Define the PDE
##############################################
##############################################

def L(U):
    # definitions
    (ux, uy, p) = split(U)
    u = as_vector((ux,uy))

    # equations
    eq1 = u - grad(p)
    eq2 = div(u)
    eq3 = curl(u)
    return (eq1,eq2,eq3)


def f():
    # right hand side equations
    rhs1 = Expression(('0.0','0.0'))
    rhs2 = Expression('-2.0*x[0]*(x[0]-1.0) - 2.0*x[1]*(x[1]-1.0)')
    rhs3 = Expression('0.0')
    return (rhs1, rhs2, rhs3)

# Compute solution on coarse mesh
##############################################
##############################################

# Create initial mesh on unit square
mesh0 = UnitSquareMesh(n_coarse, n_coarse,"crossed")

# Create Vector Space  [ux, uy, p]
V0 = VectorFunctionSpace(mesh0, "Lagrange", order, 3)

u0 = TrialFunction(V0)
v0 = TestFunction(V0)

bilinear_form = 0
for Lu,Lv in zip(L(u0),L(v0)):
    bilinear_form += inner(Lu,Lv)*dx

linear_form = 0
for rhs,Lv in  zip(f(),L(v0)):
    linear_form += inner(rhs,Lv)*dx

# Setup boundary conditions and solve
bc_ux = DirichletBC(V0.sub(0), Expression('0.0'), boundarySN)
bc_uy = DirichletBC(V0.sub(1), Expression('0.0'), boundaryWE)
bc_p = DirichletBC(V0.sub(2), Expression('0.0'), boundaryALL)
BCs = [bc_ux, bc_uy, bc_p];

U0 = Function(V0)

solve(bilinear_form == linear_form, U0, bcs=BCs, solver_parameters={"linear_solver": "cg"})#,"preconditioner":"hypre_amg"})

# Compute solution to chopped residual equation on refined meshes
##############################################
##############################################

# Initialize new mesh
mesh = UnitSquareMesh(n_coarse, n_coarse,"crossed")

for iteration in range(n_refine):

    # Refine mesh
    mesh = refine(mesh)

    # Create Chi function (this chops the right hand side)
    DG_0 = FunctionSpace(mesh, "DG", 0)
    Chi = interpolate(Chopper(home_turf), DG_0);

    # Update vectorspace and trial/test function
    V = VectorFunctionSpace(mesh, "Lagrange", order, 3)
    u = TrialFunction(V)
    v = TestFunction(V)

    # Fine grid representation of coarse grid solution
    U0_curr = interpolate(U0, V)

    # Setup chopped right hand side
    ChiRes = []
    for Lu0,rhs in zip(L(U0_curr),f()):
            ChiRes.append(Chi*(rhs - Lu0))

    # Build bilinear/linear forms
    bilinear_form = 0
    for Lu,Lv in zip(L(u),L(v)):
            bilinear_form += inner(Lu,Lv)*dx

    linear_form = 0
    for rhs,Lv in  zip(ChiRes,L(v)):
            linear_form += inner(rhs,Lv)*dx

    # Construct BCs
    bc_ux = DirichletBC(V.sub(0), Expression('0.0'), boundarySN)
    bc_uy = DirichletBC(V.sub(1), Expression('0.0'), boundaryWE)
    bc_p = DirichletBC(V.sub(2), Expression('0.0'), boundaryALL)
    BCs = [bc_ux, bc_uy, bc_p];

    # Compute solution
    dU = Function(V)
    solve(bilinear_form == linear_form, dU, bcs=BCs, solver_parameters={"linear_solver": "cg"})#,"preconditioner":"hypre_amg"})

    # Compute LSF
    LSF = 0
    for Lu,rhs in zip(L(dU),ChiRes):
            LSF += inner(Lu-rhs,Lu-rhs)*dx

    lsf = sqrt(assemble(LSF))
    print "    dofs = ", mesh.num_cells()
    print "    LSF = ", lsf
    outfile.write(str(mesh.num_cells())+' '+str(lsf)+'\n')

# plot(dU[2])
# interactive()