Mercurial > repos > public > sbplib_julia
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