changeset 1655:51f23a0b07fa feature/sbp_operators/laplace_curvilinear

Add inner product and inverse inner product for mapped grids
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2024 15:15:32 +0200
parents f4dc17cfafce
children 89456aa6fa80
files src/SbpOperators/volumeops/inner_products/inner_product.jl src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl test/SbpOperators/volumeops/inner_products/inner_product_test.jl test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl
diffstat 4 files changed, 40 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Wed Jun 26 14:41:50 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Wed Jun 26 15:15:32 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 = jacobian_determinant(g)
+    DiagonalTensor(J)∘inner_product(logicalgrid(g), stencil_set)
+end
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Wed Jun 26 14:41:50 2024 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Wed Jun 26 15:15:32 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, jacobian_determinant(g))
+    DiagonalTensor(J⁻¹)∘inner_product(logicalgrid(g), stencil_set)
+end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Jun 26 14:41:50 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Wed Jun 26 15:15:32 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	Wed Jun 26 14:41:50 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Wed Jun 26 15:15:32 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