Operator interface prototype

[MSeeker1340]

Notebook link

Updated the notations to conform with the convention discussed today. @jlperla can you take a look and confirm the notation choices? I've also expanded on the operator interface so that scalar multiplication and subtractions are now supported. The example in the last section now no longer uses field access.

As per the discussion of yesterday, currently what I'm working on is to embed coefficients into the discretizations of derivative operators along with the update_coefficients! interface, with the main focus being on the upwind operators.

@ChrisRackauckas Could use your opinion on this. I've been thinking about the update_coefficients! implementation and it occurs to me that this is not the only way to supply new coefficients to the operators. Using the current DiffEqArrayOperator as an example (this also works with the embedded coefficients approach):

>>> mu = ones(2); L = DiffEqArrayOperator(Diagonal(mu))

DiffEqArrayOperator([1.0 0.0; 0.0 1.0])

Since both Diagonal and DiffEqArrayOperator stores references instead of copying the array, we can change L directly by modifying mu:

>>> mu[1] = 100.0; L

DiffEqArrayOperator([100.0 0.0; 0.0 1.0])

So for an iterative algorithm that involves updating coefficients of a differential operator, the user can do something like

# Initialization
mu = [...]
sigma = [...]
L = DiffEqArrayOperator(Diagonal(mu)) * ... + DiffEqArrayOperator(Diagonal(sigma)) * ...
while ...
    # compute things
    # update mu and sigma
end

And note there's no need to explicitly update L because it has been handled through reference.

Now I understand this might feel unsafe due to the coefficients arrays being global variables. But as long as we package everything inside a function (which we need to do for the iterative algorithm anyway) this should be fine in my opinion.

This is not to say that we don't need update_coefficients! anymore --- the ODE solvers can still benefit from such a method. But for @jlperla's specific needs I think the direct update approach would suffice.

[ChrisRackauckas]

yup, I mentioned that here: https://github.com/JuliaDiffEq/PDERoadmap/pull/1#issuecomment-385044315 . Looks like we're on the same page.

[jlperla]

When we do move towards the implementation, I think the main change I see in the high level interface is the https://github.com/JuliaDiffEq/PDERoadmap/issues/4 style separation of domain and grid. We can discuss there separately.