What should the size of an a discretized operator from R^n -> R^n be?

[xtalax]

I've noticed that throughout the codebase, we are assuming that an AbstractDiffEqLinearOperator has a size of the form (N,M). But since they may now in general be mappings between multiple dimensions, I think this doesn't really hold, and that we should try to move away from the idea that these operators are some form of fancy matrix to some more general kind of abstraction.

D2x = CenteredDifference{1}(2, 6, dx, length(x))
D2y = CenteredDifference{2}(2, 6, dy, length(y))
D2z = CenteredDifference{3}(2, 6, dz, length(z))

L = D2x + D2y + D2z

We need a method to check that some AbstractDiffEqLinearOperator has a size that is compatible with some u. .

There is a sense in which

size(L) = ((length(x)-2, length(y)-2, length(z)-2), (length(x), length(y), length(z)))

But when L is discretized, it is concretized in to a matrix, so that it's size would be

size(Array(L)) = prod.(size(L))

The source of ambiguity

Now consider the following operator, with the intention that it will operate on some u::AbstractArray{T,3}

L = D2x+D2y

then there is a sense in which

size(L) = ((length(x)-2, length(y)-2, n-2), (length(x), length(y), n)  

where n is some arbitrary number that is as yet unknown.
But the operator has no idea that it is going to multiply a higher dimensional u. So there is also a sense in which

size(L) = ((length(x)-2, length(y)-2), (length(x), length(y)))

So the question is, what should the size of L be given as in this case?
This is also a problem when concretizing

L = D2x + D2z

since the operator doesn't know about length(y), it is still ambiguous what it's final size should be, even though we know that it will operate on a u of at least 3 dimensions.

[ChrisRackauckas]

Interesting. Maybe we just assume it's 1, and have some kind of identity derivative operator in the Y direction that can impart size to it?

[xtalax]

That or we let DerivativeOperator take size instead of length, since this will likely be known at the time of construction

[ChrisRackauckas]

That could work too.

[ajozefiak]

One thing I'll add is that the current concretization implementation works as follows:

Array(L) produces a matrix that operates on the first dimension, regardless of what the differentiation dimension of L is.

But Array(L, (x,y,z,...)), where (x,y,z,...) is a tuple of integers of length n, produces a concretization that operates on a flattened n-dimension array of shape (x,y,z,...), and which correctly operates on the dimension specified in L.