Feature/SIMD_General_Quadstrat

CHANGELOG.md

@@ -1,8 +1,12 @@
-# Changelog
+# New in 2.10.0
 
 - `solve` and `solve!`  on `GMRESSolver` allow for the defaults preconditioners to be overwritten
 - `materialize` keyword to assemble for bilinear forms allows for fine-grained control over the assembly of operators (e.g. dense representations vs H-matrix representations)
 - Restart for GMRES defaults to the maximum number of iterations
+- `assemble` of bilinear forms checks for recurring blocks and reuses them.
+
+# New in 2.9.0
+
 - Cell coloring based lock-free multi-threaded assembly for frequency domain integral operators
 - Higher order Lagrange elements (cx and c0) on segments
 - All LinearMaps can be cast into bilinear forms. Assembly of the latter trivially returns the underlying LinearMap. This enables for example the construction inverses that explicitly exploit block diagonal structure.

Project.toml

@@ -23,6 +23,7 @@ Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LiftedMaps = "d22a30c1-52ac-4762-a8c9-5838452405e0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NestedUnitRanges = "032820ab-dc03-4b49-91f4-7d58d4da98b3"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -61,10 +62,11 @@ IterativeSolvers = "0.9"
 Krylov = "0.10.1"
 LiftedMaps = "0.5.1"
 LinearMaps = "3.7 - 3.9, 3.11.2"
+LoopVectorization = "0.12.170"
 NestedUnitRanges = "0.2.3"
 OhMyThreads = "0.8.3"
-ProgressMeter = "1.11.0"
 PlotlyBase = "0.8.21"
+ProgressMeter = "1.11.0"
 Requires = "1"
 SauterSchwab3D = "0.2"
 SauterSchwabQuadrature = "2.4.0"
@@ -74,5 +76,3 @@ Suppressor = "0.2.8"
 TestItems = "0.1.1, 1"
 WiltonInts84 = "0.2.8"
 julia = "1.10"
-
-

docs/make.jl

