Taking 2D operators seriously

[clason]

First, I apologize for the presumptuous title; I couldn't think of a better one since the question is a bit broad. Please feel free to edit it to something more appropriate.

For prototyping algorithms for PDE-constrained optimization, I frequently need a simple discretization of two-dimensional elliptic partial differential operators on a uniform structured grid, for which FEniCS is overkill. (To demonstrate how simple: in Matlab I often just use A = gallery("poisson",N)...)

I am aware that I can apply such an operator via the tensor trick Au = D1*u + u*D1 and use IterativeSolvers to solve the corresponding PDEs, but for various reasons (avoiding unnecessary inexactness in linear solves, precomputing factorizations for performance, including in block operators for one-shot or Newton-type methods for solving optimality systems), I need an actual matrix representation for A (ideally, in either a sparse or banded form).

Is there any way to obtain such a matrix from DiffEqOperators? Are there any (plans for) convenience constructors?

[ChrisRackauckas]

Yes, we plan to have a lazy Kron so that way you can build the block-banded matrix, but it will take some work and I wouldn't expect it too soon.

[clason]

That is good to know (and no hurry)! In the meantime, is there some other way I can construct such a matrix (it doesn't have to be really efficient since I'm only doing it once as a setup)?

[ChrisRackauckas]

AbstractMatrix(A), Matrix(A), etc. We need to update the docs as well. Give us a a month or so to put this into shape.

[clason]

That's fair; I was in fact worried that this was a bad time for a feature (or documentation) request. But just to make sure: This should be possible right now by manually forming the Kronecker product? Because neither Matrix(D1) nor AbstractMatrix(D1) (nor kron(D1,Id) for a manually constructed identity matrix) works for me...

[ChrisRackauckas]

I thought we made those for the structured full matrices already, but we may still need to upgrade from using full (just got back from Juliacon so I need to check where we're at in more detail). Once the full/banded/etc matrix is made kron should just be the normal kron. And I wonder if there's a special kron for BlockBandedMatrices.

[clason]

Ah, got it: It's not Matrix(D1), it's now convert(Array,D1). Or sparse to directly get a sparse matrix. Then the following works:

using LinearAlgebra, SparseArrays, DiffEqOperators
N = 10
h = 1/N
D1 = sparse(DerivativeOperator{Float64}(2,2,h,N,:Dirichlet0,:Dirichlet0))
Id = sparse(I,N,N)
A = kron(D1,Id) + kron(Id,D1)

(EDIT: Although this gives me a singular matrix, so I must have made a mistake with the boundary conditions. Got it. Might be useful to others so I'm adding a brief tutorial to the /docs folder.)