Add replacement of quadrature tables with index lookup on dS integrals of symmetric quadrature rules

[jorgensd]

Avoids retabulations of the basis functions and reduced the size of the tables in the kernel.

For non-symmetric quadrature rules, this is not a valid operation, and this simplification is not applied (for instance collapsed gll on tetrahedral facets).
Inception of the idea was done in collaboration with @chrisrichardson.

The code has been generated by Claude (Sonnet 5) and has been reviewed, edited and tested by me.
The following summary is written by me, based on experiments created together with Claude.

Performance summary

Kernel runtime

As each kernel only accesses one row of the permutation table, the runtime is almost untouched by this change.
The plots below show the difference in runtime across 100 repeats.
Below we see the runtime comparison for a sweep of lagrange polynomial degrees for an order 4, 16 and 30 quadrature rule for the kernel:

 element = basix.ufl.element("Lagrange", args.cellname, args.elem_degree)
    domain = ufl.Mesh(basix.ufl.element("Lagrange", args.cellname, 1, shape=(3,)))
    space = ufl.FunctionSpace(domain, element)
    u, v = ufl.TrialFunction(space), ufl.TestFunction(space)
    dS = ufl.Measure(
        "dS",
        domain=domain,
        metadata={"quadrature_degree": args.qdegree, "quadrature_rule": args.rule},
    )
    form = ufl.inner(ufl.jump(ufl.grad(u)), ufl.jump(ufl.grad(v))) * dS

Compilation time and table sizes

As these tables now become single floating type tables + an integer table we same some space within each kernel.
However, the most noteable improvement is in the compilation time of the kernels. As we no longer have to do repeat tabulates per permutation, we same quite alot of time, as seen below.

Example code

Main branch

// Quadrature rules
static const double weights_a5e[7] = {0.1125, 0.06296959027241358, 0.06296959027241358, 0.06296959027241358, 0.06619707639425308, 0.06619707639425308, 0.06619707639425308};
// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE0_C0_D001_F_Qa5e[6][4][7][10] = {{{{....}}}};
static const double FE0_C0_D010_F_Qa5e[6][4][7][10] = {{{{....}}}};
static const double FE0_C0_D100_F_Qa5e[6][4][7][10] = {{{{}}}};
static const double FE2_C0_D100_F_Qa5e[1][1][1][4] = {{{{-1.0, 1.0, 0.0, 0.0}}}};
static const double FE2_C1_D010_F_Qa5e[1][1][1][4] = {{{{-1.0, 0.0, 1.0, 0.0}}}};
static const double FE2_C2_D001_F_Qa5e[1][1][1][4] = {{{{-1.0, 0.0, 0.0, 1.0}}}};
....
    for (int i = 0; i < 10; ++i)
    {
      temp_0[i] = fw0 * FE0_C0_D100_F_Qa5e[quadrature_permutation[0]][entity_local_index[0]][iq][i] + fw1 * FE0_C0_D010_F_Qa5e[quadrature_permutation[0]][entity_local_index[0]][iq][i] + fw2 * FE0_C0_D001_F_Qa5e[quadrature_permutation[0]][entity_local_index[0]][iq][i];
      ...
      temp_6[i] = fw6 * FE0_C0_D100_F_Qa5e[quadrature_permutation[1]][entity_local_index[1]][iq][i] + fw7 * FE0_C0_D010_F_Qa5e[quadrature_permutation[1]][entity_local_index[1]][iq][i] + fw8 * FE0_C0_D001_F_Qa5e[quadrature_permutation[1]][entity_local_index[1]][iq][i];
      ...
      }
    }
  }
  // ------------------------ 
}

}

this PR

// Precomputed values of basis functions and precomputations
// FE* dimensions: [permutation][entities][points][dofs]
static const double FE0_C0_D001_F_Qa5e[1][4][7][10] = {{{{...}}}};
static const double FE0_C0_D010_F_Qa5e[1][4][7][10] = {{{{...}}}};
static const double FE0_C0_D100_F_Qa5e[1][4][7][10] = {{{{...}}}};
static const int QPTa5e[6][7] = {{0, 1, 2, 3, 4, 5, 6},
  {0, 2, 1, 3, 5, 4, 6},
  {0, 3, 1, 2, 6, 4, 5},
  {0, 3, 2, 1, 6, 5, 4},
  {0, 2, 3, 1, 5, 6, 4},
  {0, 1, 3, 2, 4, 6, 5}};
   ....
    for (int i = 0; i < 10; ++i)
    {
      temp_0[i] = fw0 * FE0_C0_D100_F_Qa5e[0][entity_local_index[0]][QPTa5e[quadrature_permutation[0]][iq]][i] + fw1 * FE0_C0_D010_F_Qa5e[0][entity_local_index[0]][QPTa5e[quadrature_permutation[0]][iq]][i] + fw2 * FE0_C0_D001_F_Qa5e[0][entity_local_index[0]][QPTa5e[quadrature_permutation[0]][iq]][i];
     ...
     temp_6[i] = fw6 * FE0_C0_D100_F_Qa5e[0][entity_local_index[1]][QPTa5e[quadrature_permutation[1]][iq]][i] + fw7 * FE0_C0_D010_F_Qa5e[0][entity_local_index[1]][QPTa5e[quadrature_permutation[1]][iq]][i] + fw8 * FE0_C0_D001_F_Qa5e[0][entity_local_index[1]][QPTa5e[quadrature_permutation[1]][iq]][i];

[jorgensd]

Benchmark scripts can be found here:
benchmark_degree_sweep.py
benchmark_worker.py
plot_benchmark_results.py

[chrisrichardson]

It definitely seems nicer, but I guess it does add an extra layer of indirection.