--- a/Notes.md	Mon Sep 09 22:20:51 2024 +0200
+++ b/Notes.md	Tue Sep 10 20:42:41 2024 +0200
@@ -392,3 +392,38 @@
 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.
+
+
+## Tiled loops and regions in apply
+There should be easy ways to use functionalty splitting the application of a lazy array into regions and using tiled iteration. This could make the application more efficient by reducing branching and improving cache usage in the tight loop. On commit f215ac2a5c66 and before there were some early tests regarding this in a DiffOp submodule.
+
+The main ideas were:
+```julia
+function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
+    apply_region!(D, u, v, Lower, Lower)
+    apply_region!(D, u, v, Lower, Interior)
+    apply_region!(D, u, v, Lower, Upper)
+    apply_region!(D, u, v, Interior, Lower)
+    apply_region!(D, u, v, Interior, Interior)
+    apply_region!(D, u, v, Interior, Upper)
+    apply_region!(D, u, v, Upper, Lower)
+    apply_region!(D, u, v, Upper, Interior)
+    apply_region!(D, u, v, Upper, Upper)
+    return nothing
+end
+```
+
+```julia
+using TiledIteration
+function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
+    ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
+    # TODO: Pass Tilesize to function
+    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
+        for j ∈ tileaxs[2], i ∈ tileaxs[1]
+            I = ri[i,j]
+            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
+        end
+    end
+    return nothing
+end
+```

--- a/TODO.md	Mon Sep 09 22:20:51 2024 +0200
+++ b/TODO.md	Tue Sep 10 20:42:41 2024 +0200
@@ -5,7 +5,6 @@
 
 ## Coding
  - [ ] Ändra namn på variabler och funktioner så att det följer style-guide
- - [ ] Add new Laplace operator to DiffOps, probably named WaveEqOp(?!!?)
  - [ ] Create a struct that bundles the necessary Tensor operators for solving the wave equation.
  - [ ] Use `@inferred` in a lot of tests.
  - [ ] Replace `@inferred` tests with a benchmark suite that automatically tests for regressions.

--- a/docs/make.jl	Mon Sep 09 22:20:51 2024 +0200
+++ b/docs/make.jl	Tue Sep 10 20:42:41 2024 +0200
@@ -1,7 +1,6 @@
 using Documenter
 using Diffinitive
 
-using Diffinitive.DiffOps
 using Diffinitive.Grids
 using Diffinitive.LazyTensors
 using Diffinitive.RegionIndices
