Mercurial > repos > public > sbplib_julia
changeset 1888:eb70a0941cd6 allocation_testing
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 03 Feb 2023 23:02:46 +0100 |
| parents | 24590890e124 (current diff) d13b73e9f65e (diff) |
| children | c1519a9b4903 |
| files | test/Manifest.toml |
| diffstat | 33 files changed, 888 insertions(+), 682 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Thu Apr 07 17:02:51 2022 +0200 +++ b/Manifest.toml Fri Feb 03 23:02:46 2023 +0100 @@ -1,25 +1,50 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.7.1" +julia_version = "1.8.5" manifest_format = "2.0" +project_hash = "b024d6898b484706c36ee3b2a041918f3a9d2088" [[deps.Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f" +git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.3.3" +version = "3.5.0" + +[[deps.ArrayInterface]] +deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"] +git-tree-sha1 = "dedc16cbdd1d32bead4617d27572f582216ccf23" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "6.0.25" + +[[deps.ArrayInterfaceCore]] +deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d" +uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2" +version = "0.1.29" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" +[[deps.Compat]] +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.0" + [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.1+0" [[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" +[[deps.IfElse]] +git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" +uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +version = "0.1.1" + [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -29,27 +54,70 @@ [[deps.OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed" +git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.10.8" +version = "1.12.9" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" [[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" + +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.Static]] +deps = ["IfElse"] +git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.8.3" + +[[deps.SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" [[deps.TiledIteration]] -deps = ["OffsetArrays"] -git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297" +deps = ["ArrayInterface", "OffsetArrays"] +git-tree-sha1 = "1bf2bb587a7fc99fefac2ff076b18b500128e9c0" uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" -version = "0.3.1" +version = "0.4.2" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" @@ -57,3 +125,4 @@ [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.1.1+0"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/WORKFLOW.md Fri Feb 03 23:02:46 2023 +0100 @@ -0,0 +1,46 @@ +# Branch model + +The `default` branch contains stable code and the goal is that all tests should be passing on this branch. Version of the code are marked with mercurial tags. Changes in the code are developed in named branches which are closed and merged into `default` when the code is ready. During development merging `default` into the development branch is encouraged to avoid complicated conflicts. + +Branches are named using slash separated keywords. The first keyword describes the type of change being pursued in the branch. Important type keywords are + * feature + * bugfix + * refactor +Further keywords may describe where, e.g. what sub package, the change happens. The last keyword should describe the change. + +Some examples: + * refactor/grids: Branch to refactor the grids module + * bugfix/lazy_tensors/lazyfunctionarray: Branch to fix a bug in LazyFunctionArray + +## Merging a branch into `default` +The changes in a branch has been reviewed and deemed ready to merge the branch is closed and then merged. + +Before merging a development branch, `default` should be merge into the development branch to make sure the whole state of the code is reviewed and tested before it ends up on `default`. + +With the development branch active the following commands can be used to complete the merging of a development branch. +```shell +hg merge default +hg commit --close-branch -m "Close before merge" +hg update default +hg merge development/branch +hg commit -m "Merge development/branch" +``` + +# Review + +## Checklist for review +* Push and pull new changes +* Search and check TODOs +* Search and check TBDs +* Search and check REVIEWs +* Review code +* Review tests +* Review docstrings +* Render Documenter and check docstrings in browser +* Run full tests + +# Special comments +The following special comments are used: +* `# TODO: `: Something that should be done at some point. +* `# TBD: `: "To be determined", i.e a decision that has to be made. +* `# REVIEW: `: A review comment. Should only exist on development branches.
--- a/docs/Manifest.toml Thu Apr 07 17:02:51 2022 +0200 +++ b/docs/Manifest.toml Fri Feb 03 23:02:46 2023 +0100 @@ -1,7 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.7.1" +julia_version = "1.8.5" manifest_format = "2.0" +project_hash = "4f0756199bb5f6739a5f4697152617efc4e0705c" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -10,9 +11,21 @@ [[deps.Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f" +git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.3.3" +version = "3.5.0" + +[[deps.ArrayInterface]] +deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"] +git-tree-sha1 = "dedc16cbdd1d32bead4617d27572f582216ccf23" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "6.0.25" + +[[deps.ArrayInterfaceCore]] +deps = ["LinearAlgebra", "SnoopPrecompile", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "e5f08b5689b1aad068e01751889f2f615c7db36d" +uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2" +version = "0.1.29" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -20,9 +33,16 @@ [[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" +[[deps.Compat]] +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.0" + [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.1+0" [[deps.Dates]] deps = ["Printf"] @@ -30,15 +50,15 @@ [[deps.DocStringExtensions]] deps = ["LibGit2"] -git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b" +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.6" +version = "0.9.3" [[deps.Documenter]] deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "f425293f7e0acaf9144de6d731772de156676233" +git-tree-sha1 = "58fea7c536acd71f3eef6be3b21c0df5f3df88fd" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.27.10" +version = "0.27.24" [[deps.IOCapture]] deps = ["Logging", "Random"] @@ -46,15 +66,20 @@ uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" version = "0.2.2" +[[deps.IfElse]] +git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" +uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +version = "0.1.1" + [[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37" +git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.2" +version = "0.21.3" [[deps.LibGit2]] deps = ["Base64", "NetworkOptions", "Printf", "SHA"] @@ -79,22 +104,30 @@ [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" [[deps.OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed" +git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.10.8" +version = "1.12.9" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" [[deps.Parsers]] -deps = ["Dates"] -git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918" +deps = ["Dates", "SnoopPrecompile"] +git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.1.3" +version = "2.5.4" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" [[deps.Printf]] deps = ["Unicode"] @@ -110,6 +143,7 @@ [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" [[deps.Sbplib]] deps = ["TOML", "TiledIteration"] @@ -120,22 +154,47 @@ [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" + [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.Static]] +deps = ["IfElse"] +git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.8.3" + +[[deps.SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" [[deps.Test]] deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.TiledIteration]] -deps = ["OffsetArrays"] -git-tree-sha1 = "5683455224ba92ef59db72d10690690f4a8dc297" +deps = ["ArrayInterface", "OffsetArrays"] +git-tree-sha1 = "1bf2bb587a7fc99fefac2ff076b18b500128e9c0" uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" -version = "0.3.1" +version = "0.4.2" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" @@ -143,3 +202,4 @@ [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.1.1+0"
--- a/src/Grids/AbstractGrid.jl Thu Apr 07 17:02:51 2022 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -""" - AbstractGrid - -Should implement - dimension(grid::AbstractGrid) - points(grid::AbstractGrid) - -""" -abstract type AbstractGrid end -export AbstractGrid -function dimension end -function points end -export dimension, points - -""" - evalOn(g::AbstractGrid, f::Function) - -Evaluate function f on the grid g -""" -function evalOn(g::AbstractGrid, f::Function) - F(x) = f(x...) - return F.(points(g)) -end -export evalOn
--- a/src/Grids/EquidistantGrid.jl Thu Apr 07 17:02:51 2022 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,183 +0,0 @@ -export EquidistantGrid -export spacing -export inverse_spacing -export restrict -export boundary_identifiers -export boundary_grid -export refine -export coarsen - -""" - EquidistantGrid{Dim,T<:Real} <: AbstractGrid - -`Dim`-dimensional equidistant grid with coordinates of type `T`. -""" -struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid - size::NTuple{Dim, Int} - limit_lower::NTuple{Dim, T} - limit_upper::NTuple{Dim, T} - - function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T} - if any(size .<= 0) - throw(DomainError("all components of size must be postive")) - end - if any(limit_upper.-limit_lower .<= 0) - throw(DomainError("all side lengths must be postive")) - end - return new{Dim,T}(size, limit_lower, limit_upper) - end -end - -""" - EquidistantGrid(size, limit_lower, limit_upper) - -Construct an equidistant grid with corners at the coordinates `limit_lower` and -`limit_upper`. - -The length of the domain sides are given by the components of -`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined -as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative. - -The number of equidistantly spaced points in each coordinate direction are given -by the tuple `size`. -""" -function EquidistantGrid(size, limit_lower, limit_upper) - return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper) -end - -""" - EquidistantGrid{T}() - -Constructs a 0-dimensional grid. -""" -EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid - -""" - EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) - -Convenience constructor for 1D grids. -""" -function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T - return EquidistantGrid((size,),(limit_lower,),(limit_upper,)) -end - -Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T - -Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size) - -Base.size(g::EquidistantGrid) = g.size - -""" - dimension(grid::EquidistantGrid) - -The dimension of the grid. -""" -dimension(grid::EquidistantGrid{Dim}) where Dim = Dim - -""" - spacing(grid::EquidistantGrid) - -The spacing between grid points. -""" -spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1) - -""" - inverse_spacing(grid::EquidistantGrid) - -The reciprocal of the spacing between grid points. -""" -inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid) - -""" - points(grid::EquidistantGrid) - -The point of the grid as an array of tuples with the same dimension as the grid. -The points are stored as [(x1,y1), (x1,y2), … (x1,yn); - (x2,y1), (x2,y2), … (x2,yn); - ⋮ ⋮ ⋮ - (xm,y1), (xm,y2), … (xm,yn)] -""" -function points(grid::EquidistantGrid) - indices = Tuple.(CartesianIndices(grid.size)) - h = spacing(grid) - return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) -end - -""" - restrict(::EquidistantGrid, dim) - -Pick out given dimensions from the grid and return a grid for them. -""" -function restrict(grid::EquidistantGrid, dim) - size = grid.size[dim] - limit_lower = grid.limit_lower[dim] - limit_upper = grid.limit_upper[dim] - - return EquidistantGrid(size, limit_lower, limit_upper) -end - -""" - boundary_identifiers(::EquidistantGrid) - -Returns a tuple containing the boundary identifiers for the grid, stored as - (CartesianBoundary(1,Lower), - CartesianBoundary(1,Upper), - CartesianBoundary(2,Lower), - ...) -""" -boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),dimension(g)))...)...,) - - -""" - boundary_grid(grid::EquidistantGrid,id::CartesianBoundary) - boundary_grid(::EquidistantGrid{1},::CartesianBoundary{1}) - -Creates the lower-dimensional restriciton of `grid` spanned by the dimensions -orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional -grid is a zero-dimensional grid. -""" -function boundary_grid(grid::EquidistantGrid,id::CartesianBoundary) - dims = collect(1:dimension(grid)) - orth_dims = dims[dims .!= dim(id)] - if orth_dims == dims - throw(DomainError("boundary identifier not matching grid")) - end - return restrict(grid,orth_dims) -end -boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}() - - -""" - refine(grid::EquidistantGrid, r::Int) - -Refines `grid` by a factor `r`. The factor is applied to the number of -intervals which is 1 less than the size of the grid. - -See also: [`coarsen`](@ref) -""" -function refine(grid::EquidistantGrid, r::Int) - sz = size(grid) - new_sz = (sz .- 1).*r .+ 1 - return EquidistantGrid{dimension(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) -end - -""" - coarsen(grid::EquidistantGrid, r::Int) - -Coarsens `grid` by a factor `r`. The factor is applied to the number of -intervals which is 1 less than the size of the grid. If the number of -intervals are not divisible by `r` an error is raised. - -See also: [`refine`](@ref) -""" -function coarsen(grid::EquidistantGrid, r::Int) - sz = size(grid) - - if !all(n -> (n % r == 0), sz.-1) - throw(DomainError(r, "Size minus 1 must be divisible by the ratio.")) - end - - new_sz = (sz .- 1).÷r .+ 1 - - return EquidistantGrid{dimension(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) -end
--- a/src/Grids/Grids.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/Grids/Grids.jl Fri Feb 03 23:02:46 2023 +0100 @@ -2,18 +2,30 @@ using Sbplib.RegionIndices -export BoundaryIdentifier, CartesianBoundary +# Grid +export Grid +export dims +export points +export evalOn -abstract type BoundaryIdentifier end -struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end -dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim -region(::CartesianBoundary{Dim, R}) where {Dim, R} = R() +# BoundaryIdentifier +export BoundaryIdentifier +export CartesianBoundary +export dim +export region -export dim, region +# EquidistantGrid +export EquidistantGrid +export spacing +export inverse_spacing +export restrict +export boundary_identifiers +export boundary_grid +export refine +export coarsen -include("AbstractGrid.jl") -include("EquidistantGrid.jl") - -# TODO: Rename AbstractGrid to Grid and move definition here. +include("grid.jl") +include("boundary_identifier.jl") +include("equidistant_grid.jl") end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/boundary_identifier.jl Fri Feb 03 23:02:46 2023 +0100 @@ -0,0 +1,6 @@ + +abstract type BoundaryIdentifier end + +struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end +dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim +region(::CartesianBoundary{Dim, R}) where {Dim, R} = R() \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/equidistant_grid.jl Fri Feb 03 23:02:46 2023 +0100 @@ -0,0 +1,191 @@ + +""" + EquidistantGrid{Dim,T<:Real} <: Grid + +`Dim`-dimensional equidistant grid with coordinates of type `T`. +""" +struct EquidistantGrid{Dim,T<:Real} <: Grid + size::NTuple{Dim, Int} + limit_lower::NTuple{Dim, T} + limit_upper::NTuple{Dim, T} + + function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T} + if any(size .<= 0) + throw(DomainError("all components of size must be postive")) + end + if any(limit_upper.-limit_lower .<= 0) + throw(DomainError("all side lengths must be postive")) + end + return new{Dim,T}(size, limit_lower, limit_upper) + end +end + + +""" + EquidistantGrid(size, limit_lower, limit_upper) + +Construct an equidistant grid with corners at the coordinates `limit_lower` and +`limit_upper`. + +The length of the domain sides are given by the components of +`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined +as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative. + +The number of equidistantly spaced points in each coordinate direction are given +by the tuple `size`. +""" +function EquidistantGrid(size, limit_lower, limit_upper) + return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper) +end + + +""" + EquidistantGrid{T}() + +Constructs a 0-dimensional grid. +""" +EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid + + +""" + EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) + +Convenience constructor for 1D grids. +""" +function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T + return EquidistantGrid((size,),(limit_lower,),(limit_upper,)) +end + +Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T + +Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size) + +Base.size(g::EquidistantGrid) = g.size + +Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim + + + + + +""" + spacing(grid::EquidistantGrid) + +The spacing between grid points. +""" +spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1) + + +""" + inverse_spacing(grid::EquidistantGrid) + +The reciprocal of the spacing between grid points. +""" +inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid) + + +""" + points(grid::EquidistantGrid) + +The point of the grid as an array of tuples with the same dimension as the grid. +The points are stored as [(x1,y1), (x1,y2), … (x1,yn); + (x2,y1), (x2,y2), … (x2,yn); + ⋮ ⋮ ⋮ + (xm,y1), (xm,y2), … (xm,yn)] +""" +function points(grid::EquidistantGrid) + indices = Tuple.(CartesianIndices(grid.size)) + h = spacing(grid) + return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) +end + + +""" + restrict(::EquidistantGrid, dim) + +Pick out given dimensions from the grid and return a grid for them. +""" +function restrict(grid::EquidistantGrid, dim) + size = grid.size[dim] + limit_lower = grid.limit_lower[dim] + limit_upper = grid.limit_upper[dim] + + return EquidistantGrid(size, limit_lower, limit_upper) +end + + +""" + orthogonal_dims(grid::EquidistantGrid,dim) + +Returns the dimensions of grid orthogonal to that of dim. +""" +function orthogonal_dims(grid::EquidistantGrid, dim) + orth_dims = filter(i -> i != dim, dims(grid)) + if orth_dims == dims(grid) + throw(DomainError(string("dimension ",string(dim)," not matching grid"))) + end + return orth_dims +end + + +""" + boundary_identifiers(::EquidistantGrid) + +Returns a tuple containing the boundary identifiers for the grid, stored as + (CartesianBoundary(1,Lower), + CartesianBoundary(1,Upper), + CartesianBoundary(2,Lower), + ...) +""" +boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,) + + +""" + boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) + +Creates the lower-dimensional restriciton of `grid` spanned by the dimensions +orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional +grid is a zero-dimensional grid. +""" +function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) + orth_dims = orthogonal_dims(grid, dim(id)) + return restrict(grid, orth_dims) +end +boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}() + + +""" + refine(grid::EquidistantGrid, r::Int) + +Refines `grid` by a factor `r`. The factor is applied to the number of +intervals which is 1 less than the size of the grid. + +See also: [`coarsen`](@ref) +""" +function refine(grid::EquidistantGrid, r::Int) + sz = size(grid) + new_sz = (sz .- 1).*r .+ 1 + return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) +end + + +""" + coarsen(grid::EquidistantGrid, r::Int) + +Coarsens `grid` by a factor `r`. The factor is applied to the number of +intervals which is 1 less than the size of the grid. If the number of +intervals are not divisible by `r` an error is raised. + +See also: [`refine`](@ref) +""" +function coarsen(grid::EquidistantGrid, r::Int) + sz = size(grid) + + if !all(n -> (n % r == 0), sz.-1) + throw(DomainError(r, "Size minus 1 must be divisible by the ratio.")) + end + + new_sz = (sz .- 1).÷r .+ 1 + + return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/grid.jl Fri Feb 03 23:02:46 2023 +0100 @@ -0,0 +1,27 @@ +""" + Grid + +Should implement + Base.ndims(grid::Grid) + points(grid::Grid) + +""" +abstract type Grid end +function points end + +""" + dims(grid::Grid) + +A range containing the dimensions of `grid` +""" +dims(grid::Grid) = 1:ndims(grid) + +""" + evalOn(grid::Grid, f::Function) + +Evaluate function `f` on `grid` +""" +function evalOn(grid::Grid, f::Function) + F(x) = f(x...) + return F.(points(grid)) +end
--- a/src/LazyTensors/LazyTensors.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/LazyTensors/LazyTensors.jl Fri Feb 03 23:02:46 2023 +0100 @@ -3,9 +3,10 @@ export TensorApplication export TensorTranspose export TensorComposition -export DenseTensor export IdentityTensor export ScalingTensor +export DiagonalTensor +export DenseTensor export InflatedTensor export LazyOuterProduct export ⊗
--- a/src/LazyTensors/lazy_array.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/LazyTensors/lazy_array.jl Fri Feb 03 23:02:46 2023 +0100 @@ -28,7 +28,7 @@ export LazyFunctionArray function LazyFunctionArray(f::F, size::NTuple{D,Int}) where {F<:Function,D} - T = typeof(f(ones(D)...)) + T = typeof(f(ones(Int, D)...)) return LazyFunctionArray{F,T,D}(f,size) end @@ -93,17 +93,22 @@ -# NOTE: Är det knas att vi har till exempel * istället för .* ?? -# Oklart om det ens går att lösa.. +# Overload +,-,*,/ for LazyArrays +# Element wise operation for `*` and `/` are not overloaded for due to conflicts with the behavior +# of regular `*` and `/` for AbstractArrays. Use tilde versions instead. Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b +Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b + Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a +̃ b -Base.@propagate_inbounds Base.:+(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b +Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a -̃ b -Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b -Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::AbstractArray{T,D}) where {T,D} = a -̃ b +Base.@propagate_inbounds Base.:+(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a +̃ b Base.@propagate_inbounds Base.:-(a::AbstractArray{T,D}, b::LazyArray{T,D}) where {T,D} = a -̃ b -# Element wise operation for `*` and `\` are not overloaded due to conflicts with the behavior -# of regular `*` and `/` for AbstractArrays. Use tilde versions instead. +Base.@propagate_inbounds Base.:+(a::LazyArray{T,D}, b::T) where {T,D} = a +̃ b +Base.@propagate_inbounds Base.:-(a::LazyArray{T,D}, b::T) where {T,D} = a -̃ b + +Base.@propagate_inbounds Base.:+(a::T, b::LazyArray{T,D}) where {T,D} = a +̃ b +Base.@propagate_inbounds Base.:-(a::T, b::LazyArray{T,D}) where {T,D} = a -̃ b export +̃, -̃, *̃, /̃
--- a/src/LazyTensors/lazy_tensor_operations.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/LazyTensors/lazy_tensor_operations.jl Fri Feb 03 23:02:46 2023 +0100 @@ -98,7 +98,6 @@ apply_transpose(c.t2, c.t1'*v, I...) end - """ TensorComposition(tm, tmi::IdentityTensor) TensorComposition(tmi::IdentityTensor, tm) @@ -120,6 +119,8 @@ return tmi end +Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm) +Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm """ InflatedTensor{T,R,D} <: LazyTensor{T,R,D} @@ -268,6 +269,29 @@ LazyOuterProduct(tms::Vararg{LazyTensor}) = foldl(LazyOuterProduct, tms) + +""" + inflate(tm::LazyTensor, sz, dir) + +Inflate `tm` such that it gets the size `sz` in all directions except `dir`. +Here `sz[dir]` is ignored and replaced with the range and domains size of +`tm`. + +An example of when this operation is useful is when extending a one +dimensional difference operator `D` to a 2D grid of a ceratin size. In that +case we could have + +```julia +Dx = inflate(D, (10,10), 1) +Dy = inflate(D, (10,10), 2) +``` +""" +function inflate(tm::LazyTensor, sz, dir) + Is = IdentityTensor{eltype(tm)}.(sz) + parts = Base.setindex(Is, tm, dir) + return foldl(⊗, parts) +end + function check_domain_size(tm::LazyTensor, sz) if domain_size(tm) != sz throw(DomainSizeMismatch(tm,sz))
--- a/src/LazyTensors/tensor_types.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/LazyTensors/tensor_types.jl Fri Feb 03 23:02:46 2023 +0100 @@ -37,6 +37,24 @@ """ + DiagonalTensor{T,D,...} <: LazyTensor{T,D,D} + DiagonalTensor(a::AbstractArray) + +A lazy tensor with diagonal `a`. +""" +struct DiagonalTensor{T,D,AT<:AbstractArray{T,D}} <: LazyTensor{T,D,D} + diagonal::AT +end + +range_size(tm::DiagonalTensor) = size(tm.diagonal) +domain_size(tm::DiagonalTensor) = size(tm.diagonal) + + +LazyTensors.apply(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...] +LazyTensors.apply_transpose(tm::DiagonalTensor{T,D}, v::AbstractArray{<:Any,D}, I::Vararg{Any,D}) where {T,D} = tm.diagonal[I...]*v[I...] + + +""" DenseTensor{T,R,D,...}(A, range_indicies, domain_indicies) LazyTensor defined by the AbstractArray A. `range_indicies` and `domain_indicies` define which indicies of A should
--- a/src/SbpOperators/boundaryops/boundary_operator.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Fri Feb 03 23:02:46 2023 +0100 @@ -1,32 +1,3 @@ -""" - boundary_operator(grid,closure_stencil,boundary) - -Creates a boundary operator on a `Dim`-dimensional grid for the -specified `boundary`. The action of the operator is determined by `closure_stencil`. - -When `Dim=1`, the corresponding `BoundaryOperator` tensor mapping is returned. -When `Dim>1`, the `BoundaryOperator` `op` is inflated by the outer product -of `IdentityTensors` in orthogonal coordinate directions, e.g for `Dim=3`, -the boundary restriction operator in the y-direction direction is `Ix⊗op⊗Iz`. -""" -function boundary_operator(grid::EquidistantGrid, closure_stencil, boundary::CartesianBoundary) - #TODO:Check that dim(boundary) <= Dim? - - # Create 1D boundary operator - r = region(boundary) - d = dim(boundary) - op = BoundaryOperator(restrict(grid, d), closure_stencil, r) - - # Create 1D IdentityTensors for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) - Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids)) - - # Formulate the correct outer product sequence of the identity mappings and - # the boundary operator - parts = Base.setindex(Is, op, d) - return foldl(⊗, parts) -end - """ BoundaryOperator{T,R,N} <: LazyTensor{T,0,1} @@ -41,8 +12,6 @@ size::Int end -BoundaryOperator{R}(stencil::Stencil{T,N}, size::Int) where {T,R,N} = BoundaryOperator{T,R,N}(stencil, size) - """ BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region) @@ -55,6 +24,7 @@ """ closure_size(::BoundaryOperator) + The size of the closure stencil. """ closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Fri Feb 03 23:02:46 2023 +0100 @@ -8,11 +8,13 @@ On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of a `BoundaryOperator`. -See also: [`boundary_operator`](@ref). +See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref). """ function boundary_restriction(grid, closure_stencil, boundary) converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil) - return SbpOperators.boundary_operator(grid, converted_stencil, boundary) + + op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary)) + return LazyTensors.inflate(op, size(grid), dim(boundary)) end """
--- a/src/SbpOperators/boundaryops/normal_derivative.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Fri Feb 03 23:02:46 2023 +0100 @@ -8,12 +8,14 @@ On a one-dimensional `grid`, `d` is a `BoundaryOperator`. On a multi-dimensional `grid`, `d` is the inflation of a `BoundaryOperator`. -See also: [`boundary_operator`](@ref). +See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref). """ function normal_derivative(grid, closure_stencil, boundary) direction = dim(boundary) h_inv = inverse_spacing(grid)[direction] - return SbpOperators.boundary_operator(grid, scale(closure_stencil,h_inv), boundary) + + op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary)) + return LazyTensors.inflate(op, size(grid), dim(boundary)) end """
--- a/src/SbpOperators/stencil_set.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/stencil_set.jl Fri Feb 03 23:02:46 2023 +0100 @@ -14,7 +14,7 @@ """ -read_stencil_set(filename; filters) + read_stencil_set(filename; filters) Creates a `StencilSet` from a TOML file based on some key-value filters. If more than one set matches the filters an error is raised. The
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl Fri Feb 03 23:02:46 2023 +0100 @@ -7,14 +7,16 @@ `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` for the points in the closure regions. -On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the outer product of the -one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections. +On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of +a `VolumeOperator`. -See also: [`volume_operator`](@ref). +See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref). """ function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) h_inv = inverse_spacing(grid)[direction] - return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd, direction) + + D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd) + return LazyTensors.inflate(D₁, size(grid), direction) end
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl Fri Feb 03 23:02:46 2023 +0100 @@ -7,14 +7,16 @@ `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` for the points in the closure regions. -On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the outer product of the -one-dimensional operator with the `IdentityTensor`s in orthogonal coordinate dirrections. +On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of +a `VolumeOperator`. -See also: [`volume_operator`](@ref). +See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref). """ function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) h_inv = inverse_spacing(grid)[direction] - return SbpOperators.volume_operator(grid, scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even, direction) + + D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even) + return LazyTensors.inflate(D₂, size(grid), direction) end
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Fri Feb 03 23:02:46 2023 +0100 @@ -18,7 +18,7 @@ function inner_product(grid::EquidistantGrid, interior_weight, closure_weights) Hs = () - for i ∈ 1:dimension(grid) + for i ∈ dims(grid) Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights)) end
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Fri Feb 03 23:02:46 2023 +0100 @@ -15,7 +15,7 @@ function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights) H⁻¹s = () - for i ∈ 1:dimension(grid) + for i ∈ dims(grid) H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights)) end
--- a/src/SbpOperators/volumeops/laplace/laplace.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Fri Feb 03 23:02:46 2023 +0100 @@ -48,7 +48,7 @@ """ function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) - for d = 2:dimension(grid) + for d = 2:ndims(grid) Δ += second_derivative(grid, inner_stencil, closure_stencils, d) end return Δ
--- a/src/SbpOperators/volumeops/volume_operator.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/src/SbpOperators/volumeops/volume_operator.jl Fri Feb 03 23:02:46 2023 +0100 @@ -1,30 +1,6 @@ -""" - volume_operator(grid, inner_stencil, closure_stencils, parity, direction) - -Creates a volume operator on a `Dim`-dimensional grid acting along the -specified coordinate `direction`. The action of the operator is determined by -the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the -corresponding `VolumeOperator` tensor mapping is returned. When `Dim>1`, the -returned operator is the appropriate outer product of a one-dimensional -operators and `IdentityTensor`s, e.g for `Dim=3` the volume operator in the -y-direction is `I⊗op⊗I`. -""" -function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction) - #TODO: Check that direction <= Dim? - - # Create 1D volume operator in along coordinate direction - op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) - # Create 1D IdentityTensors for each coordinate direction - one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid))) - Is = IdentityTensor{eltype(grid)}.(size.(one_d_grids)) - # Formulate the correct outer product sequence of the identity mappings and - # the volume operator - parts = Base.setindex(Is, op, direction) - return foldl(⊗, parts) -end - """ VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1} + Implements a one-dimensional constant coefficients volume operator """ struct VolumeOperator{T,N,M,K} <: LazyTensor{T,1,1} @@ -59,3 +35,4 @@ r = getregion(i, closure_size(op), op.size[1]) return LazyTensors.apply(op, v, Index(i, r)) end +# TODO: Move this to LazyTensors when we have the region communication down.
--- a/test/Grids/EquidistantGrid_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,123 +0,0 @@ -using Sbplib.Grids -using Test -using Sbplib.RegionIndices - - -@testset "EquidistantGrid" begin - @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid - @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid - # constuctor - @test_throws DomainError EquidistantGrid(0,0.0,1.0) - @test_throws DomainError EquidistantGrid(1,1.0,1.0) - @test_throws DomainError EquidistantGrid(1,1.0,-1.0) - @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,)) - - @testset "Base" begin - @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64 - @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int - @test size(EquidistantGrid(4,0.0,1.0)) == (4,) - @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3) - end - - # dimension - @test dimension(EquidistantGrid(4,0.0,1.0)) == 1 - @test dimension(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2 - - # spacing - @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13 - @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13 - - # inverse_spacing - @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13 - @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13 - - # points - g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11)) - gp = points(g); - p = [(-1.,0.) (-1.,7.11/2) (-1.,7.11); - (-0.75,0.) (-0.75,7.11/2) (-0.75,7.11); - (-0.5,0.) (-0.5,7.11/2) (-0.5,7.11); - (-0.25,0.) (-0.25,7.11/2) (-0.25,7.11); - (0.,0.) (0.,7.11/2) (0.,7.11)] - for i ∈ eachindex(gp) - @test [gp[i]...] ≈ [p[i]...] atol=5e-13 - end - - # restrict - g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0)) - @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0) - @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0) - - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0) - @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0) - @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0) - @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0)) - @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) - @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0)) - @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0)) - - @testset "boundary_identifiers" begin - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(), - CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(), - CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}()) - @test boundary_identifiers(g) == bids - @inferred boundary_identifiers(g) - end - - @testset "boundary_grid" begin - @testset "1D" begin - g = EquidistantGrid(5,0.0,2.0) - (id_l, id_r) = boundary_identifiers(g) - @test boundary_grid(g,id_l) == EquidistantGrid{Float64}() - @test boundary_grid(g,id_r) == EquidistantGrid{Float64}() - @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}()) - @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}()) - end - @testset "2D" begin - g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) - (id_w, id_e, id_s, id_n) = boundary_identifiers(g) - @test boundary_grid(g,id_w) == restrict(g,2) - @test boundary_grid(g,id_e) == restrict(g,2) - @test boundary_grid(g,id_s) == restrict(g,1) - @test boundary_grid(g,id_n) == restrict(g,1) - @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) - end - @testset "3D" begin - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - (id_w, id_e, - id_s, id_n, - id_t, id_b) = boundary_identifiers(g) - @test boundary_grid(g,id_w) == restrict(g,[2,3]) - @test boundary_grid(g,id_e) == restrict(g,[2,3]) - @test boundary_grid(g,id_s) == restrict(g,[1,3]) - @test boundary_grid(g,id_n) == restrict(g,[1,3]) - @test boundary_grid(g,id_t) == restrict(g,[1,2]) - @test boundary_grid(g,id_b) == restrict(g,[1,2]) - @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) - end - end - - @testset "refine" begin - @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() - @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() - - g = EquidistantGrid((10,5),(0.,1.),(2.,3.)) - @test refine(g, 1) == g - @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.)) - @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.)) - end - - @testset "coarsen" begin - @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() - @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() - - g = EquidistantGrid((7,13),(0.,1.),(2.,3.)) - @test coarsen(g, 1) == g - @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.)) - @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) - - @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) - end -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/equidistant_grid_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -0,0 +1,125 @@ +using Sbplib.Grids +using Test +using Sbplib.RegionIndices + + +@testset "EquidistantGrid" begin + @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid + @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid + # constuctor + @test_throws DomainError EquidistantGrid(0,0.0,1.0) + @test_throws DomainError EquidistantGrid(1,1.0,1.0) + @test_throws DomainError EquidistantGrid(1,1.0,-1.0) + @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,)) + + @testset "Base" begin + @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64 + @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int + @test size(EquidistantGrid(4,0.0,1.0)) == (4,) + @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3) + @test ndims(EquidistantGrid(4,0.0,1.0)) == 1 + @test ndims(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2 + end + + @testset "spacing" begin + @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13 + @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13 + end + + @testset "inverse_spacing" begin + @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13 + @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13 + end + + @testset "points" begin + g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11)) + gp = points(g); + p = [(-1.,0.) (-1.,7.11/2) (-1.,7.11); + (-0.75,0.) (-0.75,7.11/2) (-0.75,7.11); + (-0.5,0.) (-0.5,7.11/2) (-0.5,7.11); + (-0.25,0.) (-0.25,7.11/2) (-0.25,7.11); + (0.,0.) (0.,7.11/2) (0.,7.11)] + for i ∈ eachindex(gp) + @test [gp[i]...] ≈ [p[i]...] atol=5e-13 + end + end + + @testset "restrict" begin + g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0)) + @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0) + @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0) + + g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) + @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0) + @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0) + @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0) + @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0)) + @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) + @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0)) + @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0)) + end + + @testset "boundary_identifiers" begin + g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) + bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(), + CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(), + CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}()) + @test boundary_identifiers(g) == bids + @inferred boundary_identifiers(g) + end + + @testset "boundary_grid" begin + @testset "1D" begin + g = EquidistantGrid(5,0.0,2.0) + (id_l, id_r) = boundary_identifiers(g) + @test boundary_grid(g,id_l) == EquidistantGrid{Float64}() + @test boundary_grid(g,id_r) == EquidistantGrid{Float64}() + @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}()) + @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}()) + end + @testset "2D" begin + g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) + (id_w, id_e, id_s, id_n) = boundary_identifiers(g) + @test boundary_grid(g,id_w) == restrict(g,2) + @test boundary_grid(g,id_e) == restrict(g,2) + @test boundary_grid(g,id_s) == restrict(g,1) + @test boundary_grid(g,id_n) == restrict(g,1) + @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) + end + @testset "3D" begin + g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) + (id_w, id_e, + id_s, id_n, + id_t, id_b) = boundary_identifiers(g) + @test boundary_grid(g,id_w) == restrict(g,[2,3]) + @test boundary_grid(g,id_e) == restrict(g,[2,3]) + @test boundary_grid(g,id_s) == restrict(g,[1,3]) + @test boundary_grid(g,id_n) == restrict(g,[1,3]) + @test boundary_grid(g,id_t) == restrict(g,[1,2]) + @test boundary_grid(g,id_b) == restrict(g,[1,2]) + @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) + end + end + + @testset "refine" begin + @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() + @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() + + g = EquidistantGrid((10,5),(0.,1.),(2.,3.)) + @test refine(g, 1) == g + @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.)) + @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.)) + end + + @testset "coarsen" begin + @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() + @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() + + g = EquidistantGrid((7,13),(0.,1.),(2.,3.)) + @test coarsen(g, 1) == g + @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.)) + @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) + + @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) + end +end
--- a/test/LazyTensors/lazy_array_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/LazyTensors/lazy_array_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -67,11 +67,18 @@ @test isa(v1 + v2, LazyArray) @test isa(v2 + v1, LazyArray) @test isa(v1 - v2, LazyArray) - @test isa(v2 - v1, LazyArray) + @test isa(v2 - v1, LazyArray) + @test isa(v1 + s, LazyArray) + @test isa(s + v1, LazyArray) + @test isa(v1 - s, LazyArray) + @test isa(s - v1, LazyArray) for i ∈ eachindex(v2) @test (v1 + v2)[i] == (v2 + v1)[i] == r_add_v[i] @test (v1 - v2)[i] == -(v2 - v1)[i] == r_sub_v[i] + @test (v1 + s)[i] == (s + v1)[i] == r_add_s[i] + @test (v1 - s)[i] == -(s - v1)[i] == r_sub_s[i] end + @test_throws BoundsError (v1 + v2)[4] v2 = [1., 2, 3, 4] @test_throws DimensionMismatch v1 + v2 @@ -95,4 +102,7 @@ @test_throws BoundsError LazyFunctionArray((i,j)->i*j, (3,2))[4,2] @test_throws BoundsError LazyFunctionArray((i,j)->i*j, (3,2))[2,3] + # Test that the constructor works with a restrictive function + f(x::Vararg{Int}) = sum(x) + @test LazyFunctionArray(f,(3,4)) isa LazyFunctionArray end
--- a/test/LazyTensors/lazy_tensor_operations_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -181,6 +181,14 @@ @test (Ã∘B̃*ComplexF64[1.,2.,3.,4.])[1] isa ComplexF64 @test ((Ã∘B̃)'*ComplexF64[1.,2.])[1] isa ComplexF64 + + a = 2. + v = rand(3) + @test a*Ã isa TensorComposition + @test a*Ã == Ã*a + @test range_size(a*Ã) == range_size(Ã) + @test domain_size(a*Ã) == domain_size(Ã) + @test a*Ã*v == a.*A*v end @@ -358,3 +366,17 @@ @test I1⊗Ã⊗I2 == InflatedTensor(I1, Ã, I2) end end + +@testset "inflate" begin + I = LazyTensors.inflate(IdentityTensor(),(3,4,5,6), 2) + @test I isa LazyTensor{Float64, 3,3} + @test range_size(I) == (3,5,6) + @test domain_size(I) == (3,5,6) + + @test LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 1) == InflatedTensor(IdentityTensor{Float64}(),ScalingTensor(1., (4,)),IdentityTensor(4,5,6)) + @test LazyTensors.inflate(ScalingTensor(2., (1,)),(3,4,5,6), 2) == InflatedTensor(IdentityTensor(3),ScalingTensor(2., (1,)),IdentityTensor(5,6)) + @test LazyTensors.inflate(ScalingTensor(3., (6,)),(3,4,5,6), 4) == InflatedTensor(IdentityTensor(3,4,5),ScalingTensor(3., (6,)),IdentityTensor{Float64}()) + + @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 0) + @test_throws BoundsError LazyTensors.inflate(ScalingTensor(1., (4,)),(3,4,5,6), 5) +end
--- a/test/LazyTensors/tensor_types_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/LazyTensors/tensor_types_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -1,5 +1,6 @@ using Test using Sbplib.LazyTensors +using BenchmarkTools @testset "IdentityTensor" begin @test IdentityTensor{Float64}((4,5)) isa IdentityTensor{T,2} where T @@ -58,6 +59,47 @@ @inferred (st'*v)[2,2] end +@testset "DiagonalTensor" begin + @test DiagonalTensor([1,2,3,4]) isa LazyTensor{Int,1,1} + @test DiagonalTensor([1 2 3; 4 5 6]) isa LazyTensor{Int,2,2} + @test DiagonalTensor([1. 2. 3.; 4. 5. 6.]) isa LazyTensor{Float64,2,2} + + @test range_size(DiagonalTensor([1,2,3,4])) == (4,) + @test domain_size(DiagonalTensor([1,2,3,4])) == (4,) + + @test range_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3) + @test domain_size(DiagonalTensor([1 2 3; 4 5 6])) == (2,3) + + @testset "apply size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)] + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + + @test tm*v == diag.*v + @test tm'*v == diag.*v + end + + + @testset "allocations size=$sz" for sz ∈ [(4,),(3,2),(3,4,2)] + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + + @test tm*v == diag.*v + @test tm'*v == diag.*v + end + + sz = (3,2) + diag = rand(sz...) + tm = DiagonalTensor(diag) + + v = rand(sz...) + LazyTensors.apply(tm,v, 2,1) + @test (@ballocated LazyTensors.apply($tm,$v, 2,1)) == 0 +end + @testset "DenseTensor" begin # Test a standard matrix-vector product
--- a/test/Manifest.toml Thu Apr 07 17:02:51 2022 +0200 +++ b/test/Manifest.toml Fri Feb 03 23:02:46 2023 +0100 @@ -1,10 +1,12 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.7.1" +julia_version = "1.8.5" manifest_format = "2.0" +project_hash = "23260eda65ade7d11fffed313a68520d0bc053fc" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" @@ -14,31 +16,32 @@ [[deps.BenchmarkTools]] deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] -git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93" +git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8" uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -version = "1.3.1" +version = "1.3.2" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f" +git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.14.0" +version = "1.15.7" [[deps.ChangesOfVariables]] deps = ["ChainRulesCore", "LinearAlgebra", "Test"] -git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1" +git-tree-sha1 = "844b061c104c408b24537482469400af6075aae4" uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" -version = "0.1.2" +version = "0.1.5" [[deps.Compat]] -deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4" +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.42.0" +version = "4.6.0" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.1+0" [[deps.Dates]] deps = ["Printf"] @@ -49,15 +52,11 @@ uuid = "ab62b9b5-e342-54a8-a765-a90f495de1a6" version = "1.2.0" -[[deps.DelimitedFiles]] -deps = ["Mmap"] -uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" - [[deps.DiffRules]] deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] -git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0" +git-tree-sha1 = "c5b6685d53f933c11404a3ae9822afe30d522494" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" -version = "1.10.0" +version = "1.12.2" [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] @@ -65,13 +64,17 @@ [[deps.DocStringExtensions]] deps = ["LibGit2"] -git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b" +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.6" +version = "0.9.3" [[deps.Downloads]] -deps = ["ArgTools", "LibCURL", "NetworkOptions"] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" [[deps.Glob]] git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2" @@ -84,9 +87,9 @@ [[deps.InverseFunctions]] deps = ["Test"] -git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d" +git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f" uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.3" +version = "0.1.8" [[deps.IrrationalConstants]] git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" @@ -108,10 +111,12 @@ [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" [[deps.LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "7.84.0+0" [[deps.LibGit2]] deps = ["Base64", "NetworkOptions", "Printf", "SHA"] @@ -120,6 +125,7 @@ [[deps.LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -130,9 +136,9 @@ [[deps.LogExpFunctions]] deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71" +git-tree-sha1 = "45b288af6956e67e621c5cbb2d75a261ab58300b" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.10" +version = "0.3.20" [[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" @@ -144,28 +150,34 @@ [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.0+0" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.2.1" [[deps.NaNMath]] -git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5" +deps = ["OpenLibm_jll"] +git-tree-sha1 = "a7c3d1da1189a1c2fe843a3bfa04d18d20eb3211" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.0.0" +version = "1.0.1" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -174,20 +186,21 @@ version = "0.5.5+0" [[deps.Parsers]] -deps = ["Dates"] -git-tree-sha1 = "621f4f3b4977325b9128d5fae7a8b4829a0c2222" +deps = ["Dates", "SnoopPrecompile"] +git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.2.4" +version = "2.5.4" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.8.0" [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd" +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.2.5" +version = "1.3.0" [[deps.Printf]] deps = ["Unicode"] @@ -213,13 +226,16 @@ [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -[[deps.SharedArrays]] -deps = ["Distributed", "Mmap", "Random", "Serialization"] -uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -230,9 +246,9 @@ [[deps.SpecialFunctions]] deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9" +git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.1.4" +version = "2.1.7" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -241,10 +257,12 @@ [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" [[deps.Tar]] deps = ["ArgTools", "SHA"] uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.1" [[deps.Test]] deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] @@ -258,9 +276,9 @@ [[deps.Tullio]] deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"] -git-tree-sha1 = "7830c974acc69437a3fee35dd7b510a74cbc862d" +git-tree-sha1 = "7871a39eac745697ee512a87eeff06a048a7905b" uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" -version = "0.3.3" +version = "0.3.5" [[deps.UUIDs]] deps = ["Random", "SHA"] @@ -272,15 +290,19 @@ [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.12+3" [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.1.1+0" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0"
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -6,7 +6,7 @@ using Sbplib.RegionIndices import Sbplib.SbpOperators.Stencil import Sbplib.SbpOperators.BoundaryOperator -import Sbplib.SbpOperators.boundary_operator + @testset "BoundaryOperator" begin closure_stencil = Stencil(2.,1.,3.; center = 1) @@ -14,120 +14,50 @@ g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) @testset "Constructors" begin - @testset "1D" begin - op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1]) - @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower()) - @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}()) - @test op_l isa LazyTensor{T,0,1} where T + @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T + @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T + end - op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) - @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper()) - @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}()) - @test op_r isa LazyTensor{T,0,1} where T - end - - @testset "2D" begin - e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}()) - @test e_w isa InflatedTensor - @test e_w isa LazyTensor{T,1,2} where T - end - end - op_l, op_r = boundary_operator.(Ref(g_1D), Ref(closure_stencil), boundary_identifiers(g_1D)) - op_w, op_e, op_s, op_n = boundary_operator.(Ref(g_2D), Ref(closure_stencil), boundary_identifiers(g_2D)) + op_l = BoundaryOperator(g_1D, closure_stencil, Lower()) + op_r = BoundaryOperator(g_1D, closure_stencil, Upper()) @testset "Sizes" begin - @testset "1D" begin - @test domain_size(op_l) == (11,) - @test domain_size(op_r) == (11,) - - @test range_size(op_l) == () - @test range_size(op_r) == () - end + @test domain_size(op_l) == (11,) + @test domain_size(op_r) == (11,) - @testset "2D" begin - @test domain_size(op_w) == (11,15) - @test domain_size(op_e) == (11,15) - @test domain_size(op_s) == (11,15) - @test domain_size(op_n) == (11,15) - - @test range_size(op_w) == (15,) - @test range_size(op_e) == (15,) - @test range_size(op_s) == (11,) - @test range_size(op_n) == (11,) - end + @test range_size(op_l) == () + @test range_size(op_r) == () end @testset "Application" begin - @testset "1D" begin - v = evalOn(g_1D,x->1+x^2) - u = fill(3.124) - @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] - @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] - @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] - @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] - @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] - - v = evalOn(g_1D, x->1. +x*im) - @test (op_l*v)[] isa ComplexF64 + v = evalOn(g_1D,x->1+x^2) + u = fill(3.124) + @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] + @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] + @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] + @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] + @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] - u = fill(1. +im) - @test (op_l'*u)[1] isa ComplexF64 - @test (op_l'*u)[5] isa ComplexF64 - @test (op_l'*u)[11] isa ComplexF64 - end - - @testset "2D" begin - v = rand(size(g_2D)...) - u = fill(3.124) - @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14 - @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14 - @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14 - @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14 - - - g_x = rand(size(g_2D)[1]) - g_y = rand(size(g_2D)[2]) - - G_w = zeros(Float64, size(g_2D)...) - G_w[1,:] = 2*g_y - G_w[2,:] = g_y - G_w[3,:] = 3*g_y + v = evalOn(g_1D, x->1. +x*im) + @test (op_l*v)[] isa ComplexF64 - G_e = zeros(Float64, size(g_2D)...) - G_e[end,:] = 2*g_y - G_e[end-1,:] = g_y - G_e[end-2,:] = 3*g_y - - G_s = zeros(Float64, size(g_2D)...) - G_s[:,1] = 2*g_x - G_s[:,2] = g_x - G_s[:,3] = 3*g_x - - G_n = zeros(Float64, size(g_2D)...) - G_n[:,end] = 2*g_x - G_n[:,end-1] = g_x - G_n[:,end-2] = 3*g_x + u = fill(1. +im) + @test (op_l'*u)[1] isa ComplexF64 + @test (op_l'*u)[5] isa ComplexF64 + @test (op_l'*u)[11] isa ComplexF64 - @test op_w'*g_y == G_w - @test op_e'*g_y == G_e - @test op_s'*g_x == G_s - @test op_n'*g_x == G_n - end + u = fill(3.124) + @test (op_l'*u)[Index(1,Lower)] == 2*u[] + @test (op_l'*u)[Index(2,Lower)] == u[] + @test (op_l'*u)[Index(6,Interior)] == 0 + @test (op_l'*u)[Index(10,Upper)] == 0 + @test (op_l'*u)[Index(11,Upper)] == 0 - @testset "Regions" begin - u = fill(3.124) - @test (op_l'*u)[Index(1,Lower)] == 2*u[] - @test (op_l'*u)[Index(2,Lower)] == u[] - @test (op_l'*u)[Index(6,Interior)] == 0 - @test (op_l'*u)[Index(10,Upper)] == 0 - @test (op_l'*u)[Index(11,Upper)] == 0 - - @test (op_r'*u)[Index(1,Lower)] == 0 - @test (op_r'*u)[Index(2,Lower)] == 0 - @test (op_r'*u)[Index(6,Interior)] == 0 - @test (op_r'*u)[Index(10,Upper)] == u[] - @test (op_r'*u)[Index(11,Upper)] == 2*u[] - end + @test (op_r'*u)[Index(1,Lower)] == 0 + @test (op_r'*u)[Index(2,Lower)] == 0 + @test (op_r'*u)[Index(6,Interior)] == 0 + @test (op_r'*u)[Index(10,Upper)] == u[] + @test (op_r'*u)[Index(11,Upper)] == 2*u[] end @testset "Inferred" begin
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -71,6 +71,7 @@ end @testset "Accuracy on function" begin + # 1D g = EquidistantGrid(30, 0.,1.) v = evalOn(g, x->exp(x)) @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)] @@ -79,6 +80,18 @@ @test D₁*v ≈ v rtol=tol end + + # 2D + g = EquidistantGrid((30,60), (0.,0.),(1.,2.)) + v = evalOn(g, (x,y)->exp(0.8x+1.2*y)) + @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) + Dx = first_derivative(g, stencil_set, 1) + Dy = first_derivative(g, stencil_set, 2) + + @test Dx*v ≈ 0.8v rtol=tol + @test Dy*v ≈ 1.2v rtol=tol + end end end
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -6,6 +6,8 @@ import Sbplib.SbpOperators.VolumeOperator +# TODO: Refactor these test to look more like the tests in first_derivative_test.jl. + @testset "SecondDerivative" begin operator_path = sbp_operators_path()*"standard_diagonal.toml" stencil_set = read_stencil_set(operator_path; order=4)
--- a/test/SbpOperators/volumeops/volume_operator_test.jl Thu Apr 07 17:02:51 2022 +0200 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Fri Feb 03 23:02:46 2023 +0100 @@ -7,120 +7,76 @@ import Sbplib.SbpOperators.Stencil import Sbplib.SbpOperators.VolumeOperator -import Sbplib.SbpOperators.volume_operator import Sbplib.SbpOperators.odd import Sbplib.SbpOperators.even + @testset "VolumeOperator" begin inner_stencil = CenteredStencil(1/4, 2/4, 1/4) - closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2)) - g_1D = EquidistantGrid(11,0.,1.) - g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.)) - g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.)) + closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2)) + g = EquidistantGrid(11,0.,1.) + @testset "Constructors" begin - @testset "1D" begin - op = VolumeOperator(inner_stencil,closure_stencils,(11,),even) - @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even) - @test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1) - @test op isa LazyTensor{T,1,1} where T - end - @testset "2D" begin - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2) - Ix = IdentityTensor{Float64}((11,)) - Iy = IdentityTensor{Float64}((12,)) - @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy - @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even) - @test op_x isa LazyTensor{T,2,2} where T - @test op_y isa LazyTensor{T,2,2} where T - end - @testset "3D" begin - op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2) - op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3) - Ix = IdentityTensor{Float64}((11,)) - Iy = IdentityTensor{Float64}((12,)) - Iz = IdentityTensor{Float64}((10,)) - @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz - @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz - @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even) - @test op_x isa LazyTensor{T,3,3} where T - @test op_y isa LazyTensor{T,3,3} where T - @test op_z isa LazyTensor{T,3,3} where T - end + op = VolumeOperator(inner_stencil,closure_stencils,(11,),even) + @test op == VolumeOperator(g,inner_stencil,closure_stencils,even) + @test op isa LazyTensor{T,1,1} where T end @testset "Sizes" begin - @testset "1D" begin - op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1) - @test range_size(op) == domain_size(op) == size(g_1D) - end - - @testset "2D" begin - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2) - @test range_size(op_y) == domain_size(op_y) == - range_size(op_x) == domain_size(op_x) == size(g_2D) - end - @testset "3D" begin - op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2) - op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3) - @test range_size(op_z) == domain_size(op_z) == - range_size(op_y) == domain_size(op_y) == - range_size(op_x) == domain_size(op_x) == size(g_3D) - end + op = VolumeOperator(g,inner_stencil,closure_stencils,even) + @test range_size(op) == domain_size(op) == size(g) end - op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) - op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2) - v = zeros(size(g_2D)) - Nx = size(g_2D)[1] - Ny = size(g_2D)[2] - for i = 1:Nx - v[i,:] .= i + + op_even = VolumeOperator(g, inner_stencil, closure_stencils, even) + op_odd = VolumeOperator(g, inner_stencil, closure_stencils, odd) + + N = size(g)[1] + v = rand(N) + + r_even = copy(v) + r_odd = copy(v) + + r_even[1] = (v[1] + v[2])/2 + r_odd[1] = (v[1] + v[2])/2 + + r_even[2] = 2v[1] + v[2] + r_odd[2] = 2v[1] + v[2] + + for i ∈ 3:N-2 + r_even[i] = (v[i-1] + 2v[i] + v[i+1])/4 + r_odd[i] = (v[i-1] + 2v[i] + v[i+1])/4 end - rx = copy(v) - rx[1,:] .= 1.5 - rx[Nx,:] .= (2*Nx-1)/2 - ry = copy(v) - ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny] + + r_even[N-1] = v[N-1] + 2v[N] + r_odd[N-1] = -v[N-1] - 2v[N] + + r_even[N] = (v[N-1] + v[N])/2 + r_odd[N] = -(v[N-1] + v[N])/2 + @testset "Application" begin - @test op_x*v ≈ rx rtol = 1e-14 - @test op_y*v ≈ ry rtol = 1e-14 + @test op_even*v ≈ r_even + @test op_odd*v ≈ r_odd - @test (op_x*rand(ComplexF64,size(g_2D)))[2,2] isa ComplexF64 + @test (op_even*rand(ComplexF64,size(g)))[2] isa ComplexF64 end @testset "Regions" begin - @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14 - @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14 - @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14 - @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14 - @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14 - - @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)] - @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)] + @test (op_even*v)[Index(1,Lower)] ≈ r_even[1] + @test (op_even*v)[Index(2,Lower)] ≈ r_even[2] + @test (op_even*v)[Index(6,Interior)] ≈ r_even[6] + @test (op_even*v)[Index(10,Upper)] ≈ r_even[10] + @test (op_even*v)[Index(11,Upper)] ≈ r_even[11] - @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14 - @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14 - - @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)] - @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)] + @test_throws BoundsError (op_even*v)[Index(3,Lower)] + @test_throws BoundsError (op_even*v)[Index(9,Upper)] end @testset "Inferred" begin - @test_skip @inferred apply(op_x, v,1,1) - @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower)) - @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower)) - @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower)) - @test_skip @inferred apply(op_y, v,1,1) - @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower)) - @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior)) - @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper)) + @inferred apply(op_even, v, 1) + @inferred apply(op_even, v, Index(1,Lower)) + @inferred apply(op_even, v, Index(6,Interior)) + @inferred apply(op_even, v, Index(11,Upper)) end end