Mercurial > repos > public > sbplib_julia
changeset 1355:102ebdaf7c11 feature/variable_derivatives
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 08 Feb 2023 21:21:28 +0100 |
| parents | fa0800591306 (current diff) 7ee258e5289e (diff) |
| children | 49d03d1169ef |
| files | src/Grids/EquidistantGrid.jl src/Grids/equidistant_grid.jl src/LazyTensors/lazy_tensor_operations.jl src/SbpOperators/SbpOperators.jl src/SbpOperators/boundaryops/boundary_operator.jl src/SbpOperators/stencil_set.jl src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl src/SbpOperators/volumeops/volume_operator.jl test/Grids/EquidistantGrid_test.jl test/Grids/equidistant_grid_test.jl test/SbpOperators/boundaryops/boundary_operator_test.jl |
| diffstat | 36 files changed, 862 insertions(+), 712 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Fri Feb 03 22:50:42 2023 +0100 +++ b/Manifest.toml Wed Feb 08 21:21:28 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 Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/docs/Manifest.toml Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ /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 Fri Feb 03 22:50:42 2023 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,192 +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 -# TBD: Should it be an AbstractArray? - -""" - 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 - -function Base.getindex(g::EquidistantGrid, I::Vararg{Int}) - h = spacing(g) - return g.limit_lower .+ (I.-1).*h -end - -Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...] -# TBD: Can this method be removed if `EquidistantGrid` is an AbstractArray? - -""" - 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/Grids/Grids.jl Wed Feb 08 21:21:28 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 Wed Feb 08 21:21:28 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 Wed Feb 08 21:21:28 2023 +0100 @@ -0,0 +1,195 @@ + +""" + 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 + +function Base.getindex(g::EquidistantGrid, I::Vararg{Int}) + h = spacing(g) + return g.limit_lower .+ (I.-1).*h +end + +Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...] +# TBD: Can this method be removed if `EquidistantGrid` is an AbstractArray? + +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 Wed Feb 08 21:21:28 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/lazy_array.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/LazyTensors/lazy_array.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/LazyTensors/lazy_tensor_operations.jl Wed Feb 08 21:21:28 2023 +0100 @@ -100,7 +100,6 @@ apply_transpose(c.t2, c.t1'*v, I...) end - """ TensorComposition(tm, tmi::IdentityTensor) TensorComposition(tmi::IdentityTensor, tm) @@ -122,6 +121,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} @@ -272,11 +273,20 @@ """ - inflate(tm, sz, dir) + 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`. -Inflate `tm` with identity tensors in all directions `d` for `d != dir`. +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 -# TODO: Describe when it is useful +```julia +Dx = inflate(D, (10,10), 1) +Dy = inflate(D, (10,10), 2) +``` """ function inflate(tm::LazyTensor, sz, dir) Is = IdentityTensor{eltype(tm)}.(sz)
--- a/src/LazyTensors/tuple_manipulation.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/LazyTensors/tuple_manipulation.jl Wed Feb 08 21:21:28 2023 +0100 @@ -75,7 +75,11 @@ flatten_tuple(t::Tuple) = ((flatten_tuple.(t)...)...,) # simplify? flatten_tuple(ts::Vararg) = flatten_tuple(ts) +""" + left_pad_tuple(t, val, N) +Left pad the tuple `t` to length `N` using the value `val`. +""" function left_pad_tuple(t, val, N) if N < length(t) throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements")) @@ -85,6 +89,11 @@ return (padding..., t...) end +""" + right_pad_tuple(t, val, N) + +Right pad the tuple `t` to length `N` using the value `val`. +""" function right_pad_tuple(t, val, N) if N < length(t) throw(DomainError(N, "Can't pad tuple of length $(length(t)) to $N elements"))
--- a/src/SbpOperators/SbpOperators.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/SbpOperators.jl Wed Feb 08 21:21:28 2023 +0100 @@ -20,7 +20,7 @@ export normal_derivative export first_derivative export second_derivative -export undivided_dissipation +export undivided_skewed04 using Sbplib.RegionIndices using Sbplib.LazyTensors
--- a/src/SbpOperators/boundaryops/boundary_operator.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Wed Feb 08 21:21:28 2023 +0100 @@ -1,27 +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? - - d = dim(boundary) - op = BoundaryOperator(restrict(grid, d), closure_stencil, region(boundary)) - - # 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)) - - return LazyTensors.inflate(op, size(grid), d) -end - """ BoundaryOperator{T,R,N} <: LazyTensor{T,0,1} @@ -36,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) @@ -50,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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/stencil_set.jl Wed Feb 08 21:21:28 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/dissipation.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl Wed Feb 08 21:21:28 2023 +0100 @@ -1,4 +1,16 @@ -function undivided_dissipation(g::EquidistantGrid, p, direction) +""" + undivided_skewed04(g::EquidistantGrid, p, direction) + +Create undivided difference operators approximating the `p`th derivative. The +operators do not satisfy any SBP-property and are meant to be used for +building artificial dissipation terms. + +The operators and how they are used to create accurate artifical dissipation +is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate +Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp. +57–79, Aug. 2004" +""" +function undivided_skewed04(g::EquidistantGrid, p, direction) T = eltype(g) interior_weights = T.(dissipation_interior_weights(p)) @@ -20,7 +32,7 @@ return D, Dᵀ end -undivided_dissipation(g::EquidistantGrid{1}, p) = undivided_dissipation(g, p, 1) +undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1) function dissipation_interior_weights(p) if p == 0 @@ -47,8 +59,16 @@ dissipation_lower_closure_size(weights) = midpoint(weights) - 1 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights) -dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i ), dissipation_lower_closure_size(interior_weights)) -dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights)) +function dissipation_lower_closure_stencils(interior_weights) + stencil(i) = Stencil(interior_weights..., center=i) + return ntuple(i->stencil(i), dissipation_lower_closure_size(interior_weights)) +end + +function dissipation_upper_closure_stencils(interior_weights) + center(i) = length(interior_weights) - dissipation_upper_closure_size(interior_weights) + i + stencil(i) = Stencil(interior_weights..., center=center(i)) + return ntuple(i->stencil(i), dissipation_upper_closure_size(interior_weights)) +end function dissipation_transpose_lower_closure_stencils(interior_weights) closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl Wed Feb 08 21:21:28 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/derivatives/second_derivative_variable.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl Wed Feb 08 21:21:28 2023 +0100 @@ -26,7 +26,7 @@ Δxᵢ = spacing(grid)[dir] scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2) scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2) - return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff) + return SecondDerivativeVariable{dir, ndims(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff) end function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils) @@ -61,8 +61,8 @@ end function check_coefficient(grid, coeff) - if dimension(grid) != ndims(coeff) - throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(dimension(grid))")) + if ndims(grid) != ndims(coeff) + throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(ndims(grid))")) end if size(grid) != size(coeff)
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Wed Feb 08 21:21:28 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/stencil_operator_distinct_closures.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl Wed Feb 08 21:21:28 2023 +0100 @@ -1,3 +1,17 @@ +""" + stencil_operator_distinct_closures( + grid::EquidistantGrid, + inner_stencil, + lower_closure, + upper_closure, + direction + ) + +Creates a multi-dimensional `StencilOperatorDistinctClosures` acting on grid +functions of `grid`. + +See also: [`StencilOperatorDistinctClosures`](@ref) +""" function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction) op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure) return LazyTensors.inflate(op, size(grid), direction) @@ -14,7 +28,7 @@ closures is useful for representing operators with skewed stencils like upwind operators. -See also: [`VolumeOperator`](@ref) +See also: [`VolumeOperator`](@ref), [`stencil_operator_distinct_closures`](@ref) """ struct StencilOperatorDistinctClosures{T,K,N,M,LC<:NTuple{N,Stencil{T,L}} where L, UC<:NTuple{M,Stencil{T,L}} where L} <: LazyTensor{T,1,1} inner_stencil::Stencil{T,K}
--- a/src/SbpOperators/volumeops/volume_operator.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/src/SbpOperators/volumeops/volume_operator.jl Wed Feb 08 21:21:28 2023 +0100 @@ -1,24 +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? - - op = VolumeOperator(restrict(grid, direction), inner_stencil, closure_stencils, parity) - return LazyTensors.inflate(op, size(grid), direction) -end -# TBD: Should the inflation happen here or should we remove this method and do it at the caller instead? - """ 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}
--- a/test/Grids/EquidistantGrid_test.jl Fri Feb 03 22:50:42 2023 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,136 +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 - - - @testset "getindex" begin - g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11)) - @test g[1,1] == (-1.0,0.0) - @test g[1,3] == (-1.0,7.11) - @test g[5,1] == (0.0,0.0) - @test g[5,3] == (0.0,7.11) - - @test g[4,2] == (-0.25,7.11/2) - - @test g[CartesianIndex(1,3)] == (-1.0,7.11) - 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 Wed Feb 08 21:21:28 2023 +0100 @@ -0,0 +1,137 @@ +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 "getindex" begin + g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11)) + @test g[1,1] == (-1.0,0.0) + @test g[1,3] == (-1.0,7.11) + @test g[5,1] == (0.0,0.0) + @test g[5,3] == (0.0,7.11) + + @test g[4,2] == (-0.25,7.11/2) + + @test g[CartesianIndex(1,3)] == (-1.0,7.11) + 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/test/LazyTensors/lazy_array_test.jl Wed Feb 08 21:21:28 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 Fri Feb 03 22:50:42 2023 +0100 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Wed Feb 08 21:21:28 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
--- a/test/Manifest.toml Fri Feb 03 22:50:42 2023 +0100 +++ b/test/Manifest.toml Wed Feb 08 21:21:28 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 = "85b5da0fa43588c75bb1ff986493443f821c70b7" +deps = ["Dates", "SnoopPrecompile"] +git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.2.3" +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 Fri Feb 03 22:50:42 2023 +0100 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Wed Feb 08 21:21:28 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/dissipation_test.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl Wed Feb 08 21:21:28 2023 +0100 @@ -26,19 +26,19 @@ x^k/factorial(k) end -@testset "undivided_dissipation" begin +@testset "undivided_skewed04" begin g = EquidistantGrid(20, 0., 11.) - D,Dᵀ = undivided_dissipation(g, 1) + D,Dᵀ = undivided_skewed04(g, 1) - @test D isa LazyTensor{Float64,1,1} where T - @test Dᵀ isa LazyTensor{Float64,1,1} where T + @test D isa LazyTensor{Float64,1,1} + @test Dᵀ isa LazyTensor{Float64,1,1} @testset "Accuracy conditions" begin N = 20 g = EquidistantGrid(N, 0//1,2//1) h = only(spacing(g)) @testset "D_$p" for p ∈ [1,2,3,4] - D,Dᵀ = undivided_dissipation(g, p) + D,Dᵀ = undivided_skewed04(g, p) @testset "x^$k" for k ∈ 0:p v = evalOn(g, x->monomial(x,k)) @@ -49,7 +49,7 @@ end end - @testset "tanspose equality" begin + @testset "transpose equality" begin function get_matrix(D) N = only(range_size(D)) M = only(domain_size(D)) @@ -69,7 +69,7 @@ g = EquidistantGrid(11, 0., 1.) @testset "D_$p" for p ∈ [1,2,3,4] - D,Dᵀ = undivided_dissipation(g, p) + D,Dᵀ = undivided_skewed04(g, p) D̄ = get_matrix(D) D̄ᵀ = get_matrix(Dᵀ)
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Wed Feb 08 21:21:28 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/stencil_operator_distinct_closures_test.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl Wed Feb 08 21:21:28 2023 +0100 @@ -2,7 +2,6 @@ using Sbplib.SbpOperators using Sbplib.Grids -# using Sbplib.RegionIndices using Sbplib.LazyTensors import Sbplib.SbpOperators.Stencil
--- a/test/SbpOperators/volumeops/volume_operator_test.jl Fri Feb 03 22:50:42 2023 +0100 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Wed Feb 08 21:21:28 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