Reduce Sauter-Schwab assembly allocations

Project.toml

@@ -1,6 +1,6 @@
 name = "BEAST"
 uuid = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
-version = "2.10.0"
+version = "2.11.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -67,7 +67,7 @@ ProgressMeter = "1.11.0"
 PlotlyBase = "0.8.21"
 Requires = "1"
 SauterSchwab3D = "0.2"
-SauterSchwabQuadrature = "2.4.0"
+SauterSchwabQuadrature = "2.5.0"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1, 2"
 StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1"
 Suppressor = "0.2.8"

src/BEAST.jl

@@ -230,6 +230,8 @@ include("quadrature/strategies/testrefinestrialqstrat.jl")
 include("quadrature/strategies/trialrefinestestqstrat.jl")
 include("quadrature/strategies/nonconftestbaryrefoftrialqstrat.jl")
 include("quadrature/strategies/timedomain/nothingqstrat.jl")
+include("quadrature/quadraturebuffer.jl")
+include("quadrature/quadrulecallbacks.jl")
 
 
 include("excitation.jl")

src/coloring.jl

@@ -41,7 +41,7 @@ function GraphsColoring.conflicts(
     end
 
     for i in eachindex(conflictindices)
-        conflictindices[i] = unique(conflictindices[i])
+        unique!(conflictindices[i])
     end
 
     return eachindex(elements),

src/decoupled/dpops.jl

@@ -20,9 +20,11 @@ end
 function momintegrals!(out, op::CurlSingleLayerDP3D,
     test_local_space::RTRefSpace, tptr, test_triangular_element,
     trial_local_space::LagrangeRefSpace, bptr, trial_triangular_element,
-    qrule::SauterSchwabStrategy)
+    qrule::SauterSchwabStrategy, qbuffer)
 
