Mercurial > repos > public > sbplib_julia
changeset 1785:96f8cad255b4 feature/sbp_operators/laplace_curvilinear
Merge feature/grids/manifolds
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 16 Sep 2024 10:33:47 +0200 |
| parents | f3d7e2d7a43f (diff) c5070edd0ebb (current diff) |
| children | 1f42944d4a72 |
| files | Manifest.toml Project.toml src/Grids/Grids.jl |
| diffstat | 16 files changed, 402 insertions(+), 46 deletions(-) [+] |
line wrap: on
line diff
diff -r c5070edd0ebb -r 96f8cad255b4 ext/SbplibPlotsExt.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/SbplibPlotsExt.jl Mon Sep 16 10:33:47 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 + +
diff -r c5070edd0ebb -r 96f8cad255b4 ext/SbplibSparseArrayKitExt.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/SbplibSparseArrayKitExt.jl Mon Sep 16 10:33:47 2024 +0200 @@ -0,0 +1,14 @@ +module SbplibSparseArrayKitExt + +using Sbplib +using Sbplib.LazyTensors + +using SparseArrayKit +using Tokens + +function SparseArrayKit.SparseArray(t::LazyTensor) + v = ArrayToken(:v, domain_size(t)...) + return Tokens._to_tensor(t*v, range_size(t), domain_size(t)) +end + +end
diff -r c5070edd0ebb -r 96f8cad255b4 ext/SbplibSparseArraysExt.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/SbplibSparseArraysExt.jl Mon Sep 16 10:33:47 2024 +0200 @@ -0,0 +1,17 @@ +module SbplibSparseArraysExt + +using Sbplib +using Sbplib.LazyTensors + +using SparseArrays +using Tokens + +function SparseArrays.sparse(t::LazyTensor) + v = ArrayToken(:v, prod(domain_size(t))) + + v̄ = reshape(v,domain_size(t)...) + tv = reshape(t*v̄, :) + return Tokens._to_matrix(tv, prod(range_size(t)), prod(domain_size(t))) +end + +end
diff -r c5070edd0ebb -r 96f8cad255b4 src/Grids/Grids.jl --- a/src/Grids/Grids.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/Grids/Grids.jl Mon Sep 16 10:33:47 2024 +0200 @@ -70,6 +70,8 @@ export metric_tensor export metric_tensor_inverse +abstract type BoundaryIdentifier end + include("manifolds.jl") include("grid.jl") include("tensor_grid.jl")
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/SbpOperators.jl --- a/src/SbpOperators/SbpOperators.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/SbpOperators.jl Mon Sep 16 10:33:47 2024 +0200 @@ -43,6 +43,7 @@ using Diffinitive.RegionIndices using Diffinitive.LazyTensors using Diffinitive.Grids +using LinearAlgebra # Includes include("stencil.jl")
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/boundaryops/boundary_restriction.jl --- a/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Mon Sep 16 10:33:47 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
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/boundaryops/normal_derivative.jl --- a/src/SbpOperators/boundaryops/normal_derivative.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Mon Sep 16 10:33:47 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
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl --- a/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Mon Sep 16 10:33:47 2024 +0200 @@ -140,6 +140,7 @@ ## x-direction function apply_lower(op::SecondDerivativeVariable{1}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[:,j] c̃ = @view op.coefficient[:,j] @@ -147,6 +148,7 @@ end function apply_interior(op::SecondDerivativeVariable{1}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[:,j] c̃ = @view op.coefficient[:,j] @@ -154,6 +156,7 @@ end function apply_upper(op::SecondDerivativeVariable{1}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[:,j] c̃ = @view op.coefficient[:,j] @@ -165,6 +168,7 @@ ## y-direction function apply_lower(op::SecondDerivativeVariable{2}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[i,:] c̃ = @view op.coefficient[i,:] @@ -172,6 +176,7 @@ end function apply_interior(op::SecondDerivativeVariable{2}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[i,:] c̃ = @view op.coefficient[i,:] @@ -179,6 +184,7 @@ end function apply_upper(op::SecondDerivativeVariable{2}, v, i, j) + Base.@constprop :aggressive ṽ = @view v[i,:] c̃ = @view op.coefficient[i,:]
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/volumeops/inner_products/inner_product.jl --- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Mon Sep 16 10:33:47 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
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl --- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Sep 16 10:33:47 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
diff -r c5070edd0ebb -r 96f8cad255b4 src/SbpOperators/volumeops/laplace/laplace.jl --- a/src/SbpOperators/volumeops/laplace/laplace.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Mon Sep 16 10:33:47 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)
diff -r c5070edd0ebb -r 96f8cad255b4 test/SbpOperators/boundaryops/boundary_restriction_test.jl --- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl Mon Sep 16 10:33:47 2024 +0200 @@ -6,6 +6,8 @@ using Diffinitive.RegionIndices using Diffinitive.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
diff -r c5070edd0ebb -r 96f8cad255b4 test/SbpOperators/boundaryops/normal_derivative_test.jl --- a/test/SbpOperators/boundaryops/normal_derivative_test.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl Mon Sep 16 10:33:47 2024 +0200 @@ -6,54 +6,130 @@ using Diffinitive.RegionIndices import Diffinitive.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, LowerBoundary()) - @test d_l == normal_derivative(g_1D, stencil_set, LowerBoundary()) - @test d_l isa BoundaryOperator{T,LowerBoundary} 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,LowerBoundary}()) - d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,UpperBoundary}()) - Ix = IdentityTensor{Float64}((size(g_2D)[1],)) - Iy = IdentityTensor{Float64}((size(g_2D)[2],)) - d_l = normal_derivative(g_2D.grids[1], stencil_set, LowerBoundary()) - d_r = normal_derivative(g_2D.grids[2], stencil_set, UpperBoundary()) - @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}()) - @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, LowerBoundary()) + @test d_l == normal_derivative(g_1D, stencil_set, LowerBoundary()) + @test d_l isa BoundaryOperator{T,LowerBoundary} 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,LowerBoundary}()) + d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,UpperBoundary}()) + Ix = IdentityTensor{Float64}((size(g_2D)[1],)) + Iy = IdentityTensor{Float64}((size(g_2D)[2],)) + d_l = normal_derivative(g_2D.grids[1], stencil_set, LowerBoundary()) + d_r = normal_derivative(g_2D.grids[2], stencil_set, UpperBoundary()) + @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}()) + @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
diff -r c5070edd0ebb -r 96f8cad255b4 test/SbpOperators/volumeops/inner_products/inner_product_test.jl --- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Mon Sep 16 10:33:47 2024 +0200 @@ -6,6 +6,9 @@ import Diffinitive.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
diff -r c5070edd0ebb -r 96f8cad255b4 test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl --- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon Sep 16 10:33:47 2024 +0200 @@ -6,6 +6,8 @@ import Diffinitive.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
diff -r c5070edd0ebb -r 96f8cad255b4 test/SbpOperators/volumeops/laplace/laplace_test.jl --- a/test/SbpOperators/volumeops/laplace/laplace_test.jl Mon Sep 16 09:53:29 2024 +0200 +++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl Mon Sep 16 10:33:47 2024 +0200 @@ -4,6 +4,8 @@ using Diffinitive.Grids using Diffinitive.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