def METHOD_NAME(solve_type):
    tape = Tape()
    set_working_tape(tape)
    # Defining the domain, 100 points from 0 to 1
    mesh = IntervalMesh(100, 0, 1)
    # Defining function space, test and trial functions
    V = FunctionSpace(mesh, "CG", 1)
    u = TrialFunction(V)
    u_ = Function(V)
    v = TestFunction(V)
    # Dirichlet boundary conditions
    bc_left = DirichletBC(V, 1, 1)
    bc_right = DirichletBC(V, 2, 2)
    bc = [bc_left, bc_right]
    # Some variables
    T = 0.5
    dt = 0.1
    f = Function(V)
    f.vector()[:] = 1
    u_1 = Function(V)
    u_1.vector()[:] = 1
    control = Control(u_1)
    a = u_1 * u * v * dx + dt * f * inner(grad(u), grad(v)) * dx
    L = u_1 * v * dx
    if solve_type == "LVS":
        problem = LinearVariationalProblem(a, L, u_, bcs=bc, constant_jacobian=False)
        solver = LinearVariationalSolver(problem)
    # Time loop
    t = dt
    while t <= T:
        if solve_type == "LVS":
            solver.solve()
        else:
            solve(a == L, u_, bc)
        u_1.assign(u_)
        t += dt
    J = assemble(u_1 ** 2 * dx)
    Jhat = ReducedFunctional(J, control)
    h = Function(V)
    h.vector()[:] = rand(h.dof_dset.size)
    u_1.tlm_value = h
    tape.evaluate_tlm()
    assert (taylor_test(Jhat, control.tape_value(), h, dJdm=J.block_variable.tlm_value) > 1.9)