Document working vs WIP features

The current documentation can be misleading regarding which features are ready vs in-progress. I'm not sure if there's a standard practice for this, but if we can summarize the current state in this issue I'll update the README.

Relatedly, we might consider adding a header to the README that warns users that this library is still a WIP, and not to depend on it for anything critical. Again, not sure if there's already a standard practice for that.

So far, I'm aware that the following are broken / not implemented:

• UpwindDifference
• generally anything for 3+ dimentions

Let's make that list complete and specific.

This will hopefully avoid surprising users, as well as reduce the flow of "yes, we know" issue submissions.

[ChrisRackauckas]

I think this is just too broad to be useful. We should open issues on every specific thing that we run into, and of course there will always be uncertainty somewhere, but this is too broad. This is definitely not a core SciML library yet, which is why we don't directly advertise it anywhere, but the basics should be working according to the tests. But yes, there are edges we need to soften.

UpwindDifference

I think it needs to be more specific. There's tests on things like its action, turning into an array, etc. so most things should be working. Here's an entire set of tests covering what we thought was everything (https://github.com/SciML/DiffEqOperators.jl/blob/73044151e6ba55f3348e0f5b59e8bea91f865923/test/upwind_operators_interface.jl), but it looks like we missed something. Banded concretization seems to be a specific issue covered in #240 . Is there anything else?

In fact, that set of tests covers banded concretizations, so something peculiar is going on with the example or the test is faulty.

generally anything for 3+ dimensions

Here's 1000 lines of tests on 3 dimensional actions:

DiffEqOperators.jl/test/2D_3D_fast_multiplication.jl

Lines 722 to 1647 in 7304415

	@testset "3D Multiplication with identical bpc and dx = dy = dz = 1.0" begin
	N = 100
	M = zeros(N+2,N+2,N+2)
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
	end
	end
	end
	# Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = dy = dz = 1.0
	# Lx3 and Lx4 have the same number of boundary points
	Lx3 = CenteredDifference{1}(3,4,1.0,N)
	Lx4 = CenteredDifference{1}(4,4,1.0,N)
	A = Lx3 + Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)+(Lx4*M))
	# Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = dy = dz = 1.0, no coefficient
	# Ly3 and Ly4 have the same number of boundary points
	Ly3 = CenteredDifference{2}(3,4,1.0,N)
	Ly4 = CenteredDifference{2}(4,4,1.0,N)
	A = Ly3 + Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)+(Ly4*M))
	# Test a single axis, multiple operators: (Lzzz + Lzzzz)*M, dx = dy = dz = 1.0, no coefficient
	# Lz3 and Lz4 have the same number of boundary points
	Lz3 = CenteredDifference{3}(3,4,1.0,N)
	Lz4 = CenteredDifference{3}(4,4,1.0,N)
	A = Lz3 + Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz3*M)+(Lz4*M))
	# Test (Lxxxx + Lyyyy)*M, dx = 1.0, no coefficient, two boundary points on each axis
	M_temp = zeros(N,N,N+2)
	A = Lx4 + Ly4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Lx3 + Ly3
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:])
	# Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
	A = Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test (Lxxxx + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
	M_temp = zeros(N,N+2,N)
	A = Lx4 + Lz4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test (Lxxx + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Lx3 + Lz3
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N])
	# Test multiple operators on both axis: (Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
	A = Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test (Lyyyy + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
	M_temp = zeros(N+2,N,N)
	A = Ly4 + Lz4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test (Lyyy + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Ly3 + Lz3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N])
	# Test multiple operators on both axis: (Lyyy + Lzzz + Lyyyy + Lzzzz)*M, no coefficient
	A = Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test a single operator on each axis: (Lx3 + Ly3 + Lz3)*M, no coefficient
	A = Lx3 + Ly3 + Lz3
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N])
	# Test a single operator on each axis: (Lx4 + Ly4 + Lz4)*M, no coefficient
	A = Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
	# Test multiple operators on each axis: (Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4)*M, no coefficient
	A = Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
	end
	@show "3D BPC"
	@testset "3D Multiplication with differing bpc and dx = dy = dz = 1.0" begin
	N = 100
	M = zeros(N+2,N+2,N+2)
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = CenteredDifference{1}(2,2,1.0,N)
	# Lx3 has 1 boundary point
	Lx3 = CenteredDifference{1}(3,3,1.0,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,1.0,N)
	# Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
	A = Lx2+Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx4*M))
	# Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
	A += Lx3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))
	# Ly2 has 0 boundary points
	Ly2 = CenteredDifference{2}(2,2,1.0,N)
	# Ly3 has 1 boundary point
	Ly3 = CenteredDifference{2}(3,3,1.0,N)
	# Ly4 has 2 boundary points
	Ly4 = CenteredDifference{2}(4,4,1.0,N)
	# Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
	A = Ly2+Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly4*M))
	# Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
	A += Ly3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))
	# Lz2 has 0 boundary points
	Lz2 = CenteredDifference{3}(2,2,1.0,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,1.0,N)
	# Lz4 has 2 boundary points
	Lz4 = CenteredDifference{3}(4,4,1.0,N)
	# Test a single axis, multiple operators: (Lzy+Lzzzz)*M, dz = 1.0
	A = Lz2+Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz4*M))
	# Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M, dz = 1.0
	A += Lz3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))
	# Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
	A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
	A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
	A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
	end
	@testset "3D Multiplication with identical bpc and dx = 0.5, dy = 0.25, dz = 0.05, no coefficients" begin
	N = 100
	M = zeros(N+2,N+2,N+2)
	dx = 0.5
	dy = 0.25
	dz = 0.05
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
	end
	end
	end
	# Test a single axis, multiple operators: (Lxxx + Lxxxx)*M
	# Lx3 and Lx4 have the same number of boundary points
	Lx3 = CenteredDifference{1}(3,4,dx,N)
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	A = Lx3 + Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)+(Lx4*M))
	# Test a single axis, multiple operators: (Lyyy + Lyyyy)*M
	# Ly3 and Ly4 have the same number of boundary points
	Ly3 = CenteredDifference{2}(3,4,dy,N)
	Ly4 = CenteredDifference{2}(4,4,dy,N)
	A = Ly3 + Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)+(Ly4*M))
	# Test a single axis, multiple operators: (Lzzz + Lzzzz)*M
	# Lz3 and Lz4 have the same number of boundary points
	Lz3 = CenteredDifference{3}(3,4,dz,N)
	Lz4 = CenteredDifference{3}(4,4,dz,N)
	A = Lz3 + Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz3*M)+(Lz4*M))
	# Test (Lxxxx + Lyyyy)*M, two boundary points on each axis
	M_temp = zeros(N,N,N+2)
	A = Lx4 + Ly4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Lx3 + Ly3
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:])
	# Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M
	A = Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test (Lxxxx + Lzzzz)*M, two boundary points on each axis
	M_temp = zeros(N,N+2,N)
	A = Lx4 + Lz4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test (Lxxx + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Lx3 + Lz3
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N])
	# Test multiple operators on both axis: (Lxxx + Lzzz + Lxxxx + Lzzzz)*M
	A = Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test (Lyyyy + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
	M_temp = zeros(N+2,N,N)
	A = Ly4 + Lz4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test (Lyyy + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
	A = Ly3 + Lz3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N])
	# Test multiple operators on both axis: (Lyyy + Lzzz + Lyyyy + Lzzzz)*M
	A = Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test a single operator on each axis: (Lx3 + Ly3 + Lz3)*M
	A = Lx3 + Ly3 + Lz3
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N])
	# Test a single operator on each axis: (Lx4 + Ly4 + Lz4)*M
	A = Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
	# Test multiple operators on each axis: (Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4)*M
	A = Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])
	end
	@testset "3D Multiplication with differing bpc and dx = 0.5, dy = 0.25, dz = 0.05, no coefficients" begin
	N = 100
	M = zeros(N+2,N+2,N+2)
	dx = 0.5
	dy = 0.25
	dz = 0.05
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = CenteredDifference{1}(2,2,dx,N)
	# Lx3 has 1 boundary point
	Lx3 = CenteredDifference{1}(3,3,dx,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	# Test a single axis, multiple operators: (Lxx+Lxxxx)*M
	A = Lx2+Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx4*M))
	# Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M
	A += Lx3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))
	# Ly2 has 0 boundary points
	Ly2 = CenteredDifference{2}(2,2,dy,N)
	# Ly3 has 1 boundary point
	Ly3 = CenteredDifference{2}(3,3,dy,N)
	# Ly4 has 2 boundary points
	Ly4 = CenteredDifference{2}(4,4,dy,N)
	# Test a single axis, multiple operators: (Lyy+Lyyyy)*M
	A = Ly2+Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly4*M))
	# Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M
	A += Ly3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))
	# Lz2 has 0 boundary points
	Lz2 = CenteredDifference{3}(2,2,dz,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,dz,N)
	# Lz4 has 2 boundary points
	Lz4 = CenteredDifference{3}(4,4,dz,N)
	# Test a single axis, multiple operators: (Lzy+Lzzzz)*M
	A = Lz2+Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz4*M))
	# Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M
	A += Lz3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))
	# Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M
	A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
	A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
	A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
	end
	@show "3D with coefficients"
	# This is the same test set as the above test set with the addition of coefficients
	@testset "3D Multiplication with coefficients" begin
	N = 100
	M = zeros(N+2,N+2,N+2)
	dx = 0.5
	dy = 0.25
	dz = 0.05
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = 1.234*CenteredDifference{1}(2,2,dx,N)
	# Lx3 has 1 boundary point
	Lx3 = 2.456*CenteredDifference{1}(3,3,dx,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	# Test a single axis, multiple operators: (Lxx+Lxxxx)*M
	A = Lx2+Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx4*M))
	# Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M
	A += Lx3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))
	# Ly2 has 0 boundary points
	Ly2 = CenteredDifference{2}(2,2,dy,N)
	# Ly3 has 1 boundary point
	Ly3 = 4.014*CenteredDifference{2}(3,3,dy,N)
	# Ly4 has 2 boundary points
	Ly4 = 1.49*CenteredDifference{2}(4,4,dy,N)
	# Test a single axis, multiple operators: (Lyy+Lyyyy)*M
	A = Ly2+Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly4*M))
	# Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M
	A += Ly3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))
	# Lz2 has 0 boundary points
	Lz2 = 1.546*CenteredDifference{3}(2,2,dz,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,dz,N)
	# Lz4 has 2 boundary points
	Lz4 = 0.55*CenteredDifference{3}(4,4,dz,N)
	# Test a single axis, multiple operators: (Lzy+Lzzzz)*M
	A = Lz2+Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz4*M))
	# Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M
	A += Lz3
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))
	# Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M
	A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
	A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
	A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
	end
	@show "3D irregular"
	@testset "x, y, and z are all irregular grids" begin
	N = 100
	dx = cumsum(rand(N+2))
	dy = cumsum(rand(N+2))
	dz = cumsum(rand(N+2))
	M = zeros(N+2,N+2,N+2)
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx[i])+sin(dy[j]) + exp(dz[k])
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = CenteredDifference{1}(2,2,dx,N)
	# Lx3 has 1 boundary point
	Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	# Ly2 has 0 boundary points
	Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
	# Ly3 has 1 boundary point
	Ly3 = CenteredDifference{2}(3,3,dy,N)
	# Ly4 has 2 boundary points
	Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)
	# Lz2 has 0 boundary points
	Lz2 = CenteredDifference{3}(2,2,dz,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,dz,N)
	# Lz4 has 2 boundary points
	Lz4 = CenteredDifference{3}(4,4,dz,N)
	# Test composition of all first-dimension operators
	A = Lx2+Lx3+Lx4
	M_temp = zeros(N,N+2,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)
	# Test composition of all second-dimension operators
	A = Ly2+Ly3+Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)
	# Test composition of all third-dimension operators
	A = Lz2+Lz3+Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ (Lz2*M + Lz3*M + Lz4*M)
	# Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
	A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	# Numerical errors accumulating for this test case, more so than the other tests
	@test_broken M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])
	# Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
	A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])
	# Test composition of all operators
	A = Lx2+Lx3+Lx4+Ly2+Ly3+Ly4+Lz2+Lz3+Lz4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]
	end
	@testset "irregular x grid and regular y and z grids (dy = 0.25, dz = 0.05)" begin
	N = 100
	dy = 0.25
	dz = 0.05
	dy_vec = 0.25:0.25:(N+2)*0.25
	dz_vec = 0.05:0.05:(N+2)*0.05
	dx = cumsum(rand(N+2))
	M = zeros(N+2,N+2,N+2)
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx[i])+sin(dy*j)+exp(dz*k)
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = CenteredDifference{1}(2,2,dx,N)
	# Lx3 has 1 boundary point
	Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	# Ly2 has 0 boundary points
	Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
	# Ly3 has 1 boundary point
	Ly3 = CenteredDifference{2}(3,3,dy,N)
	# Ly4 has 2 boundary points
	Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)
	# Lz2 has 0 boundary points
	Lz2 = 1.4*CenteredDifference{3}(2,2,dz,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,dz,N)
	# Lz4 has 2 boundary points
	Lz4 = CenteredDifference{3}(4,4,dz,N)
	# Test that composition of both x and y operators works
	A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Lx3*M)[1:N,2:N+1,:]+(Lx4*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Ly3*M)[2:N+1,1:N,:]+(Ly4*M)[2:N+1,1:N,:])
	# Test that composition of both x and z operators works
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4
	M_temp = zeros(N,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lx3*M)[1:N,:,2:N+1]+(Lx4*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lz3*M)[2:N+1,:,1:N]+(Lz4*M)[2:N+1,:,1:N])
	# Test that composition of x, y, and z operators works
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4 + Ly2 + Ly3 + Ly4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N])
	# Last case where we now have some `irregular-grid` operators operating on the
	# regular-spaced axis y
	# These operators are operating on the regular grid y, but are constructed as though
	# they were irregular grid operators. Hence we test if we can seperate irregular and
	# regular gird operators on the same axis
	# Ly2 has 0 boundary points
	_Ly2 = CenteredDifference{2}(2,2,dy_vec,N)
	# Lz3 has 1 boundary point
	_Lz3 = CenteredDifference{3}(3,3,dz_vec,N)
	# Ly4 has 2 boundary points
	_Ly4 = 12.1*CenteredDifference{2}(4,4,dy_vec,N)
	A += _Ly2 + _Lz3 + _Ly4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N] + (_Ly2*M)[2:N+1,1:N,2:N+1]+(_Lz3*M)[2:N+1,2:N+1,1:N]+(_Ly4*M)[2:N+1,1:N,2:N+1])
	end
	@testset "irregular y grid and regular x and z gris (dx = 0.25, dz = 0.05)" begin
	N = 100
	dx = 0.25
	dz = 0.05
	dx_vec = 0.25:0.25:(N+2)*0.25
	dz_vec = 0.05:0.05:(N+2)*0.05
	dy = cumsum(rand(N+2))
	M = zeros(N+2,N+2,N+2)
	for i in 1:N+2
	for j in 1:N+2
	for k in 1:N+2
	M[i,j,k] = cos(dx*i)+sin(dy[j])+exp(dz*k)
	end
	end
	end
	# Lx2 has 0 boundary points
	Lx2 = CenteredDifference{1}(2,2,dx,N)
	# Lx3 has 1 boundary point
	Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
	# Lx4 has 2 boundary points
	Lx4 = CenteredDifference{1}(4,4,dx,N)
	# Ly2 has 0 boundary points
	Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
	# Ly3 has 1 boundary point
	Ly3 = CenteredDifference{2}(3,3,dy,N)
	# Ly4 has 2 boundary points
	Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)
	# Lz2 has 0 boundary points
	Lz2 = 1.4*CenteredDifference{3}(2,2,dz,N)
	# Lz3 has 1 boundary point
	Lz3 = CenteredDifference{3}(3,3,dz,N)
	# Lz4 has 2 boundary points
	Lz4 = CenteredDifference{3}(4,4,dz,N)
	# Test that composition of both x and y operators works
	A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
	M_temp = zeros(N,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Lx3*M)[1:N,2:N+1,:]+(Lx4*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Ly3*M)[2:N+1,1:N,:]+(Ly4*M)[2:N+1,1:N,:])
	# Test that composition of both y and z operators works
	A = Ly2 + Ly3 + Ly4 + Lz2 + Lz3 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	# Need to figure out why this test is exploding
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	############################################################################
	# Tests to isolate the above problem
	############################################################################
	A = Ly2 + Ly3 + Ly4
	M_temp = zeros(N+2,N,N+2)
	mul!(M_temp, A, M)
	@test M_temp ≈ Ly2*M + Ly3*M + Ly4*M
	### Test the addition of Lz2, Lz3, Lz4 seperately
	A = Ly2 + Ly3 + Ly4 + Lz2
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N])
	A = Ly2 + Ly3 + Ly4 + Lz3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N])
	A = Ly2 + Ly3 + Ly4 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz4*M)[:,2:N+1,1:N])
	###
	A = Ly2 + Ly3 + Ly4 + Lz2 + Lz3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])
	A = Ly2 + Ly3 + Ly4 + Lz2 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	A = Ly2 + Ly3 + Ly4 + Lz3 + Lz4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	###
	A = Ly2 + Ly3 + Ly4 + Lz2 + Lz2
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz2*M)[:,2:N+1,1:N])
	# It appears that multiple z operators with some y operators is causing the issues
	###
	A = Lz2 + Lz3 + Lz4
	M_temp = zeros(N+2,N+2,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ Lz2*M + Lz3*M + Lz4*M
	A = Lz2 + Lz3 + Lz4 + Ly2
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	A = Lz2 + Lz3 + Lz4 + Ly3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	# It seems that the y paddign could be the issue
	A = Lz2 + Lz3 + Lz4 + Ly4
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	# It seems that the y paddign could be the issue
	A = Lz2 + Lz3 + Lz4 + Ly4 + Ly3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	# It seems that the y paddign could be the issue
	A = Lz2 + Lz3 + Lz4 + Ly4 + Ly2
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
	# It seems that the y paddign could be the issue
	###
	A = Lz2 + Lz3 +Ly3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test_broken M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])
	A = Lz3 +Ly3
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N])
	A = Lz2 + Lz3 +Ly2
	M_temp = zeros(N+2,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])
	# It seems that the padding of y is forcing multiple z operators to fail
	############################################################################
	############################################################################
	# Test that composition of x, y, and z operators works
	A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4 + Ly2 + Ly3 + Ly4
	M_temp = zeros(N,N,N)
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N])
	# Last case where we now have some `irregular-grid` operators operating on the
	# regular-spaced axis x and z
	# These operators are operating on the regular grid x and z, but are constructed as though
	# they were irregular grid operators. Hence we test if we can seperate irregular and
	# regular gird operators on the same axis
	# Lx2 has 0 boundary points
	_Lx2 = CenteredDifference{1}(2,2,dx_vec,N)
	# Lz3 has 1 boundary point
	_Lz3 = CenteredDifference{3}(3,3,dz_vec,N)
	# Lx4 has 2 boundary points
	_Lx4 = 12.1*CenteredDifference{1}(4,4,dx_vec,N)
	A += _Lx2 + _Lz3 + _Lx4
	mul!(M_temp, A, M)
	@test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
	+ (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N] + (_Lx2*M)[1:N,2:N+1,2:N+1]+(_Lz3*M)[2:N+1,2:N+1,1:N]+(_Lx4*M)[1:N,2:N+1,2:N+1])
	end

. Regular grids seem fine. Irregular grids have some remaining issues. 3D BCs are 🤷 though.