diff --git a/CHANGELOG.md b/CHANGELOG.md index 32b9e660..23574579 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,20 @@ KernelInterpolation.jl follows the interpretation of used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.3 from v0.2.x + +#### Added + +- General floating point support ([#121]). + +#### Changed + +- The functions `random_hypersphere` and `random_hypersphere_boundary` not require a `Tuple` for + the argument `center`. Before, e.g., a `Vector` was allowed ([#121]). +- The element type of `NodeSet`s will now always be converted to a floating point type, i.e., also when + integer values are passed. This is more consistent for an interpolation framework makes many things easier. + A similar approach is also used in the Meshes.jl/CoordRefSystems.jl ecosystem ([#121]). + ## Changes in the v0.2 lifecycle #### Added diff --git a/Project.toml b/Project.toml index f15923ab..9bc56b97 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" KernelInterpolationMeshesExt = "Meshes" [compat] -DiffEqCallbacks = "3, 4" +DiffEqCallbacks = "4" ForwardDiff = "0.10.36" LinearAlgebra = "1" Meshes = "0.52.1, 0.53" @@ -36,11 +36,11 @@ Printf = "1" Random = "1" ReadVTK = "0.2" RecipesBase = "1.3.4" -Reexport = "1.2" -SciMLBase = "2.56" +Reexport = "1.2.2" +SciMLBase = "2.78" SimpleUnPack = "1.1" SpecialFunctions = "2" -StaticArrays = "1.9" +StaticArrays = "1.9.7" TimerOutputs = "0.5.23" TrixiBase = "0.1.3" TypedPolynomials = "0.4.1" diff --git a/examples/interpolation/interpolation_2d_sphere.jl b/examples/interpolation/interpolation_2d_sphere.jl index 88f3257c..3c4fe06a 100644 --- a/examples/interpolation/interpolation_2d_sphere.jl +++ b/examples/interpolation/interpolation_2d_sphere.jl @@ -5,7 +5,7 @@ using Plots f(x) = x[1] * x[2] r = 2.0 -center = [-1.0, 2.0] +center = (-1.0, 2.0) n = 40 nodeset = random_hypersphere(n, r, center) nodeset_boundary = random_hypersphere_boundary(20, r, center) diff --git a/src/KernelInterpolation.jl b/src/KernelInterpolation.jl index 2afc5fb8..813d588c 100644 --- a/src/KernelInterpolation.jl +++ b/src/KernelInterpolation.jl @@ -13,7 +13,7 @@ module KernelInterpolation using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect using ForwardDiff: ForwardDiff -using LinearAlgebra: Symmetric, I, norm, tr, muladd, dot, diagind +using LinearAlgebra: Symmetric, I, norm, tr, dot, diagind using Printf: @sprintf using Random: Random using ReadVTK: VTKFile, get_points, get_point_data, get_data diff --git a/src/basis.jl b/src/basis.jl index 874d9e16..45c4bdf8 100644 --- a/src/basis.jl +++ b/src/basis.jl @@ -59,14 +59,14 @@ The basis functions are given by where `K` is the kernel and `x_j` are the nodes in `centers`. """ -struct StandardBasis{Kernel} <: AbstractBasis - centers::NodeSet +struct StandardBasis{Dim, RealT, Kernel} <: AbstractBasis + centers::NodeSet{Dim, RealT} kernel::Kernel function StandardBasis(centers::NodeSet, kernel::Kernel) where {Kernel} if dim(kernel) != dim(centers) throw(DimensionMismatch("The dimension of the kernel and the centers must be the same")) end - new{typeof(kernel)}(centers, kernel) + new{dim(centers), eltype(centers), typeof(kernel)}(centers, kernel) end end @@ -85,36 +85,35 @@ already includes polynomial augmentation of degree `m` defaulting to `order(kern which means that the [`kernel_matrix`](@ref) of this basis is the identity matrix making it suitable for interpolation. Since the basis already includes polynomials no additional polynomial augmentation is needed for interpolation with this basis. """ -struct LagrangeBasis{Kernel, I <: AbstractInterpolation, Monomials, PolyVars} <: +struct LagrangeBasis{Dim, RealT, Kernel, I <: AbstractInterpolation, Monomials, PolyVars} <: AbstractBasis - centers::NodeSet + centers::NodeSet{Dim, RealT} kernel::Kernel basis_functions::Vector{I} ps::Monomials xx::PolyVars - function LagrangeBasis(centers::NodeSet, kernel::Kernel; - m = order(kernel)) where {Kernel} + function LagrangeBasis(centers::NodeSet, kernel::AbstractKernel; + m = order(kernel)) if dim(kernel) != dim(centers) throw(DimensionMismatch("The dimension of the kernel and the centers must be the same")) end + RealT = eltype(centers) K = length(centers) - values = zeros(K) - values[1] = 1.0 + values = zeros(RealT, K) + values[1] = one(RealT) b = interpolate(centers, values, kernel; m = m) basis_functions = Vector{typeof(b)}(undef, K) basis_functions[1] = b for i in 2:K - values[i - 1] = 0.0 - values[i] = 1.0 + values[i - 1] = zero(RealT) + values[i] = one(RealT) basis_functions[i] = interpolate(centers, values, kernel; m = m) end # All basis functions have same polynomials ps = first(basis_functions).ps xx = first(basis_functions).xx - new{typeof(kernel), eltype(basis_functions), typeof(ps), typeof(xx)}(centers, - kernel, - basis_functions, - ps, xx) + new{dim(centers), eltype(centers), typeof(kernel), eltype(basis_functions), + typeof(ps), typeof(xx)}(centers, kernel, basis_functions, ps, xx) end end diff --git a/src/discretization.jl b/src/discretization.jl index 17dfef1f..106492fc 100644 --- a/src/discretization.jl +++ b/src/discretization.jl @@ -176,7 +176,8 @@ function rhs!(dc, c, semi, t) boundary_condition.(Ref(t), nodeset_boundary)] end # dc = -pde_boundary_matrix * c + rhs_vector - @trixi_timeit timer() "muladd" dc[:]=muladd(pde_boundary_matrix, -c, rhs_vector) + @trixi_timeit timer() "muladd" dc[:]=Base.muladd(pde_boundary_matrix, -c, + rhs_vector) end return nothing end diff --git a/src/interpolation.jl b/src/interpolation.jl index c1507326..8a9ca721 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -159,7 +159,7 @@ Otherwise, `nodeset` is set to `centers(basis)` or `centers`. A regularization can be applied to the kernel matrix using the `regularization` argument, cf. [`regularize!`](@ref). """ function interpolate(basis::AbstractBasis, values::Vector{RealT}, - nodeset::NodeSet{Dim} = centers(basis); + nodeset::NodeSet{Dim, RealT} = centers(basis); m = order(basis), regularization = NoRegularization()) where {Dim, RealT} @assert dim(basis) == Dim @@ -174,27 +174,29 @@ function interpolate(basis::AbstractBasis, values::Vector{RealT}, else system_matrix = least_squares_matrix(basis, nodeset, ps, regularization) end - b = [values; zeros(q)] + b = [values; zeros(RealT, q)] c = system_matrix \ b - return Interpolation(basis, nodeset, c, system_matrix, ps, xx) + return Interpolation{typeof(basis), dim(basis), eltype(nodeset), typeof(system_matrix), + typeof(ps), typeof(xx)}(basis, nodeset, c, system_matrix, ps, xx) end function interpolate(centers::NodeSet{Dim, RealT}, nodeset::NodeSet{Dim, RealT}, - values::AbstractVector{RealT}, kernel = GaussKernel{Dim}(); + values::AbstractVector{RealT}, kernel = GaussKernel{Dim, RealT}(); kwargs...) where {Dim, RealT} interpolate(StandardBasis(centers, kernel), values, nodeset; kwargs...) end function interpolate(centers::NodeSet{Dim, RealT}, - values::AbstractVector{RealT}, kernel = GaussKernel{Dim}(); + values::AbstractVector{RealT}, + kernel = GaussKernel{Dim}(; shape_parameter = RealT(1.0)); kwargs...) where {Dim, RealT} interpolate(StandardBasis(centers, kernel), values; kwargs...) end # Evaluate interpolant function (itp::Interpolation)(x) - s = 0 bas = basis(itp) c = kernel_coefficients(itp) + s = zero(eltype(x)) for j in eachindex(c) s += c[j] * bas[j](x) end diff --git a/src/kernel_matrices.jl b/src/kernel_matrices.jl index 4fc630d3..a075373b 100644 --- a/src/kernel_matrices.jl +++ b/src/kernel_matrices.jl @@ -76,7 +76,7 @@ function interpolation_matrix(basis::AbstractBasis, ps, regularize!(k_matrix, regularization) p_matrix = polynomial_matrix(centers(basis), ps) system_matrix = [k_matrix p_matrix - p_matrix' zeros(q, q)] + p_matrix' zeros(eltype(k_matrix), q, q)] return Symmetric(system_matrix) end @@ -112,7 +112,7 @@ function least_squares_matrix(basis::AbstractBasis, nodeset::NodeSet, ps, p_matrix1 = polynomial_matrix(nodeset, ps) p_matrix2 = polynomial_matrix(centers(basis), ps) system_matrix = [k_matrix p_matrix1 - p_matrix2' zeros(q, q)] + p_matrix2' zeros(eltype(k_matrix), q, q)] return system_matrix end diff --git a/src/kernels/radialsymmetric_kernel.jl b/src/kernels/radialsymmetric_kernel.jl index 336f2b9a..2092cc6c 100644 --- a/src/kernels/radialsymmetric_kernel.jl +++ b/src/kernels/radialsymmetric_kernel.jl @@ -278,21 +278,21 @@ function Base.show(io::IO, kernel::WendlandKernel{Dim}) where {Dim} kernel.shape_parameter, ", d = ", kernel.d, ")") end -function phi(kernel::WendlandKernel, r::Real) +function phi(kernel::WendlandKernel, r::RealT) where {RealT <: Real} a_r = kernel.shape_parameter * r if a_r >= 1 - return 0.0 + return RealT(0.0) end - l = floor(kernel.d / 2) + kernel.k + 1 + l = floor(Int, kernel.d / 2) + kernel.k + 1 if kernel.k == 0 return (1 - a_r)^l elseif kernel.k == 1 return (1 - a_r)^(l + 1) * ((l + 1) * a_r + 1) elseif kernel.k == 2 - return 1 / 3 * (1 - a_r)^(l + 2) * + return 1 // 3 * (1 - a_r)^(l + 2) * ((l^2 + 4 * l + 3) * a_r^2 + (3 * l + 6) * a_r + 3) elseif kernel.k == 3 - return 1 / 15 * (1 - a_r)^(l + 3) * + return 1 // 15 * (1 - a_r)^(l + 3) * ((l^3 + 9 * l^2 + 23 * l + 15) * a_r^3 + (6 * l^2 + 36 * l + 45) * a_r^2 + (15 * l + 45) * a_r + 15) end @@ -346,10 +346,10 @@ function Base.show(io::IO, kernel::WuKernel{Dim}) where {Dim} ", shape_parameter = ", kernel.shape_parameter, ")") end -function phi(kernel::WuKernel, r::Real) +function phi(kernel::WuKernel, r::RealT) where {RealT <: Real} a_r = kernel.shape_parameter * r if a_r >= 1 - return 0.0 + return RealT(0.0) end if kernel.l == 0 # k = 0 @@ -358,29 +358,29 @@ function phi(kernel::WuKernel, r::Real) if kernel.k == 0 return (1 - a_r)^3 * (a_r^2 + 3 * a_r + 1) elseif kernel.k == 1 - return 1 / 2 * (1 - a_r)^2 * (a_r + 2) + return 1 // 2 * (1 - a_r)^2 * (a_r + 2) end elseif kernel.l == 2 if kernel.k == 0 return (1 - a_r)^5 * (a_r^4 + 5 * a_r^3 + 9 * a_r^2 + 5 * a_r + 1) elseif kernel.k == 1 - return 1 / 4 * (1 - a_r)^4 * (3 * a_r^3 + 12 * a_r^2 + 16 * a_r + 4) + return 1 // 4 * (1 - a_r)^4 * (3 * a_r^3 + 12 * a_r^2 + 16 * a_r + 4) elseif kernel.k == 2 - return 1 / 8 * (1 - a_r)^3 * (3 * a_r^2 + 9 * a_r + 8) + return 1 // 8 * (1 - a_r)^3 * (3 * a_r^2 + 9 * a_r + 8) end elseif kernel.l == 3 if kernel.k == 0 - return 1 / 5 * (1 - a_r)^7 * + return 1 // 5 * (1 - a_r)^7 * (5 * a_r^6 + 35 * a_r^5 + 101 * a_r^4 + 147 * a_r^3 + 101 * a_r^2 + 35 * a_r + 5) elseif kernel.k == 1 - return 1 / 6 * (1 - a_r)^6 * + return 1 // 6 * (1 - a_r)^6 * (5 * a_r^5 + 30 * a_r^4 + 72 * a_r^3 + 82 * a_r^2 + 36 * a_r + 6) elseif kernel.k == 2 - return 1 / 8 * (1 - a_r)^5 * + return 1 // 8 * (1 - a_r)^5 * (5 * a_r^4 + 25 * a_r^3 + 48 * a_r^2 + 40 * a_r + 8) elseif kernel.k == 3 - return 1 / 16 * (1 - a_r)^4 * (5 * a_r^3 + 20 * a_r^2 + 29 * a_r + 16) + return 1 // 16 * (1 - a_r)^4 * (5 * a_r^3 + 20 * a_r^2 + 29 * a_r + 16) end end end @@ -544,8 +544,8 @@ function Base.show(io::IO, kernel::Matern32Kernel{Dim}) where {Dim} kernel.shape_parameter, ")") end -function phi(kernel::Matern32Kernel, r::Real) - y = sqrt(3) * kernel.shape_parameter * r +function phi(kernel::Matern32Kernel, r::RealT) where {RealT <: Real} + y = RealT(sqrt(3)) * kernel.shape_parameter * r return (1 + y) * exp(-y) end order(::Matern32Kernel) = 0 @@ -580,9 +580,9 @@ function Base.show(io::IO, kernel::Matern52Kernel{Dim}) where {Dim} print(io, "Matern52Kernel{", Dim, "}(shape_parameter = ", kernel.shape_parameter, ")") end -function phi(kernel::Matern52Kernel, r::Real) - y = sqrt(5) * kernel.shape_parameter * r - return 1 / 3 * (3 + 3 * y + y^2) * exp(-y) +function phi(kernel::Matern52Kernel, r::RealT) where {RealT <: Real} + y = RealT(sqrt(5)) * kernel.shape_parameter * r + return 1 // 3 * (3 + 3 * y + y^2) * exp(-y) end order(::Matern52Kernel) = 0 @@ -616,8 +616,8 @@ function Base.show(io::IO, kernel::Matern72Kernel{Dim}) where {Dim} print(io, "Matern72Kernel{", Dim, "}(shape_parameter = ", kernel.shape_parameter, ")") end -function phi(kernel::Matern72Kernel, r::Real) - y = sqrt(7) * kernel.shape_parameter * r +function phi(kernel::Matern72Kernel, r::RealT) where {RealT <: Real} + y = RealT(sqrt(7)) * kernel.shape_parameter * r return (1 + y + 6 * y^2 / 15 + y^3 / 15) * exp(-y) end order(::Matern72Kernel) = 0 diff --git a/src/kernels/special_kernel.jl b/src/kernels/special_kernel.jl index c5fdb53e..5f2dffe4 100644 --- a/src/kernels/special_kernel.jl +++ b/src/kernels/special_kernel.jl @@ -57,7 +57,7 @@ end function (kernel::ProductKernel)(x, y) @assert length(x) == length(y) - res = 1.0 + res = eltype(x)(1.0) for k in kernel.kernels res *= k(x, y) end diff --git a/src/nodes.jl b/src/nodes.jl index f8fcef0a..245b5981 100644 --- a/src/nodes.jl +++ b/src/nodes.jl @@ -5,7 +5,7 @@ Set of interpolation nodes. """ mutable struct NodeSet{Dim, RealT} nodes::Vector{MVector{Dim, RealT}} - q::Float64 + q::RealT end function Base.show(io::IO, nodeset::NodeSet) @@ -36,27 +36,30 @@ end # Constructors function NodeSet(nodes::Vector{MVector{Dim, RealT}}) where {Dim, RealT} q = separation_distance(nodes) - NodeSet{Dim, RealT}(nodes, q) + # Convert nodes to floats by design + NodeSet{Dim, float(RealT)}(nodes, q) end function NodeSet(nodes::Vector{SVector{Dim, RealT}}) where {Dim, RealT} - nodes = MVector.(nodes) - q = separation_distance(nodes) - NodeSet{Dim, RealT}(nodes, q) + NodeSet(MVector.(nodes)) end function NodeSet(nodes::AbstractVector{Vector{RealT}}) where {RealT} n = length(nodes) @assert n > 0 ndims = length(nodes[1]) @assert ndims > 0 - data = [MVector{ndims, RealT}(nodes[i]) for i in 1:n] + data = Vector{MVector{ndims, RealT}}(undef, n) + for i in 1:n + data[i] = MVector{ndims, RealT}(nodes[i]) + end return NodeSet(data) end function NodeSet(nodes::AbstractMatrix{RealT}) where {RealT} - NodeSet(Vector{eltype(nodes)}[eachrow(nodes)...]) + d = size(nodes, 2) + NodeSet(MVector{d, RealT}.(eachrow(nodes))) end function NodeSet(nodes::AbstractVector{RealT}) where {RealT} @assert length(nodes) > 0 - return NodeSet([[node] for node in nodes]) + return NodeSet(MVector{1}.([[node] for node in nodes])) end NodeSet(nodeset::NodeSet) = nodeset @@ -137,7 +140,7 @@ function Base.setindex!(nodeset::NodeSet{Dim, RealT}, v::Vector{RealT}, # could be done more efficiently update_separation_distance!(nodeset) end -function Base.setindex!(nodeset::NodeSet{1, RealT}, v::RealT, i::Int) where {RealT} +function Base.setindex!(nodeset::NodeSet{1, RealT}, v, i::Int) where {RealT} @assert dim(nodeset) == 1 nodeset.nodes[i] = [v] # update separation distance of nodeset because it possibly changed @@ -232,21 +235,21 @@ function random_hypercube(n, x_min = 0.0, x_max = 1.0; kwargs...) random_hypercube(Random.default_rng(), n, x_min, x_max; kwargs...) end -function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::Real = 0.0, - x_max::Real = 1.0; dim = 1) - nodes = x_min .+ (x_max - x_min) .* rand(rng, n, dim) +function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::RealT = 0.0, + x_max::RealT = 1.0; dim = 1) where {RealT} + nodes = x_min .+ (x_max - x_min) .* rand(rng, float(RealT), n, dim) return NodeSet(nodes) end -function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::NTuple{Dim}, - x_max::NTuple{Dim}; - dim = Dim) where {Dim} +function random_hypercube(rng::Random.AbstractRNG, n::Int, x_min::NTuple{Dim, RealT}, + x_max::NTuple{Dim, RealT}; + dim = Dim) where {Dim, RealT} @assert dim == Dim - nodes = rand(rng, n, dim) + nodes = rand(rng, float(RealT), n, dim) for i in 1:dim nodes[:, i] = x_min[i] .+ (x_max[i] - x_min[i]) .* view(nodes, :, i) end - return NodeSet(nodes) + return NodeSet(SVector{Dim}.(eachrow(nodes))) end """ @@ -262,8 +265,8 @@ function random_hypercube_boundary(n, x_min = 0.0, x_max = 1.0; kwargs...) random_hypercube_boundary(Random.default_rng(), n, x_min, x_max; kwargs...) end -function random_hypercube_boundary(rng::Random.AbstractRNG, n::Int, x_min::Real = 0.0, - x_max::Real = 1.0; dim = 1) +function random_hypercube_boundary(rng::Random.AbstractRNG, n::Int, x_min::RealT = 0.0, + x_max::RealT = 1.0; dim = 1) where {RealT} random_hypercube_boundary(rng, n, ntuple(_ -> x_min, dim), ntuple(_ -> x_max, dim)) end @@ -286,7 +289,7 @@ function random_hypercube_boundary(rng::Random.AbstractRNG, n::Int, x_min::NTupl x_max::NTuple{Dim}; dim = Dim) where {Dim} @assert dim == Dim - if dim == 1 && n >= 2 + if Dim == 1 && n >= 2 @warn "For one dimension the boundary of the hypercube consists only of 2 points" return NodeSet([x_min[1], x_max[1]]) end @@ -326,13 +329,13 @@ function homogeneous_hypercube(n::Int, x_min::NTuple{Dim}, x_max::NTuple{Dim}; end function homogeneous_hypercube(n::NTuple{Dim, Int}, - x_min::NTuple{Dim} = ntuple(_ -> 0.0, Dim), - x_max::NTuple{Dim} = ntuple(_ -> 1.0, Dim); - dim = Dim) where {Dim} + x_min::NTuple{Dim, RealT} = ntuple(_ -> 0.0, Dim), + x_max::NTuple{Dim, RealT} = ntuple(_ -> 1.0, Dim); + dim = Dim) where {Dim, RealT} @assert Dim == dim - nodes = Vector{MVector{Dim, Float64}}(undef, prod(n)) + nodes = Vector{MVector{Dim, float(RealT)}}(undef, prod(n)) for (i, indices) in enumerate(Iterators.product(ntuple(j -> 1:n[j], Dim)...)) - node = MVector{Dim, Float64}(undef) + node = MVector{Dim, float(RealT)}(undef) for j in 1:dim node[j] = x_min[j] + (x_max[j] - x_min[j]) * (indices[j] - 1) / (n[j] - 1) end @@ -379,15 +382,15 @@ end # TODO: Is there a better way to create these `NodeSet`s? function homogeneous_hypercube_boundary(n::NTuple{Dim}, - x_min::NTuple{Dim} = ntuple(_ -> 0.0, Dim), - x_max::NTuple{Dim} = ntuple(_ -> 1.0, Dim); - dim = Dim) where {Dim} - if dim == 1 && n[1] >= 2 + x_min::NTuple{Dim, RealT} = ntuple(_ -> 0.0, Dim), + x_max::NTuple{Dim, RealT} = ntuple(_ -> 1.0, Dim); + dim = Dim) where {Dim, RealT} + if Dim == 1 && n[1] >= 2 # @warn "For one dimension the boundary of the hypercube consists only of 2 points" return NodeSet([x_min[1], x_max[1]]) end @assert Dim == dim - nodes = Vector{MVector{Dim, Float64}}(undef, number_of_nodes(n, Dim)) + nodes = Vector{MVector{Dim, float(RealT)}}(undef, number_of_nodes(n, Dim)) local i = 1 # Left side is like homogeneous hypercube in `dim - 1` hypercube for indices in Iterators.product(ntuple(j -> 1:n[j + 1], Dim - 1)...) @@ -421,10 +424,10 @@ function homogeneous_hypercube_boundary(n::NTuple{Dim}, end """ - random_hypersphere([rng], n, r = 1.0, center = zeros(dim); [dim]) + random_hypersphere([rng], n, r = 1.0, center = Tuple(zeros(dim)); [dim]) Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` inside a hypersphere with -radius `r` around the center `center`. +radius `r` around the center `center`, which is given as a tuple. If `dim` is not given explicitly, it is inferred by the length of `center` if possible. Optionally, pass a random number generator `rng`. """ @@ -436,25 +439,28 @@ function random_hypersphere(n, r, center; kwargs...) random_hypersphere(Random.default_rng(), n, r, center; kwargs...) end -function random_hypersphere(rng::Random.AbstractRNG, n, r = 1.0; dim = 2) - random_hypersphere(rng, n, r, zeros(dim)) +function random_hypersphere(rng::Random.AbstractRNG, n, r::RealT = 1.0; + dim = 2) where {RealT} + random_hypersphere(rng, n, r, Tuple(zeros(float(RealT), dim))) end -function random_hypersphere(rng::Random.AbstractRNG, n, r, - center; dim = length(center)) - @assert length(center) == dim - nodes = randn(rng, n, dim) +function random_hypersphere(rng::Random.AbstractRNG, n, r::RealT, + center::NTuple{Dim, RealT}; dim = Dim) where {Dim, RealT} + @assert Dim == dim + nodes = randn(rng, float(RealT), n, dim) for i in 1:n - nodes[i, :] .= center .+ r .* nodes[i, :] ./ norm(nodes[i, :]) * rand(rng)^(1 / dim) + nodes[i, :] .= center .+ + r .* nodes[i, :] ./ norm(nodes[i, :]) * + rand(rng, float(RealT))^(1 / dim) end - return NodeSet(nodes) + return NodeSet(SVector{Dim}.(eachrow(nodes))) end """ - random_hypersphere_boundary([rng], n, r = 1.0, center = zeros(dim); [dim]) + random_hypersphere_boundary([rng], n, r = 1.0, center = Tuple(zeros(dim)); [dim]) Create a [`NodeSet`](@ref) with `n` random nodes each of dimension `dim` at the boundary of a -hypersphere with radius `r` around the center `center`. +hypersphere with radius `r` around the center `center`, which is given as a tuple. If `dim` is not given explicitly, it is inferred by the length of `center` if possible. Optionally, pass a random number generator `rng`. """ @@ -466,21 +472,22 @@ function random_hypersphere_boundary(n, r, center; kwargs...) random_hypersphere_boundary(Random.default_rng(), n, r, center; kwargs...) end -function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r = 1.0; dim = 2) - random_hypersphere_boundary(rng, n, r, zeros(dim)) +function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r::RealT = 1.0; + dim = 2) where {RealT} + random_hypersphere_boundary(rng, n, r, Tuple(zeros(float(RealT), dim))) end -function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r, - center; - dim = length(center)) - @assert length(center) == dim - if dim == 1 && n >= 2 +function random_hypersphere_boundary(rng::Random.AbstractRNG, n, r::RealT, + center::NTuple{Dim, RealT}; + dim = Dim) where {Dim, RealT} + @assert Dim == dim + if Dim == 1 && n >= 2 @warn "For one dimension the boundary of the hypersphere consists only of 2 points" return NodeSet([-r, r]) end - nodes = randn(rng, n, dim) + nodes = randn(rng, float(RealT), n, dim) for i in 1:n nodes[i, :] .= center .+ r .* nodes[i, :] ./ norm(nodes[i, :]) end - return NodeSet(nodes) + return NodeSet(SVector{Dim}.(eachrow(nodes))) end diff --git a/src/util.jl b/src/util.jl index 1714cccc..02ef552e 100644 --- a/src/util.jl +++ b/src/util.jl @@ -42,4 +42,4 @@ end # Create `d` polyvars from `TypedPolynomials.jl`, don't use `@polyvars` because of # https://github.com/JuliaAlgebra/TypedPolynomials.jl/issues/51, instead use the # workaround from there -polyvars(d) = ntuple(i -> Variable{Symbol("x[", i, "]")}(), d) +polyvars(d) = ntuple(i -> Variable{Symbol("x[", i, "]")}(), Val(d)) diff --git a/test/Project.toml b/test/Project.toml index 9b4ab7c2..8fc62337 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,15 +16,15 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] Aqua = "0.8.3" -ExplicitImports = "1.0.1" +ExplicitImports = "1.6" LinearAlgebra = "1" Meshes = "0.52.1, 0.53" -OrdinaryDiffEqNonlinearSolve = "1.3" -OrdinaryDiffEqRosenbrock = "1.3" +OrdinaryDiffEqNonlinearSolve = "1.6" +OrdinaryDiffEqRosenbrock = "1.9" Plots = "1.25.11" QuasiMonteCarlo = "0.3.1" Random = "1" -StaticArrays = "1.9" +StaticArrays = "1.9.7" Test = "1" TestItemRunner = "1" TestItems = "1" diff --git a/test/test_aqua.jl b/test/test_aqua.jl index b7f4103d..003f0e68 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -1,8 +1,14 @@ @testitem "Aqua.jl" begin import Aqua - using ExplicitImports: check_no_implicit_imports, check_no_stale_explicit_imports - Aqua.test_all(ambiguities = false, KernelInterpolation) + using ExplicitImports: check_no_implicit_imports, check_no_stale_explicit_imports, + check_all_explicit_imports_via_owners, + check_all_qualified_accesses_via_owners, + check_no_self_qualified_accesses + # We need `unbound_args = false` due to https://github.com/JuliaTesting/Aqua.jl/issues/265#issuecomment-2173168334 + Aqua.test_all(KernelInterpolation, ambiguities = false, unbound_args = false) @test isnothing(check_no_implicit_imports(KernelInterpolation)) @test isnothing(check_no_stale_explicit_imports(KernelInterpolation; ignore = (:trixi_include,))) + @test isnothing(check_all_qualified_accesses_via_owners(KernelInterpolation)) + @test isnothing(check_no_self_qualified_accesses(KernelInterpolation)) end diff --git a/test/test_unit.jl b/test/test_unit.jl index 2eae6287..02cbb73d 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -283,17 +283,18 @@ end @test nodeset3[1] == [1.0, 2.0] @test_nowarn nodeset3[1] = MVector{2}([2.0, 3.0]) @test_nowarn nodeset3[1] = [2.0, 3.0] - nodeset4 = @test_nowarn similar(nodeset1, Int64) - @test nodeset4 isa NodeSet{2, Int64} + nodeset4 = @test_nowarn similar(nodeset1, Float32) + @test nodeset4 isa NodeSet{2, Float32} nodeset5 = @test_nowarn similar(nodeset1, 10) @test nodeset5 isa NodeSet{2, Float64} @test length(nodeset5) == 10 - nodeset6 = @test_nowarn similar(nodeset1, Int64, 10) - @test nodeset6 isa NodeSet{2, Int64} + nodeset6 = @test_nowarn similar(nodeset1, Float32, 10) + @test nodeset6 isa NodeSet{2, Float32} @test length(nodeset6) == 10 + # Integer nodes should be converted to float by design nodeset7 = @test_nowarn NodeSet(1:4) - @test eltype(nodeset7) == Int64 + @test eltype(nodeset7) == Float64 @test dim(nodeset7) == 1 @test length(nodeset7) == 4 @test size(nodeset7) == (4, 1) @@ -571,7 +572,7 @@ end end r = 2.0 - center = [-1.0, 3.0, 2.0, -pi] + center = (-1.0, 3.0, 2.0, -pi) nodeset13 = @test_nowarn random_hypersphere(50, r, center) @test nodeset13 isa NodeSet{4, Float64} for node in nodeset13 @@ -1106,6 +1107,69 @@ end @test isapprox(titp(t, x), u2(t, x, pde), atol = 0.12) end +@testitem "Different floating point types" setup=[Setup] begin + # Special nodesets + @test eltype(@inferred random_hypercube(10, (0.5f0, 0.5f0), (1.0f0, 1.0f0))) == Float32 + @test eltype(@inferred random_hypercube_boundary(10, (0.5f0, 0.5f0), (1.0f0, 1.0f0))) == + Float32 + @test eltype(@inferred homogeneous_hypercube(10, (0.5f0, 0.5f0), (1.0f0, 1.0f0))) == + Float32 + @test eltype(@inferred homogeneous_hypercube_boundary(10, (0.5f0, 0.5f0), + (1.0f0, 1.0f0))) == Float32 + @test eltype(@inferred random_hypersphere(10, 1.0f0, (1.0f0, 1.0f0))) == Float32 + @test eltype(@inferred random_hypersphere_boundary(10, 1.0f0, (1.0f0, 1.0f0))) == + Float32 + + # Interpolation with `StandardBasis` + centers = NodeSet(Float32[0.0 0.0 + 1.0 0.0 + 0.0 1.0 + 1.0 1.0]) + @test eltype(centers) == Float32 + kernel = MultiquadricKernel{dim(centers)}(; shape_parameter = 0.5f0) + f(x) = x[1] + x[2] + ff = f.(centers) + itp = @test_nowarn interpolate(centers, ff, kernel) + @test eltype(coefficients(itp)) == Float32 + @test eltype(system_matrix(itp)) == Float32 + @test typeof(@inferred itp([0.5f0, 0.5f0])) == Float32 + + # Interpolation with `LagrangeBasis` + basis = @test_nowarn LagrangeBasis(centers, kernel) + K = kernel_matrix(basis) + @test eltype(K) == Float32 + nodes = NodeSet(Float32[0.0 0.0 + 1.0 0.0 + 0.5 0.5 + 0.0 1.0 + 1.0 1.0]) + ff = f.(nodes) + itp = @test_nowarn interpolate(basis, ff, nodes) + @test eltype(coefficients(itp)) == Float32 + @test eltype(system_matrix(itp)) == Float32 + @test typeof(@inferred itp([0.5f0, 0.5f0])) == Float32 + + # Solving stationary PDE + nodeset_inner = NodeSet(Float32[0.25 0.25 + 0.75 0.25 + 0.25 0.75 + 0.75 0.75]) + u1(x) = x[1] * (x[1] - 1) + (x[2] - 1) * x[2] + f1(x, equations) = -4.0f0 # -Δu + nodeset_boundary = NodeSet(Float32[0.0 0.0 + 1.0 0.0 + 0.0 1.0 + 1.0 1.0]) + g1(x) = u1(x) + kernel = WendlandKernel{2}(3; shape_parameter = 0.5f0) + pde = PoissonEquation(f1) + sd = SpatialDiscretization(pde, nodeset_inner, g1, nodeset_boundary, kernel) + itp = @test_nowarn solve_stationary(sd) + @test eltype(coefficients(itp)) == Float32 + @test eltype(system_matrix(itp)) == Float32 + @test typeof(@inferred itp([0.5f0, 0.5f0])) == Float32 +end + @testitem "Callbacks" setup=[Setup, AdditionalImports] begin # AliveCallback alive_callback = AliveCallback(dt = 0.1)