Mercurial > repos > public > sbplib_julia

changeset 676:bf48761c1345 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/boundary_quads
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 06 Feb 2021 14:03:14 +0100
parents 2d56a53a1646 (current diff) 538ccbaeb1f8 (diff)
children 011863b3f24c
files
diffstat 2 files changed, 103 insertions(+), 35 deletions(-) [+]
line wrap: on
line diff
--- a/src/SbpOperators/volumeops/quadratures/quadrature.jl	Sun Jan 31 22:19:53 2021 +0100
+++ b/src/SbpOperators/volumeops/quadratures/quadrature.jl	Sat Feb 06 14:03:14 2021 +0100
@@ -1,36 +1,69 @@
 """
-    Quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils)
+    quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils)
+    quadrature(grid::EquidistantGrid, closure_stencils)
 
 Creates the quadrature operator `H` as a `TensorMapping`
 
-The quadrature approximates the integral operator on the grid using
+`H` approximiates the integral operator on `grid` the using the stencil
 `inner_stencil` in the interior and a set of stencils `closure_stencils`
-for the points in the closure regions.
+for the points in the closure regions. If `inner_stencil` is omitted a central
+interior stencil with weight 1 is used.
 
 On a one-dimensional `grid`, `H` is a `VolumeOperator`. On a multi-dimensional
 `grid`, `H` is the outer product of the 1-dimensional quadrature operators in
 each coordinate direction. Also see the documentation of
 `SbpOperators.volume_operator(...)` for more details.
 """
-function Quadrature(grid::EquidistantGrid{Dim}, inner_stencil, closure_stencils) where Dim
+function quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils) where Dim
     h = spacing(grid)
     H = SbpOperators.volume_operator(grid, scale(inner_stencil,h[1]), scale.(closure_stencils,h[1]), even, 1)
-    for i ∈ 2:Dim
+    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ᵢ
     end
     return H
 end
-export Quadrature
+export quadrature
+
+function quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}) where {M,T}
+    inner_stencil = Stencil(Tuple{T}(1),center=1)
+    return quadrature(grid, inner_stencil, closure_stencils)
+end
 
 """
-    DiagonalQuadrature(grid::EquidistantGrid, closure_stencils)
+    boundary_quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils, id::CartesianBoundary)
+    boundary_quadrature(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, id)
+    boundary_quadrature(grid::EquidistantGrid, closure_stencils, id)
 
-Creates the quadrature operator with the inner stencil 1/h and 1-element sized
-closure stencils (i.e the operator is diagonal)
+Creates the lower-dimensional quadrature operator associated with the boundary
+of `grid` specified by `id`. The quadrature operator is defined on the grid
+spanned by the dimensions orthogonal to the boundary coordinate direction.
+If the dimension of `grid` is 1, then the boundary quadrature is the 0-dimensional
+`IdentityMapping`. If `inner_stencil` is omitted a central interior stencil with
+weight 1 is used.
 """
-function DiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {M,T}
+function boundary_quadrature(grid::EquidistantGrid, inner_stencil, closure_stencils, id::CartesianBoundary)
+    return quadrature(orthogonal_grid(grid,dim(id)),inner_stencil,closure_stencils)
+end
+export boundary_quadrature
+
+function boundary_quadrature(grid::EquidistantGrid{1}, inner_stencil::Stencil{T}, closure_stencils::NTuple{M,Stencil{T}}, id::CartesianBoundary{1}) where {M,T}
+    return IdentityMapping{T}()
+end
+
+function boundary_quadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T}}, id::CartesianBoundary) where {M,T}
     inner_stencil = Stencil(Tuple{T}(1),center=1)
-    return Quadrature(grid, inner_stencil, closure_stencils)
+    return boundary_quadrature(grid,inner_stencil,closure_stencils,id)
 end
-export DiagonalQuadrature
+
+"""
+    orthogonal_grid(grid,dim)
+
+Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
+orthogonal to `dim`.
+"""
+function orthogonal_grid(grid,dim)
+    dims = collect(1:dimension(grid))
+    orth_dims = dims[dims .!= dim]
+    return restrict(grid,orth_dims)
+end
--- a/test/testSbpOperators.jl	Sun Jan 31 22:19:53 2021 +0100
+++ b/test/testSbpOperators.jl	Sat Feb 06 14:03:14 2021 +0100
@@ -290,7 +290,7 @@
                 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
             end
         end
