Mercurial > repos > public > sbplib_julia

changeset 1594:d68d02dd882f feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Sat, 25 May 2024 16:07:10 -0700
parents 615eeb6e662e (current diff) efe1fc4cb6b0 (diff)
children 611ae2308aa1
files Notes.md src/LazyTensors/lazy_tensor_operations.jl test/SbpOperators/volumeops/laplace/laplace_test.jl
diffstat 30 files changed, 383 insertions(+), 160 deletions(-) [+]
line wrap: on
line diff
--- a/Manifest.toml	Tue Jan 23 20:48:25 2024 +0100
+++ b/Manifest.toml	Sat May 25 16:07:10 2024 -0700
@@ -1,30 +1,32 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.2"
 manifest_format = "2.0"
 project_hash = "a36735c53cfa4453f39635046eeaa47a4ea1231b"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0"
+git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.1"
+version = "4.0.4"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
     AdaptStaticArraysExt = "StaticArrays"
 
 [[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
+deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.7.0"
+version = "7.9.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
     ArrayInterfaceTrackerExt = "Tracker"
 
@@ -32,7 +34,9 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
@@ -41,9 +45,9 @@
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b"
+git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.12.0"
+version = "4.14.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -52,7 +56,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -82,19 +86,19 @@
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -124,9 +128,9 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
+git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.8"
+version = "0.8.10"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
@@ -141,9 +145,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
+git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.1"
+version = "1.9.3"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/Notes.md	Tue Jan 23 20:48:25 2024 +0100
+++ b/Notes.md	Sat May 25 16:07:10 2024 -0700
@@ -204,18 +204,41 @@
 
 The interaction of the map methods with the probable design of multiblock functions involving nested indecies complicate the picture slightly. It's clear at the time of writing how this would work with `Base.map`. Perhaps we want to implement our own versions of both eager and lazy map.
 
+
+### 2024-04
+MappedArrays.jl provides a simple array type and function like the description
+of LazyMapping above. One option is to remove `eval_on` completely and rely on
+destructuring arguments if handling the function input as a vector is
+undesirable.
+
+If we can let multi-block grids be iterators over grid points we could even
+handle those by specialized implementation of `map` and `mappedarray`.
+
 ## Multiblock implementation
 We want multiblock things to work very similarly to regular one block things.
 
 ### Grid functions
-Should probably support a nested indexing so that we first have an index for subgrid and then an index for nodes on that grid. E.g `g[1,2][2,3]` or `g[3][43,21]`.
+Should probably support a nested indexing so that we first have an index for
+subgrid and then an index for nodes on that grid. E.g `g[1,2][2,3]` or
+`g[3][43,21]`.
 
-We could also possibly provide a combined indexing style `g[1,2,3,4]` where the first group of indices are for the subgrid and the remaining are for the nodes.
+We could also possibly provide a combined indexing style `g[1,2,3,4]` where
+the first group of indices are for the subgrid and the remaining are for the
+nodes.
 
-We should make sure the underlying buffer for gridfunctions are continuously stored and are easy to convert to, so that interaction with for example DifferentialEquations is simple and without much boilerplate.
+We should make sure the underlying buffer for grid functions are continuously
+stored and are easy to convert to, so that interaction with for example
+DifferentialEquations is simple and without much boilerplate.
 
 #### `map` and `collect` and nested indexing
-We need to make sure `collect`, `map` and a potential lazy map work correctly through the nested indexing.
+We need to make sure `collect`, `map` and a potential lazy map work correctly
+through the nested indexing. Also see notes on `eval_on` above.
+
+Possibly this can be achieved by providing special nested indexing but not
+adhering to an array interface at the top level, instead being implemented as
+an iterator over the grid points. A custom trait can let map and other methods
+know the shape (or structure) of the nesting so that they can efficiently
+allocate result arrays.
 
 ### Tensor applications
 Should behave as grid functions
@@ -343,3 +366,11 @@
 * This would impact how tensor grid works.
 * To make things homogenous maybe these index collections should be used for the more simple grids too.
 * The function `to_indices` from Base could be useful to implement for `IndexCollection`
+
+
+## Stencil application pipeline
+We should make sure that `@inbounds` and `Base.@propagate_inbounds` are
+applied correctly throughout the stack. When testing the performance of
+stencil application on the bugfix/sbp_operators/stencil_return_type branch
+there seemed to be some strange results where such errors could be the
+culprit.
--- a/benchmark/Manifest.toml	Tue Jan 23 20:48:25 2024 +0100
+++ b/benchmark/Manifest.toml	Sat May 25 16:07:10 2024 -0700
@@ -1,19 +1,19 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.2"
 manifest_format = "2.0"
 project_hash = "25bba7b4a00465d5a2b00b589eb10e3301c31f2a"
 
 [[deps.AbstractTrees]]
-git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.4.4"
+version = "0.4.5"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0"
+git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.1"
+version = "4.0.4"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -24,16 +24,18 @@
 version = "1.1.1"
 
 [[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
+deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.7.0"
+version = "7.9.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
     ArrayInterfaceTrackerExt = "Tracker"
 
@@ -41,7 +43,9 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
@@ -53,15 +57,15 @@
 
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "f1f03a9fa24271160ed7e73051fba3c1a759b53f"
+git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.4.0"
+version = "1.5.0"
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b"
+git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.12.0"
+version = "4.14.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -70,7 +74,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.DataAPI]]
 git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
@@ -192,7 +196,7 @@
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OrderedCollections]]
 git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
@@ -218,15 +222,15 @@
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -279,9 +283,9 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
+git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.8"
+version = "0.8.10"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
@@ -296,9 +300,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
+git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.1"
+version = "1.9.3"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/benchmark/benchmark_laplace.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/benchmark/benchmark_laplace.jl	Sat May 25 16:07:10 2024 -0700
@@ -10,7 +10,7 @@
 
 function benchmark_const_coeff_1d(;N = 100, order = 4)
     stencil_set = read_stencil_set(operator_path; order=order)
-    g = equidistant_grid(N, 0., 1.)
+    g = equidistant_grid(0., 1., N)
     D = second_derivative(g, stencil_set)
     u = rand(size(g)...)
     u_xx = rand(size(g)...)
@@ -25,7 +25,7 @@
 
 function benchmark_var_coeff_1d(;N = 100, order = 4)
     stencil_set = read_stencil_set(operator_path; order=order)
-    g = equidistant_grid(N, 0., 1.)
+    g = equidistant_grid(0., 1., N)
     c = rand(size(g)...)
     c_lz = eval_on(g, x -> 0.5)
     D = second_derivative_variable(g, c, stencil_set)
@@ -49,7 +49,7 @@
 
 function benchmark_const_coeff_2d(;N = 100, order = 4)
     stencil_set = read_stencil_set(operator_path; order=order)
-    g = equidistant_grid((N,N), (0.,0.,),(1.,1.))
+    g = equidistant_grid((0.,0.,),(1.,1.), N, N)
     D = Laplace(g, stencil_set)
     u = rand(size(g)...)
     u_xx = rand(size(g)...)
@@ -71,7 +71,7 @@
 
 function benchmark_var_coeff_2d(;N = 100, order = 4)
     stencil_set = read_stencil_set(operator_path; order=order)
-    g = equidistant_grid((N,N), (0.,0.,),(1.,1.))
+    g = equidistant_grid((0.,0.,),(1.,1.), N, N)
     c = rand(size(g)...)
     c_lz = eval_on(g, x-> 0.5)
     D = second_derivative_variable(g, c, stencil_set, 1) + second_derivative_variable(g, c, stencil_set, 2)
@@ -217,4 +217,4 @@
     for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,end-clz_sz+1:end]
         u_xx[I] = tm[Index{Upper}(I[1]),Index{Upper}(I[2])]
     end
-end
\ No newline at end of file
+end
--- a/benchmark/benchmark_utils.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/benchmark/benchmark_utils.jl	Sat May 25 16:07:10 2024 -0700
@@ -23,7 +23,7 @@
 For control over what happens to the benchmark result datastructure see the
 different methods of [`run_benchmark`](@ref)
 """
-function main(;rev=nothing, target=nothing, baseline=nothing , kwargs...)
+function main(;rev=nothing, target=nothing, baseline=nothing, name=nothing, kwargs...)
     if !isnothing(rev)
         r = run_benchmark(rev; kwargs...)
     elseif !isnothing(baseline)
@@ -37,7 +37,7 @@
         r = run_benchmark(;kwargs...)
     end
 
-    file_path = write_result_html(r)
+    file_path = write_result_html(r; name)
     open_in_default_browser(file_path)
 end
 
@@ -137,9 +137,14 @@
     Mustache.render(io, template, Dict("title"=>dt, "content"=>content))
 end
 
-function write_result_html(r)
+function write_result_html(r; name=nothing)
     dt = Dates.format(PkgBenchmark.date(r), "yyyy-mm-dd HHMMSS")
-    file_path = joinpath(results_dir, dt*".html")
+
+    if isnothing(name)
+        file_path = joinpath(results_dir, dt*".html")
+    else
+        file_path = joinpath(results_dir, dt*" "*name*".html")
+    end
 
     open(file_path, "w") do io
         write_result_html(io, r)
--- a/benchmark/benchmarks.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/benchmark/benchmarks.jl	Sat May 25 16:07:10 2024 -0700
@@ -15,9 +15,9 @@
 ll(d) = ntuple(i->0., d)
 lu(d) = ntuple(i->1., d)
 
-g1 = equidistant_grid(sz(1)[1],ll(1)[1],lu(1)[1])
-g2 = equidistant_grid(sz(2),ll(2),lu(2))
-g3 = equidistant_grid(sz(3),ll(3),lu(3))
+g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)[1])
+g2 = equidistant_grid(ll(2), lu(2), sz(2))
+g3 = equidistant_grid(ll(3), lu(3), sz(3))
 
 v1 = rand(sz(1)...)
 v2 = rand(sz(2)...)
--- a/src/Grids/Grids.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/Grids/Grids.jl	Sat May 25 16:07:10 2024 -0700
@@ -13,10 +13,11 @@
 export boundary_indices
 export boundary_identifiers
 export boundary_grid
+export coarsen
+export refine
 export eval_on
-export coarsen
-export getcomponent
-export refine
+export componentview
+export ArrayComponentView
 
 export BoundaryIdentifier
 export TensorGridBoundary
--- a/src/Grids/equidistant_grid.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/Grids/equidistant_grid.jl	Sat May 25 16:07:10 2024 -0700
@@ -88,7 +88,7 @@
 
 
 """
-    equidistant_grid(size::Dims, limit_lower, limit_upper)
+    equidistant_grid(limit_lower, limit_upper, dims...)
 
 Construct an equidistant grid with corners at the coordinates `limit_lower` and
 `limit_upper`.
@@ -99,24 +99,27 @@
 of the grid are not allowed to be negative.
 
 The number of equispaced points in each coordinate direction are given
-by the tuple `size`.
+by the tuple `dims`.
 
-Note: If `limit_lower` and `limit_upper` are integers and `size` would allow a
+Note: If `limit_lower` and `limit_upper` are integers and `dims` would allow a
 completely integer grid, `equidistant_grid` will still return a floating point
-grid. This simlifies the implementation and avoids certain surprise
-behaviours.
+grid. This simplifies the implementation and avoids certain surprise
+behaviors.
 """
-function equidistant_grid(size::Dims, limit_lower, limit_upper)
-    gs = map(equidistant_grid, size, limit_lower, limit_upper)
+function equidistant_grid(limit_lower, limit_upper, dims::Vararg{Int})
+    if !(length(limit_lower) == length(limit_upper) == length(dims))
+        throw(ArgumentError("All arguments must be of the same length"))
+    end
+    gs = map(equidistant_grid, limit_lower, limit_upper, dims)
     return TensorGrid(gs...)
 end
 
 """
-    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
+    equidistant_grid(limit_lower::T, limit_upper::T, size::Int)
 
 Constructs a 1D equidistant grid.
 """
-function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+function equidistant_grid(limit_lower::T, limit_upper::T, size::Int) where T
     if any(size .<= 0)
         throw(DomainError("size must be postive"))
     end
--- a/src/Grids/grid.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/Grids/grid.jl	Sat May 25 16:07:10 2024 -0700
@@ -27,18 +27,60 @@
 """
     coordinate_size(g)
 
-The lenght of the coordinate vector of `Grid` `g`.
+The length of the coordinate vector of `Grid` `g`.
 """
 coordinate_size(::Type{<:Grid{T}}) where T = _ncomponents(T)
 coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?!
 
 """
-    component_type(g)
+    component_type(gf)
+
+The type of the components of the elements of `gf`. I.e if `gf` is a vector
+valued grid function, `component_view(gf)` is the element type of the vectors
+at each grid point.
+
+# Examples
+```julia-repl
+julia> component_type([[1,2], [2,3], [3,4]])
+Int64
+```
+"""
+component_type(T::Type) = eltype(eltype(T))
+component_type(t) = component_type(typeof(t))
+
+"""
+    componentview(gf, component_index...)
+
+A view of `gf` with only the components specified by `component_index...`.
 
-The type of the components of the coordinate vector of `Grid` `g`.
+# Examples
+```julia-repl
+julia> componentview([[1,2], [2,3], [3,4]],2)
+3-element ArrayComponentView{Int64, Vector{Int64}, 1, Vector{Vector{Int64}}, Tuple{Int64}}:
+ 2
+ 3
+ 4
+```
 """
-component_type(::Type{<:Grid{T}}) where T = eltype(T)
-component_type(g::Grid) = component_type(typeof(g))
+componentview(gf, component_index...) = ArrayComponentView(gf, component_index)
+
+struct ArrayComponentView{CT,T,D,AT <: AbstractArray{T,D}, IT} <: AbstractArray{CT,D}
+    v::AT
+    component_index::IT
+
+    function ArrayComponentView(v, component_index)
+        CT = typeof(first(v)[component_index...])
+        return new{CT, eltype(v), ndims(v), typeof(v), typeof(component_index)}(v,component_index)
+    end
+end
+
+Base.size(cv) = size(cv.v)
+Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...]
+Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...]
+IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)
+
+# TODO: Implement `setindex!`?
+# TODO: Implement a more general ComponentView that can handle non-AbstractArrays.
 
 """
     refine(g::Grid, r)
@@ -99,6 +141,8 @@
     if hasmethod(f, (Any,))
         return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g))
     else
