def METHOD_NAME(backend=PSYDAC_BACKEND_PYTHON):
    # ... abstract model
    V = ScalarFunctionSpace('V', domain)
    x,y = domain.coordinates
    F = element_of(V,'F')
    v = element_of(V, 'v')
    u = element_of(V, 'u')
    expr = dot(grad(v), grad(u))
    a = BilinearForm((v,u), expr)
    expr = 2*pi**2*sin(pi*x)*sin(pi*y)*v
    l = LinearForm(v, expr)
    error = F -sin(pi*x)*sin(pi*y)
    l2norm = Norm(error, domain, kind='l2', name='u')
    h1norm = Norm(error, domain, kind='h1', name='u')
    # ...
    domain_h = discretize(domain, ncells=(2**8,2**8))
    # ...
    # ... discrete spaces
    Vh = discretize(V, domain_h, degree=(3,3))
    # ...
    Timing = namedtuple('Timing', ['kind', 'value'])
    # ...
    ah = discretize(a, domain_h, [Vh, Vh], backend=backend)
    tb = time.time()
    Ah = ah.assemble()
    te = time.time()
    matrix_timing = Timing('matrix', te-tb)
    # ...
    # ...
    lh = discretize(l, domain_h, Vh, backend=backend)
    tb = time.time()
    Lh = lh.assemble()
    te = time.time()
    rhs_timing = Timing('rhs', te-tb)
    # ...
    # ... coeff of phi are 0
    phi = FemField( Vh, 'phi' )
    # ...
    # ...
    l2norm_h = discretize(l2norm, domain_h, Vh, backend=backend)
    tb = time.time()
    error = l2norm_h.assemble(F=phi)
    te = time.time()
    l2norm_timing = Timing('l2norm', te-tb)
    # ...
    # ...
    if mpi_rank == 0:
        print_timing([matrix_timing, rhs_timing, l2norm_timing], backend)
    # ...