Mercurial > repos > public > sbplib_julia
changeset 642:f4a16b403487 feature/volume_and_boundary_operators
Implement the inverse quadrature operator as a volume operator and update tests.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Mon, 04 Jan 2021 17:17:40 +0100 |
| parents | 5e50e9815732 |
| children | 0928bbc3ee8b |
| files | src/SbpOperators/SbpOperators.jl src/SbpOperators/volumeops/quadratures/inverse_diagonal_quadrature.jl src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl test/testSbpOperators.jl |
| diffstat | 4 files changed, 121 insertions(+), 152 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Mon Jan 04 17:16:04 2021 +0100 +++ b/src/SbpOperators/SbpOperators.jl Mon Jan 04 17:17:40 2021 +0100 @@ -12,7 +12,7 @@ include("volumeops/derivatives/secondderivative.jl") include("volumeops/laplace/laplace.jl") include("volumeops/quadratures/quadrature.jl") -include("volumeops/quadratures/inverse_diagonal_quadrature.jl") +include("volumeops/quadratures/inverse_quadrature.jl") include("boundaryops/boundary_operator.jl") include("boundaryops/boundary_restriction.jl") include("boundaryops/normal_derivative.jl")
--- a/src/SbpOperators/volumeops/quadratures/inverse_diagonal_quadrature.jl Mon Jan 04 17:16:04 2021 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -""" -inverse_diagonal_quadrature(g,quadrature_closure) - -Constructs the inverse `Hi` of a `DiagonalQuadrature` on a grid of `Dim` dimensions as -a `TensorMapping`. The one-dimensional operator is a `InverseDiagonalQuadrature`, while -the multi-dimensional operator is the outer-product of the one-dimensional operators -in each coordinate direction. -""" -function inverse_diagonal_quadrature(g::EquidistantGrid{Dim}, quadrature_closure) where Dim - Hi = InverseDiagonalQuadrature(restrict(g,1), quadrature_closure) - for i ∈ 2:Dim - Hi = Hi⊗InverseDiagonalQuadrature(restrict(g,i), quadrature_closure) - end - return Hi -end -export inverse_diagonal_quadrature - - -""" - InverseDiagonalQuadrature{T,M} <: TensorMapping{T,1,1} - -Implements the inverse of a one-dimensional `DiagonalQuadrature` as a `TensorMapping` -The operator is defined by the reciprocal of the quadrature interval length `h_inv`, the -reciprocal of the quadrature closure weights `closure` and the number of quadrature intervals `size`. The -interior stencil has the weight 1. -""" -struct InverseDiagonalQuadrature{T<:Real,M} <: TensorMapping{T,1,1} - h_inv::T - closure::NTuple{M,T} - size::Tuple{Int} -end -export InverseDiagonalQuadrature - -""" - InverseDiagonalQuadrature(g, quadrature_closure) - -Constructs the `InverseDiagonalQuadrature` on the `EquidistantGrid` `g` with -closure given by the reciprocal of `quadrature_closure`. -""" -function InverseDiagonalQuadrature(g::EquidistantGrid{1}, quadrature_closure) - return InverseDiagonalQuadrature(inverse_spacing(g)[1], 1 ./ quadrature_closure, size(g)) -end - -""" - domain_size(Hi::InverseDiagonalQuadrature) - -The size of an object in the range of `Hi` -""" -LazyTensors.range_size(Hi::InverseDiagonalQuadrature) = Hi.size - -""" - domain_size(Hi::InverseDiagonalQuadrature) - -The size of an object in the domain of `Hi` -""" -LazyTensors.domain_size(Hi::InverseDiagonalQuadrature) = Hi.size - -""" - apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T -Implements the application `(Hi*v)[i]` an `Index{R}` where `R` is one of the regions -`Lower`,`Interior`,`Upper`. If `i` is another type of index (e.g an `Int`) it will first -be converted to an `Index{R}`. -""" -function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Lower}) where T - return @inbounds Hi.h_inv*Hi.closure[Int(i)]*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Upper}) where T - N = length(v); - return @inbounds Hi.h_inv*Hi.closure[N-Int(i)+1]*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i::Index{Interior}) where T - return @inbounds Hi.h_inv*v[Int(i)] -end - -function LazyTensors.apply(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T - N = length(v); - r = getregion(i, closure_size(Hi), N) - return LazyTensors.apply(Hi, v, Index(i, r)) -end - -LazyTensors.apply_transpose(Hi::InverseDiagonalQuadrature{T}, v::AbstractVector{T}, i) where T = LazyTensors.apply(Hi,v,i) - -""" - closure_size(Hi) -Returns the size of the closure stencil of a InverseDiagonalQuadrature `Hi`. -""" -closure_size(Hi::InverseDiagonalQuadrature{T,M}) where {T,M} = M
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/quadratures/inverse_quadrature.jl Mon Jan 04 17:17:40 2021 +0100 @@ -0,0 +1,41 @@ + +""" + InverseQuadrature(grid::EquidistantGrid, inv_inner_stencil, inv_closure_stencils) + +Creates the inverse `H⁻¹` of the quadrature operator as a `TensorMapping` + +The inverse quadrature approximates the integral operator on the grid using +`inv_inner_stencil` in the interior and a set of stencils `inv_closure_stencils` +for the points in the closure regions. + +On a one-dimensional `grid`, `H⁻¹` is a `VolumeOperator`. On a multi-dimensional +`grid`, `H` is the outer product of the 1-dimensional inverse quadrature operators in +each coordinate direction. Also see the documentation of +`SbpOperators.volume_operator(...)` for more details. +""" +function InverseQuadrature(grid::EquidistantGrid{Dim}, inv_inner_stencil, inv_closure_stencils) where Dim + 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:Dim + Hᵢ⁻¹ = SbpOperators.volume_operator(grid,scale(inv_inner_stencil,h⁻¹[i]),scale.(inv_closure_stencils,h⁻¹[i]),even,i) + H⁻¹ = H⁻¹∘Hᵢ⁻¹ + end + return H⁻¹ +end +export InverseQuadrature + +""" + InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils) + +Creates the inverse of the diagonal quadrature operator defined by the inner stencil +1/h and a set of 1-element closure stencils in `closure_stencils`. Note that +the closure stencils are those of the quadrature operator (and not the inverse). +""" +function InverseDiagonalQuadrature(grid::EquidistantGrid, closure_stencils::NTuple{M,Stencil{T,1}}) where {T,M} + inv_inner_stencil = Stencil(Tuple{T}(1),center=1) + inv_closure_stencils = reciprocal_stencil.(closure_stencils) + return InverseQuadrature(grid, inv_inner_stencil, inv_closure_stencils) +end +export InverseDiagonalQuadrature + +reciprocal_stencil(s::Stencil{T}) where T = Stencil(s.range,one(T)./s.weights)
--- a/test/testSbpOperators.jl Mon Jan 04 17:16:04 2021 +0100 +++ b/test/testSbpOperators.jl Mon Jan 04 17:17:40 2021 +0100 @@ -489,68 +489,85 @@ end end -# @testset "InverseDiagonalQuadrature" begin -# op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) -# Lx = π/2. -# Ly = Float64(π) -# g_1D = EquidistantGrid(77, 0.0, Lx) -# g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) -# @testset "Constructors" begin -# # 1D -# Hi_x = InverseDiagonalQuadrature(inverse_spacing(g_1D)[1], 1. ./ op.quadratureClosure, size(g_1D)); -# @test Hi_x == InverseDiagonalQuadrature(g_1D,op.quadratureClosure) -# @test Hi_x == inverse_diagonal_quadrature(g_1D,op.quadratureClosure) -# @test Hi_x isa TensorMapping{T,1,1} where T -# @test Hi_x' isa TensorMapping{T,1,1} where T -# -# # 2D -# Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) -# @test Hi_xy isa TensorMappingComposition -# @test Hi_xy isa TensorMapping{T,2,2} where T -# @test Hi_xy' isa TensorMapping{T,2,2} where T -# end -# -# @testset "Sizes" begin -# # 1D -# Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) -# @test domain_size(Hi_x) == size(g_1D) -# @test range_size(Hi_x) == size(g_1D) -# # 2D -# Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) -# @test domain_size(Hi_xy) == size(g_2D) -# @test range_size(Hi_xy) == size(g_2D) -# end -# -# @testset "Application" begin -# # 1D -# H_x = diagonal_quadrature(g_1D,op.quadratureClosure) -# Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) -# v_1D = evalOn(g_1D,x->sin(x)) -# u_1D = evalOn(g_1D,x->x^3-x^2+1) -# @test Hi_x*H_x*v_1D ≈ v_1D rtol = 1e-15 -# @test Hi_x*H_x*u_1D ≈ u_1D rtol = 1e-15 -# @test Hi_x*v_1D == Hi_x'*v_1D -# # 2D -# H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) -# Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) -# v_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y)) -# u_2D = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) -# @test Hi_xy*H_xy*v_2D ≈ v_2D rtol = 1e-15 -# @test Hi_xy*H_xy*u_2D ≈ u_2D rtol = 1e-15 -# @test Hi_xy*v_2D ≈ Hi_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason? -# end -# -# @testset "Inferred" begin -# Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) -# Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) -# v_1D = ones(Float64, size(g_1D)) -# v_2D = ones(Float64, size(g_2D)) -# @inferred Hi_x*v_1D -# @inferred Hi_x'*v_1D -# @inferred Hi_xy*v_2D -# @inferred Hi_xy'*v_2D -# end -# end +@testset "InverseDiagonalQuadrature" begin + Lx = π/2. + Ly = Float64(π) + g_1D = EquidistantGrid(77, 0.0, Lx) + g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + @testset "Constructors" begin + op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) + @testset "1D" begin + Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure); + inner_stencil = Stencil((1.,),center=1) + closures = () + for i = 1:length(op.quadratureClosure) + closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights)) + end + @test Hi == InverseQuadrature(g_1D,inner_stencil,closures) + @test Hi isa TensorMapping{T,1,1} where T + end + @testset "2D" begin + Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) + Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure) + Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure) + @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) + @testset "1D" begin + Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) + @test domain_size(Hi) == size(g_1D) + @test range_size(Hi) == size(g_1D) + end + @testset "2D" begin + Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) + @test domain_size(Hi) == size(g_2D) + @test range_size(Hi) == size(g_2D) + end + end + + @testset "Accuracy" begin + @testset "1D" begin + 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 = DiagonalQuadrature(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) + Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) + @test Hi*H*v ≈ v rtol = 1e-15 + @test Hi*H*u ≈ u rtol = 1e-15 + end + end + @testset "2D" begin + 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 = DiagonalQuadrature(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) + Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) + @test Hi*H*v ≈ v rtol = 1e-15 + @test Hi*H*u ≈ u rtol = 1e-15 + end + end + end +end @testset "BoundaryOperator" begin closure_stencil = Stencil((0,2), (2.,1.,3.))