+        # TBD This branch can be removed if we accept the trade off that we define f with the syntax f((x,y)) instead if we don't want to handle the vector in the body of f. (Add an example in the docs)
+        # Also see Notes.md
         return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g))
     end
 end
--- a/src/LazyTensors/lazy_array.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/LazyTensors/lazy_array.jl	Sat May 25 16:07:10 2024 -0700
@@ -1,7 +1,7 @@
 """
     LazyArray{T,D} <: AbstractArray{T,D}
 
-Array which is calcualted lazily when indexing.
+Array which is calculated lazily when indexing.
 
 A subtype of `LazyArray` will use lazy version of `+`, `-`, `*`, `/`.
 """
@@ -42,7 +42,7 @@
 
 """
     LazyElementwiseOperation{T,D,Op} <: LazyArray{T,D}
-Struct allowing for lazy evaluation of elementwise operations on `AbstractArray`s.
+Struct allowing for lazy evaluation of element-wise operations on `AbstractArray`s.
 
 A `LazyElementwiseOperation` contains two arrays together with an operation.
 The operations are carried out when the `LazyElementwiseOperation` is indexed.
--- a/src/LazyTensors/lazy_tensor_operations.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Sat May 25 16:07:10 2024 -0700
@@ -5,7 +5,7 @@
 
 Allows the result of a `LazyTensor` applied to a vector to be treated as an `AbstractArray`.
 With a mapping `m` and a vector `v` the TensorApplication object can be created by `m*v`.
-The actual result will be calcualted when indexing into `m*v`.
+The actual result will be calculated when indexing into `m*v`.
 """
 struct TensorApplication{T,R,D, TM<:LazyTensor{<:Any,R,D}, AA<:AbstractArray{<:Any,D}} <: LazyArray{T,R}
     t::TM
