diff --git a/examples/PDEs/laplace_2d_annulus.jl b/examples/PDEs/laplace_2d_annulus.jl new file mode 100644 index 0000000..b8575b6 --- /dev/null +++ b/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) diff --git a/examples/PDEs/poisson_2d_basic.jl b/examples/PDEs/poisson_2d_basic.jl index b691fed..d385076 100644 --- a/examples/PDEs/poisson_2d_basic.jl +++ b/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) diff --git a/test/Project.toml b/test/Project.toml index 31de9f0..b6f8a51 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/test_examples_pde.jl b/test/test_examples_pde.jl index 58a676e..6fb6931 100644 --- a/test/test_examples_pde.jl +++ b/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, diff --git a/test/test_util.jl b/test/test_util.jl index 5e4632e..53f772c 100644 --- a/test/test_util.jl +++ b/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)