def METHOD_NAME():
    mesh = UnitSquareMesh(3, 3)
    V_h = VectorElement("CG", mesh.ufl_cell(), 2)
    Q_h = FiniteElement("CG", mesh.ufl_cell(), 1)
    W = FunctionSpace(mesh, V_h * Q_h)
    M = FunctionSpace(mesh, "CG", 2)
    m = Function(M)
    m.vector()[:] = 2.0
    expr = Expression(("m*m*m", "m*m"), m=m, degree=1)
    deriv = Expression(("3*m*m", "2*m"), m=m, tlm_input=1.0, degree=1)
    deriv.user_defined_derivatives = {m: Expression(("6*m", "2"), m=m, degree=1)}
    expr.user_defined_derivatives = {m: deriv}
    noslip = DirichletBC(W.sub(0), (0, 0),
                         "on_boundary && (x[1] >= 0.9 || x[1] < 0.1)")
    inflow = DirichletBC(W.sub(0), expr, "on_boundary && x[0] <= 0.1 && x[1] < 0.9 && x[1] > 0.1")
    bcs = [inflow, noslip]
    s = Function(W, name="State")
    v, q = TestFunctions(W)
    x = TrialFunction(W)
    u, p = split(x)
    nu = Constant(1)
    a = (nu * inner(grad(u), grad(v)) * dx
         - inner(p, div(v)) * dx
         - inner(q, div(u)) * dx
         )
    L = inner(Constant((0.0, 0.0)), v) * dx
    A, b = assemble_system(a, L, bcs)
    solve(A, s.vector(), b)
    u, p = split(s)
    J = assemble(u**2*dx)
    Jhat = ReducedFunctional(J, Control(m))
    h = Function(M)
    h.vector()[:] = rand(M.dim())
    dJdm = Jhat.derivative()._ad_dot(h)
    Hm = Jhat.hessian(h)._ad_dot(h)
    assert taylor_test(Jhat, m, h, dJdm=dJdm, Hm=Hm) > 2.9