@@ -102,7 +102,7 @@
     TensorComposition(tm, tmi::IdentityTensor)
     TensorComposition(tmi::IdentityTensor, tm)
 
-Composes a `Tensormapping` `tm` with an `IdentityTensor` `tmi`, by returning `tm`
+Composes a `LazyTensor` `tm` with an `IdentityTensor` `tmi`, by returning `tm`
 """
 function TensorComposition(tm::LazyTensor{T,R,D}, tmi::IdentityTensor{T,D}) where {T,R,D}
     @boundscheck check_domain_size(tm, range_size(tmi))
@@ -126,7 +126,7 @@
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}
 
-An inflated `LazyTensor` with dimensions added before and afer its actual dimensions.
+An inflated `LazyTensor` with dimensions added before and after its actual dimensions.
 """
 struct InflatedTensor{T,R,D,D_before,R_middle,D_middle,D_after, TM<:LazyTensor{T,R_middle,D_middle}} <: LazyTensor{T,R,D}
     before::IdentityTensor{T,D_before}
@@ -220,7 +220,7 @@
 @doc raw"""
     LazyOuterProduct(tms...)
 
-Creates a `TensorComposition` for the outerproduct of `tms...`.
+Creates a `TensorComposition` for the outer product of `tms...`.
 This is done by separating the outer product into regular products of outer products involving only identity mappings and one non-identity mapping.
 
 First let
@@ -279,7 +279,7 @@
 `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
+dimensional difference operator `D` to a 2D grid of a certain size. In that
 case we could have
 
 ```julia
--- a/src/LazyTensors/tensor_types.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/LazyTensors/tensor_types.jl	Sat May 25 16:07:10 2024 -0700
@@ -1,7 +1,7 @@
 """
     IdentityTensor{T,D} <: LazyTensor{T,D,D}
 
