diff --git a/docs/src/examples/graded_mesh.jl b/docs/src/examples/graded_mesh.jl index 54598e05..6e6dd69c 100644 --- a/docs/src/examples/graded_mesh.jl +++ b/docs/src/examples/graded_mesh.jl @@ -1,26 +1,34 @@ using Inti using StaticArrays using Meshes -using GLMakie +using CairoMakie using IterativeSolvers -using FMMLIB2D +using HMatrices using LinearAlgebra -## Parameters +TEST = true +PLOT = false # can be quite slow at large k + +## Frequency k = 10π λ = 2π / k -meshsize = λ / 10 -TEST = false -PLOT = true # can be quite slow at large k -## Create the geometry as an equilateral +## Tolerance +tol = 1e-12 + +## Target points to test solution +target = [5 * SVector(cos(θ), sin(θ)) for θ in 0:2π/100:2π] +vals_pointsource = Dict() +vals_planewave = Dict() + +## Geometry definition Γ = let - P = 2 + P = 3 χ = Inti.kress_change_of_variables_periodic(P) + # χ = identity p1 = SVector(0.0, 0.0) p2 = SVector(1.0, 0.0) p3 = SVector(0.5, sqrt(3) / 2) - l1 = Inti.parametric_curve( t -> p1 + χ(t) * (p2 - p1) + 0.0 * SVector(0, sin(2π * t)), 0.0, @@ -31,73 +39,135 @@ PLOT = true # can be quite slow at large k Inti.Domain(l1, l2, l3) end -msh = Inti.meshgen(Γ; meshsize) -Q = Inti.Quadrature(msh; qorder = 3) +ppw = [1, 2, 4, 8, 16, 32, 64, 128, 256] -## +## Loop +for n in ppw + println("Computing with $n points-per-wavelength") + meshsize = λ / n + ## Create the geometry as an equilateral triangle -## Create integral operators -pde = Inti.Helmholtz(; dim = 2, k) -S, D = Inti.single_double_layer(; - pde, - target = Q, - source = Q, - compression = (method = :fmm, tol = 1e-6), - correction = (method = :dim,), -) -L = I / 2 + D - im * k * S + msh = Inti.meshgen(Γ; meshsize) + Q = Inti.Quadrature(msh; qorder = 2) + + # Make sure all quadrature nodes are different and weights nonzero (floating + # point paranoia) + qcoords = [q.coords for q in Q] + @assert length(unique!(qcoords)) == length(Q) + @assert all(q -> q.weight > 0, Q) + + println("|--- $(length(Q)) dofs") + + ## Create integral operators + pde = Inti.Helmholtz(; dim = 2, k) + S, D = Inti.single_double_layer(; + pde, + target = Q, + source = Q, + compression = (method = :hmatrix, tol), + # compression = (method = :none,), + # correction = ( + # method = :adaptive, + # tol = tol, + # maxdist = 5 * meshsize, + # maxsplit = 100_000, + # ), + correction = (method = :dim,), + ) + L = I / 2 + D - im * k * S + + 𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) -𝒮, 𝒟 = Inti.single_double_layer_potential(; pde, source = Q) + ## Test with a source inside for validation + if TEST + G = Inti.SingleLayerKernel(pde) + uₑ = x -> G(x, SVector(0.4, 0.4)) + rhs = map(q -> uₑ(q.coords), Q) + σ, hist = gmres(L, rhs; restart = 1000, maxiter = 1000, abstol = tol, log = true) + println("|--- converged in $(hist.iters)") + uₕ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) + vals_pointsource[n] = map(uₕ, target) + end -## Test with a source inside for validation -if TEST - G = Inti.SingleLayerKernel(pde) - uₑ = x -> G(x, SVector(0.4, 0.4)) - rhs = map(q -> uₑ(q.coords), Q) - σ = gmres(L, rhs; restart = 100, maxiter = 100, abstol = 1e-8, verbose = true) + ## Test with a planewave + θ = 0 * π / 4 + uᵢ = x -> exp(im * k * dot(x, SVector(cos(θ), sin(θ)))) + rhs = map(q -> -uᵢ(q.coords), Q) + σ, hist = gmres(L, rhs; restart = 1000, maxiter = 1000, abstol = tol, log = true) + println("|--- converged in $(hist.iters)") uₕ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) - xₜ = SVector(2, 0) - @show uₕ(xₜ), uₑ(xₜ) -end + vals_planewave[n] = map(uₕ, target) + ## Visualize a solution + if PLOT && n > 8 + PLOT = false + xx = -1:λ/10:2 + yy = -1:λ/10:2 + target = [SVector(x, y) for x in xx, y in yy] -## Test with a planewave -θ = 0 * π / 4 -uᵢ = x -> exp(im * k * dot(x, SVector(cos(θ), sin(θ)))) -rhs = map(q -> -uᵢ(q.coords), Q) + Spot, Dpot = Inti.single_double_layer(; + pde, + target = target, + source = Q, + compression = (method = :fmm, tol = 1e-4), + correction = (method = :none, target_location = :outside, maxdist = meshsize), + ) -σ = gmres(L, rhs; restart = 1000, maxiter = 1000, abstol = 1e-8, verbose = true) -uₕ = x -> 𝒟[σ](x) - im * k * 𝒮[σ](x) -@info "Done solving..." + vals = Dpot * σ - im * k * Spot * σ + vals = map( + i -> Inti.isinside(target[i], Q) ? NaN + im * NaN : vals[i] + uᵢ(target[i]), + 1:length(target), + ) + vals = reshape(vals, length(xx), length(yy)) -## Visualize a solution -if PLOT - xx = -1:λ/10:2 - yy = -1:λ/10:2 - target = [SVector(x, y) for x in xx, y in yy] + heatmap( + xx, + yy, + real.(vals); + colorrange = (-1, 1), + colormap = :inferno, + axis = (aspect = DataAspect(),), + ) + viz!(msh[Γ]; showsegments = true) + Makie.current_figure() + end +end - Spot, Dpot = Inti.single_double_layer(; - pde, - target = target, - source = Q, - compression = (method = :fmm, tol = 1e-4), - correction = (method = :none, target_location = :outside, maxdist = meshsize), - ) +## +nmax = ppw[end] +er_planewave = Float64[] +er_pointsource = Float64[] - vals = Dpot * σ - im * k * Spot * σ - vals = map( - i -> Inti.isinside(target[i], Q) ? NaN + im * NaN : vals[i] + uᵢ(target[i]), - 1:length(target), - ) - vals = reshape(vals, length(xx), length(yy)) - - heatmap( - xx, - yy, - real.(vals); - colorrange = (-1, 1), - colormap = :inferno, - axis = (aspect = DataAspect(),), - ) - viz!(msh[Γ]; showsegments = true) - Makie.current_figure() +for n in ppw + n == nmax && continue + push!(er_planewave, norm(vals_planewave[n] - vals_planewave[nmax], Inf)) + push!(er_pointsource, norm(vals_pointsource[n] - vals_pointsource[nmax], Inf)) end + +fig = Figure() +ax1 = Axis( + fig[1, 1]; + yscale = log10, + xscale = log2, + xlabel = "points-per-wavelength", + ylabel = "|u - uref|∞", +) +scatterlines!(ax1, ppw[1:end-1], er_pointsource; label = "error pointsource") +refslope = 3 +refvals = (1 ./ ppw[1:end-1] .^ refslope) * ppw[end-1]^refslope * er_pointsource[end] +scatterlines!(ax1, ppw[1:end-1], refvals; label = "slope $refslope") +axislegend() + +ax2 = Axis( + fig[1, 2]; + yscale = log10, + xscale = log2, + xlabel = "points-per-wavelength", + ylabel = "|u - uref|∞", +) +scatterlines!(ax2, ppw[1:end-1], er_planewave; label = "error planewave") +refslope = 3 +refvals = (1 ./ ppw[1:end-1] .^ refslope) * ppw[end-1]^refslope * er_planewave[end] +scatterlines!(ax2, ppw[1:end-1], refvals; label = "slope $refslope") +axislegend() + +fig