Unable to successfully call PETSc function 'MatSetValuesLocal'.

Hi All,
Could you please help with this error?
The attached code is modified from this: http://www.sciencedirect.com/science/article/pii/S0098300416302746
Thanks a lot!

from dolfin import *
import numpy

def gravity(u):
    '''
    right hand side of Eq. 7
    '''
    Ra = 1.0
    val = as_vector([0.0, Ra*u])
    return val

def solver(nx, order, t, dt, T):
    '''
    nx: mesh density
    order: approximation order of pressure field
    t: simulation start time
    dt: time step
    T: stop time
    '''
    mesh =UnitSquareMesh(nx, nx)

    #W = MixedFunctionSpace([DGO, BDM, DG1])
    P1 = FiniteElement("DG", mesh.ufl_cell(), order)
    P2 = FiniteElement("BDM", mesh.ufl_cell(), order+1)
    P3 = FiniteElement("DG", mesh.ufl_cell(), order+1)

    DGO = FunctionSpace(mesh, P1)
    # function space for Velocity
    BDM = FunctionSpace(mesh, P2)
    # function space for Temperature
    DG1 = FunctionSpace(mesh, P3)

    WE = MixedElement([P1,P2,P3])
    W = FunctionSpace(mesh, WE)

    #source terms
    f_1 = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element = DGO.ufl_element())
    f_2 = Expression(('0.0','(-sin(x[0]) - sin(x[1]))*cos(t)'), t=t, element = BDM.ufl_element())
    f_3 = Expression('-(sin(x[0]) +sin(x[1]))*sin(t) + (-cos(x[0])*cos(x[0]) + sin(x[1])*cos(x[1]))*cos(t)*cos(t) + (sin(x[0])+sin(x[1]))*cos(t)', t=t, element = DG1.ufl_element())

    # exact solutions
    p_exact = Expression('(sin(x[0]) + cos(x[1]))*cos(t)', t=t, element = DGO.ufl_element()) # pressure
    uu_exact = Expression(('-cos(x[0])*cos(t)','sin(x[1])*cos(t)'), t=t, element = BDM.ufl_element()) # velocity
    u_exact = Expression('(sin(x[0]) + sin(x[1]))*cos(t)', t=t, element = DG1.ufl_element()) # tempreture

    #(p, uu, u) = TrialFunctions(W)  # store the solution for current step
    #(p0, uu0, u0) = TrialFunctions(W) # store the solution for previous step
    w = Function(W)
    w0 = Function(W)

    # ititialize solution according to exact solution
    assign(w.sub(0), interpolate(p_exact, DGO))
    assign(w.sub(1), interpolate(uu_exact, BDM))
    assign(w.sub(2), interpolate(u_exact, DG1))
    assign(w0.sub(0), interpolate(p_exact, DGO))
    assign(w0.sub(1), interpolate(uu_exact, BDM))
    assign(w0.sub(2), interpolate(u_exact, DG1))

    # split w into p (pressure), uu (velocity) and u (temperature)
    (p, uu, u) = w.split()
    (p0, uu0, u0) = w0.split()

    # define test functions
    (q, vv, v) = TestFunctions(W)

    n = FacetNormal(mesh)

    # penalty weights
    alpha = Constant(5000.0)
    gamma = Constant(10000.0)

    # cell size
    h = CellSize(mesh)

    # Eq. 20, flow
    F_flo = nabla_div(uu)*q*dx - f_1*q*dx
    #F_flo = inner(nabla_div(uu),q)*dx - inner(f_1,q)*dx

    # Eq. 21, velocity
    F_vel = (dot(uu,vv) - div(vv)*p - inner(gravity(u),vv) - inner(f_2,vv))*dx + dot(n,vv)*p_exact*ds

    # un = un on the outflow facet, otherwise 0
    un = (dot(uu,n) +abs(dot(uu,n))) / 2.0

    # Eq. 26 , uu is not divergence-free here
    a_a = dot(grad(v), -uu*u)*dx - u*v*nabla_div(uu)*dx - f_3*v*dx + dot(jump(v), un('+')*u('+') - un('-')*u('-'))*dS + dot(v,un*u)*ds

    # Eq. 28
    a_d = dot(grad(v),grad(u))*dx - dot(avg(grad(v)),jump(u,n))*dS - dot(jump(v,n), avg(grad(u)))*dS\
        - dot(grad(v),n*u)*ds - dot(v*n, grad(u))*ds\
        + (alpha('+')/h('+'))*dot(jump(v,n),jump(u,n))*dS + (gamma/h)*v*u*ds\
        + u_exact*dot(grad(v),n)*ds - (gamma/h)*u_exact*v*ds

    a_tim = (u-u0)/dt*v*dx

    #Eq. 29
    F_a_d = a_tim +a_a +a_d

    #Eq. 32
    F = F_flo + F_vel + F_a_d

    # Time-stepping loop
    while t<T:
        t += dt
        print "nx=%d order=%d t=%f" % (nx, order, t)
        # update time in time-dependent expressions
        f_1.t = t; f_2.t = t; f_3.t = t
        u_exact.t = t; p_exact.t = t; uu_exact.t = t
        # solve the system and store solution in w
        solve(F==0, w)
        # update w0 with w
        w0.vector()[:] = w.vector()
        # save solution in .xml file
        File('Solution.xml') << w


solver(10, 1, 0, 0.01, 0.1)

nx=10 order=1 t=0.010000
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.704e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
File "hydro_thermo.py", line 125, in
solver(10, 1, 0, 0.01, 0.1)
File "hydro_thermo.py", line 118, in solver
solve(F==0, w)
File "/Users/geomechanics/opt/fenics-dev/g2jtvjzvjksg/master/lib/python2.7/site-packages/dolfin/fem/solving.py", line 300, in solve
_solve_varproblem(*args, **kwargs)
File "/Users/geomechanics/opt/fenics-dev/g2jtvjzvjksg/master/lib/python2.7/site-packages/dolfin/fem/solving.py", line 349, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /Users/geomechanics/dev/fenics-dev/master/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0

*** DOLFIN version: 2017.1.0.dev0
*** Git changeset: ec44213345f8397f914ec6f4b124c9d5998116f3
*** -------------------------------------------------------------------------

Thanks!