Mercurial > repos > public > sbplib_julia

changeset 817:a1d556611e3c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge refactor/inner_product_recursive
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 13 Jan 2022 12:29:11 +0100
parents 181c3c93463d (current diff) f64c0bf9c1eb (diff)
children c7af0d04efee 3c1dd7692797
files
diffstat 4 files changed, 54 insertions(+), 62 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl	Thu Jan 13 12:29:11 2022 +0100
@@ -6,8 +6,7 @@
 
 `inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)` creates
 `H` on `grid` the using a set of stencils `closure_stencils` for the points in
-the closure regions and the stencil and `inner_stencil` in the interior. If
-`inner_stencil` is omitted a central interior stencil with weight 1 is used.
+the closure regions and the stencil and `inner_stencil` in the interior.
 
 On a 1-dimensional `grid`, `H` is a `VolumeOperator`. On a N-dimensional
 `grid`, `H` is the outer product of the 1-dimensional inner product operators in
@@ -15,15 +14,21 @@
 `SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
 `H` is a 0-dimensional `IdentityMapping`.
 """
-function inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil = CenteredStencil(one(eltype(grid))))
-    h = spacing(grid)
-    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:dimension(grid)
-        Hᵢ = SbpOperators.volume_operator(grid, scale(inner_stencil,h[i]), scale.(closure_stencils,h[i]), even, i)
-        H = H∘Hᵢ
+function inner_product(grid::EquidistantGrid, closure_stencils, inner_stencil)
+    Hs = ()
+
+    for i ∈ 1:dimension(grid)
+        Hs = (Hs..., inner_product(restrict(grid, i), closure_stencils, inner_stencil))
     end
-    return H
+
+    return foldl(⊗, Hs)
 end
 export inner_product
 
+function inner_product(grid::EquidistantGrid{1}, closure_stencils, inner_stencil)
+    h = spacing(grid)
+    H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
+    return H
+end
+
 inner_product(grid::EquidistantGrid{0}, closure_stencils, inner_stencil) = IdentityMapping{eltype(grid)}()
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Thu Jan 13 12:29:11 2022 +0100
@@ -1,6 +1,5 @@
 """
     inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-    inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})
 
 Creates the inverse inner product operator `H⁻¹` as a `TensorMapping` on an
 equidistant grid. `H⁻¹` is defined implicitly by `H⁻¹∘H = I`, where
@@ -8,12 +7,7 @@
 
 `inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)`
 constructs `H⁻¹` using a set of stencils `inv_closure_stencils` for the points
-in the closure regions and the stencil `inv_inner_stencil` in the interior. If
-`inv_closure_stencils` is omitted, a central interior stencil with weight 1 is used.
-
-`inverse_inner_product(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}})`
-constructs a diagonal inverse inner product operator where `closure_stencils` are the
-closure stencils of `H` (not `H⁻¹`!).
+in the closure regions and the stencil `inv_inner_stencil` in the interior.
 
 On a 1-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a N-dimensional
 `grid`, `H⁻¹` is the outer product of the 1-dimensional inverse inner product
@@ -21,23 +15,23 @@
 `SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
 `H⁻¹` is a 0-dimensional `IdentityMapping`.
 """
-function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil = CenteredStencil(one(eltype(grid))))
+function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil)
+    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)
     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
-
 reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Thu Jan 13 12:29:11 2022 +0100
@@ -16,20 +16,19 @@
     @testset "inner_product" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "0D" begin
-            H = inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
+            H = inner_product(EquidistantGrid{Float64}(), op.quadratureClosure, CenteredStencil(1.))
             @test H == IdentityMapping{Float64}()
             @test H isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            H = inner_product(g_1D,op.quadratureClosure)
-            inner_stencil = CenteredStencil(1.)
-            @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
+            H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+            @test H == inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
             @test H isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            H = inner_product(g_2D,op.quadratureClosure)
-            H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
-            H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
+            H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+            H_x = inner_product(restrict( g_2D,1),op.quadratureClosure, CenteredStencil(1.))
+            H_y = inner_product(restrict( g_2D,2),op.quadratureClosure, CenteredStencil(1.))
             @test H == H_x⊗H_y
             @test H isa TensorMapping{T,2,2} where T
         end
@@ -38,12 +37,12 @@
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = inner_product(g_1D,op.quadratureClosure)
+            H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = inner_product(g_2D,op.quadratureClosure)
+            H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -60,7 +59,7 @@
 
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -69,7 +68,7 @@
 
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -83,13 +82,13 @@
             u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-4
             end
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = inner_product(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Dec 21 15:34:51 2021 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Thu Jan 13 12:29:11 2022 +0100
@@ -14,24 +14,18 @@
     @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), CenteredStencil(1.))
             @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), CenteredStencil(1.));
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
-            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
-            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure, CenteredStencil(1.))
+            Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure, CenteredStencil(1.))
+            Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure, CenteredStencil(1.))
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
@@ -40,12 +34,12 @@
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_1D,op.quadratureClosure, CenteredStencil(1.))
             @test domain_size(Hi) == size(g_1D)
             @test range_size(Hi) == size(g_1D)
         end
         @testset "2D" begin
-            Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+            Hi = inverse_inner_product(g_2D,op.quadratureClosure, CenteredStencil(1.))
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -57,15 +51,15 @@
             u = evalOn(g_1D,x->x^3-x^2+1)
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_1D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+                Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
                 @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)
-                Hi = inverse_inner_product(g_1D,op.quadratureClosure)
+                H = inner_product(g_1D, op.quadratureClosure, CenteredStencil(1.))
+                Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -75,15 +69,15 @@
             u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = inner_product(g_2D,op.quadratureClosure)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+                Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
                 @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)
-                Hi = inverse_inner_product(g_2D,op.quadratureClosure)
+                H = inner_product(g_2D, op.quadratureClosure, CenteredStencil(1.))
+                Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end