Mercurial > repos > public > sbplib_julia

changeset 1679:529b533a1dbb feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/sbp_operators/laplace_curvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 30 Jun 2024 10:57:05 +0200
parents de6300bd36cc (diff) 13a7a4ff49e3 (current diff)
children b30db2ea34ed
files Project.toml src/Grids/Grids.jl
diffstat 14 files changed, 255 insertions(+), 48 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Sun Jun 30 10:50:44 2024 +0200
+++ b/Manifest.toml	Sun Jun 30 10:57:05 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/Project.toml	Sun Jun 30 10:57:05 2024 +0200
@@ -4,6 +4,7 @@
 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"
--- a/src/Grids/Grids.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/Grids/Grids.jl	Sun Jun 30 10:57:05 2024 +0200
@@ -4,6 +4,31 @@
 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
 
 export ParameterSpace
 export HyperBox
--- a/src/SbpOperators/SbpOperators.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Sun Jun 30 10:57:05 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/boundaryops/boundary_restriction.jl	Sun Jun 30 10:57:05 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/boundaryops/normal_derivative.jl	Sun Jun 30 10:57:05 2024 +0200
@@ -28,3 +28,30 @@
     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⁻¹ = geometric_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
+
+    # 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Sun Jun 30 10:57:05 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, geometric_tensor(g))
+    DiagonalTensor(J)∘inner_product(logicalgrid(g), stencil_set)
+end
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sun Jun 30 10:57:05 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, geometric_tensor(g))
+    DiagonalTensor(J⁻¹)∘inner_product(logicalgrid(g), stencil_set)
+end
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 30 10:57:05 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, geometric_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/SbpOperators/boundaryops/boundary_restriction_test.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sun Jun 30 10:57:05 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sun Jun 30 10:57:05 2024 +0200
@@ -6,54 +6,99 @@
 using Sbplib.RegionIndices
 import Sbplib.SbpOperators.BoundaryOperator
 
+using StaticArrays
+
 @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)
+
+        b_w, b_e, b_s, b_n = boundary_identifiers(mg)
+
+        @test_broken normal_derivative(mg, stencil_set, b_w) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_e) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_s) isa LazyTensor{<:Any, 1, 2}
+        @test_broken normal_derivative(mg, stencil_set, b_n) isa LazyTensor{<:Any, 1, 2}
+
+        @testset "Accuracy" begin
+            v = map(x̄ -> NaN, mg)
+            dₙv = map(x̄ -> NaN, mg)
 
-            @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 "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(mg), Ref(stencil_set), boundary_identifiers(mg))
+
+                @test_broken d_w*v ≈ dₙv atol = 1e-13
+                @test_broken d_e*v ≈ dₙv atol = 1e-13
+                @test_broken d_s*v ≈ dₙv atol = 1e-13
+                @test_broken d_n*v ≈ dₙv 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
+            @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(mg), Ref(stencil_set), boundary_identifiers(mg))
+
+                @test_broken d_w*v ≈ dₙv atol = 1e-13
+                @test_broken d_e*v ≈ dₙv atol = 1e-13
+                @test_broken d_s*v ≈ dₙv atol = 1e-13
+                @test_broken d_n*v ≈ dₙv atol = 1e-13
+            end
         end
     end
 end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sun Jun 30 10:57:05 2024 +0200
@@ -6,6 +6,8 @@
 
 import Sbplib.SbpOperators.ConstantInteriorScalingOperator
 
+using StaticArrays
+
 @testset "Diagonal-stencil inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
@@ -94,4 +96,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 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/inner_products/inverse_inner_product_test.jl	Sun Jun 30 10:50:44 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sun Jun 30 10:57:05 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	Sun Jun 30 10:50:44 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sun Jun 30 10:57:05 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, 30,30)
+
+        @test laplace(g, stencil_set) isa LazyTensor{<:Any,2,2}
+
+        f((x,y)) = sin(4(x + y))
+        Δf((x,y)) = -16sin(4(x + y))
+        gf = map(f,g)
+
+        Δ = laplace(g, stencil_set)
+
+        @test collect(Δ*gf) isa Array{<:Any,2}
+        @test Δ*gf ≈ map(Δf, g)
+    end
 end
 
 @testset "sat_tensors" begin