Mercurial > repos > public > sbplib_julia
changeset 1700:87d603499fc3 feature/sbp_operators/laplace_curvilinear
Merge feature/lazy_tensors/sparse_conversions
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 03 Sep 2024 08:38:27 +0200 |
| parents | 3eb4e584425c (diff) 3e9c3986930d (current diff) |
| children | 3684db043add |
| files | Project.toml |
| diffstat | 23 files changed, 1137 insertions(+), 51 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Mon Sep 02 15:35:54 2024 +0200 +++ b/Manifest.toml Tue Sep 03 08:38:27 2024 +0200 @@ -2,7 +2,7 @@ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "a36735c53cfa4453f39635046eeaa47a4ea1231b" +project_hash = "32fac879810099260f177c27318d3f26de4a00cc" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"]
--- a/Project.toml Mon Sep 02 15:35:54 2024 +0200 +++ b/Project.toml Tue Sep 03 08:38:27 2024 +0200 @@ -4,17 +4,20 @@ version = "0.1.1" [deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" [weakdeps] +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Tokens = "040c2ec2-8d69-4aca-bf03-7d3a7092f2f6" [extensions] +SbplibPlotsExt = "Plots" SbplibMakieExt = "Makie" SbplibSparseArrayKitExt = ["SparseArrayKit", "Tokens"] SbplibSparseArraysExt = ["SparseArrays", "Tokens"]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/SbplibPlotsExt.jl Tue Sep 03 08:38:27 2024 +0200 @@ -0,0 +1,72 @@ +module SbplibPlotsExt + +using Sbplib.Grids +using Plots + +@recipe f(::Type{<:Grid}, g::Grid) = map(Tuple,g)[:] + +@recipe function f(c::Chart{2,<:Rectangle}, n=5, m=n; draw_border=true, bordercolor=1) + Ξ = parameterspace(c) + ξs = range(limits(Ξ,1)..., n) + ηs = range(limits(Ξ,2)..., m) + + label := false + seriescolor --> 2 + for ξ ∈ ξs + @series adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1)) + end + + for η ∈ ηs + @series adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2)) + end + + if ~draw_border + return + end + + for ξ ∈ limits(Ξ,1) + @series begin + linewidth --> 3 + seriescolor := bordercolor + adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1)) + end + end + + for η ∈ limits(Ξ,2) + @series begin + linewidth --> 3 + seriescolor := bordercolor + adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2)) + end + end +end + +function adapted_curve_grid(g, minmax) + t1, _ = PlotUtils.adapted_grid(t->g(t)[1], minmax) + t2, _ = PlotUtils.adapted_grid(t->g(t)[2], minmax) + + ts = sort(vcat(t1,t2)) + + x = map(ts) do t + g(t)[1] + end + y = map(ts) do t + g(t)[2] + end + + return x, y +end + +# get_axis_limits(plt, :x) + + +# ReicpesPipline/src/user_recipe.jl +# @recipe function f(f::FuncOrFuncs{F}) where {F<:Function} + +# @recipe function f(f::Function, xmin::Number, xmax::Number) + +# _scaled_adapted_grid(f, xscale, yscale, xmin, xmax) + +end + +
--- a/src/Grids/Grids.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/Grids/Grids.jl Tue Sep 03 08:38:27 2024 +0200 @@ -1,8 +1,34 @@ +# TODO: Double check that the interfaces for indexing and iterating are fully implemented and tested for all grids. module Grids using Sbplib.RegionIndices using Sbplib.LazyTensors using StaticArrays +using LinearAlgebra + +export ParameterSpace +export HyperBox +export Simplex +export Interval +export Rectangle +export Box +export Triangle +export Tetrahedron + +export limits +export unitinterval +export unitsquare +export unitcube +export unithyperbox + +export verticies +export unittriangle +export unittetrahedron +export unitsimplex + +export Chart +export ConcreteChart +export parameterspace # Grid export Grid @@ -19,6 +45,7 @@ export eval_on export componentview export ArrayComponentView +export normal export BoundaryIdentifier export TensorGridBoundary @@ -33,11 +60,22 @@ export equidistant_grid +# MappedGrid +export MappedGrid +export jacobian +export logicalgrid +export mapped_grid +export jacobian_determinant +export metric_tensor +export metric_tensor_inverse + abstract type BoundaryIdentifier end +include("manifolds.jl") include("grid.jl") include("tensor_grid.jl") include("equidistant_grid.jl") include("zero_dim_grid.jl") +include("mapped_grid.jl") end # module
--- a/src/Grids/equidistant_grid.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/Grids/equidistant_grid.jl Tue Sep 03 08:38:27 2024 +0200 @@ -117,21 +117,36 @@ end """ - equidistant_grid(limit_lower::T, limit_upper::T, size::Int) + equidistant_grid(limit_lower::Number, limit_upper::Number, size::Int) Constructs a 1D equidistant grid. """ -function equidistant_grid(limit_lower::T, limit_upper::T, size::Int) where T - if any(size .<= 0) +function equidistant_grid(limit_lower::Number, limit_upper::Number, size::Int) + if size <= 0 throw(DomainError("size must be postive")) end - if any(limit_upper.-limit_lower .<= 0) + if limit_upper-limit_lower <= 0 throw(DomainError("side length must be postive")) end + return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead? end + +equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(hb.a, hb.b, dims...) +# TODO: One dimensional grids shouldn't have vector eltype right?, Change here or in HyperBox? + +function equidistant_grid(c::Chart, dims::Vararg{Int}) + lg = equidistant_grid(parameterspace(c), dims...) + return MappedGrid( + lg, + map(c,lg), + map(ξ->jacobian(c, ξ), lg), + ) +end + + CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
--- a/src/Grids/grid.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/Grids/grid.jl Tue Sep 03 08:38:27 2024 +0200 @@ -115,6 +115,8 @@ """ function boundary_identifiers end +# TBD: Boundary identifiers for charts and atlases? + """ boundary_grid(g::Grid, id::BoundaryIdentifier)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/manifolds.jl Tue Sep 03 08:38:27 2024 +0200 @@ -0,0 +1,207 @@ +""" + ParameterSpace{D} + +A space of parameters of dimension `D`. Used with `Chart` to indicate which +parameters are valid for that chart. + +Common parameter spaces are created using the functions unit sized spaces +* `unitinterval` +* `unitrectangle` +* `unitbox` +* `unittriangle` +* `unittetrahedron` +* `unithyperbox` +* `unitsimplex` + +See also: [`Interval`](@ref), [`Rectangle`](@ref), [`Box`](@ref), +[`Triangle`](@ref), [`Tetrahedron`](@ref), [`HyperBox`](@ref), +[`Simplex`](@ref), +""" +abstract type ParameterSpace{D} end +Base.ndims(::ParameterSpace{D}) where D = D +# TBD: Should implement domain_dim? + +struct HyperBox{T,D} <: ParameterSpace{D} + a::SVector{D,T} + b::SVector{D,T} +end + +function HyperBox(a,b) + T = SVector{length(a)} + HyperBox(convert(T,a), convert(T,b)) +end + +Interval{T} = HyperBox{T,1} +Rectangle{T} = HyperBox{T,2} +Box{T} = HyperBox{T,3} + +limits(box::HyperBox, d) = (box.a[d], box.b[d]) +limits(box::HyperBox) = (box.a, box.b) + +unitinterval(T=Float64) = unithyperbox(T,1) +unitsquare(T=Float64) = unithyperbox(T,2) +unitcube(T=Float64) = unithyperbox(T,3) +unithyperbox(T, D) = HyperBox((@SVector zeros(T,D)), (@SVector ones(T,D))) +unithyperbox(D) = unithyperbox(Float64,D) + + +struct Simplex{T,D,NV} <: ParameterSpace{D} + verticies::NTuple{NV,SVector{D,T}} +end + +Simplex(verticies::Vararg{AbstractArray}) = Simplex(Tuple(SVector(v...) for v ∈ verticies)) + +verticies(s::Simplex) = s.verticies + +Triangle{T} = Simplex{T,2} +Tetrahedron{T} = Simplex{T,3} + +unittriangle(T=Float64) = unitsimplex(T,2) +unittetrahedron(T=Float64) = unitsimplex(T,3) +function unitsimplex(T,D) + z = @SVector zeros(T,D) + unitelement = one(eltype(z)) + verticies = ntuple(i->setindex(z, unitelement, i), D) + return Simplex((z,verticies...)) +end +unitsimplex(D) = unitsimplex(Float64, D) + +""" + Chart{D} + +A parametrized description of a manifold or part of a manifold. +""" +struct Chart{D, PST<:ParameterSpace{D}, MT} + mapping::MT + parameterspace::PST +end + +domain_dim(::Chart{D}) where D = D +(c::Chart)(ξ) = c.mapping(ξ) +parameterspace(c::Chart) = c.parameterspace + +""" + jacobian(c::Chart, ξ) + +The jacobian of the mapping evaluated at `ξ`. This defers to the +implementation of `jacobian` for the mapping itself. If no implementation is +available one can easily be specified for either the mapping function or the +chart itself. +```julia +c = Chart(f, ps) +jacobian(f::typeof(f), ξ) = f′(ξ) +``` +or +```julia +c = Chart(f, ps) +jacobian(c::typeof(c),ξ) = f′(ξ) +``` +which will both allow calling `jacobian(c,ξ)`. +""" +jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ) +# TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian? + + +# TBD: Should Charts, parameterspaces have boundary names? + +""" + Atlas + +A collection of charts and their connections. +Should implement methods for `charts` and +""" +abstract type Atlas end + +""" + charts(::Atlas) + +The colloction of charts in the atlas. +""" +function charts end + +""" + connections + +TBD: What exactly should this return? + +""" + +struct CartesianAtlas <: Atlas + charts::Matrix{Chart} +end + +charts(a::CartesianAtlas) = a.charts + +struct UnstructuredAtlas <: Atlas + charts::Vector{Chart} + connections +end + +charts(a::UnstructuredAtlas) = a.charts + + +### +# Geometry +### + +abstract type Curve end +abstract type Surface end + + +struct Line{PT} <: Curve + p::PT + tangent::PT +end + +(c::Line)(s) = c.p + s*c.tangent + + +struct LineSegment{PT} <: Curve + a::PT + b::PT +end + +(c::LineSegment)(s) = (1-s)*c.a + s*c.b + + +function linesegments(ps...) + return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1] +end + + +function polygon_edges(ps...) + n = length(ps) + return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(Ps)] +end + +struct Circle{T,PT} <: Curve + c::PT + r::T +end + +(c::Circle)(θ) = c.c + r*@SVector[cos(Θ), sin(Θ)] + +struct TransfiniteInterpolationSurface{T1,T2,T3,T4} <: Surface + c₁::T1 + c₂::T2 + c₃::T3 + c₄::T4 +end + +function (s::TransfiniteInterpolationSurface)(u,v) + c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄ + P₀₀ = c₁(0) + P₁₀ = c₂(0) + P₁₁ = c₃(0) + P₀₁ = c₄(0) + return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - ( + (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁ + ) +end + +function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray) + s(ξ̄...) +end + +# TODO: Implement jacobian() for the different mapping helpers +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/mapped_grid.jl Tue Sep 03 08:38:27 2024 +0200 @@ -0,0 +1,142 @@ +struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D} + logicalgrid::GT + physicalcoordinates::CT + jacobian::JT +end + +jacobian(g::MappedGrid) = g.jacobian +logicalgrid(g::MappedGrid) = g.logicalgrid + + +# Indexing interface +Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...] +Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid) + +Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d) +Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d) +# TODO: axes(...)? + +# Iteration interface + +Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates) +Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state) + +Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}() +Base.length(g::MappedGrid) = length(g.logicalgrid) +Base.size(g::MappedGrid) = size(g.logicalgrid) +Base.size(g::MappedGrid, d) = size(g.logicalgrid, d) + +boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logicalgrid) +boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, id) + +function boundary_grid(g::MappedGrid, id::TensorGridBoundary) + b_indices = boundary_indices(g.logicalgrid, id) + + # Calculate indices of needed jacobian components + D = ndims(g) + all_indices = SVector{D}(1:D) + free_variable_indices = deleteat(all_indices, grid_id(id)) + jacobian_components = (:, free_variable_indices) + + # Create grid function for boundary grid jacobian + boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...) + boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...] + + return MappedGrid( + boundary_grid(g.logicalgrid, id), + boundary_physicalcoordinates, + boundary_jacobian, + ) +end + +# TBD: refine and coarsen could be implemented once we have a simple manifold implementation. +# Before we do, we should consider the overhead of including such a field in the mapped grid struct. + +function mapped_grid(x, J, size...) + D = length(size) + lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...) + return MappedGrid( + lg, + map(x,lg), + map(J,lg), + ) +end +# TODO: Delete this function + + +function jacobian_determinant(g::MappedGrid) + return map(jacobian(g)) do ∂x∂ξ + det(∂x∂ξ) + end +end + +function metric_tensor(g::MappedGrid) + return map(jacobian(g)) do ∂x∂ξ + ∂x∂ξ'*∂x∂ξ + end +end + +function metric_tensor_inverse(g::MappedGrid) + return map(jacobian(g)) do ∂x∂ξ + inv(∂x∂ξ'*∂x∂ξ) + end +end + +function min_spacing(g::MappedGrid{T,1} where T) + n, = size(g) + + ms = Inf + for i ∈ 1:n-1 + ms = min(ms, norm(g[i+1]-g[i])) + end + + return ms +end + +function min_spacing(g::MappedGrid{T,2} where T) + n, m = size(g) + + ms = Inf + for i ∈ 1:n-1, j ∈ 1:m-1 # loop over each cell of the grid + + ms = min( + ms, + norm(g[i+1,j]-g[i,j]), + norm(g[i,j+1]-g[i,j]), + + norm(g[i+1,j]-g[i+1,j+1]), + norm(g[i,j+1]-g[i+1,j+1]), + + norm(g[i+1,j+1]-g[i,j]), + norm(g[i+1,j]-g[i,j+1]), + ) + # NOTE: This could be optimized to avoid checking all interior edges twice. + end + + return ms +end + +""" + normal(g::MappedGrid, boundary) + +The outward pointing normal as a grid function on the boundary +""" +function normal(g::MappedGrid, boundary) + b_indices = boundary_indices(g, boundary) + σ =_boundary_sign(component_type(g), boundary) + return map(jacobian(g)[b_indices...]) do ∂x∂ξ + ∂ξ∂x = inv(∂x∂ξ) + k = grid_id(boundary) + σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:]) + end +end + +function _boundary_sign(T, boundary) + if boundary_id(boundary) == Upper() + return one(T) + elseif boundary_id(boundary) == Lower() + return -one(T) + else + throw(ArgumentError("The boundary identifier must be either `Lower()` or `Upper()`")) + end +end
--- a/src/Grids/tensor_grid.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/Grids/tensor_grid.jl Tue Sep 03 08:38:27 2024 +0200 @@ -1,3 +1,5 @@ +# TODO: Check this file and other grids for duplicate implementation of general methods implemented for Grid + """ TensorGrid{T,D} <: Grid{T,D}
--- a/src/SbpOperators/SbpOperators.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/SbpOperators.jl Tue Sep 03 08:38:27 2024 +0200 @@ -43,6 +43,7 @@ using Sbplib.RegionIndices using Sbplib.LazyTensors using Sbplib.Grids +using LinearAlgebra # Includes include("stencil.jl")
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Tue Sep 03 08:38:27 2024 +0200 @@ -25,3 +25,7 @@ converted_stencil = convert(Stencil{eltype(g)}, closure_stencil) return BoundaryOperator(g, converted_stencil, boundary) end + +function boundary_restriction(g::MappedGrid, stencil_set::StencilSet, boundary) + return boundary_restriction(logicalgrid(g), stencil_set, boundary) +end
--- a/src/SbpOperators/boundaryops/normal_derivative.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Tue Sep 03 08:38:27 2024 +0200 @@ -28,3 +28,36 @@ scaled_stencil = scale(closure_stencil,h_inv) return BoundaryOperator(g, scaled_stencil, boundary) end + +function normal_derivative(g::MappedGrid, stencil_set::StencilSet, boundary) + k = grid_id(boundary) + b_indices = boundary_indices(g, boundary) + + # Compute the weights for the logical derivatives + g⁻¹ = metric_tensor_inverse(g) + α = map(CartesianIndices(g⁻¹)[b_indices...]) do I # TODO: Fix iterator here + gᵏⁱ = g⁻¹[I][k,:] + gᵏᵏ = g⁻¹[I][k,k] + + gᵏⁱ./sqrt(gᵏᵏ) + end + + σ = ScalingTensor( + Grids._boundary_sign(component_type(g), boundary), + size(boundary_grid(g,boundary)), + ) + + + # Assemble difference operator + mapreduce(+,1:ndims(g)) do i + if i == k + ∂_ξᵢ = normal_derivative(logicalgrid(g), stencil_set, boundary) + else + e = boundary_restriction(logicalgrid(g), stencil_set, boundary) + ∂_ξᵢ = σ ∘ e ∘ first_derivative(logicalgrid(g), stencil_set, i) + end + + αᵢ = componentview(α,i) + DiagonalTensor(αᵢ) ∘ ∂_ξᵢ + end +end
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Tue Sep 03 08:38:27 2024 +0200 @@ -50,3 +50,8 @@ """ inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}() + +function inner_product(g::MappedGrid, stencil_set) + J = map(sqrt∘det, metric_tensor(g)) + DiagonalTensor(J)∘inner_product(logicalgrid(g), stencil_set) +end
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Tue Sep 03 08:38:27 2024 +0200 @@ -49,3 +49,8 @@ Implemented to simplify 1D code for SBP operators. """ inverse_inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}() + +function inverse_inner_product(g::MappedGrid, stencil_set) + J⁻¹ = map(inv∘sqrt∘det, metric_tensor(g)) + DiagonalTensor(J⁻¹)∘inverse_inner_product(logicalgrid(g), stencil_set) +end
--- a/src/SbpOperators/volumeops/laplace/laplace.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Tue Sep 03 08:38:27 2024 +0200 @@ -51,8 +51,31 @@ end return Δ end + laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set) +function laplace(grid::MappedGrid, stencil_set) + J = jacobian_determinant(grid) + J⁻¹ = DiagonalTensor(map(inv, J)) + + Jg = map(*, J, metric_tensor_inverse(grid)) + lg = logicalgrid(grid) + + return mapreduce(+, CartesianIndices(first(Jg))) do I + i, j = I[1], I[2] + Jgⁱʲ = componentview(Jg, i, j) + + if i == j + J⁻¹∘second_derivative_variable(lg, Jgⁱʲ, stencil_set, i) + else + Dᵢ = first_derivative(lg, stencil_set, i) + Dⱼ = first_derivative(lg, stencil_set, j) + J⁻¹∘Dᵢ∘DiagonalTensor(Jgⁱʲ)∘Dⱼ + end + end +end + + """ sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
--- a/test/Grids/equidistant_grid_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/Grids/equidistant_grid_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -2,6 +2,7 @@ using Test using Sbplib.RegionIndices using Sbplib.LazyTensors +using StaticArrays @testset "EquidistantGrid" begin @@ -112,6 +113,7 @@ @testset "equidistant_grid" begin @test equidistant_grid(0.0,1.0, 4) isa EquidistantGrid @test equidistant_grid((0.0,0.0),(8.0,5.0), 4, 3) isa TensorGrid + @test equidistant_grid((0.0,),(8.0,), 4) isa TensorGrid # constuctor @test_throws DomainError equidistant_grid(0.0, 1.0, 0) @@ -150,6 +152,23 @@ @test [gp[i]...] ≈ [p[i]...] atol=5e-13 end end + + + @testset "equidistant_grid(::ParameterSpace)" begin + ps = HyperBox((0,0),(2,1)) + + @test equidistant_grid(ps, 3,4) == equidistant_grid((0,0), (2,1), 3,4) + end + + + @testset "equidistant_grid(::Chart)" begin + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ, 3η] + end + Grids.jacobian(c::typeof(c), ξ̄) = @SMatrix[2 0; 0 3] + + @test equidistant_grid(c, 5, 4) isa Grid + end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/manifolds_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -0,0 +1,63 @@ +using Test + +using Sbplib.Grids +using Sbplib.RegionIndices +using Sbplib.LazyTensors + +# using StaticArrays + +@testset "ParameterSpace" begin + @test ndims(HyperBox([1,1], [2,2])) == 2 + @test ndims(unittetrahedron()) == 3 +end + +@testset "HyperBox" begin + @test HyperBox{<:Any, 2} <: ParameterSpace{2} + @test HyperBox([1,1], [2,2]) isa HyperBox{Int, 2} + + @test limits(HyperBox([1,2], [3,4])) == ([1,2], [3,4]) + @test limits(HyperBox([1,2], [3,4]), 1) == (1,3) + @test limits(HyperBox([1,2], [3,4]), 2) == (2,4) + + @test unitinterval() isa HyperBox{Float64,1} + @test limits(unitinterval()) == ([0], [1]) + + @test unitinterval(Int) isa HyperBox{Int,1} + @test limits(unitinterval(Int)) == ([0], [1]) + + @test unitsquare() isa HyperBox{Float64,2} + @test limits(unitsquare()) == ([0,0],[1,1]) + + @test unitcube() isa HyperBox{Float64,3} + @test limits(unitcube()) == ([0,0,0],[1,1,1]) + + @test unithyperbox(4) isa HyperBox{Float64,4} + @test limits(unithyperbox(4)) == ([0,0,0,0],[1,1,1,1]) +end + +@testset "Simplex" begin + @test Simplex{<:Any, 3} <: ParameterSpace{3} + @test Simplex([1,2], [3,4]) isa Simplex{Int, 2} + @test Simplex([1,2,3], [4,5,6],[1,1,1]) isa Simplex{Int, 3} + + @test verticies(Simplex([1,2], [3,4])) == ([1,2], [3,4]) + + @test unittriangle() isa Simplex{Float64,2} + @test verticies(unittriangle()) == ([0,0], [1,0], [0,1]) + + @test unittetrahedron() isa Simplex{Float64,3} + @test verticies(unittetrahedron()) == ([0,0,0], [1,0,0], [0,1,0],[0,0,1]) + + @test unitsimplex(4) isa Simplex{Float64,4} +end + +@testset "Chart" begin + c = Chart(x->2x, unitsquare()) + @test c isa Chart{2} + @test c([3,2]) == [6,4] + @test parameterspace(c) == unitsquare() +end + +@testset "Atlas" begin + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/mapped_grid_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -0,0 +1,278 @@ +using Sbplib.Grids +using Sbplib.RegionIndices +using Test +using StaticArrays +using LinearAlgebra + + +_skew_mapping(a,b) = (ξ̄->ξ̄[1]*a + ξ̄[2]*b, ξ̄->[a b]) + +function _partially_curved_mapping() + x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] + J((ξ, η)) = @SMatrix[ + 1 0; + η*(2ξ-1) 1+ξ*(ξ-1); + ] + + return x̄, J +end + +function _fully_curved_mapping() + x̄((ξ, η)) = @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + J((ξ, η)) = @SMatrix[ + 2 1-2η; + (2+η)*ξ 3+1/2*ξ^2; + ] + + return x̄, J +end + +@testset "MappedGrid" begin + lg = equidistant_grid((0,0), (1,1), 11, 11) # TODO: Change dims of the grid to be different + x̄ = map(ξ̄ -> 2ξ̄, lg) + J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) + mg = MappedGrid(lg, x̄, J) + + # TODO: Test constructor for different dims of range and domain for the coordinates + # TODO: Test constructor with different type than TensorGrid. a dummy type? + + @test_broken false # @test_throws ArgumentError("Sizes must match") MappedGrid(lg, map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)) + + + @test mg isa Grid{SVector{2, Float64},2} + + @test jacobian(mg) isa Array{<:AbstractMatrix} + @test logicalgrid(mg) isa Grid + + @testset "Indexing Interface" begin + mg = MappedGrid(lg, x̄, J) + @test mg[1,1] == [0.0, 0.0] + @test mg[4,2] == [0.6, 0.2] + @test mg[6,10] == [1., 1.8] + + @test mg[begin, begin] == [0.0, 0.0] + @test mg[end,end] == [2.0, 2.0] + @test mg[begin,end] == [0., 2.] + + @test eachindex(mg) == CartesianIndices((11,11)) + + @testset "cartesian indexing" begin + cases = [ + (1,1) , + (3,5) , + (10,6), + (1,1) , + (3,2) , + ] + + @testset "i = $is" for (lg, is) ∈ cases + @test mg[CartesianIndex(is...)] == mg[is...] + end + end + + @testset "eachindex" begin + @test eachindex(mg) == CartesianIndices((11,11)) + end + + @testset "firstindex" begin + @test firstindex(mg, 1) == 1 + @test firstindex(mg, 2) == 1 + end + + @testset "lastindex" begin + @test lastindex(mg, 1) == 11 + @test lastindex(mg, 2) == 11 + end + end + # TODO: Test with different types of logical grids + + @testset "Iterator interface" begin + sg = MappedGrid( + equidistant_grid((0,0), (1,1), 15, 11), + map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11) + ) + + @test eltype(mg) == SVector{2,Float64} + @test eltype(sg) == SVector{3,Float64} + + @test eltype(typeof(mg)) == SVector{2,Float64} + @test eltype(typeof(sg)) == SVector{3,Float64} + + @test size(mg) == (11,11) + @test size(sg) == (15,11) + + @test size(mg,2) == 11 + @test size(sg,2) == 11 + + @test length(mg) == 121 + @test length(sg) == 165 + + @test Base.IteratorSize(mg) == Base.HasShape{2}() + @test Base.IteratorSize(typeof(mg)) == Base.HasShape{2}() + + @test Base.IteratorSize(sg) == Base.HasShape{2}() + @test Base.IteratorSize(typeof(sg)) == Base.HasShape{2}() + + element, state = iterate(mg) + @test element == lg[1,1].*2 + element, _ = iterate(mg, state) + @test element == lg[2,1].*2 + + element, state = iterate(sg) + @test element == sg.physicalcoordinates[1,1] + element, _ = iterate(sg, state) + @test element == sg.physicalcoordinates[2,1] + + @test collect(mg) == 2 .* lg + end + + @testset "Base" begin + @test ndims(mg) == 2 + end + + @testset "boundary_identifiers" begin + @test boundary_identifiers(mg) == boundary_identifiers(lg) + end + + @testset "boundary_indices" begin + @test boundary_indices(mg, CartesianBoundary{1,Lower}()) == boundary_indices(lg,CartesianBoundary{1,Lower}()) + @test boundary_indices(mg, CartesianBoundary{2,Lower}()) == boundary_indices(lg,CartesianBoundary{2,Lower}()) + @test boundary_indices(mg, CartesianBoundary{1,Upper}()) == boundary_indices(lg,CartesianBoundary{1,Upper}()) + end + + @testset "boundary_grid" begin + x̄, J = _partially_curved_mapping() + mg = mapped_grid(x̄, J, 10, 11) + J1((ξ, η)) = @SMatrix[ + 1 ; + η*(2ξ-1); + ] + J2((ξ, η)) = @SMatrix[ + 0; + 1+ξ*(ξ-1); + ] + + function test_boundary_grid(mg, bId, Jb) + bg = boundary_grid(mg, bId) + + lg = logicalgrid(mg) + expected_bg = MappedGrid( + boundary_grid(lg, bId), + map(x̄, boundary_grid(lg, bId)), + map(Jb, boundary_grid(lg, bId)), + ) + + @testset let bId=bId, bg=bg, expected_bg=expected_bg + @test collect(bg) == collect(expected_bg) + @test logicalgrid(bg) == logicalgrid(expected_bg) + @test jacobian(bg) == jacobian(expected_bg) + # TODO: Implement equality of a curvilinear grid and simlify the above + end + end + + @testset test_boundary_grid(mg, TensorGridBoundary{1, Lower}(), J2) + @testset test_boundary_grid(mg, TensorGridBoundary{1, Upper}(), J2) + @testset test_boundary_grid(mg, TensorGridBoundary{2, Lower}(), J1) + @testset test_boundary_grid(mg, TensorGridBoundary{2, Upper}(), J1) + end +end + +@testset "mapped_grid" begin + x̄, J = _partially_curved_mapping() + mg = mapped_grid(x̄, J, 10, 11) + @test mg isa MappedGrid{SVector{2,Float64}, 2} + + lg = equidistant_grid((0,0), (1,1), 10, 11) + @test logicalgrid(mg) == lg + @test collect(mg) == map(x̄, lg) +end + +@testset "jacobian_determinant" begin + @test_broken false +end + +@testset "metric_tensor" begin + @test_broken false +end + +@testset "metric_tensor_inverse" begin + @test_broken false +end + +@testset "min_spacing" begin + let g = mapped_grid(identity, x->@SMatrix[1], 11) + @test min_spacing(g) ≈ 0.1 + end + + let g = mapped_grid(x->x+x.^2/2, x->@SMatrix[1 .+ x], 11) + @test min_spacing(g) ≈ 0.105 + end + + let g = mapped_grid(x->x + x.*(1 .- x)/2, x->@SMatrix[1.5 .- x], 11) + @test min_spacing(g) ≈ 0.055 + end + + let g = mapped_grid(identity, x->@SMatrix[1 0; 0 1], 11,11) + @test min_spacing(g) ≈ 0.1 + end + + let g = mapped_grid(identity, x->@SMatrix[1 0; 0 1], 11,21) + @test min_spacing(g) ≈ 0.05 + end + + + @testset let a = @SVector[1,0], b = @SVector[1,1]/√2 + g = mapped_grid(_skew_mapping(a,b)...,11,11) + + @test min_spacing(g) ≈ 0.1*norm(b-a) + end + + @testset let a = @SVector[1,0], b = @SVector[-1,1]/√2 + g = mapped_grid(_skew_mapping(a,b)...,11,11) + + @test min_spacing(g) ≈ 0.1*norm(a+b) + end +end + +@testset "normal" begin + g = mapped_grid(_partially_curved_mapping()...,10, 11) + + @test normal(g, CartesianBoundary{1,Lower}()) == fill(@SVector[-1,0], 11) + @test normal(g, CartesianBoundary{1,Upper}()) == fill(@SVector[1,0], 11) + @test normal(g, CartesianBoundary{2,Lower}()) == fill(@SVector[0,-1], 10) + @test normal(g, CartesianBoundary{2,Upper}()) ≈ map(boundary_grid(g,CartesianBoundary{2,Upper}())|>logicalgrid) do ξ̄ + α = 1-2ξ̄[1] + @SVector[α,1]/√(α^2 + 1) + end + + g = mapped_grid(_fully_curved_mapping()...,5,4) + + unit(v) = v/norm(v) + @testset let bId = CartesianBoundary{1,Lower}() + lbg = boundary_grid(logicalgrid(g), bId) + @test normal(g, bId) ≈ map(lbg) do (ξ, η) + -unit(@SVector[1/2, η/3-1/6]) + end + end + + @testset let bId = CartesianBoundary{1,Upper}() + lbg = boundary_grid(logicalgrid(g), bId) + @test normal(g, bId) ≈ map(lbg) do (ξ, η) + unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2)) + end + end + + @testset let bId = CartesianBoundary{2,Lower}() + lbg = boundary_grid(logicalgrid(g), bId) + @test normal(g, bId) ≈ map(lbg) do (ξ, η) + -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ)) + end + end + + @testset let bId = CartesianBoundary{2,Upper}() + lbg = boundary_grid(logicalgrid(g), bId) + @test normal(g, bId) ≈ map(lbg) do (ξ, η) + unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ)) + end + end +end
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -6,6 +6,8 @@ using Sbplib.RegionIndices using Sbplib.SbpOperators: BoundaryOperator, Stencil +using StaticArrays + @testset "boundary_restriction" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) e_closure = parse_stencil(stencil_set["e"]["closure"]) @@ -33,7 +35,7 @@ end @testset "Application" begin - @testset "1D" begin + @testset "EquidistantGrid" begin e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D)) v = eval_on(g_1D,x->1+x^2) u = fill(3.124) @@ -43,7 +45,7 @@ @test (e_r*v)[1] == v[end] end - @testset "2D" begin + @testset "TensorGrid" begin e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) v = rand(11, 15) u = fill(3.124) @@ -53,5 +55,22 @@ @test e_s*v == v[:,1] @test e_n*v == v[:,end] end + + @testset "MappedGrid" begin + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + end + Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2] + + mg = equidistant_grid(c, 10,13) + + e_w, e_e, e_s, e_n = boundary_restriction.(Ref(mg), Ref(stencil_set), boundary_identifiers(mg)) + v = rand(10, 13) + + @test e_w*v == v[1,:] + @test e_e*v == v[end,:] + @test e_s*v == v[:,1] + @test e_n*v == v[:,end] + end end end
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -6,54 +6,130 @@ using Sbplib.RegionIndices import Sbplib.SbpOperators.BoundaryOperator +using StaticArrays +using LinearAlgebra + @testset "normal_derivative" begin - g_1D = equidistant_grid(0.0, 1.0, 11) - g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 12) - @testset "normal_derivative" begin - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - @testset "1D" begin - d_l = normal_derivative(g_1D, stencil_set, Lower()) - @test d_l == normal_derivative(g_1D, stencil_set, Lower()) - @test d_l isa BoundaryOperator{T,Lower} where T - @test d_l isa LazyTensor{T,0,1} where T - end - @testset "2D" begin - d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) - d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}()) - Ix = IdentityTensor{Float64}((size(g_2D)[1],)) - Iy = IdentityTensor{Float64}((size(g_2D)[2],)) - d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower()) - d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper()) - @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) - @test d_w == d_l⊗Iy - @test d_n == Ix⊗d_r - @test d_w isa LazyTensor{T,1,2} where T - @test d_n isa LazyTensor{T,1,2} where T + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) + + @testset "EquidistantGrid" begin + g_1D = equidistant_grid(0.0, 1.0, 11) + + d_l = normal_derivative(g_1D, stencil_set, Lower()) + @test d_l == normal_derivative(g_1D, stencil_set, Lower()) + @test d_l isa BoundaryOperator{T,Lower} where T + @test d_l isa LazyTensor{T,0,1} where T + end + + @testset "TensorGrid" begin + g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 12) + d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) + d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}()) + Ix = IdentityTensor{Float64}((size(g_2D)[1],)) + Iy = IdentityTensor{Float64}((size(g_2D)[2],)) + d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower()) + d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper()) + @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) + @test d_w == d_l⊗Iy + @test d_n == Ix⊗d_r + @test d_w isa LazyTensor{T,1,2} where T + @test d_n isa LazyTensor{T,1,2} where T + + @testset "Accuracy" begin + v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y) + v∂x = eval_on(g_2D, (x,y)-> 2*x + y) + v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x) + # TODO: Test for higher order polynomials? + @testset "2nd order" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) + d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) + + @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 + @test d_e*v ≈ v∂x[end,:] atol = 1e-13 + @test d_s*v ≈ -v∂y[:,1] atol = 1e-13 + @test d_n*v ≈ v∂y[:,end] atol = 1e-13 + end + + @testset "4th order" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) + d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) + + @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 + @test d_e*v ≈ v∂x[end,:] atol = 1e-13 + @test d_s*v ≈ -v∂y[:,1] atol = 1e-13 + @test d_n*v ≈ v∂y[:,end] atol = 1e-13 + end end end - @testset "Accuracy" begin - v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y) - v∂x = eval_on(g_2D, (x,y)-> 2*x + y) - v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x) - # TODO: Test for higher order polynomials? - @testset "2nd order" begin - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) + + @testset "MappedGrid" begin + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + end + Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2] + mg = equidistant_grid(c, 10,13) + - @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 - @test d_e*v ≈ v∂x[end,:] atol = 1e-13 - @test d_s*v ≈ -v∂y[:,1] atol = 1e-13 - @test d_n*v ≈ v∂y[:,end] atol = 1e-13 + # x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] + # J((ξ, η)) = @SMatrix[ + # 1 0; + # η*(2ξ-1) 1+ξ*(ξ-1); + # ] + # mg = mapped_grid(x̄, J, 20, 21) + + + # x̄((ξ, η)) = @SVector[ξ,η] + # J((ξ, η)) = @SMatrix[ + # 1 0; + # 0 1; + # ] + # mg = mapped_grid(identity, J, 10, 11) + + for bid ∈ boundary_identifiers(mg) + @testset let bid=bid + @test normal_derivative(mg, stencil_set, bid) isa LazyTensor{<:Any, 1, 2} + end end - @testset "4th order" begin - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) - - @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 - @test d_e*v ≈ v∂x[end,:] atol = 1e-13 - @test d_s*v ≈ -v∂y[:,1] atol = 1e-13 - @test d_n*v ≈ v∂y[:,end] atol = 1e-13 + @testset "Consistency" begin + v = map(identity, mg) + + @testset "4nd order" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) + + for bid ∈ boundary_identifiers(mg) + @testset let bid=bid + d = normal_derivative(mg, stencil_set, bid) + @test d*v ≈ normal(mg, bid) rtol=1e-13 + end + end + end + end + + @testset "Accuracy" begin + v = function(x̄) + sin(norm(x̄+@SVector[1,1])) + end + ∇v = function(x̄) + ȳ = x̄+@SVector[1,1] + cos(norm(ȳ))*(ȳ/norm(ȳ)) + end + + mg = equidistant_grid(c, 80,80) + v̄ = map(v, mg) + + @testset for (order, atol) ∈ [(2,4e-2),(4,2e-3)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=order) + + @testset for bId ∈ boundary_identifiers(mg) + ∂ₙv = map(boundary_grid(mg,bId),normal(mg,bId)) do x̄,n̂ + n̂⋅∇v(x̄) + end + + dₙ = normal_derivative(mg, stencil_set, bId) + @test dₙ*v̄ ≈ ∂ₙv atol=atol + end + end end end end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -6,6 +6,9 @@ import Sbplib.SbpOperators.ConstantInteriorScalingOperator +using StaticArrays +using LinearAlgebra + @testset "Diagonal-stencil inner_product" begin Lx = π/2. Ly = Float64(π) @@ -94,4 +97,43 @@ end end end + + @testset "MappedGrid" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + end + Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2] + + mg = equidistant_grid(c, 10,13) + + @test inner_product(mg, stencil_set) isa LazyTensor{<:Any, 2,2} + + @testset "Accuracy" begin + v = function(x̄) + log(norm(x̄-@SVector[.5, .5]))/2π + log(norm(x̄-@SVector[1.5, 3]))/2π + end + ∇v = function(x̄) + ∇log(ȳ) = ȳ/(ȳ⋅ȳ) + ∇log(x̄-@SVector[.5, .5])/2π + ∇log(x̄-@SVector[1.5, 3])/2π + end + + mg = equidistant_grid(c, 80,80) + v̄ = map(v, mg) + + @testset for (order, atol) ∈ [(2,1e-3),(4,1e-7)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=order) + + @test sum(boundary_identifiers(mg)) do bId + ∂ₙv = map(boundary_grid(mg,bId),normal(mg,bId)) do x̄,n̂ + n̂⋅∇v(x̄) + end + Hᵧ = inner_product(boundary_grid(mg,bId), stencil_set) + sum(Hᵧ*∂ₙv) + end ≈ 2 atol=atol + + end + end + @test_broken false # Test that it calculates the right thing + end end
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -6,6 +6,8 @@ import Sbplib.SbpOperators.ConstantInteriorScalingOperator +using StaticArrays + @testset "Diagonal-stencil inverse_inner_product" begin Lx = π/2. Ly = Float64(π) @@ -82,4 +84,17 @@ end end end + + @testset "MappedGrid" begin + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + end + Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2] + + mg = equidistant_grid(c, 10,13) + + @test inverse_inner_product(mg, stencil_set) isa LazyTensor{<:Any, 2,2} + @test_broken false # Test that it calculates the right thing + end end
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl Mon Sep 02 15:35:54 2024 +0200 +++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl Tue Sep 03 08:38:27 2024 +0200 @@ -4,6 +4,8 @@ using Sbplib.Grids using Sbplib.LazyTensors +using StaticArrays + @testset "Laplace" begin # Default stencils (4th order) operator_path = sbp_operators_path()*"standard_diagonal.toml" @@ -72,12 +74,12 @@ g_1D = equidistant_grid(0.0, 1., 101) g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52) - @testset "1D" begin + @testset "EquidistantGrid" begin Δ = laplace(g_1D, stencil_set) @test Δ == second_derivative(g_1D, stencil_set) @test Δ isa LazyTensor{Float64,1,1} end - @testset "3D" begin + @testset "TensorGrid" begin Δ = laplace(g_3D, stencil_set) @test Δ isa LazyTensor{Float64,3,3} Dxx = second_derivative(g_3D, stencil_set, 1) @@ -86,6 +88,26 @@ @test Δ == Dxx + Dyy + Dzz @test Δ isa LazyTensor{Float64,3,3} end + + @testset "MappedGrid" begin + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ + η*(1-η), 3η+(1+η/2)*ξ^2] + end + Grids.jacobian(c::typeof(c), (ξ,η)) = @SMatrix[2 1-2η; (2+η)*ξ 3+ξ^2/2] + + g = equidistant_grid(c, 60,60) + + @test laplace(g, stencil_set) isa LazyTensor{<:Any,2,2} + + f((x,y)) = sin(4(x + y)) + Δf((x,y)) = -32sin(4(x + y)) + gf = map(f,g) + + Δ = laplace(g, stencil_set) + + @test collect(Δ*gf) isa Array{<:Any,2} + @test Δ*gf ≈ map(Δf, g) rtol=2e-2 + end end @testset "sat_tensors" begin