-    I, J, K, L = SauterSchwabQuadrature.reorder(
+    sbuffer = sauterschwab_buffer(qbuffer, qrule)
+    I, J, K, L = sbuffer.I, sbuffer.J, sbuffer.K, sbuffer.L
+    sauterschwab_reorder!(I, J, K, L,
         test_triangular_element.vertices,
         trial_triangular_element.vertices, qrule)
 

src/integralop.jl

@@ -97,7 +97,7 @@ function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store;
     else
         quadstrat
     end
-    
+
     qd = quaddata(biop, tshapes, bshapes, test_elements, bsis_elements, qs)
     zlocal = zeros(scalartype(biop, tfs, bfs), 2num_tshapes, 2num_bshapes)
     # @show "after" qs
@@ -142,7 +142,6 @@ end
     @test BlockArrays.blocksizes(M) == [(n1,n1) (n1,n2); (n2,n1) (n2,n2)]
 end
 
-
 function assemblechunk_body!(biop, test_space, trial_space,
     test_elements, test_element_ptrs, test_assembly_data, active_test_els,
     trial_elements, trial_element_ptrs, trial_assembly_data, active_trial_els,
@@ -156,6 +155,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
         @local begin
             zlocal = zeros(scalartype(biop, test_space, trial_space), num_tshapes, num_bshapes)
             tadjq = Vector{eltype(trial_assembly_data.data)}(undef, size(trial_assembly_data.data,1))
+            qbuffer = quadraturebuffer(quadstrat)
         end
         P = active_test_els[p]
         tcell = test_elements[P]
@@ -166,10 +166,11 @@ function assemblechunk_body!(biop, test_space, trial_space,
             bptr = trial_element_ptrs[Q]
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, refspace(test_space), refspace(trial_space),
+            apply = ApplyMomintegrals(zlocal, biop,
+                test_space, tptr, tcell, trial_space, bptr, bcell,
+                qbuffer)
+            quadrule(apply, biop, refspace(test_space), refspace(trial_space),
                 P, tcell, Q, bcell, qd, quadstrat)
-            momintegrals!(zlocal, biop,
-                test_space,  tptr, tcell, trial_space, bptr, bcell, qrule)
             for j in 1 : num_bshapes
                 tadjq .= @view trial_assembly_data.data[:,j,q]
                 for i in 1 : num_tshapes
@@ -184,6 +185,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
 end end end end end end end
 
 
+
 # function assemblechunk_body_colored!(biop,
 #         test_space, testad, testelementids,
 #         trial_space, trialad, trialelementids,
@@ -202,7 +204,7 @@ end end end end end end end
 #         @set scheduler = scheduler
 #         @local zlocal = zeros(scalartype(biop, test_space, trial_space), num_tshapes, num_bshapes)
 #         tcell, tptr = test_elements[p], test_cell_ptrs[p]
-    
+
 #         for q in trialelementids
 #             bcell, bptr = trial_elements[q], trial_cell_ptrs[q]
 #             fill!(zlocal, 0)
@@ -432,7 +434,7 @@ end
 #                             m′ = get(test_id_in_blk, m, 0)
 #                             m′ == 0 && continue
 #                             store(a*zlocal[i,j]*b, m′, n′)
-#     end end end end end end 
+#     end end end end end end
 #     # put!(zlocals, zlocal)
 # end
 
@@ -629,6 +631,7 @@ function assemblerow_body!(biop,
     zlocal, quadrature_data, store; quadstrat)
 
     test_function = test_functions.fns[1]
+    qbuffer = quadraturebuffer(quadstrat)
     for shape in test_function
         p = shape.cellid
         i = shape.refid
@@ -637,11 +640,12 @@ function assemblerow_body!(biop,
         for (q,bcell) in enumerate(trial_elements)
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(zlocal, biop,
+            apply = ApplyMomintegrals(zlocal, biop,
                 test_functions, nothing, tcell,
                 trial_functions, nothing, bcell,
-                qrule)
+                qbuffer)
+            quadrule(apply, biop, test_shapes, trial_shapes,
+                p, tcell, q, bcell, quadrature_data, quadstrat)
 
             for j in 1:size(zlocal,2)
                 for (n,b) in trial_assembly_data[q,j]
@@ -682,6 +686,7 @@ function assemblecol_body!(biop,
     zlocal, quadrature_data, store; quadstrat)
 
     trial_function = trial_functions.fns[1]
+    qbuffer = quadraturebuffer(quadstrat)
     for shape in trial_function
         q = shape.cellid
         j = shape.refid
@@ -691,15 +696,14 @@ function assemblecol_body!(biop,
         for (p,tcell) in enumerate(test_elements)
 
             fill!(zlocal, 0)
-            qrule = quadrule(biop, test_shapes, trial_shapes, p, tcell, q, bcell, quadrature_data, quadstrat)
-            momintegrals!(zlocal, biop,
+            apply = ApplyMomintegrals(zlocal, biop,
                 test_functions, nothing, tcell,
-                trial_functions, nothing, bcell, qrule)
+                trial_functions, nothing, bcell,
+                qbuffer)
+            quadrule(apply, biop, test_shapes, trial_shapes,
+                p, tcell, q, bcell, quadrature_data, quadstrat)
 
             for i in 1:size(zlocal,1)
                 for (m,a) in test_assembly_data[p,i]
                     store(a*zlocal[i,j]*b, m, 1)
 end end end end end
-
-
-

src/maxwell/bogaertints.jl

@@ -1,7 +1,7 @@
 function momintegrals!(z, op::MWSingleLayer3D,
     g::RTRefSpace, tptr, t,
     f::RTRefSpace, bptr, s,
-    strat::BogaertStrategy)
+    strat::BogaertStrategy, qbuffer)
 
     T, GG  = GetIntegrals(t, s, op.gamma, strat)
 
@@ -58,7 +58,7 @@ end
 function momintegrals!(z, op::MWDoubleLayer3D,
     g::RTRefSpace, tptr, τ,
     f::RTRefSpace, bptr, σ,
-    strat::BogaertStrategy)
+    strat::BogaertStrategy, qbuffer)
 
     # Get the primitives
     r = τ.vertices

src/maxwell/sauterschwabints_bdm_rt.jl

@@ -60,9 +60,11 @@ end
 
 function momintegrals!(op::MWDoubleLayer3D,
     test_local_space::BDMRefSpace, trial_local_space::RTRefSpace,
-    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
+    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy, qbuffer)
 
-    I, J, K, L = SauterSchwabQuadrature.reorder(
+    sbuffer = sauterschwab_buffer(qbuffer, strat)
+    I, J, K, L = sbuffer.I, sbuffer.J, sbuffer.K, sbuffer.L
+    sauterschwab_reorder!(I, J, K, L,
         test_triangular_element.vertices,
         trial_triangular_element.vertices, strat)
 

src/maxwell/wiltonints.jl

@@ -5,6 +5,8 @@ mutable struct WiltonSERule{P,Q} <: SingularityExtractionRule
     regularpart_quadrule::Q
 end
 
+quadraturebuffer(qrule::WiltonSERule) = quadraturebuffer(qrule.regularpart_quadrule)
+
 function innerintegrals!(op::MWSingleLayer3DSng, p, g, f, t, s, z,
         strat::WiltonSERule, dx)
 

src/operator.jl

@@ -246,12 +246,13 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
 
     testelements, testad, trialelements, trialad, qdata, _ = assembleblock_primer(
         operator, testfunctions, trialfunctions; quadstrat=quadstrat)
-    
+
     # testad = (testelements, testad, 1:length(testelements))
-    # trialad = (trialelements, trialad, 1:length(trialelements)) 
+    # trialad = (trialelements, trialad, 1:length(trialelements))
 
     testelementcolors = color(testfunctions, scheduler; addata=(testelements, testad, 1:length(testelements)))
     trialelementcolors = [1:length(trialelements)]
+    # trialelementcolors = color(trialfunctions, scheduler; addata=(trialelements, trialad, 1:length(trialelements)))
 
     qs = if CompScienceMeshes.refines(geometry(testfunctions), geometry(trialfunctions))
         TestRefinesTrialQStrat(quadstrat)
@@ -265,7 +266,7 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
 
     for (i, coloredtestelements) in enumerate(testelementcolors)
         testad1 = AssemblyData(testad.data[:,:,coloredtestelements])
-        for (j, coloredtrialelements) in enumerate(trialelementcolors)            
+        for (j, coloredtrialelements) in enumerate(trialelementcolors)
             # assemblechunk_body_colored!(operator,
             #     testfunctions, testad, coloredtestelements,
             #     trialfunctions, trialad, coloredtrialelements,
@@ -276,7 +277,7 @@ function assemble!(operator::Operator, testfunctions::Space, trialfunctions::Spa
                 trialelements, eachindex(trialelements), trialad, coloredtrialelements,
                 qdata, nothing, store; quadstrat=qs, scheduler)
             next!(pbar; step = length(testelementcolors[i]) * length(trialelementcolors[j]))
-    end end 
+    end end
     finish!(pbar)
 end
 
@@ -347,7 +348,7 @@ function assemble!(op::AbstractOperator, tfs::DirectProductSpace, bfs::DirectPro
     store, threading=Threading{:multi};
     quadstrat=defaultquadstrat(op, tfs[1], bfs[1]),
     kwargs...)
-    
+
     I = Int[0]
     for s in tfs.factors push!(I, last(I) + numfunctions(s)) end
     for (i,s) in enumerate(tfs.factors)
@@ -377,7 +378,7 @@ function assemble!(op::BlockDiagonalOperator, U::DirectProductSpace, V::DirectPr
     store, threading=Threading{:multi};
     quadstrat = defaultquadstrat(op, U, V),
     kwargs...)
-    
+
     @assert length(U.factors) == length(V.factors)
     I = Int[0]; for u in U.factors push!(I, last(I) + numfunctions(u)) end
     J = Int[0]; for v in V.factors push!(J, last(J) + numfunctions(v)) end
@@ -404,7 +405,7 @@ function assemble!(op::BlockFullOperators, U::DirectProductSpace, V::DirectProdu
     store, threading;
     quadstrat = defaultquadstrat(op, U, V),
     kwargs...)
-    
+
     I = Int[0]; for u in U.factors push!(I, last(I) + numfunctions(u)) end
     J = Int[0]; for v in V.factors push!(J, last(J) + numfunctions(v)) end
 

src/operators/quasilocalops.jl

@@ -70,6 +70,7 @@ function assemblechunk!(op::QuasiLocalOperator, tfs::Space, bfs::Space, store321
 
     qd = quaddata(op, trefs, brefs, tels, bels, quadstrat)
     zlocal = zeros(T, num_trefs, num_brefs)
+    qbuffer = quadraturebuffer(quadstrat)
     tree = elementstree(bels, 1.1)
 
     δ = oprange(op)
@@ -93,10 +94,12 @@ function assemblechunk!(op::QuasiLocalOperator, tfs::Space, bfs::Space, store321
                 @assert q <= size(bad.data, 3)
 
                 fill!(zlocal, 0)
-                qrule = quadrule(op, trefs, brefs, p, tcell, q, bcell, qd, quadstrat)
-                momintegrals!(zlocal, op,
+                apply = ApplyMomintegrals(zlocal, op,
                     tfs, tptr, tcell,
-                    bfs, bptr, bcell, qrule)
+                    bfs, bptr, bcell,
+                    qbuffer)
+                quadrule(apply, op, trefs, brefs,
+                    p, tcell, q, bcell, qd, quadstrat)
 
                 for j in 1 : length(bad[q])
                     for i in 1 : length(tad[p])

src/quadrature/SauterSchwabQuadrature1D.jl

@@ -6,9 +6,9 @@ export SauterSchwabStrategy1D
 export CommonEdge, CommonVertex
 
 # functions
-export sauterschwab_parameterized1D, _MRWrules, reorder
+export sauterschwab_parameterized1D, _MRWrules, reorder, reorder!
 
 # -------- included files
 include("doublesauterschwabint.jl")
 
-end
+end

src/quadrature/commonfaceoverlappingedgeqstrat.jl

@@ -10,17 +10,23 @@ end
 function quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd,
     qs::CommonFaceOverlappingEdgeQStrat)
 
+    return quadrule(ReturnQuadrule(), a, 𝒳, 𝒴, i, τ, j, σ, qd, qs)
+end
+
+function quadrule(f::QuadruleCallback, a, 𝒳, 𝒴, i, τ, j, σ, qd,
+    qs::CommonFaceOverlappingEdgeQStrat)
+
     if CompScienceMeshes.overlap(τ, σ)
-        return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+        return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
     end
 
     for (i,λ) in pairs(faces(τ))
         for (j,μ) in pairs(faces(σ))
             if CompScienceMeshes.overlap(λ, μ)
-                return NonConformingTouchQRule(qs.conforming_qstrat, i, j)
+                return f(NonConformingTouchQRule(qs.conforming_qstrat, i, j))
     end end end
 
     # Either positive distance, common face, or common vertex, which can
     # be handled directly by the parent quadrature strategy
-    return quadrule(a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)    
-end
+    return quadrule(f, a, 𝒳, 𝒴, i, τ, j, σ, qd, qs.conforming_qstrat)
+end

src/quadrature/doublenumints.jl

@@ -10,19 +10,19 @@ momintegrals!(biop, tshs, bshs, tcell, bcell, interactions, strat)
 Function for the computation of moment integrals using simple double quadrature.
 """
 function momintegrals!(biop,
-    tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule)
+    tshs, bshs, tcell, bcell, z, strat::DoubleQuadRule, qbuffer)
 
     igd = Integrand(biop, tshs, bshs, tcell, bcell)
 
     womps = strat.outer_quad_points
     wimps = strat.inner_quad_points
-    
+
     for womp in womps
         tgeo = womp.point
         tvals = womp.value
         M = length(tvals)
         jx = womp.weight
-        
+
         for wimp in wimps
             bgeo = wimp.point
             bvals = wimp.value
@@ -43,4 +43,4 @@ function momintegrals!(biop,
 end
 
 
-_TransposedStrat(a::DoubleQuadRule) = a
+_TransposedStrat(a::DoubleQuadRule) = a