Construct differentiation matrix without sympy

[pbrubeck]

On simplex elements, if the degree is high enough, the main bottleneck is the construction of the differentiation matrix. We were doing this with sympy and that is extremely slow.

I replaced the tabulation of first and second derivatives with the derivative of the recurrence relation from (Kirby, 2010), and I am evaluating dmats lazily only if we ever need to tabulate derivatives.

This PR depends on #46

Here are the runtimes for a single call to Lagrange(cell, degree) before and after this PR:

Before:

Lagrange(triangle,  1) setup time 0.074 seconds
Lagrange(triangle,  2) setup time 0.059 seconds
Lagrange(triangle,  3) setup time 0.207 seconds
Lagrange(triangle,  4) setup time 0.473 seconds
Lagrange(triangle,  5) setup time 1.023 seconds
Lagrange(triangle,  6) setup time 1.935 seconds
Lagrange(triangle,  7) setup time 3.615 seconds
Lagrange(triangle,  8) setup time 6.519 seconds
Lagrange(triangle,  9) setup time 12.053 seconds

Lagrange(tetrahedron,  1) setup time 0.024 seconds
Lagrange(tetrahedron,  2) setup time 0.142 seconds
Lagrange(tetrahedron,  3) setup time 0.688 seconds
Lagrange(tetrahedron,  4) setup time 1.892 seconds
Lagrange(tetrahedron,  5) setup time 4.614 seconds
Lagrange(tetrahedron,  6) setup time 10.868 seconds

After:

Lagrange(triangle,  1): setup time 0.001 seconds
Lagrange(triangle,  2): setup time 0.001 seconds
Lagrange(triangle,  3): setup time 0.001 seconds
Lagrange(triangle,  4): setup time 0.001 seconds
Lagrange(triangle,  5): setup time 0.001 seconds
Lagrange(triangle,  6): setup time 0.001 seconds
Lagrange(triangle,  7): setup time 0.002 seconds
Lagrange(triangle,  8): setup time 0.002 seconds
Lagrange(triangle,  9): setup time 0.003 seconds
Lagrange(triangle, 10): setup time 0.004 seconds
Lagrange(triangle, 20): setup time 0.035 seconds
Lagrange(triangle, 30): setup time 0.150 seconds
Lagrange(triangle, 40): setup time 0.472 seconds
Lagrange(triangle, 50): setup time 1.163 seconds
Lagrange(triangle, 60): setup time 2.574 seconds
Lagrange(triangle, 70): setup time 4.916 seconds
Lagrange(triangle, 80): setup time 9.069 seconds
Lagrange(triangle, 90): setup time 16.213 seconds

Lagrange(tetrahedron,  1): setup time 0.002 seconds
Lagrange(tetrahedron,  2): setup time 0.001 seconds
Lagrange(tetrahedron,  3): setup time 0.001 seconds
Lagrange(tetrahedron,  4): setup time 0.002 seconds
Lagrange(tetrahedron,  5): setup time 0.004 seconds
Lagrange(tetrahedron,  6): setup time 0.007 seconds
Lagrange(tetrahedron,  7): setup time 0.013 seconds
Lagrange(tetrahedron,  8): setup time 0.023 seconds
Lagrange(tetrahedron,  9): setup time 0.037 seconds
Lagrange(tetrahedron, 10): setup time 0.060 seconds
Lagrange(tetrahedron, 15): setup time 0.456 seconds
Lagrange(tetrahedron, 20): setup time 2.262 seconds
Lagrange(tetrahedron, 25): setup time 8.945 seconds
Lagrange(tetrahedron, 30): setup time 30.647 seconds

[rckirby]

I support this endeavor and will review the code more carefully. A bit of FIAT history that is relevant to the discussion

• Originally, I had all the classic Duffy business and numerical dmats in FIAT, much as you are proposing. After writing a paper reworking the recurrence relations to avoid the corner singularity (https://dl.acm.org/doi/pdf/10.1145/1644001.1644006), I redid a lot of it, not yet symbolically.
• In a fit of laziness, I then used AD (but picked the Wrong Tool and created an awkward dependency) to get derivatives instead of using the derivative rules I worked out.
• To get around that (and for other reasons that @dham might recall), Somebody went through and stripped out the numerical dmat evaluation that was there (and you've reconstructed) in favor of sympy.
• I support going back to a numerical approach for speed reasons, but being able to get the basis functions as polynomials might be helpful, so making sure that the reworked evaluation rules let you put in x, y from sympy would be helpful.
  -- Note: if you avoid coordinate singularities, the recurrence relations are obviously polynomials and don't rely on any cancellation. It's trivial for sympy to DTRT here, but it's not clear whether it would have problems cancelling Jacobians symbolically if you go through Duffy.
• The code still seems to tabulate derivatives via _tabulate_dpts in sympy, but does the dmats numerically? I'm really afraid that we're going to make the whole thing very confusing. We need to decide what it should do and then make it do that.