Mercurial > repos > public > sbplib_julia

changeset 1557:9113f437431d feature/grids/curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 13 Apr 2024 23:52:40 +0200
parents 69790e9d1652 (current diff) ec5e7926c37b (diff)
children 43e6acbefdd1 5d32ecb98db8
files src/Grids/equidistant_grid.jl src/Grids/grid.jl
diffstat 23 files changed, 225 insertions(+), 109 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Tue Apr 09 15:26:49 2024 +0200
+++ b/Notes.md	Sat Apr 13 23:52:40 2024 +0200
@@ -186,18 +186,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
@@ -325,3 +348,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 Apr 09 15:26:49 2024 +0200
+++ b/benchmark/Manifest.toml	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/benchmark/benchmark_laplace.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/benchmark/benchmark_utils.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/benchmark/benchmarks.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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/equidistant_grid.jl	Tue Apr 09 15:26:49 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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,25 +99,24 @@
 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 simplifies the implementation and avoids certain surprise
 behaviors.
 """
-# TODO: Change signature to `equidistant_grid(limit_lower,limit_upper, size...)
-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})
+    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 Apr 09 15:26:49 2024 +0200
+++ b/src/Grids/grid.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -142,6 +142,7 @@
         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/SbpOperators/stencil.jl	Tue Apr 09 15:26:49 2024 +0200
+++ b/src/SbpOperators/stencil.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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/test/Grids/equidistant_grid_test.jl	Tue Apr 09 15:26:49 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -105,34 +105,34 @@
 
 
 @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)
 
     @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 Apr 09 15:26:49 2024 +0200
+++ b/test/Grids/grid_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -47,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))
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Tue Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/stencil_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/second_derivative_variable_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -8,8 +8,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
@@ -69,8 +69,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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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 Apr 09 15:26:49 2024 +0200
+++ b/test/SbpOperators/volumeops/volume_operator_test.jl	Sat Apr 13 23:52:40 2024 +0200
@@ -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)