Add Laplace example on annulus (#52)

* add Laplace example on annulus

* format

examples/PDEs/laplace_2d_annulus.jl

@@ -0,0 +1,47 @@
+using KernelInterpolation
+using QuasiMonteCarlo: sample, HaltonSample
+using LinearAlgebra: norm
+using Plots
+
+# right-hand-side of Poisson equation is zero -> Laplace equation
+f(x, equations) = 0.0
+pde = PoissonEquation(f)
+
+# Create inner nodes of annulus-like domain
+function create_annulus_inner(n, inner_radius = 1.0, outer_radius = 2.0)
+    nodes_matrix = sample(n, [-outer_radius, -outer_radius], [outer_radius, outer_radius],
+                          HaltonSample())
+    nodes_halton = NodeSet(nodes_matrix')
+    # Remove nodes inside unit circle
+    to_delete = []
+    for i in eachindex(nodes_halton)
+        if norm(nodes_halton[i]) > outer_radius || norm(nodes_halton[i]) < inner_radius
+            push!(to_delete, i)
+        end
+    end
+    deleteat!(nodes_halton, to_delete)
+    return nodes_halton
+end
+
+inner_radius = 1.0
+outer_radius = 2.0
+nodeset_inner = create_annulus_inner(300, inner_radius, outer_radius)
+# boundary nodes are the nodes on circles with radius 1.0 and 2.0
+nodes_boundary_outer = random_hypersphere_boundary(80, outer_radius; dim = 2)
+nodes_boundary_inner = random_hypersphere_boundary(20, inner_radius; dim = 2)
+nodeset_boundary = merge(nodes_boundary_outer, nodes_boundary_inner)
+# Dirichlet boundary condition (here 0.0 at inner boundary and sin(x[1]) at outer boundary)
+g(x) = isapprox(norm(x), outer_radius) ? cos(5.0 * acos(0.5 * x[1])) : 0.0
+
+kernel = WendlandKernel{2}(3, shape_parameter = 0.3)
+sd = SpatialDiscretization(pde, nodeset_inner, g, nodeset_boundary, kernel)
+itp = solve_stationary(sd)
+
+many_nodes_inner = create_annulus_inner(800)
+many_nodes_boundary_outer = random_hypersphere_boundary(160, outer_radius; dim = 2)
+many_nodes_boundary_inner = random_hypersphere_boundary(40, inner_radius; dim = 2)
+many_nodes = merge(many_nodes_inner, many_nodes_boundary_outer, many_nodes_boundary_inner)
+
+OUT = "out"
+ispath(OUT) || mkpath(OUT)
+vtk_save(joinpath(OUT, "laplace_2d_annulus"), many_nodes, itp)

examples/PDEs/poisson_2d_basic.jl

@@ -2,11 +2,11 @@ using KernelInterpolation
 using Plots
 
 # right-hand-side of Poisson equation
-f(x, equations) = 5 / 4 * pi^2 * sin(pi * x[1]) * cos(pi * x[2] / 2)
+f(x, equations) = 5 / 4 * pi^2 * sinpi(x[1]) * cospi(x[2] / 2)
 pde = PoissonEquation(f)
 
 # analytical solution of equation
-u(x, equations) = sin(pi * x[1]) * cos(pi * x[2] / 2)
+u(x, equations) = sinpi(x[1]) * cospi(x[2] / 2)
 
 n = 10
 nodeset_inner = homogeneous_hypercube(n, (0.1, 0.1), (0.9, 0.9); dim = 2)

test/Project.toml

@@ -4,6 +4,7 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/test_examples_pde.jl

@@ -17,6 +17,12 @@ EXAMPLES_DIR = joinpath(examples_dir(), "PDEs")
                               pde_test=true)
     end
 
+    @ki_testset "laplace_2d_annulus.jl" begin
+        # No analytical solution available (don't compare l2 and linf norms)
+        @test_include_example(joinpath(EXAMPLES_DIR, "laplace_2d_annulus.jl"),
+                              pde_test=true, atol=1e-11)
+    end
+
     @ki_testset "anisotropic_elliptic_2d_basic.jl" begin
         @test_include_example(joinpath(EXAMPLES_DIR, "anisotropic_elliptic_2d_basic.jl"),
                               l2=0.3680733486177618, linf=0.09329545570900866,

test/test_util.jl

@@ -28,6 +28,8 @@ macro test_include_example(example, args...)
         # if present, compare l2 and linf against reference values
         if !isnothing($l2) || !isnothing($linf)
             if !$pde_test
+                # interpolation test
+                # assumes `many_nodes` and `values` are defined in the example
                 values_test = itp.(nodeset)
                 # Check interpolation at interpolation nodes
                 @test isapprox(norm(values .- values_test, Inf), 0;
@@ -39,6 +41,9 @@ macro test_include_example(example, args...)
                 @test isapprox(norm(many_values .- many_values_test, Inf), $linf;
                                atol = $atol, rtol = $rtol)
             else
+                # PDE test
+                # assumes `many_nodes`, `nodes_inner` and `nodeset_boundary` are defined in the example
+                # if `u` is defined, it is used to compare the solution (analytical solution or initial condition) using the l2 and linf norms
                 if pde isa KernelInterpolation.AbstractStationaryEquation
                     rhs_values = KernelInterpolation.rhs(nodeset_inner, pde)
                     for i in eachindex(nodeset_inner)
@@ -49,13 +54,17 @@ macro test_include_example(example, args...)
                     # Because of some namespace issues
                     itp2 = itp
 
-                    many_values = u.(many_nodes, Ref(pde))
+                    if @isdefined u
+                        many_values = u.(many_nodes, Ref(pde))
+                    end
                 elseif pde isa KernelInterpolation.AbstractTimeDependentEquation
                     t = last(tspan)
                     values_boundary = g.(t, nodeset_boundary)
                     itp2 = titp(t)
 
-                    many_values = u.(Ref(t), many_nodes, Ref(pde))
+                    if @isdefined u
+                        many_values = u.(Ref(t), many_nodes, Ref(pde))
+                    end
                 else
                     error("Unknown PDE type")
                 end
@@ -64,11 +73,13 @@ macro test_include_example(example, args...)
                     @test isapprox(itp2(node), value, atol = $atol, rtol = $rtol)
                 end
 
-                many_values_test = itp2.(many_nodes)
-                @test isapprox(norm(many_values .- many_values_test), $l2;
-                               atol = $atol, rtol = $rtol)
-                @test isapprox(norm(many_values .- many_values_test, Inf), $linf;
-                               atol = $atol, rtol = $rtol)
+                if @isdefined many_values
+                    many_values_test = itp2.(many_nodes)
+                    @test isapprox(norm(many_values .- many_values_test), $l2;
+                                   atol = $atol, rtol = $rtol)
+                    @test isapprox(norm(many_values .- many_values_test, Inf), $linf;
+                                   atol = $atol, rtol = $rtol)
+                end
             end
         end
         println("═"^100)