Mercurial > repos > public > sbplib_julia
changeset 786:937e19326795 refactor/inner_product_recursive
Refactor inverse_inner_product
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Mon, 19 Jul 2021 22:05:09 +0200 |
parents | 21d6bc4a6d65 |
children | 961c36f5c60f |
files | src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl |
diffstat | 2 files changed, 18 insertions(+), 24 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Jul 19 21:16:14 2021 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon Jul 19 22:05:09 2021 +0200 @@ -22,22 +22,22 @@ `H⁻¹` is a 0-dimensional `IdentityMapping`. """ function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid)))) + H⁻¹s = () + + for i ∈ 1:dimension(grid) + H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), inv_closure_stencils, inv_inner_stencil)) + end + + return foldl(⊗, H⁻¹s) +end + +function inverse_inner_product(grid::EquidistantGrid{1}, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid)))) h⁻¹ = inverse_spacing(grid) - H⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[1]),scale.(inv_closure_stencils,h⁻¹[1]),even,1) - for i ∈ 2:dimension(grid) - Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i) - H⁻¹ = H⁻¹∘Hᵢ⁻¹ - end + H⁻¹ = SbpOperators.volume_operator(grid, scale(inv_inner_stencil, h⁻¹[1]), scale.(inv_closure_stencils, h⁻¹[1]),even,1) return H⁻¹ end export inverse_inner_product -inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil) = IdentityMapping{eltype(grid)}() - -function inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T} - inv_closure_stencils = reciprocal_stencil.(closure_stencils) - inv_inner_stencil = CenteredStencil(one(T)) - return inverse_inner_product(grid, inv_closure_stencils, inv_inner_stencil) -end +inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid)))) = IdentityMapping{eltype(grid)}() reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon Jul 19 21:16:14 2021 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon Jul 19 22:05:09 2021 +0200 @@ -14,18 +14,12 @@ @testset "inverse_inner_product" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) @testset "0D" begin - Hi = inverse_inner_product(EquidistantGrid{Float64}(),op.quadratureClosure) + Hi = inverse_inner_product(EquidistantGrid{Float64}(),SbpOperators.reciprocal_stencil.(op.quadratureClosure)) @test Hi == IdentityMapping{Float64}() @test Hi isa TensorMapping{T,0,0} where T end @testset "1D" begin - Hi = inverse_inner_product(g_1D, op.quadratureClosure); - inner_stencil = CenteredStencil(1.) - closures = () - for i = 1:length(op.quadratureClosure) - closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights)) - end - @test Hi == inverse_inner_product(g_1D,closures,inner_stencil) + Hi = inverse_inner_product(g_1D, SbpOperators.reciprocal_stencil.(op.quadratureClosure)); @test Hi isa TensorMapping{T,1,1} where T end @testset "2D" begin @@ -58,14 +52,14 @@ @testset "2nd order" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.)) - Hi = inverse_inner_product(g_1D,op.quadratureClosure) + Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(op.quadratureClosure)) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end @testset "4th order" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.)) - Hi = inverse_inner_product(g_1D,op.quadratureClosure) + Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(op.quadratureClosure)) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end @@ -76,14 +70,14 @@ @testset "2nd order" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.)) - Hi = inverse_inner_product(g_2D,op.quadratureClosure) + Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(op.quadratureClosure)) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end @testset "4th order" begin op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.)) - Hi = inverse_inner_product(g_2D,op.quadratureClosure) + Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(op.quadratureClosure)) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end