@@ -33,7 +32,6 @@
     "matrix_and_tensor_representations.md",
     "Submodules" => [
         "submodules/grids.md",
-        "submodules/diff_ops.md",
         "submodules/lazy_tensors.md",
         "submodules/region_indices.md",
         "submodules/sbp_operators.md",

--- a/docs/src/submodules/diff_ops.md	Mon Sep 09 22:20:51 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-# DiffOps
-
-## Contents
-```@contents
-Pages = ["diff_ops.md"]
-```
-
-## Index
-```@index
-Pages = ["diff_ops.md"]
-```
-
-## Public interface
-```@autodocs
-Modules = [DiffOps]
-Private = false # Hide unexported objects
-```
-
-## Internal interface
-```@autodocs
-Modules = [DiffOps]
-Public = false # Hide exported objects
-```

--- a/src/DiffOps/DiffOps.jl	Mon Sep 09 22:20:51 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,103 +0,0 @@
-module DiffOps
-
-using Diffinitive.RegionIndices
-using Diffinitive.SbpOperators
-using Diffinitive.Grids
-using Diffinitive.LazyTensors
-
-"""
-    DiffOp
-
-Supertype of differential operator discretisations.
-The action of the DiffOp is defined in the method
-    apply(D::DiffOp, v::AbstractVector, I...)
-"""
-abstract type DiffOp end
-
-function apply end
-
-function matrixRepresentation(D::DiffOp)
-    error("not implemented")
-end
-
-abstract type DiffOpCartesian{Dim} <: DiffOp end
-
-# DiffOp must have a grid of dimension Dim!!!
-function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
-    for I ∈ eachindex(D.grid)
-        u[I] = apply(D, v, I)
-    end
-
-    return nothing
-end
-export apply!
-
-function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region!(D, u, v, Lower, Lower)
-    apply_region!(D, u, v, Lower, Interior)
-    apply_region!(D, u, v, Lower, Upper)
-    apply_region!(D, u, v, Interior, Lower)
-    apply_region!(D, u, v, Interior, Interior)
-    apply_region!(D, u, v, Interior, Upper)
-    apply_region!(D, u, v, Upper, Lower)
-    apply_region!(D, u, v, Upper, Interior)
-    apply_region!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-# Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
-function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    for I ∈ regionindices(D.grid.size, closuresize(D.op), (r1,r2))
-        @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
-        @inbounds u[I] = apply(D, v, indextuple)
-    end
-    return nothing
-end
-export apply_region!
-
-function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
-    apply_region_tiled!(D, u, v, Lower, Lower)
-    apply_region_tiled!(D, u, v, Lower, Interior)
-    apply_region_tiled!(D, u, v, Lower, Upper)
-    apply_region_tiled!(D, u, v, Interior, Lower)
-    apply_region_tiled!(D, u, v, Interior, Interior)
-    apply_region_tiled!(D, u, v, Interior, Upper)
-    apply_region_tiled!(D, u, v, Upper, Lower)
-    apply_region_tiled!(D, u, v, Upper, Interior)
-    apply_region_tiled!(D, u, v, Upper, Upper)
-    return nothing
-end
-
-using TiledIteration
-function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
-    ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
-    # TODO: Pass Tilesize to function
-    for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
-        for j ∈ tileaxs[2], i ∈ tileaxs[1]
-            I = ri[i,j]
-            u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
-        end
-    end
-    return nothing
-end
-export apply_region_tiled!
-
-function apply(D::DiffOp, v::AbstractVector)::AbstractVector
-    u = zeros(eltype(v), size(v))
-    apply!(D,v,u)
-    return u
-end
-
-# TODO: This conflicts with LazyTensors. Shouldn't DiffOps be LazyTensorOperators and use that apply?
-# export apply
-
-
-"""
-    BoundaryCondition
-A BoundaryCondition should implement the method
-    sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
-"""
-abstract type BoundaryCondition end
-
-
-end # module

--- a/src/Diffinitive.jl	Mon Sep 09 22:20:51 2024 +0200
+++ b/src/Diffinitive.jl	Tue Sep 10 20:42:41 2024 +0200
@@ -4,12 +4,10 @@
 include("LazyTensors/LazyTensors.jl")
 include("Grids/Grids.jl")
 include("SbpOperators/SbpOperators.jl")
-include("DiffOps/DiffOps.jl")
 
 export RegionIndices
 export LazyTensors
 export Grids
 export SbpOperators
-export DiffOps
 
 end

--- a/test/DiffOps/DiffOps_test.jl	Mon Sep 09 22:20:51 2024 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,195 +0,0 @@
-using Test
-using Diffinitive.DiffOps
-using Diffinitive.Grids
-using Diffinitive.SbpOperators
-using Diffinitive.RegionIndices
-using Diffinitive.LazyTensors
-
-#
-# @testset "BoundaryValue" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
-#
-#     e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
-#     e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
-#     e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
-#     e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = zeros(Float64, 4, 5)
-#     v[:,5] = [1, 2, 3,4]
-#     v[:,4] = [1, 2, 3,4]
-#     v[:,3] = [4, 5, 6, 7]
-#     v[:,2] = [7, 8, 9, 10]
-#     v[:,1] = [10, 11, 12, 13]
-#
-#     @test e_w  isa LazyTensor{T,2,1} where T
-#     @test e_w' isa LazyTensor{T,1,2} where T
-#
-#     @test domain_size(e_w, (3,2)) == (2,)
-#     @test domain_size(e_e, (3,2)) == (2,)
-#     @test domain_size(e_s, (3,2)) == (3,)
-#     @test domain_size(e_n, (3,2)) == (3,)
-#
-#     @test size(e_w'*v) == (5,)
-#     @test size(e_e'*v) == (5,)
-#     @test size(e_s'*v) == (4,)
-#     @test size(e_n'*v) == (4,)
-#
-#     @test collect(e_w'*v) == [10,7,4,1.0,1]
-#     @test collect(e_e'*v) == [13,10,7,4,4.0]
-#     @test collect(e_s'*v) == [10,11,12,13.0]
-#     @test collect(e_n'*v) == [1,2,3,4.0]
-#
-#     g_x = [1,2,3,4.0]
-#     g_y = [5,4,3,2,1.0]
-#
-#     G_w = zeros(Float64, (4,5))
-#     G_w[1,:] = g_y
-#
-#     G_e = zeros(Float64, (4,5))
-#     G_e[4,:] = g_y
-#
-#     G_s = zeros(Float64, (4,5))
-#     G_s[:,1] = g_x
-#
-#     G_n = zeros(Float64, (4,5))
-#     G_n[:,5] = g_x
-#
-#     @test size(e_w*g_y) == (UnknownDim,5)
-#     @test size(e_e*g_y) == (UnknownDim,5)
-#     @test size(e_s*g_x) == (4,UnknownDim)
-#     @test size(e_n*g_x) == (4,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken collect(e_w*g_y) == G_w
-#     @test_broken collect(e_e*g_y) == G_e
-#     @test_broken collect(e_s*g_x) == G_s
-#     @test_broken collect(e_n*g_x) == G_n
-# end
-#
-# @testset "NormalDerivative" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
-#
-#     d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
-#     d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
-#     d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
-#     d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
-#
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#     v∂x = evalOn(g, (x,y)-> 2*x + y)
-#     v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
-#
-#     @test d_w  isa LazyTensor{T,2,1} where T
-#     @test d_w' isa LazyTensor{T,1,2} where T
-#
-#     @test domain_size(d_w, (3,2)) == (2,)
-#     @test domain_size(d_e, (3,2)) == (2,)
-#     @test domain_size(d_s, (3,2)) == (3,)
-#     @test domain_size(d_n, (3,2)) == (3,)
-#
-#     @test size(d_w'*v) == (6,)
-#     @test size(d_e'*v) == (6,)
-#     @test size(d_s'*v) == (5,)
-#     @test size(d_n'*v) == (5,)
-#
-#     @test collect(d_w'*v) ≈ v∂x[1,:]
-#     @test collect(d_e'*v) ≈ v∂x[5,:]
-#     @test collect(d_s'*v) ≈ v∂y[:,1]
-#     @test collect(d_n'*v) ≈ v∂y[:,6]
-#
-#
-#     d_x_l = zeros(Float64, 5)
-#     d_x_u = zeros(Float64, 5)
-#     for i ∈ eachindex(d_x_l)
-#         d_x_l[i] = op.dClosure[i-1]
-#         d_x_u[i] = -op.dClosure[length(d_x_u)-i]
-#     end
-#
-#     d_y_l = zeros(Float64, 6)
-#     d_y_u = zeros(Float64, 6)
-#     for i ∈ eachindex(d_y_l)
-#         d_y_l[i] = op.dClosure[i-1]
-#         d_y_u[i] = -op.dClosure[length(d_y_u)-i]
-#     end
-#
-#     function prod_matrix(x,y)
-#         G = zeros(Float64, length(x), length(y))
-#         for I ∈ CartesianIndices(G)
-#             G[I] = x[I[1]]*y[I[2]]
-#         end
-#
-#         return G
-#     end
-#
-#     g_x = [1,2,3,4.0,5]
-#     g_y = [5,4,3,2,1.0,11]
-#
-#     G_w = prod_matrix(d_x_l, g_y)
-#     G_e = prod_matrix(d_x_u, g_y)
-#     G_s = prod_matrix(g_x, d_y_l)
-#     G_n = prod_matrix(g_x, d_y_u)
-#
-#
-#     @test size(d_w*g_y) == (UnknownDim,6)
-#     @test size(d_e*g_y) == (UnknownDim,6)
-#     @test size(d_s*g_x) == (5,UnknownDim)
-#     @test size(d_n*g_x) == (5,UnknownDim)
-#
-#     # These tests should be moved to where they are possible (i.e we know what the grid should be)
-#     @test_broken collect(d_w*g_y) ≈ G_w
-#     @test_broken collect(d_e*g_y) ≈ G_e
-#     @test_broken collect(d_s*g_x) ≈ G_s
-#     @test_broken collect(d_n*g_x) ≈ G_n
-# end
-#
-# @testset "BoundaryQuadrature" begin
-#     op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-#     g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
-#
-#     H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
-#     H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
-#     H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
-#     H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
-#
-#     v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-#
-#     function get_quadrature(N)
-#         qc = op.quadratureClosure
-#         q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
-#         @assert length(q) == N
-#         return q
-#     end
-#
-#     v_w = v[1,:]
-#     v_e = v[10,:]
-#     v_s = v[:,1]
-#     v_n = v[:,11]
-#
-#     q_x = spacing(g)[1].*get_quadrature(10)
-#     q_y = spacing(g)[2].*get_quadrature(11)
-#
-#     @test H_w isa TensorOperator{T,1} where T
-#
-#     @test domain_size(H_w, (3,)) == (3,)
-#     @test domain_size(H_n, (3,)) == (3,)
-#
-#     @test range_size(H_w, (3,)) == (3,)
-#     @test range_size(H_n, (3,)) == (3,)
-#
-#     @test size(H_w*v_w) == (11,)
-#     @test size(H_e*v_e) == (11,)
-#     @test size(H_s*v_s) == (10,)
-#     @test size(H_n*v_n) == (10,)
-#
-#     @test collect(H_w*v_w) ≈ q_y.*v_w
-#     @test collect(H_e*v_e) ≈ q_y.*v_e
-#     @test collect(H_s*v_s) ≈ q_x.*v_s
-#     @test collect(H_n*v_n) ≈ q_x.*v_n
-#
-#     @test collect(H_w'*v_w) == collect(H_w'*v_w)
-#     @test collect(H_e'*v_e) == collect(H_e'*v_e)
-#     @test collect(H_s'*v_s) == collect(H_s'*v_s)
-#     @test collect(H_n'*v_n) == collect(H_n'*v_n)
-# end