speed up periodic boundary coupling

[pjaap]

By precomputing suitable search areas on the target grid boundary region.

Compute for each face the bounding box, transfer it to the opposite side and check intersections with the bounding boxes of the faces there.

A not so trivial benchmark:

xgrid = simplexgrid(0:0.05:1.0, 0:0.05:1.0, 0:0.05:1.0)
FES = FESpace{H1P2{2, 3}}(xgrid)
function give_opposite!(y, x)
    y .= x
    y[1] = 1 - x[1]
    return nothing
end

Then,
@time A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8, heuristic_search = true)
takes about 31 seconds,

while the current method
@time A = get_periodic_coupling_matrix(FES, xgrid, 4, 2, give_opposite!, sparsity_tol = 1.0e-8, heuristic_search = false)
takes about 5 minutes!

The switch heuristic_search should probably be removed when this improvement is accepted.

Remark: the allocations are crazy (independent of the switch): 31.423613 seconds (36.58 M allocations: 31.542 GiB, 5.39% gc time) ... we should look into this.

[pjaap]

@chmerdon what do you think? Seems to cause no problems so far, so I would remove the heuristic_search entry and merge.

[chmerdon]

I put the wrap of the point evaluator outside the loop which improved the situation further:
14.768451 seconds (32.95 M allocations: 21.125 GiB, 6.16% gc time, 5.17% compilation time)
There are still a lot of allocations in that main loop that we should try to bring down...

[pjaap]

Nice. We should still tackle the allocation problem.
I did some deep-dive into this an most allocation happen in
edgemoments::Array{T, 2} = zeros(T, nmoments, nitems)
(https://github.com/WIAS-PDELib/ExtendableFEMBase.jl/blob/master/src/interpolations.jl#L249)

[chmerdon]

Yes, so far the interpolate was expected to be used for all elements at once, here it is used in a loop to interpolate on each single boundary face. Maybe we can reduce the number of interpolate calls somehow by treating several boundary faces at once, that have non overlapping supports on the other side? Or we need to find a way to rewrite the interpolate code, such that variables like edgemoments are reused...

[chmerdon]

Got rid of some more minor allocations that brought also a little speedup... I am now at
11.737367 seconds (16.93 M allocations: 20.625 GiB, 7.81% gc time)

[pjaap]

Great. I removed the heuristic_search switch, since we agree that anything else than true is not suitable.

[chmerdon]

I accidently ran another Formatter, before I remember to run the pre-commit, maybe that caused it. At least I don't see a need for these commas and white spaces myself ^^

[pjaap]

There shall not be another formatter than Runic 😆