@@ -27,6 +27,8 @@ makedocs(;
             "Setting the Quadrature Strategy" => "manual/quadstrat.md",
             "Custom Quadrature Rules" => "manual/quadrule.md",
             "Custom Operators" => "manual/customop.md",
+            "Nested Quadrature Strategies" => "manual/nestedquadstrat.md",
+            "SIMDDoubleNum Quadrature Strategy" => "manual/SIMDquadstrat.md",
             "Application Examples"=>Any[
                 "Time-Harmonic"=>Any[
                     "EFIE"=>"manual/examplesTH/efie.md",

docs/src/manual/SIMDquadstrat.md

@@ -0,0 +1,22 @@
+# SIMDDoubleNum Quadstrat
+This quadrature strategy uses the SIMD properties of the hardware to accelerate the double-num interactions. It should be used in combination with a quadrature strategy for the near interacting triangels, for example with the sacuterschwab quadrature.
+The SIMD quadtarture is implemented for the Maxwell3D.singlelayer operator assembled in a raviart-thomas basis on triangles and for the composed operator structures at any basis embedded in a 3D space.
+The Maxwell3D.doublelayer operator is assembled by casting the structure to the correpsonding compsoed operator.
+Example
+```julia
+    op = Maxwell3D.singlelayer(wavenumber=1.0)
+    assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6))
+    op = Maxwell3D.doublelayer(wavenumber=1.0)
+    assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6))
+    op = Maxwell3D.doublelayer(wavenumber=1.0)
+    assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6))
+    op = BEAST.CompDoubleInt(x->x, BEAST.Dot(),BEAST.HH3DGreen(1.0*im),BEAST.STimesV(), x->x)
+    assemble(op,X,X,quadstrat = BEAST.SIMDDoubleNumSauterQstrat(6,6,6,6,6,6))
+```
+
+Note: The operators used in the operator fields of the CompDoubleInt should be of the specific types shown below
+vector scalar multiplication: BEAST.VTimesS
+scalar vector multiplication: BEAST.STimesV
+scalar scalar multiplication: BEAST.STimesS
+dot product: BEAST.Dot
+cross product: BEAST.Cross

docs/src/manual/nestedquadstrat.md

@@ -0,0 +1,13 @@
+# Nested Quadrature Strategies
+
+Some quadrature strategies are only useble for a subest of all simplex-simplex interactions. If two simplicces do fullfill the right conditions this quadrature strategie returns a specified quadrature rule. If not, it calls the nested quadrature strategy on the simplex-simplex pair. This nested quadrature strategy is stored in the nested_strat field that a nested quadrature strategy must have.
+An example of a nested quadstrat is given by
+```julia
+struct SauterQStrat{NestedStrat,S} <: NestedQuadStrat
+    nested_strat::NestedStrat
+    sauter_schwab_common_tetr::S
+    sauter_schwab_common_face::S
+    sauter_schwab_common_edge::S
+    sauter_schwab_common_vert::S
+end
+```

docs/src/manual/quadstrat.md

@@ -118,4 +118,6 @@ Z3 = assemble(𝑇, RT, RT; quadstrat=myquadstrat2)
 
 Best practice is too return `BEAST.defaultquadstrat(op, testnfs, trialfns)` by default to ensure that all operators are supported, also those for which no explicit overwrite is specified.
 
+# Quadrature cache
 
+The `quadcache` function is called for each thread seperately and ats writable working memory to the quadrature data.

examples/pmchwt.jl

@@ -59,12 +59,25 @@ h = (n × H) × n
 @hilbertspace k l
 
 α, α′ = 1/η, 1/η′
-pmchwt = @discretise(
-    (η*T+η′*T′)[k,j] -      (K+K′)[k,m] +
-         (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l],
-    j∈X, m∈X, k∈X, l∈X)
+# pmchwt = @discretise(
+#     (η*T+η′*T′)[k,j] -      (K+K′)[k,m] +
+#          (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l],
+#     j∈X, m∈X, k∈X, l∈X)
 
-u = solve(pmchwt)
+a = (
+    η*T[k,j] + η′*T′[k,j] - K[k,m] - K′[k,m]
+    + K[l,j] + K′[l,j]    + α*T[l,m] + α′*T′[l,m])
+l = -e[k] - h[l]
+
+𝕏 = X × X
+A = assemble(a, 𝕏, 𝕏; threading=:cellcoloring)
+b = assemble(l, 𝕏)
+
+A⁻¹ = BEAST.GMRESSolver(A, reltol=1e-5, maxiter=1000)
+u = A⁻¹ * b
+
+# error()
+# u = solve(pmchwt)
 
 #preconditioner
 #=

src/BEAST.jl

@@ -340,6 +340,7 @@ include("composedoperators/displacementmesh.jl")
 include("composedoperators/potentials.jl")
 include("composedoperators/trace.jl")
 include("composedoperators/analytic_excitation.jl")
+include("quadrature/simddoublenum.jl")
 
 const x̂ = point(1, 0, 0)
 const ŷ = point(0, 1, 0)

src/bases/local/bdm3dlocal.jl

@@ -1,6 +1,7 @@
 struct BDM3DRefSpace{T} <: RefSpace{T} end
 
 numfunctions(f::BDM3DRefSpace, ch::CompScienceMeshes.ReferenceSimplex{3}) = 12
+numfunctions(f::Type{<:BDM3DRefSpace}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 12
 
 function (f::BDM3DRefSpace)(p)
 

src/bases/local/bdmlocal.jl

@@ -26,6 +26,7 @@ end
 divergence(ref::BDMRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))]
 
 numfunctions(x::BDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
+numfunctions(x::Type{<:BDMRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 6
 
 const _vert_perms_bdm = [
     (1,2,3),

src/bases/local/gwpdivlocal.jl

@@ -6,6 +6,10 @@ function numfunctions(x::GWPDivRefSpace{<:Any,D},
     dom::CompScienceMeshes.ReferenceSimplex{2}) where {D}
     (D+1)*(D+3)
 end
+function numfunctions(x::Type{<:GWPDivRefSpace{<:Any,D}},
+    dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {D}
+    (D+1)*(D+3)
+end
 function dimtype(x::GWPDivRefSpace{<:Any,D},
     dom::CompScienceMeshes.ReferenceSimplex{2}) where {D}
     Val{(D+1)*(D+3)}

src/bases/local/gwplocal.jl

@@ -6,6 +6,10 @@ function numfunctions(x::GWPCurlRefSpace{<:Any,D},
         dom::CompScienceMeshes.ReferenceSimplex{2}) where {D}
         (D+1)*(D+3)
 end
+function numfunctions(x::Type{GWPCurlRefSpace{<:Any,D}},
+        dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {D}
+        (D+1)*(D+3)
+end
 function dimtype(x::GWPCurlRefSpace{<:Any,D},
     dom::CompScienceMeshes.ReferenceSimplex{2}) where {D}
     Val{(D+1)*(D+3)}

src/bases/local/laglocal.jl

@@ -9,6 +9,14 @@ numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{
 numfunctions(s::LagrangeRefSpace{T,2,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 6
 numfunctions(s::LagrangeRefSpace{T,Dg}, ch::CompScienceMeshes.ReferenceSimplex{D}) where {T,Dg,D} = binomial(D+Dg,Dg)
 
+
+
+numfunctions(s::Type{<:LagrangeRefSpace{T,D,2}}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{1}}) where {T,D} = D+1
+numfunctions(s::LagrangeRefSpace{T,0,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 1
+numfunctions(s::LagrangeRefSpace{T,1,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 3
+numfunctions(s::LagrangeRefSpace{T,2,3}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) where {T} = 6
+numfunctions(s::LagrangeRefSpace{T,Dg}, ch::Type{<:CompScienceMeshes.ReferenceSimplex{D}}) where {T,Dg,D} = binomial(D+Dg,Dg)
+
 # valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
 #         SVector{numfunctions(ref), Tuple{T,T}}
 valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T

src/bases/local/localcomposedbasis.jl

@@ -1,6 +1,8 @@
 abstract type _LocalBasisOperations{T} <: RefSpace{T} end
 numfunctions(a::_LocalBasisOperations) = coalesce(numfunctions(a.el1) , numfunctions(a.el2))
 numfunctions(a::_LocalBasisOperations,simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp))
+numfunctions(a::Type{<:_LocalBasisOperations}) = coalesce(numfunctions(a.el1) , numfunctions(a.el2))
+numfunctions(a::Type{<:_LocalBasisOperations},simp) = coalesce(numfunctions(a.el1,simp) , numfunctions(a.el2,simp))
 # struct TriangleSupport end
 # struct TetraSupport end
 

src/bases/local/ncrossbdmlocal.jl

@@ -20,3 +20,4 @@ function (f::NCrossBDMRefSpace{T})(p) where T
 end
 
 numfunctions(x::NCrossBDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
+numfunctions(x::Type{<:NCrossBDMRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 6

src/bases/local/nd2local.jl

@@ -27,4 +27,5 @@ function (ϕ::ND2RefSpace)(nbd)
 
 end
 
-numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8
+numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8
+numfunctions(x::Type{<:ND2RefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 8

src/bases/local/ndlcclocal.jl

@@ -59,6 +59,7 @@ function (ϕ::NDLCCRefSpace)(ndlc)
 end
 
 numfunctions(x::NDLCCRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 6
+numfunctions(x::Type{<:NDLCCRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 6
 
 #check orientation
 function curl(ref::NDLCCRefSpace, sh, el)

src/bases/local/ndlcdlocal.jl

@@ -31,6 +31,7 @@ function (ϕ::NDLCDRefSpace)(ndlc)
 end
 
 numfunctions(x::NDLCDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 4
+numfunctions(x::Type{<:NDLCDRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{3}}) = 4
 
 function ntrace(x::NDLCDRefSpace, el, q, fc)
     t = zeros(scalartype(x),1,4)

src/bases/local/ndlocal.jl

@@ -33,6 +33,7 @@ function (ϕ::NDRefSpace)(nbd)
 end
 
 numfunctions(x::NDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3
+numfunctions(x::Type{<:NDRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 3
 
 function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T
 

src/bases/local/rt2local.jl

@@ -34,6 +34,7 @@ function (f::RT2RefSpace)(p)
 end
 
 numfunctions(x::RT2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8
+numfunctions(x::Type{<:RT2RefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 8
 
 function interpolate(fields, interpolant::BEAST.RT2RefSpace, chart)
 

src/bases/local/rtlocal.jl

@@ -27,6 +27,7 @@ function (ϕ::RTRefSpace)(mp)
 end
 
 numfunctions(x::RTRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3
+numfunctions(x::Type{<:RTRefSpace}, dom::Type{<:CompScienceMeshes.ReferenceSimplex{2}}) = 3
 
 divergence(ref::RTRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/volume(el))]
 

src/bases/local/rtqlocal.jl

@@ -23,6 +23,7 @@ function (ϕ::RTQuadRefSpace{T})(p) where {T}
 end
 
 function numfunctions(ϕ::RTQuadRefSpace, dom::CompScienceMeshes.RefQuadrilateral) 4 end
+function numfunctions(ϕ::Type{<:RTQuadRefSpace}, dom::Type{<:CompScienceMeshes.RefQuadrilateral}) 4 end
 
 function interpolate(fields, interpolant::RTQuadRefSpace{T}, chart) where {T}
 

src/helmholtz2d/hh2dops.jl

@@ -236,12 +236,12 @@ defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumSauterQstrat(30,3
 
 function quaddata(op::HelmholtzOperator2D,
     test_local_space::RefSpace, trial_local_space::RefSpace,
-    test_charts, trial_charts, qs::DoubleNumSauterQstrat)
+    test_charts, trial_charts, qs::SauterQStrat)
 
     T = coordtype(test_charts[1])
 
-    tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
-    bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
+    # tqd = quadpoints(test_local_space,  test_charts,  (qs.outer_rule,))
+    # bqd = quadpoints(trial_local_space, trial_charts, (qs.inner_rule,))
 
     leg = (
       convert.(NTuple{2,T},_legendre(qs.sauter_schwab_common_vert,0,1)),
@@ -255,22 +255,20 @@ function quaddata(op::HelmholtzOperator2D,
      convert.(NTuple{2,T},BEAST.SauterSchwabQuadrature1D._MRWrules(qs.sauter_schwab_common_face,0,1)),
     )
 
-    return (tpoints=tqd, bpoints=bqd, gausslegendre=leg, marokhlinwandura=mrw)
+    return (gausslegendre=leg, marokhlinwandura=mrw,
+    nestedqd = quaddata(op, test_local_space, trial_local_space, test_charts, trial_charts, qs.nested_strat))
 end
 
 function quadrule(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace,
     i, τ::CompScienceMeshes.Simplex{<:Any, 1},
     j, σ::CompScienceMeshes.Simplex{<:Any, 1},
-    qd, qs::DoubleNumSauterQstrat)
+    qd, qs::SauterQStrat)
 
     hits = _numhits(τ, σ)
     @assert hits <= 2
 
     hits == 2 && return BEAST.SauterSchwabQuadrature1D.CommonEdge(qd.marokhlinwandura[2], qd.gausslegendre[2])
     hits == 1 && return BEAST.SauterSchwabQuadrature1D.CommonVertex(qd.marokhlinwandura[1], qd.gausslegendre[1])
 
-    return DoubleQuadRule(
-        qd.tpoints[1,i],
-        qd.bpoints[1,j],
-    )
+    return quadrule(op, g, f, i, τ, j, σ, qd.nestedqd, qs.nested_strat)
 end

src/helmholtz3d/nitsche.jl

@@ -4,28 +4,28 @@ mutable struct NitscheHH3{T} <: MaxwellOperator3D{T,T}
     gamma::T
 end
 
-defaultquadstrat(::NitscheHH3, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumWiltonSauterQStrat(10,8,10,8,3,3,3,3)
+defaultquadstrat(::NitscheHH3, ::LagrangeRefSpace, ::LagrangeRefSpace) = DoubleNumQStrat(10,8)
 
-function quaddata(operator::NitscheHH3,
-    localtestbasis::LagrangeRefSpace,
-    localtrialbasis::LagrangeRefSpace,
-    testelements, trialelements, qs::DoubleNumWiltonSauterQStrat)
+# function quaddata(operator::NitscheHH3,
+#     localtestbasis::LagrangeRefSpace,
+#     localtrialbasis::LagrangeRefSpace,
+#     testelements, trialelements, qs::SauterQStrat)
 
-  tqd = quadpoints(localtestbasis,  testelements,  (qs.outer_rule_far,))
-  bqd = quadpoints(x -> localtrialbasis(x), trialelements, (qs.inner_rule_far,))
+# #   tqd = quadpoints(localtestbasis,  testelements,  (qs.outer_rule_far,))
+# #   bqd = quadpoints(localtrialbasis, trialelements, (qs.inner_rule_far,))
 
-  #return QuadData(tqd, bqd)
-  return (tpoints=tqd, bpoints=bqd)
-end
+#   #return QuadData(tqd, bqd)
+#   return (nestedqd = quaddata(operator, localtestbasis, localtrialbasis, testelements, trialelements, qs.nested_strat),)
+# end
 
-function quadrule(op::NitscheHH3, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd,
-        qs::DoubleNumWiltonSauterQStrat)
+# function quadrule(op::NitscheHH3, g::LagrangeRefSpace, f::LagrangeRefSpace, i, τ, j, σ, qd,
+#         qs::DoubleNumWiltonSauterQStrat)
 
-    DoubleQuadRule(
-        qd.tpoints[1,i],
-        qd.bpoints[1,j]
-    )
-end
+#     DoubleQuadRule(
+#         qd.tpoints[1,i],
+#         qd.bpoints[1,j]
+#     )
+# end
 
 
 struct KernelValsMaxwell3D{T,U,P,Q}

src/integralop.jl

@@ -141,7 +141,8 @@ end
     @test BlockArrays.blocksize(M) == (2,2)
     @test BlockArrays.blocksizes(M) == [(n1,n1) (n1,n2); (n2,n1) (n2,n2)]
 end
-
+# quadrule(biop, tshapes, bshapes, p, tcell, q, bcell, qd, qdcache, quadstrat) =
+#     quadrule(biop, tshapes, bshapes, p, tcell, q, bcell, qd, quadstrat)
 
 function assemblechunk_body!(biop, test_space, trial_space,
     test_elements, test_element_ptrs, test_assembly_data, active_test_els,
@@ -154,6 +155,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
     @tasks for p in eachindex(active_test_els)
         @set scheduler = scheduler
         @local begin
+            qd = quadcache(biop, refspace(test_space), refspace(trial_space), test_elements, trial_elements, qd, quadstrat)
             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))
         end
@@ -167,7 +169,7 @@ function assemblechunk_body!(biop, test_space, trial_space,
 
             fill!(zlocal, 0)
             qrule = quadrule(biop, refspace(test_space), refspace(trial_space),
-                P, tcell, Q, bcell, qd, quadstrat)
+                P, tcell, Q, bcell, qd,quadstrat)
             momintegrals!(zlocal, biop,
                 test_space,  tptr, tcell, trial_space, bptr, bcell, qrule)
             for j in 1 : num_bshapes