def METHOD_NAME(mesh, iterations=100):
    # test function space
    v = FunctionSpace(mesh, "DG", 1)
    m = VectorFunctionSpace(mesh, "CG", 1)
    # advecting velocity
    if m.shape == (1, ):
        u0 = as_vector([1])
    else:
        u0 = as_vector([1, 0])
    u = Function(m).interpolate(u0)
    # advection problem
    dt = 1. / iterations
    phi = TestFunction(v)
    D = TrialFunction(v)
    n = FacetNormal(mesh)
    un = 0.5 * (dot(u, n) + abs(dot(u, n)))  # upwind value
    a_mass = inner(D, phi) * dx
    a_int = inner(-u * D, grad(phi)) * dx
    a_flux = inner(un('+') * D('+') - un('-') * D('-'), jump(phi)) * dS
    arhs = a_mass - dt * (a_int + a_flux)
    dD1 = Function(v)
    D1 = Function(v)
    x = SpatialCoordinate(mesh)
    # Initial Conditions
    D0 = conditional(real(x[0]) < 0.5, 1., 0.)
    D = Function(v).interpolate(D0)
    D1.assign(D)
    t = 0.0
    T = iterations * dt
    problem = LinearVariationalProblem(a_mass, action(arhs, D1), dD1)
    solver = LinearVariationalSolver(problem, solver_parameters={'ksp_type': 'cg'})
    # Make slope limiter
    limiter = VertexBasedLimiter(v)
    limiter.apply(D)
    while t < (T - dt / 2):
        D1.assign(D)
        limiter.apply(D1)
        solver.solve()
        D1.assign(dD1)
        limiter.apply(D1)
        solver.solve()
        D1.assign(0.75 * D + 0.25 * dD1)
        limiter.apply(D1)
        solver.solve()
        D.assign((1.0 / 3.0) * D + (2.0 / 3.0) * dD1)
        limiter.apply(D1)
        t += dt
    assert np.max(u.dat.data_ro) <= 1.0, "Failed by exceeding max values"
    assert np.min(u.dat.data_ro) >= 0.0, "Failed by exceeding min values"