-The lazy identity LazyTensor for a given size. Usefull for building up higher dimensional tensor mappings from lower
+The lazy identity LazyTensor for a given size. Useful for building up higher dimensional tensor mappings from lower
 dimensional ones through outer products. Also used in the Implementation for InflatedTensor.
 """
 struct IdentityTensor{T,D} <: LazyTensor{T,D,D}
@@ -57,8 +57,8 @@
 """
     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
-be considerd the range and domain of the LazyTensor. Each set of indices must be ordered in ascending order.
+LazyTensor defined by the AbstractArray A. `range_indicies` and `domain_indicies` define which indices of A should
+be considered the range and domain of the LazyTensor. Each set of indices must be ordered in ascending order.
 
 For instance, if A is a m x n matrix, and range_size = (1,), domain_size = (2,), then the DenseTensor performs the
 standard matrix-vector product on vectors of size n.
--- a/src/LazyTensors/tuple_manipulation.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/LazyTensors/tuple_manipulation.jl	Sat May 25 16:07:10 2024 -0700
@@ -3,7 +3,7 @@
 
 Splits the multi-index `I` into two parts. One part which is expected to be
 used as a view, and one which is expected to be used as an index.
-Eg.
+E.g.
 ```julia-repl
 julia> LazyTensors.split_index(1, 3, 2, 1, (1,2,3,4)...)
 ((1, Colon(), Colon(), Colon(), 4), (2, 3))
@@ -33,7 +33,7 @@
     split_tuple(t, szs)
 
 Split the tuple `t` into a set of tuples of the sizes given in `szs`.
-`sum(szs)` should equal `lenght(t)`.
+`sum(szs)` should equal `length(t)`.
 
 E.g
 ```julia-repl
