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