[RFC] Add flag skipassemblefars to HMatrix function to build only near field blocks

.gitignore

@@ -6,3 +6,4 @@
 /docs/build/
 lcov.info
 *.vscode
+/dev/

example/skipfars.jl

@@ -0,0 +1,96 @@
+using LinearAlgebra
+using CompScienceMeshes
+using BEAST
+using ParallelKMeans
+using H2Trees
+using AdaptiveCrossApproximation
+using Krylov
+using PlotlyJS
+
+struct BellCurveSingleLayerIOP{T,C} <: BEAST.IntegralOperator
+    width::T
+    gamma::C
+end
+
+function (igd::BEAST.Integrand{<:BellCurveSingleLayerIOP})(p,q,f,g)
+    α = igd.operator.width
+    γ = igd.operator.gamma
+    γ⁻¹ = 1 / γ
+
+    x = CompScienceMeshes.cartesian(p)
+    y = CompScienceMeshes.cartesian(q)
+    R = LinearAlgebra.norm(x-y)
+    G = exp(-(R/α)^2) * exp(-γ * R) / (4 * π * R) #* 4 / (sqrt(π) * α)
+
+    BEAST._integrands(f,g) do fi,gj
+        -γ * dot(fi.value, gj.value) * G - γ⁻¹ * dot(fi.divergence, gj.divergence) * G 
+    end
+end
+
+function BEAST.scalartype(::BellCurveSingleLayerIOP{T,C}) where {T,C}
+    promote_type(T, C)
+end
+
+
+
+Γ = meshsphere(1.0, 0.04);
+X = raviartthomas(Γ);
+Y = buffachristiansen(Γ);
+@show numfunctions(X)
+@show numfunctions(Y)
+
+κ = 1.0
+
+E = Maxwell3D.planewave(; direction=ẑ, polarization=x̂, wavenumber=κ)
+e = (n × E) × n
+rhs = BEAST.assemble(e, X)
+
+a = Maxwell3D.singlelayer(wavenumber=κ);
+d = BEAST.NCross();
+b1 = Maxwell3D.singlelayer(wavenumber=κ);
+
+A = BEAST.assemble(a,X,X; threading=:cellcoloring)
+A⁻¹ = BEAST.GMRESSolver(A; abstol=1e-4, reltol=1e-4, maxiter=1000)
+u0, ch0 = BEAST.solve(A⁻¹, rhs)
+
+D = BEAST.assemble(d,X,Y; threading=:cellcoloring)
+B1 = BEAST.assemble(b1,Y,Y; threading=:cellcoloring)
+
+D⁻¹ = BEAST.GMRESSolver(D; reltol=1e-10, maxiter=1000, verbose=false);
+D⁻ᵀ = BEAST.GMRESSolver(D'; reltol=1e-10, maxiter=1000, verbose=false);
+P1 = D⁻ᵀ * B1 * D⁻¹
+A⁻¹ = BEAST.GMRESSolver(A; reltol=1e-4, maxiter=1000, left_preconditioner=P1)
+u1, ch1 = BEAST.solve(A⁻¹, rhs)
+
+alpha = 0.1
+b2 = BellCurveSingleLayerIOP(alpha, κ); 
+@time B2 = BEAST.assemble(b2,Y,Y; threading=:cellcoloring)
+P2 = D⁻ᵀ * B2 * D⁻¹
+A⁻¹ = BEAST.GMRESSolver(A; reltol=1e-4, maxiter=1000, left_preconditioner=P2)
+u2, ch2 = BEAST.solve(A⁻¹, rhs)
+
+# error()
+
+function isnear(treea::H2Trees.BoundingBallTree, treeb::H2Trees.BoundingBallTree, nodea::Int, nodeb::Int)
+    ths = H2Trees.radius(treea, nodea)
+    shs = H2Trees.radius(treeb, nodeb)
+    dist = norm(H2Trees.center(treea, nodea) - H2Trees.center(treeb, nodeb)) - (ths + shs)
+    dist < 4 * alpha
+end
+
+ttree = KMeansTree(X.pos, 2; minvalues=100)
+tree = BlockTree(ttree, ttree)
+
+
+@time B3 = HMatrix(b2, Y, Y, tree;
+    tol=1e-4,
+    maxrank=40,
+    isnear=isnear,
+    spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
+    skipassemblefars=true
+)
+P3 = D⁻ᵀ * B3 * D⁻¹
+A⁻¹ = BEAST.GMRESSolver(A; reltol=1e-4, maxiter=1000, left_preconditioner=P3)
+u3, ch3 = BEAST.solve(A⁻¹, rhs)
+
+@show ch0 ch1 ch2 ch3

src/hmatrix/abstracthmatrix.jl

@@ -72,6 +72,7 @@ Assemble a hierarchical matrix approximation of an operator on test and trial sp
   - `nearmatrixdata`: optional data passed to near-field assembly
   - `farmatrixdata`: optional data passed to far-field assembly
   - `scheduler`: thread scheduler used for assembly
+  - `skipassemblefars`: if true, skip the assembly of far-field blocks
 
 # Returns
 
@@ -99,6 +100,7 @@ function HMatrix(
     nearmatrixdata=defaultmatrixdata(operator, testspace, trialspace),
     farmatrixdata=defaultfarmatrixdata(operator, testspace, trialspace),
     scheduler=DynamicScheduler(),
+    skipassemblefars=false
 )
     spaceordering(tree, testspace, trialspace)
 
@@ -113,18 +115,20 @@ function HMatrix(
         scheduler=scheduler,
     )
 
-    fars = assemblefars(
-        operator,
-        testspace,
-        trialspace,
-        tree,
-        spaceordering;
-        maxrank=maxrank,
-        compressor=compressor,
-        isnear=isnear,
-        matrixdata=farmatrixdata,
-        scheduler=scheduler,
-    )
+    fars = skipassemblefars ?
+        BlockSparseMatrix[] :
+        assemblefars(
+            operator,
+            testspace,
+            trialspace,
+            tree,
+            spaceordering;
+            maxrank=maxrank,
+            compressor=compressor,
+            isnear=isnear,
+            matrixdata=farmatrixdata,
+            scheduler=scheduler,
+        )
 
     return HMatrix{eltype(nears)}(nears, fars, (length(testspace), length(trialspace)))
 end

test/runtests.jl

@@ -16,6 +16,7 @@ using AdaptiveCrossApproximation
     include("test_acabeast.jl")
 
     include("test_hmatrix.jl")
+    include("test_skipfars.jl")
 end
 
 @testitem "Code quality (Aqua.jl)" begin

test/test_skipfars.jl

@@ -0,0 +1,44 @@
+using AdaptiveCrossApproximation, Test
+
+using LinearAlgebra
+using CompScienceMeshes
+using BEAST
+using ParallelKMeans
+using H2Trees
+
+@testset "HMatrix: skipfars" begin
+
+    m0 = meshsphere(radius=1.0, h=0.45)
+    m1 = CompScienceMeshes.translate(m0, [0.0, 0.0, -10.0])
+    m2 = CompScienceMeshes.translate(m0, [0.0, 0.0, +10.0])
+    m = CompScienceMeshes.weld(m1, m2)
+
+    X = raviartthomas(m)
+    n = div(numfunctions(X), 2)
+
+    a = Maxwell3D.singlelayer(; gamma=2.0);
+    ttree = KMeansTree(X.pos, 2; minvalues=100)
+    tree = BlockTree(ttree, ttree)
+
+    function isnear(treea::H2Trees.BoundingBallTree, treeb::H2Trees.BoundingBallTree, nodea::Int, nodeb::Int)
+        ths = H2Trees.radius(treea, nodea)
+        shs = H2Trees.radius(treeb, nodeb)
+        dist = norm(H2Trees.center(treea, nodea) - H2Trees.center(treeb, nodeb)) - (ths + shs)
+        dist < 5 * alpha
+    end
+
+    A = HMatrix(a, X, X, tree;
+        tol=1e-3,
+        maxrank=40,
+        isnear=AdaptiveCrossApproximation.isnear(),
+        spaceordering=AdaptiveCrossApproximation.PreserveSpaceOrder(),
+        skipassemblefars=true
+    )
+
+    I1 = findall(getindex.(X.pos, 3) .> 0.0)
+    I2 = findall(getindex.(X.pos, 3) .< 0.0)
+    u = ones(2n)
+    u[I1] .= 0
+    v = A * u
+    @test norm(v[I1]) ≈ 0 atol=1e-8
+end