diff --git a/.gitignore b/.gitignore index c3e81f4..3b0d2a5 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ /docs/build/ lcov.info *.vscode +/dev/ \ No newline at end of file diff --git a/example/skipfars.jl b/example/skipfars.jl new file mode 100644 index 0000000..f75b2d4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/hmatrix/abstracthmatrix.jl b/src/hmatrix/abstracthmatrix.jl index c8911d6..8fcfe7b 100644 --- a/src/hmatrix/abstracthmatrix.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 30a62ec..c9cd83a 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_skipfars.jl b/test/test_skipfars.jl new file mode 100644 index 0000000..77b8816 --- /dev/null +++ b/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