Reduce Sauter-Schwab assembly allocations

[djukic14]

This PR reduces allocations in BEAST’s Sauter-Schwab assembly code.

The main idea is to reuse small temporary arrays instead of allocating them again and again during assembly. It also changes the quadrature dispatch so that quadrule can directly call momintegrals!, instead of first returning a quadrature rule that is passed on later.

This PR depends on the accompanying SauterSchwabQuadrature.jl PR krcools/SauterSchwabQuadrature.jl#10, which adds the in-place reorder! function.

Changes

• Add quadraturebuffer(quadstrat) for reusable temporary arrays.
• Use these buffers for the I, J, K, and L arrays needed by reorder!.
• Add callback-based quadrature dispatch with QuadruleCallback, ApplyMomintegrals, and ApplyLocalMomintegrals.
• Let quadrule directly apply momintegrals! in the assembly hot path.
• Pass qbuffer explicitly through momintegrals!.
• Use SauterSchwabQuadrature.reorder! instead of allocating in reorder.
• Update tests for the explicit qbuffer argument.
• Bump BEAST to 2.11.0.
• Require SauterSchwabQuadrature = "2.5.0".

Motivation

Profiling showed many small allocations during singular and near-singular assembly. A large part of these came from Sauter-Schwab reordering and from passing many different quadrature rule types through the assembly code.

Reusing buffers and calling momintegrals! directly from quadrule reduces these allocations in the assembly.

[krcools]

The main thing I need convincing of is the apply -> quadrule idea. Why is having momintegrals inside of quadrule better? Surely the dynamic dispatch has to happen somewhere, right?

Dispatch now happens in quadrulecallbacks.jl:30, is that correct?

If we go this way, we need to come up with better names to make the code at least somewhat readable. And a page in the docs to explain the ideas here. Maybe:

apply -> integrands
quadrule(apply, op, ...., qdata, qstrat) -> integrate!(integrands, qstrat, qdata, qbuffer, zlocal)
(f::ApplyLocalMomintegrals)(qrule) -> integrate!(integrands, qrule, qbuffer, zlocal)

This is a bit longer, but:

• splits immutable and mutable arguments
• splits the conceptual (integrate integrands), and the implementation details (optimisation through buffer reuse)
• sleeping better knowing I understand the code.

I would also insist on removing the type ApplyMomintegrals from the assembly loop. Given the overlap in fields of ApplyMomintegrals and the arguments to quadrule, maybe it can be removed altogether. So maybe just:

integrate!(op,
    test_functions, test_cellptr, test_chart,
    trial_functions, trial_cellptr, trial_chart,
    qstrat, qdata, qbuffer, zlocal)

integrate!(op,
    test_functions, test_cellptr, test_chart,
    trial_functions, trial_cellptr, trial_chart,
    qrule, qbuffer, zlocal)

User want a new strategy: overload the first and a corresponding pair (quaddata, quadbuffer). Wants a new quadrature rule: overload the second.

ReturnQuadrule has quadrule behave as before: without introducing calls to the appropriate momintegrals!. Where is this needed? This adds to the mental load when traversing the code. I appreciate that quadrule selection and quadrule selection and integration share all code, but could this not be done in another way?

function quadrule_selection(integrand, qstrat, qbuff, qdata)
    if numhits == 3
        return SSCommonFace(....)
    ...
end

function integrate!(integrand, qstrat, qbuff, qdata, zlocal)
    qr = quadrule_selection(integrand, qstrat, qbuff, qdata)
    return integrate!(integrand, qr, qbuff, zlocal)
end

I suppose this would introduce the allocations you try to combat...