Apply Sum Factorization Again

[jowezarek]

Issue

The sum factorization (SF) algorithm is currently implemented for 3D bilinear forms defined on the interior of a domain (when not using openmp).
Psydac decides whether to use the SF algorithm by checking the dimension of the domain, whether openmp is used to parallelize the code, and by checking whether kernel_expr is an instance of sympde.expr.evaluation.DomainExpression.
The latter does not work in practice, as apparently in most cases (always?) kernel_expr is a tuple.

Fix

I changed this code snippet in api/discretization.py from

    if isinstance(a, sym_BilinearForm):
        # Currently sum factorization can only be used for interior domains
        from sympde.expr.evaluation import DomainExpression
        is_interior_expr = isinstance(kernel_expr, DomainExpression) # <- This check must be changed

        if kwargs.pop('sum_factorization') and is_interior_expr:
            return DiscreteBilinearForm_SF(a, kernel_expr, *args, **kwargs)
        else:
            return DiscreteBilinearForm(a, kernel_expr, *args, **kwargs)

to

    if isinstance(a, sym_BilinearForm):
        # Currently sum factorization can only be used for interior domains
        from sympde.expr.evaluation import DomainExpression
        if isinstance(kernel_expr, tuple):
            is_interior_expr = isinstance(kernel_expr[0], DomainExpression)
        else:
            is_interior_expr = isinstance(kernel_expr, DomainExpression)

        if kwargs.pop('sum_factorization') and is_interior_expr:
            return DiscreteBilinearForm_SF(a, kernel_expr, *args, **kwargs)
        else:
            return DiscreteBilinearForm(a, kernel_expr, *args, **kwargs)

[yguclu]

Good catch for finding this bug!

I have a couple of comments. Also, please update the file CHANGELOG.md.

[yguclu]

I am not sure that we should report all of these numbers. In particular, the timings seem to oscillate quite a bit. Maybe one should run these benchmarks a few times and report the average value?

[yguclu]

In practice, is kernel_expr anything else than a tuple? Could it be a list?