Segmentation fault, but debug step by step error gone!

Question

Thanks in advance!
The question is simple, please help me.
I think I find a bug if I am right.

When I compute errornorm inside of my_adaptive_solver(), it's OK!
But I get segmentation fault when computing errornorm out of my_adaptive_solver().

The source code is below:

from dolfin import *

def my_adaptive_solver(u_e, f, mesh):
    V = FunctionSpace(mesh, "Lagrange", 1)

    bc = DirichletBC(V, u_e, DomainBoundary())

    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v))*dx()
    L = f*v*dx

    # Define function for the solution
    u = Function(V)
    # Define goal functional (quantity of interest)
    M = u*dx()
    tol = 1.e-3

    max_iter = 10
    problem = LinearVariationalProblem(a, L, u, bc)
    solver = AdaptiveLinearVariationalSolver(problem, M)
    solver.parameters['max_iterations'] = max_iter
    solver.solve(tol)

    u_h = []; u_h.append(u.root_node())
    num_dofs = []; num_dofs.append(u.root_node().function_space().dim())
    temp = u.root_node()
    while temp.has_child():
        u_h.append(temp.child())
        num_dofs.append(temp.child().function_space().dim())
        temp = temp.child()

    return u_h, num_dofs


def compute_convergence_rate(u_e, u_h_list, num_dofs_list):
    '''
    u_e is the exact solution
    u_h is a list of finite element solution.
    num_dofs is a list of number of dofs, using it to get mesh size.
    '''
    from math import log as ln
    rate = []
    geo_dim = 2
    norm_type = 'H1'

    h = [num_dofs_list[i]**(-1.0/geo_dim) for i in range(len(u_h_list))]
    error = [errornorm(u_e, u_h_list[i], norm_type=norm_type) for i in range(len(u_h_list))]
    h.sort(reverse=True); error.sort(reverse=True)

    for i in range(len(u_h_list) - 1):
        rate.append(ln(error[i]/error[i+1]) / ln(h[i]/h[i+1]))

    return h, rate


if __name__ == '__main__':
    import pdb
    pdb.set_trace()

    mesh = UnitSquareMesh(4, 4)
    u_e = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
    f = Constant(-6.0)

    u_h_list, num_dofs_list = my_adaptive_solver(u_e, f, mesh)

    h, rate = compute_convergence_rate(u_e, u_h_list, num_dofs_list)

    print h
    print rate

    interactive()