Mercurial > repos > public > sbplib_julia

changeset 797:92fafe5980dd operator_storage_array_of_table

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update inverse_inner_product to use the new storage
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 25 Jul 2021 14:58:50 +0200
parents f682e4fe3cef
children 997d6e641bf0
files src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl
diffstat 2 files changed, 44 insertions(+), 41 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sun Jul 25 14:39:41 2021 +0200
+++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl	Sun Jul 25 14:58:50 2021 +0200
@@ -1,37 +1,28 @@
 """
-    inverse_inner_product(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils)
-
-Creates the inverse inner product operator `H⁻¹` as a `TensorMapping` on an
-equidistant grid. `H⁻¹` is defined implicitly by `H⁻¹∘H = I`, where
-`H` is the corresponding inner product operator and `I` is the `IdentityMapping`.
+    inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
 
-`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.
+Constructs the inverse inner product operator `H⁻¹` as a `TensorMapping` using the weights of `H`, `interior_weight`, `closure_weights`. `H⁻¹` is inverse of the inner product operator `H`. The weights are the
 
-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
-operators in each coordinate direction. Also see the documentation of
-`SbpOperators.volume_operator(...)` for more details. On a 0-dimensional `grid`,
+On a 1-dimensional grid, `H⁻¹` is a `ConstantInteriorScalingOperator`. On an N-dimensional
+grid, `H⁻¹` is the outer product of the 1-dimensional inverse inner product
+operators for each coordinate direction. On a 0-dimensional `grid`,
 `H⁻¹` is a 0-dimensional `IdentityMapping`.
 """
-function inverse_inner_product(grid::EquidistantGrid, inv_closure_stencils, inv_inner_stencil)
+function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights)
     H⁻¹s = ()
 
     for i ∈ 1:dimension(grid)
-        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), inv_closure_stencils, inv_inner_stencil))
+        H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights))
     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)
+function inverse_inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights)
+    h⁻¹ = inverse_spacing(grid)[1]
+    H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights)
     return H⁻¹
 end
 export inverse_inner_product
 
-inverse_inner_product(grid::EquidistantGrid{0}, inv_closure_stencils, inv_inner_stencil) = IdentityMapping{eltype(grid)}()
-
-reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
+inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityMapping{eltype(grid)}()
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sun Jul 25 14:39:41 2021 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sun Jul 25 14:58:50 2021 +0200
@@ -12,34 +12,38 @@
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
     @testset "inverse_inner_product" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+        quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
         @testset "0D" begin
-            Hi = inverse_inner_product(EquidistantGrid{Float64}(),SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+            Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure)
             @test Hi == IdentityMapping{Float64}()
             @test Hi isa TensorMapping{T,0,0} where T
         end
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D, SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.));
+            Hi = inverse_inner_product(g_1D,  quadrature_interior, quadrature_closure)
             @test Hi isa TensorMapping{T,1,1} where T
         end
         @testset "2D" begin
-            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.))
+            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
+            Hi_x = inverse_inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure)
+            Hi_y = inverse_inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure)
             @test Hi == Hi_x⊗Hi_y
             @test Hi isa TensorMapping{T,2,2} where T
         end
     end
 
     @testset "Sizes" begin
-        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+        quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
         @testset "1D" begin
-            Hi = inverse_inner_product(g_1D,op.quadratureClosure, CenteredStencil(1.))
+            Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
             @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, CenteredStencil(1.))
+            Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
             @test domain_size(Hi) == size(g_2D)
             @test range_size(Hi) == size(g_2D)
         end
@@ -50,16 +54,20 @@
             v = evalOn(g_1D,x->sin(x))
             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, CenteredStencil(1.))
-                Hi = inverse_inner_product(g_1D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
                 @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,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_1D, quadrature_interior, quadrature_closure)
+                Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end
@@ -68,16 +76,20 @@
             v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
             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, CenteredStencil(1.))
-                Hi = inverse_inner_product(g_2D,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
+                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @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,SbpOperators.reciprocal_stencil.(op.quadratureClosure), CenteredStencil(1.))
+                stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+                quadrature_interior = parse_rational(stencil_set["H"]["inner"])
+                quadrature_closure = parse_rational.(stencil_set["H"]["closure"])
+                H = inner_product(g_2D, quadrature_interior, quadrature_closure)
+                Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
             end