--- a/src/SbpOperators/stencil.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/SbpOperators/stencil.jl	Sat May 25 16:07:10 2024 -0700
@@ -69,22 +69,29 @@
 end
 
 Base.@propagate_inbounds @inline function apply_stencil(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(promote_type(eltype(s),eltype(v)))
-    @simd for k ∈ 1:length(s)
-        w += s.weights[k]*v[i + s.range[k]]
+    return sum(enumerate(s.weights)) do (k,w) #TBD: Which optimizations are needed here?
+        w*v[i + @inbounds s.range[k]]
     end
-
-    return w
 end
 
 Base.@propagate_inbounds @inline function apply_stencil_backwards(s::Stencil, v::AbstractVector, i::Int)
-    w = zero(promote_type(eltype(s),eltype(v)))
-    @simd for k ∈ length(s):-1:1
-        w += s.weights[k]*v[i - s.range[k]]
+    return sum(enumerate(s.weights)) do (k,w) #TBD: Which optimizations are needed here?
+        w*v[i - @inbounds s.range[k]]
     end
-    return w
 end
 
+# There are many options for the implementation of `apply_stencil` and
+# `apply_stencil_backwards`. Some alternatives were tried on the branch
+# bugfix/sbp_operators/stencil_return_type and can be found at the following
+# revision:
+#
+# * 237b980ffb91 (baseline)
+# * a72bab15228e (mapreduce)
+# * ffd735354d54 (multiplication)
+# * b5abd5191f2c (promote_op)
+# * 8d56846185fc (return_type)
+#
+
 function left_pad(s::Stencil, N)
     weights = LazyTensors.left_pad_tuple(s.weights, zero(eltype(s)), N)
     range = (first(s.range) - (N - length(s.weights))):last(s.range)
--- a/src/SbpOperators/stencil_set.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/src/SbpOperators/stencil_set.jl	Sat May 25 16:07:10 2024 -0700
@@ -5,7 +5,7 @@
     StencilSet
 
 A `StencilSet` contains a set of associated stencils. The stencils
-are are stored in a table, and can be accesed by indexing into the `StencilSet`.
+are are stored in a table, and can be accessed by indexing into the `StencilSet`.
 """
 struct StencilSet
     table
@@ -21,7 +21,7 @@
 table of the `StencilSet` is a parsed TOML intended for functions like
 `parse_scalar` and `parse_stencil`.
 
-The `StencilSet` table is not parsed beyond the inital TOML parse. To get usable
+The `StencilSet` table is not parsed beyond the initial TOML parse. To get usable
 stencils use the `parse_stencil` functions on the fields of the stencil set.
 
 The reason for this is that since stencil sets are intended to be very
--- a/test/Grids/equidistant_grid_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -105,34 +105,36 @@
 
 
 @testset "equidistant_grid" begin
-    @test equidistant_grid(4,0.0,1.0) isa EquidistantGrid
-    @test equidistant_grid((4,3),(0.0,0.0),(8.0,5.0)) isa TensorGrid
+    @test equidistant_grid(0.0,1.0, 4) isa EquidistantGrid
+    @test equidistant_grid((0.0,0.0),(8.0,5.0), 4, 3) isa TensorGrid
 
     # constuctor
-    @test_throws DomainError equidistant_grid(0,0.0,1.0)
-    @test_throws DomainError equidistant_grid(1,1.0,1.0)
-    @test_throws DomainError equidistant_grid(1,1.0,-1.0)
+    @test_throws DomainError equidistant_grid(0.0, 1.0, 0)
+    @test_throws DomainError equidistant_grid(1.0, 1.0, 1)
+    @test_throws DomainError equidistant_grid(1.0, -1.0, 1)
 
-    @test_throws DomainError equidistant_grid((0,0),(0.0,0.0),(1.0,1.0))
-    @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(1.0,1.0))
-    @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(-1.0,-1.0))
+    @test_throws DomainError equidistant_grid((0.0,0.0),(1.0,1.0), 0, 0)
+    @test_throws DomainError equidistant_grid((1.0,1.0),(1.0,1.0), 1, 1)
+    @test_throws DomainError equidistant_grid((1.0,1.0),(-1.0,-1.0), 1, 1)
+
+    @test_throws ArgumentError equidistant_grid((0.0,),(8.0,5.0), 4, 3, 4)
 
     @testset "Base" begin
-        @test eltype(equidistant_grid(4,0.0,1.0)) == Float64
-        @test eltype(equidistant_grid((4,3),(0,0),(1,3))) <: AbstractVector{Float64}
+        @test eltype(equidistant_grid(0.0, 1.0, 4)) == Float64
+        @test eltype(equidistant_grid((0,0),(1,3), 4, 3)) <: AbstractVector{Float64}
 
-        @test size(equidistant_grid(4,0.0,1.0)) == (4,)
-        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+        @test size(equidistant_grid(0.0, 1.0, 4)) == (4,)
+        @test size(equidistant_grid((0.0,0.0), (2.0,1.0), 5, 3)) == (5,3)
 
-        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0)),1) == 5
-        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0)),2) == 3
+        @test size(equidistant_grid((0.0,0.0), (2.0,1.0), 5, 3), 1) == 5
+        @test size(equidistant_grid((0.0,0.0), (2.0,1.0), 5, 3), 2) == 3
 
-        @test ndims(equidistant_grid(4,0.0,1.0)) == 1
-        @test ndims(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == 2
+        @test ndims(equidistant_grid(0.0, 1.0, 4)) == 1
+        @test ndims(equidistant_grid((0.0,0.0), (2.0,1.0), 5, 3)) == 2
     end
 
     @testset "getindex" begin
-        g = equidistant_grid((5,3), (-1.0,0.0), (0.0,7.11))
+        g = equidistant_grid((-1.0,0.0), (0.0,7.11), 5, 3)
         gp = collect(g);
         p = [(-1.,0.)      (-1.,7.11/2)   (-1.,7.11);
             (-0.75,0.)    (-0.75,7.11/2) (-0.75,7.11);
--- a/test/Grids/grid_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/Grids/grid_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -15,19 +15,21 @@
     @test coordinate_size(DummyGrid{SVector{3,Float64}, 2}()) == 3
 
     @test coordinate_size(DummyGrid{SVector{3,Float64}, 2}) == 3
+end
 
-    @testset "component_type" begin
-        @test component_type(DummyGrid{Int,1}()) == Int
-        @test component_type(DummyGrid{Float64,1}()) == Float64
-        @test component_type(DummyGrid{Rational,1}()) == Rational
+@testset "component_type" begin
+    @test component_type(DummyGrid{Int,1}()) == Int
+    @test component_type(DummyGrid{Float64,1}()) == Float64
+    @test component_type(DummyGrid{Rational,1}()) == Rational
 
-        @test component_type(DummyGrid{SVector{3,Int},2}()) == Int
-        @test component_type(DummyGrid{SVector{2,Float64},3}()) == Float64
-        @test component_type(DummyGrid{SVector{4,Rational},4}()) == Rational
+    @test component_type(DummyGrid{SVector{3,Int},2}()) == Int
+    @test component_type(DummyGrid{SVector{2,Float64},3}()) == Float64
+    @test component_type(DummyGrid{SVector{4,Rational},4}()) == Rational
 
-        @test component_type(DummyGrid{Float64,1}) == Float64
-        @test component_type(DummyGrid{SVector{2,Float64},3}) == Float64
-    end
+    @test component_type(DummyGrid{Float64,1}) == Float64
+    @test component_type(DummyGrid{SVector{2,Float64},3}) == Float64
+
+    @test component_type(fill(@SVector[1,2], 4,2)) == Int
 end
 
 @testset "eval_on" begin
@@ -45,7 +47,7 @@
     @test eval_on(EquidistantGrid(range(0,1,length=4)), x->2x) == 2 .* range(0,1,length=4)
 
 
-    g = equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))
+    g = equidistant_grid((0.0,0.0), (2.0,1.0), 5, 3)
 
     @test eval_on(g, x̄ -> 0.) isa LazyArray
     @test eval_on(g, x̄ -> 0.) == fill(0., (5,3))
@@ -63,6 +65,57 @@
     @test eval_on(g, f) == map(x̄->f(x̄...), g)
 end
 
+@testset "componentview" begin
+    v = [@SMatrix[1 3; 2 4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+    @test componentview(v, 1, 1) isa AbstractArray
+    @test componentview(v, 1, :) isa AbstractArray
+
+    A = @SMatrix[
+            1 4 7;
+            2 5 8;
+            3 6 9;
+        ]
+    v = [A .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 2:3, 1:2) isa AbstractArray
+
+    # Correctness of the result is tested in ArrayComponentView
+end
+
+@testset "ArrayComponentView" begin
+    v = [@SMatrix[1 3; 2 4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+    @testset "==" begin
+        @test ArrayComponentView(v, (1,1)) == ArrayComponentView(v, (1,1))
+        @test ArrayComponentView(v, (1,1)) == ArrayComponentView(copy(v), (1,1))
+        @test ArrayComponentView(v, (1,1)) == [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4] == ArrayComponentView(v, (1,1))
+    end
+
+    @testset "components" begin
+        v = [@SMatrix[1 3; 2 4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+        @test ArrayComponentView(v, (1, 1))  == [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (1, 2))  == [3 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (2, 1))  == [2 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+        @test ArrayComponentView(v, (1, :))  == [@SVector[1,3] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (2, :))  == [@SVector[2,4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (:, 1))  == [@SVector[1,2] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (:, 2))  == [@SVector[3,4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+
+        A = @SMatrix[
+            1 4 7;
+            2 5 8;
+            3 6 9;
+        ]
+        v = [A .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (1:2, 1:2)) == [@SMatrix[1 4;2 5] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+        @test ArrayComponentView(v, (2:3, 1:2)) == [@SMatrix[2 5;3 6] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    end
+end
+
 @testset "_ncomponents" begin
     @test Grids._ncomponents(Int) == 1
     @test Grids._ncomponents(Float64) == 1
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -9,8 +9,8 @@
 @testset "boundary_restriction" begin
 	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
 	e_closure = parse_stencil(stencil_set["e"]["closure"])
-    g_1D = equidistant_grid(11, 0.0, 1.0)
-    g_2D = equidistant_grid((11,15), (0.0, 0.0), (1.0,1.0))
+    g_1D = equidistant_grid(0.0, 1.0, 11)
+    g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15)
 
     @testset "boundary_restriction" begin
         @testset "1D" begin
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -7,8 +7,8 @@
 import Sbplib.SbpOperators.BoundaryOperator
 
 @testset "normal_derivative" begin
-    g_1D = equidistant_grid(11, 0.0, 1.0)
-    g_2D = equidistant_grid((11,12), (0.0, 0.0), (1.0,1.0))
+    g_1D = equidistant_grid(0.0, 1.0, 11)
+    g_2D = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 12)
     @testset "normal_derivative" begin
     	stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "1D" begin
--- a/test/SbpOperators/stencil_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/stencil_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -1,8 +1,10 @@
 using Test
 using Sbplib.SbpOperators
+using StaticArrays
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.NestedStencil
 import Sbplib.SbpOperators.scale
+import Sbplib.SbpOperators: apply_stencil, apply_stencil_backwards
 
 @testset "Stencil" begin
     s = Stencil(-2:2, (1.,2.,2.,3.,4.))
@@ -44,6 +46,52 @@
         @test promote(Stencil(1,1;center=1), Stencil(2.,2.;center=2)) == (Stencil(1.,1.;center=1), Stencil(2.,2.;center=2))
     end
 
+    @testset "apply_stencil" begin
+        v = [1, 2, 4, 8, 16, 32, 64, 128]
+        s = Stencil(1,2,3,4, center = 2)
+        @test apply_stencil(s,v, 2) == v[1] + 2*v[2] + 3*v[3] + 4*v[4]
+        @test apply_stencil(s,v, 4) == v[3] + 2*v[4] + 3*v[5] + 4*v[6]
+        @test apply_stencil_backwards(s,v, 3) == 4*v[1] + 3*v[2] + 2*v[3] + 1*v[4]
+        @test apply_stencil_backwards(s,v, 7) == 4*v[5] + 3*v[6] + 2*v[7] + 1*v[8]
+        @test apply_stencil(s,v, 2) isa Int
+        @test apply_stencil_backwards(s,v, 7) isa Int
+
+        v = [1, 2, 4, 8, 16, 32, 64, 128]
+        s = Stencil(1.,2.,3.,4., center = 2)
+        @test apply_stencil(s,v, 4) == v[3] + 2. *v[4] + 3. *v[5] + 4. *v[6]
+        @test apply_stencil_backwards(s,v, 7) == 4. *v[5] + 3. *v[6] + 2. *v[7] + v[8]
+        @test apply_stencil(s,v, 2) isa Float64
+        @test apply_stencil_backwards(s,v, 7) isa Float64
+
+        v = [1., 2., 4., 8., 16., 32., 64., 128.]
+        s = Stencil(1,2,3,4, center = 2)
+        @test apply_stencil(s,v, 2) == v[1] + 2*v[2] + 3*v[3] + 4*v[4]
+        @test apply_stencil_backwards(s,v, 3) == 4*v[1] + 3*v[2] + 2*v[3] + 1*v[4]
+        @test apply_stencil(s,v, 2) isa Float64
+        @test apply_stencil_backwards(s,v, 3) isa Float64
+
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        s = Stencil(1,2, center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7, 10]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11, 14]
+        @test apply_stencil(s,v,1) isa SVector{2, Int}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Int}
+
+        v = [@SVector[1., 2.], @SVector[3., 4.], @SVector[5., 6.], @SVector[7., 8.]]
+        s = Stencil(1,2, center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7., 10.]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11., 14.]
+        @test apply_stencil(s,v,1) isa SVector{2, Float64}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Float64}
+
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        s = Stencil(1.,2., center = 1)
+        @test apply_stencil(s,v,1) == @SVector[7., 10.]
+        @test apply_stencil_backwards(s,v,3) == @SVector[11., 14.]
+        @test apply_stencil(s,v,1) isa SVector{2, Float64}
+        @test apply_stencil_backwards(s,v,3) isa SVector{2, Float64}
+    end
+
     @testset "type stability" begin
         s_int = CenteredStencil(1,2,3)
         s_float = CenteredStencil(1.,2.,3.)
@@ -155,6 +203,27 @@
 
         @test SbpOperators.apply_stencil(ns, c, v, 4) == 5*7 + 11*11 + 6*13
         @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -3*3 - 7*5 - 4*7
+
+        # Different types in vector and stencil
+        ns = NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-2.,2.), center=2)
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) isa Stencil{Float64, 3}
+        @test SbpOperators.apply_inner_stencils(ns, c, 4) == Stencil(4.,9.,10.; center=2)
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) isa Stencil{Float64, 3}
+        @test SbpOperators.apply_inner_stencils_backwards(ns, c, 4) == Stencil(-5.,-9.,-8.; center=2)
+
+        @test SbpOperators.apply_stencil(ns, c, v, 4) isa Float64
+        @test SbpOperators.apply_stencil(ns, c, v, 4) == 193.
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) isa Float64
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 4) == -158.
+
+        # Arrays of vectors
+        ns = NestedStencil((-1.,1.,0.),(-1.,0.,1.),(0.,-2.,2.), center=2)
+        c = [  1,  3,  6, 10]
+        v = [@SVector[1, 2], @SVector[3, 4], @SVector[5, 6], @SVector[7, 8]]
+        @test SbpOperators.apply_stencil(ns, c, v, 2) isa SVector{2,Float64}
+        @test SbpOperators.apply_stencil(ns, c, v, 2) == 2v[1] + 5v[2] + 6v[3]
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 2) isa SVector{2,Float64}
+        @test SbpOperators.apply_stencil_backwards(ns, c, v, 2) == -4v[1] - 5v[2] - 3v[3]
     end
 
     @testset "type stability" begin
--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -33,7 +33,7 @@
     @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3)
 
     @testset "Grid constructor" begin
-        g = equidistant_grid(11, 0., 2.)
+        g = equidistant_grid(0., 2., 11)
         @test ConstantInteriorScalingOperator(g, 3., (.1,.2)) isa ConstantInteriorScalingOperator{Float64}
     end
 end
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -27,7 +27,7 @@
 end
 
 @testset "undivided_skewed04" begin
-    g = equidistant_grid(20, 0., 11.)
+    g = equidistant_grid(0., 11., 20)
     D,Dᵀ = undivided_skewed04(g, 1)
 
     @test D isa LazyTensor{Float64,1,1}
@@ -35,7 +35,7 @@
 
      @testset "Accuracy conditions" begin
         N = 20
-        g = equidistant_grid(N, 0//1,2//1)
+        g = equidistant_grid(0//1, 2//1, N)
         h = only(spacing(g))
         @testset "D_$p" for p ∈ [1,2,3,4]
             D,Dᵀ = undivided_skewed04(g, p)
@@ -67,7 +67,7 @@
             return Dmat
         end
 
-        g = equidistant_grid(11, 0., 1.)
+        g = equidistant_grid(0., 1., 11)
         @testset "D_$p" for p ∈ [1,2,3,4]
             D,Dᵀ = undivided_skewed04(g, p)
 
@@ -80,7 +80,7 @@
 
     @testset "2D" begin
         N = 20
-        g = equidistant_grid((N,2N), (0,0), (2,1))
+        g = equidistant_grid((0,0), (2,1), N, 2N)
         h = spacing.(g.grids)
 
         D,Dᵀ = undivided_skewed04(g, 3, 2)
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -24,8 +24,8 @@
     @testset "Constructors" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
 
-        g₁ = equidistant_grid(11, 0., 1.)
-        g₂ = equidistant_grid((11,14), (0.,1.), (1.,3.))
+        g₁ = equidistant_grid(0., 1., 11)
+        g₂ = equidistant_grid((0.,1.), (1.,3.), 11, 14)
         
         @test first_derivative(g₁, stencil_set) isa LazyTensor{Float64,1,1}
         @test first_derivative(g₂, stencil_set, 2) isa LazyTensor{Float64,2,2}
@@ -38,7 +38,7 @@
 
     @testset "Accuracy conditions" begin
         N = 20
-        g = equidistant_grid(N, 0//1,2//1)
+        g = equidistant_grid(0//1, 2//1, N)
         @testset for order ∈ [2,4]
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
             D₁ = first_derivative(g, stencil_set)
@@ -68,7 +68,7 @@
 
     @testset "Accuracy on function" begin
         @testset "1D" begin
-            g = equidistant_grid(30, 0.,1.)
+            g = equidistant_grid(0., 1., 30)
             v = eval_on(g, x->exp(x))
             @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)]
                 stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
@@ -79,7 +79,7 @@
         end
 
         @testset "2D" begin
-            g = equidistant_grid((30,60), (0.,0.),(1.,2.))
+            g = equidistant_grid((0.,0.),(1.,2.), 30, 60)
             v = eval_on(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)
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -15,8 +15,8 @@
     closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
     Lx = 3.5
     Ly = 3.
-    g_1D = equidistant_grid(121, 0.0, Lx)
-    g_2D = equidistant_grid((121,123), (0.0, 0.0), (Lx, Ly))
+    g_1D = equidistant_grid(0.0, Lx, 121)
+    g_2D = equidistant_grid((0.0, 0.0), (Lx, Ly), 121, 123)
 
     @testset "Constructors" begin
         @testset "1D" begin
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -12,7 +12,7 @@
     stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
 
     @testset "1D" begin
-        g = equidistant_grid(11, 0., 1.)
+        g = equidistant_grid(0., 1., 11)
         c = [  1.,  3.,  6., 10., 15., 21., 28., 36., 45., 55., 66.]
 
         @testset "checking c" begin
@@ -27,7 +27,7 @@
 
         @testset "application" begin
             function apply_to_functions(; v, c)
-                g = equidistant_grid(11, 0., 10.) # h = 1
+                g = equidistant_grid(0., 10., 11) # h = 1
                 c̄ = eval_on(g,c)
                 v̄ = eval_on(g,v)
 
@@ -44,12 +44,12 @@
     end
 
     @testset "2D" begin
-        g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1
+        g = equidistant_grid((0.,0.), (10.,8.), 11, 9) # h = 1
         c = eval_on(g, (x,y)->x+y)
 
         @testset "application" begin
             function apply_to_functions(dir; v, c)
-                g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1
+                g = equidistant_grid((0.,0.), (10.,8.), 11, 9) # h = 1
                 c̄ = eval_on(g,c)
                 v̄ = eval_on(g,v)
 
@@ -89,7 +89,7 @@
                 Dxv(x,y) = cos(x)*exp(x) - (exp(x) + exp(1.5 - 1.5y))*sin(x)
                 Dyv(x,y) = -1.5(1.5exp(x) + 1.5exp(1.5 - 1.5y))*cos(1.5 - 1.5y) - 2.25exp(1.5 - 1.5y)*sin(1.5 - 1.5y)
 
-                g₁ = equidistant_grid((60,67), (0.,0.), (1.,2.))
+                g₁ = equidistant_grid((0.,0.), (1.,2.), 60, 67)
                 g₂ = refine(g₁,2)
 
                 c̄₁ = eval_on(g₁, c)
@@ -155,7 +155,7 @@
         @testset "application" begin
 
             function apply_to_functions(; v, c)
-                g = equidistant_grid(11, 0., 10.) # h = 1
+                g = equidistant_grid(0., 10., 11) # h = 1
                 c̄ = eval_on(g,c)
                 v̄ = eval_on(g,v)
 
@@ -171,7 +171,7 @@
         end
 
         @testset "type stability" begin
-            g = equidistant_grid(11, 0., 10.) # h = 1
+            g = equidistant_grid(0., 10., 11) # h = 1
             c̄ = eval_on(g,x-> -1)
             v̄ = eval_on(g,x->1.)
 
@@ -185,7 +185,7 @@
     end
 
     @testset "2D" begin
-        g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1
+        g = equidistant_grid((0.,0.), (10.,8.), 11, 9) # h = 1
         c = eval_on(g, (x,y)->x+y)
         @testset "Constructors" begin
             @test SecondDerivativeVariable(c, interior_stencil, closure_stencils, 1) isa LazyTensor
@@ -210,7 +210,7 @@
 
         @testset "application" begin
             function apply_to_functions(dir; v, c)
-                g = equidistant_grid((11,9), (0.,0.), (10.,8.)) # h = 1
+                g = equidistant_grid((0.,0.), (10.,8.), 11, 9) # h = 1
                 c̄ = eval_on(g,c)
                 v̄ = eval_on(g,v)
 
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -10,9 +10,9 @@
     Lx = π/2.
     Ly = Float64(π)
     Lz = 1.
-    g_1D = equidistant_grid(77, 0.0, Lx)
-    g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly))
-    g_3D = equidistant_grid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
+    g_1D = equidistant_grid(0.0, Lx, 77)
+    g_2D = equidistant_grid((0.0, 0.0), (Lx,Ly), 77, 66)
+    g_3D = equidistant_grid((0.0, 0.0, 0.0), (Lx,Ly,Lz), 10, 10, 10)
     @testset "inner_product" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "0D" begin
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -9,8 +9,8 @@
 @testset "Diagonal-stencil inverse_inner_product" begin
     Lx = π/2.
     Ly = Float64(π)
-    g_1D = equidistant_grid(77, 0.0, Lx)
-    g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly))
+    g_1D = equidistant_grid(0.0, Lx, 77)
+    g_2D = equidistant_grid((0.0, 0.0), (Lx,Ly), 77, 66)
     @testset "inverse_inner_product" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4)
         @testset "0D" begin
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -9,8 +9,8 @@
     # Default stencils (4th order)
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)
-    g_1D = equidistant_grid(101, 0.0, 1.)
-    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+    g_1D = equidistant_grid(0.0, 1., 101)
+    g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
     @testset "Constructors" begin
         @testset "1D" begin
@@ -70,8 +70,8 @@
 @testset "laplace" begin
     operator_path = sbp_operators_path()*"standard_diagonal.toml"
     stencil_set = read_stencil_set(operator_path; order=4)
-    g_1D = equidistant_grid(101, 0.0, 1.)
-    g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
+    g_1D = equidistant_grid(0.0, 1., 101)
+    g_3D = equidistant_grid((0.0, -1.0, 0.0), (1., 1., 1.), 51, 101, 52)
 
     @testset "1D" begin
         Δ = laplace(g_1D, stencil_set)
--- a/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -8,7 +8,7 @@
 import Sbplib.SbpOperators.StencilOperatorDistinctClosures
 
 @testset "StencilOperatorDistinctClosures" begin
-    g = equidistant_grid(11, 0., 1.)
+    g = equidistant_grid(0., 1., 11)
 
     lower_closure = (
         Stencil(-1,1, center=1),
--- a/test/SbpOperators/volumeops/volume_operator_test.jl	Tue Jan 23 20:48:25 2024 +0100
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Sat May 25 16:07:10 2024 -0700
@@ -14,7 +14,7 @@
 @testset "VolumeOperator" begin
     inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
     closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2))
-    g = equidistant_grid(11,0.,1.)
+    g = equidistant_grid(0.,1., 11)
 
     @testset "Constructors" begin
         op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)