-        
+
         @testset "2D" begin
             l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
             binomials = ()
@@ -390,38 +390,73 @@
     end
 end
 
-@testset "DiagonalQuadrature" begin
+@testset "Quadrature diagonal" begin
     Lx = π/2.
     Ly = Float64(π)
+    Lz = 1.
     g_1D = EquidistantGrid(77, 0.0, Lx)
     g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
     integral(H,v) = sum(H*v)
-    @testset "Constructors" begin
+    @testset "quadrature" begin
+        op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
+        @testset "1D" begin
+            H = quadrature(g_1D,op.quadratureClosure)
+            inner_stencil = Stencil((1.,),center=1)
+            @test H == quadrature(g_1D,inner_stencil,op.quadratureClosure)
+            @test H isa TensorMapping{T,1,1} where T
+        end
+        @testset "2D" begin
+            H = quadrature(g_2D,op.quadratureClosure)
+            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
+            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
+            @test H == H_x⊗H_y
+            @test H isa TensorMapping{T,2,2} where T
+        end
+    end
+
+    @testset "boundary_quadrature" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
-            inner_stencil = Stencil((1.,),center=1)
-            @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure)
-            @test H isa TensorMapping{T,1,1} where T
+            (id_l, id_r) = boundary_identifiers(g_1D)
+            @test boundary_quadrature(g_1D,op.quadratureClosure,id_l) == IdentityMapping{Float64}()
+            @test boundary_quadrature(g_1D,op.quadratureClosure,id_r) == IdentityMapping{Float64}()
+
+        end
+        @testset "2D" begin
+            (id_w, id_e, id_s, id_n) = boundary_identifiers(g_2D)
+            H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
+            H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
+            @test boundary_quadrature(g_2D,op.quadratureClosure,id_w) == H_y
+            @test boundary_quadrature(g_2D,op.quadratureClosure,id_e) == H_y
+            @test boundary_quadrature(g_2D,op.quadratureClosure,id_s) == H_x
+            @test boundary_quadrature(g_2D,op.quadratureClosure,id_n) == H_x
         end
-        @testset "1D" begin
-            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
-            H_x = DiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
-            H_y = DiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
-            @test H == H_x⊗H_y
-            @test H isa TensorMapping{T,2,2} where T
+        @testset "3D" begin
+            (id_w, id_e,
+             id_s, id_n,
+             id_t, id_b) = boundary_identifiers(g_3D)
+            H_xy = quadrature(restrict(g_3D,[1,2]),op.quadratureClosure)
+            H_xz = quadrature(restrict(g_3D,[1,3]),op.quadratureClosure)
+            H_yz = quadrature(restrict(g_3D,[2,3]),op.quadratureClosure)
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_w) == H_yz
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_e) == H_yz
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_s) == H_xz
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_n) == H_xz
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_t) == H_xy
+            @test boundary_quadrature(g_3D,op.quadratureClosure,id_b) == H_xy
         end
     end
 
     @testset "Sizes" begin
         op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
-            H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+            H = quadrature(g_1D,op.quadratureClosure)
             @test domain_size(H) == size(g_1D)
             @test range_size(H) == size(g_1D)
         end
         @testset "2D" begin
-            H = DiagonalQuadrature(g_2D,op.quadratureClosure)
+            H = quadrature(g_2D,op.quadratureClosure)
             @test domain_size(H) == size(g_2D)
             @test range_size(H) == size(g_2D)
         end
@@ -438,7 +473,7 @@
 
             @testset "2nd order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 for i = 1:2
                     @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
                 end
@@ -447,7 +482,7 @@
 
             @testset "4th order" begin
                 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-                H = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 for i = 1:4
                     @test integral(H,v[i]) ≈ v[i+1][end] -  v[i+1][1] rtol = 1e-14
                 end
@@ -461,13 +496,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 = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 @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 = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
                 @test integral(H,u) ≈ π rtol = 1e-8
             end
@@ -521,14 +556,14 @@
             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 = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_1D,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 = DiagonalQuadrature(g_1D,op.quadratureClosure)
+                H = quadrature(g_1D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15
@@ -539,14 +574,14 @@
             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 = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_2D,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 = DiagonalQuadrature(g_2D,op.quadratureClosure)
+                H = quadrature(g_2D,op.quadratureClosure)
                 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
                 @test Hi*H*v ≈ v rtol = 1e-15
                 @test Hi*H*u ≈ u rtol = 1e-15