Discussion on design of domain and grid interface

[jlperla]

This issue isolates and continues discussions in https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260 with @dlfivefifty and @ChrisRackauckas on a joint design of the interface for domains, etc. between the ApproxFun, DiffEqOperators, and ModelingToolkit libraries. We hope to implement it in the prototype: https://github.com/JuliaDiffEq/PDERoadmap/pull/1

Also see https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/test/domains.jl

For a starting point in one dimension,

d = Interval(1.0, 10.0) #Continuous domain
x = RegularGrid(d, N)

Note: I think it is better to give points than dx, but I am open. The reason is that if you give a dx then you never know how well things line up, or how many points it ends up creating.

When things move to an irregular grid,

d = Interval(1.0, 10.0) #Continuous domain
grid_points = 1.0:0.01:10.0 #This is regular, but you get the point.
x = IrregularGrid(d, grid_points)

For the product of a discrete domain with 2 values and continuous domain (which comes up a lot in economics applications),

d = DiscreteDomain(2) * Interval(0.0,5.0) #Or tensor product?

Or, when we later want to enable naming of types for the discrete domain (e.g. for a continuous-time markov chain of switching between L and H productivity in a model)

d1 = Interval(1.0, 2.0)
d2 = DiscreteDomain([:L, :H])
d = d1 ⊗ d2 #Is the tensor product the right one to use here?

Or, instead, does it make sense to use an enumeration? I really don't like how scoping works with enumerations in Julia, which limits the clarity of the code.

[ChrisRackauckas]

To get @MSeeker1340 up to speed, this is for being able to define the operators without the discretization. DiffEqOperators is all about how do you build the computational operator to compute, but there's a higher level question we should be asking as to, how can I just symbolically/mathematically specify some differential operator and have the correct computational operators be generated? At this level, it's not a finite difference, FEM, spectral, etc. discretization: just some operator on a domain which would make it easy to switch between forms (if we get it right).

ModelingToolkit.jl is an intermediate representation (IR) for models. You can think of it almost as a mini-CAS that DSLs can write into and compute on. The variables of this IR can hold domains (https://github.com/JuliaDiffEq/ModelingToolkit.jl/blob/master/src/variables.jl#L10) and there are differential operators. We can modify it at will to make the IR have what we need to generate operators.

However, I think we should get some operators working first before we talk about auto-generating. In some sense, this is the "FEniCS" part and we need dolphin first.

[MSeeker1340]

At this level, it's not a finite difference, FEM, spectral, etc. discretization: just some operator on a domain which would make it easy to switch between forms (if we get it right).

Well this is basically symbolic math right? Maybe we can get some inspiration from Mathematica? (I don't have much experience with it though).

The type system of lazy operator composition/combination already has a strong resemblance to a parse tree. So I wouldn't be surprised if we can write something that translates an abstract operator defined in ModelingToolkit to a concrete one that shares the same tree structure.

But yes I agree we should stick to simple finite-difference realizations of the operators at the moment, while maybe think a bit about incorporating ApproxFun.

[ChrisRackauckas]

Maybe we can get some inspiration from Mathematica?

We're building this with help from REDUCE.

But yes I agree we should stick to simple finite-difference realizations of the operators at the moment, while maybe think a bit about incorporating ApproxFun.

Yes, first things first.

[dlfivefifty]

The tree-vs-type hierarchy question is important to get right. In ApproxFun I overused templated variables, which puts a limit on how complicated the operator tree can be. However, building operators should be fast. It’s not clear to me what the right balance is.

[ChrisRackauckas]

The tree-vs-type hierarchy question is important to get right. In ApproxFun I overused templated variables, which puts a limit on how complicated the operator tree can be. However, building operators should be fast. It’s not clear to me what the right balance is.

We did it a little wrong once. https://github.com/JuliaDiffEq/ModelingToolkit.jl/issues/62 . I think the key is to separate "compile time" and "runtime" here. The resulting operator can be fully typed, just the computation to get it may not be inferred. That's fine because building the operator would either be at the top scope or above a function barrier, i.e. building it is a form of compilation so like how the Julia AST itself is untyped, the specification of the operator should be fully untyped but builds strongly typed code. We'll at least give it a try.

[jlperla]

I understand the design concerns down the road when we might want to mess around with the IR, but does that change the basic types for writing the domains and grids?

In particular, isn't the design of https://github.com/JuliaDiffEq/PDERoadmap/issues/4#issue-318138562 style interfaces largely just going to help with dispatch and the user's API?

[ChrisRackauckas]

I understand the design concerns down the road when we might want to mess around with the IR, but does that change the basic types for writing the domains and grids?

Nope it shouldn't. It just carries those along but doesn't strictly use them. We can use them to do things like generate ODEs in local coordinates on known manifolds, but in terms of the IR its mostly for storing information.

[jlperla]

@dlfivefifty @ChrisRackauckas What are the thoughts on using Interval(...) here, exactly as in ApproxFun? vs. 0.0:1.0 vs. 0.0 .. 1.0?

I was experimenting with instead using the return from the range built into Julia (i.e. 1.0:10.0) but

julia> 1.0:10.0 |> typeof
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
julia> 1.0:10.0
1.0:1.0:10.0

As much as it would be great in notation to support defining intervals with 0.0:1.0, I do not think this usage scenario fits the Range concept, which is inherently iterable.

The better fit is probably https://github.com/JuliaIntervals/IntervalRootFinding.jl What about using the .. operator?:

julia> using IntervalArithmetic

julia> 1.0..10.0 |> typeof
IntervalArithmetic.Interval{Float64}

julia> 1.0..10.0
[1, 10]

Are there other precedents for using this notation for intervals?

[dlfivefifty]

a .. b is the agreed upon convention for intervals, see IntervalSets.jl

Abusing a UnitRange is a bad idea.

[dlfivefifty]

(It looks like you aren’t aware that ApproxFun also supports a .. b )

[jlperla]

Perfect! a .. b as sugar to create the Interval it is. Though this tells me that ApproxFun and IntervalArithmetic are incompatible, and explains the error I saw. Should some sort of shared interval base library be established?

(Even in the longrun, with my ADL suggestion in https://discourse.julialang.org/t/function-name-conflict-adl-function-merging/10335 wouldn't help here, since the operator .. needs to work on two floats in both cases)

For the grids, are there any other libraries that would be nice to maintain compatibility with? For https://github.com/JuliaMath/Interpolations.jl it seems that they use the StepRangeLen type directly for the grid, but I wonder if we are better off maintaining a connection from a grid to a domain.

[jlperla]

Also https://github.com/sisl/GridInterpolations.jl has its own grid objects

[dlfivefifty]

The whole point of IntervalSets.jl is to be a shared package for a .. b .

[jlperla]

Done. What about grids? Any movement towards a coordinated library for those types?

[dlfivefifty]

Not on METADATA but we’ve been working on https://github.com/daanhb/Domains.jl

Unfortunately it was triggering a compile-time Julia bug when used in ApproxFun, which put a hold on development. We also decided “domain” should be an interface rather than an abstract type, and anything that supports in can be treated as a domain.

[jlperla]

@dlfivefifty What are your thoughts on domains with discrete numbers of elements? Did you have that in mind with Domains.jl? In particular, having them "named"? I don't think using a enumeration is a good idea due to namespace clashes, but symbols work well.

@ChrisRackauckas Thoughts on getting uniform grids up and running for the prototype in DiffEqOperators.jl ? I think that IntervalSets.jl as a dependency is an easy addition? A temporary type to hold grids may be reasonable until things we are able to coordinate on a central grid type.

[dlfivefifty]

Making domain an interface has the benefit that, since 1 in 1 works, numbers can be interpreted as points and vectors can be interpreted as grids. So in 1D I’d just use linspace. For 2D one would want an array of SVector I think.

[jlperla]

What is the connection between domains, intervals, etc.? It sounds like you are saying that the domain is the discretized grid that the function is defined on. Is that correct?

If so, then what do we need for the generic algorithms to work:

• in is necessary, as you say
• Some testing of equality is needed?
  • The reason is ultimately to be able to ensure that operators we are composing are defined on the same domains, and potentially expand to the superset of the domain
• Don't we need iteration as well for generic algorithms to be implemented on the grids themselves? If so, then https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-iteration-1 would also be a requirement of a domain?

Where I start to get confused in the interface, is how this would all work for algorithms with multiple dimensions. I see how the in works, but how do you write generic code? I think this becomes especially important with designing the operator interface, because I presume the iteration is used to go through the state with the update_coefficients! design.

[dlfivefifty]

In my usage domain is a possibly infinite set. So all of these would implement the domain interface:

a..b 
a:b
1
[1,2,3]
Set([1,2,3])

The first is infinite (an interval) the rest are discrete.

Iteration is not possible with infinite sets but is for any of the discrete ones.

To generalise to higher dimensions, I interpret the “ambient space” of the domain as all elements of the type, and would use something like SVector{2,T} to represent 2D points.

@daanhb May have some thoughts.

[jlperla]

OK, then for the "named" discrete sets, vectors of symbols would work perfectly.

[:L, :H]

Thinking about your domain, the generality makes complete sense for spectral methods. I guess the question with finite-difference methods is whether there is need to have separate generic types to for the continuous domain and the discretization of the continuous domain. That was the assumption I had with my suggested

N = 20
d = 1.0..10.0
x = RegularGrid(d, N)

But, is the d really needed to keep around, or is the x enough? For example, we could just go

x = 1.0:0.1:10.0 #For a regular grid
x = linspace(1.0, 10.0, 20) #Generates the same type
x_irregular = [1.0, 1.1, 1.5, 2.0, 5.0, 10.0] #Irregular grid

One worry in that approach would be whether you can have specialized code for regular vs. irregular grids. But I think it works completely fine, since the type of the range and the vector are different and it could dispatch on StepRangeLen. That is

julia> x = linspace(0.0, 1.0, 10)
0.0:0.1111111111111111:1.0

julia> typeof(x) <: StepRangeLen
true

julia> typeof(collect(x)) <: StepRangeLen
false

julia> step(x)
0.1111111111111111

[ChrisRackauckas]

Yes yes yes and yes.

Domains as an interface with in, iteration, array interface, and SVectors or scalars. SVector is important so that way (u,p,t,x) is universal for N-dimensional. We're only talking about grids now, continuous domains are for a higher level.

This isn't pertinent right now though. We can add grids to the operators quite easily once we have that setup. It will be important when defining the irregular grids, but it seems that interface works with arrays of scalars / SVector so it's fine.

[jlperla]

We're only talking about grids now, continuous domains are for a higher level. OK, then no need to have the separation of the continuous domain and the discretization of it.

Some specific quetsion:

• So, lets say we are composing some operators. Would we then have then keep the "grid" with each of them (which could be used to ensure that the composed operators are defined on the same grid type). That is,

  x = linspace(0.0, 1.0, 100) #100 node grid
  sigma_2 = 1.1
  mu = -0.2
  L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)

  But what if there are the functions for the coefficients? Is the idea that the "functions" would conform to whatever the "grid"?

  x = linspace(0.0, 1.0, 100) #100 node grid
  sigma_2(u,p,t,x) = 1.1 * x^2
  mu(u,p,t,x) = -0.2 * x
  L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)

• Alternatively, if we want the grid associated with every operator type, there could be a wrapper, as we had discussed,

  
  x = linspace(0.0, 1.0, 100) #100 node grid
  sigma_2_tilde(u,p,t,x) = 1.1 * x^2
  mu_tilde(u,p,t,x) = -0.2 * x
  sigma_2 = DiffEqFunction(sigma_2_tilde, x)
  mu= DiffEqFunction(mu_tilde, x)

L = mu DriftOperator(x, :forwards) +sigma_2 DiffusionOperator(x)

> SVector is important so that way (u,p,t,x) is universal for N-dimensional.
I am not sure I understand you here.  Are you saying that `x` in the interface would be a `SVector`?  I would think a tuple is more appropriate, especially when there are discrete states in `x`?  (Yes, I know `SVector` is implemented as a tuple, but it is a tuple where everything is of the same type).

[jlperla]

Although I would love it if for multiple dimensions we could use a named tuple. In particular, the canonical case in economics is a discrete i in {L,H} state which determines the drift of a continuous state, y. For example, if the drift is 0.1 * y if i = :L and 0.5 * y if i=:H, I would love to be able to write this as something like:

julia> x = @NT(y=0.1, i=:L)
(y = 0.1, i = :L)

julia> μ(u,p,t,x) = x[:i]==:L ? 0.1*x[:y] : 0.5*x[:y]
μ (generic function with 1 method)

julia> μ(0,0,0,x)
0.010000000000000002

julia> μ(0,0,0,@NT(y=0.5, i=:H))
0.25

(where the @NT is obviously dropped in v0.7

[ChrisRackauckas]

Although I would love it if for multiple dimensions we could use a named tuple. In particular, the canonical case in economics is a discrete i in {L,H} state which determines the drift of a continuous state, y. For example, if the drift is 0.1 y if i = :L and 0.5 y if i=:H, I would love to be able to write this as something like:

NamedTuples satisfy the interface. You can do that if you want.

But what if there are the functions for the coefficients? Is the idea that the "functions" would conform to whatever the "grid"?

Yes, it would just be on the user side the x is of the right dimension. It would be like an FEM node listing.

[jlperla]

NamedTuples satisfy the interface. You can do that if you want.

Is it that simple? Presumably the interface of what is passed into the x of the function has something to do with the types of the dimensions of the grid? Forgetting about the named tuples for a second, it seems like in the construction of the tuple itself the underlying code in the library is going to have to take the Cartesian product of all of the elements of the types?

It seems like if I want to have an operator which works on the product of a continuous and discrete grid, I would have to put it in a tuple, or something like that:

domain = ( linspace(0.0, 1.0, 10), [:L, :H] ) #A tuple of the domain

...or add syntactic sugar for the cartesian product. Then generic code calling the functions in the library could construct a tuple to call the underlying functions of type Tuple{eltype(domain[1]), eltype(domain[2])}

And in order to have names, wouldn't we need to have those named associated with the underlying domain. That is, something like AxisArrays?

i_domain = Axis{:i}([:L, :H])
y_domain= Axis{:y}(1.0:0.1:2.0)

That way, a named tuple could be created to call the user functions by getting the names of each axis with AxisArrays.axisname(i_domain) etc. and then the type with eltype(i_domain), etc.

[dlfivefifty]

Non-lazy Cartesian products of grids would be SVector.(x, y’). Lazy could be done with BroadcastArray currently in InfiniteArrays.jl, or via a special type.

[jlperla]

OK, but

julia> x = collect(linspace(0.0, 1.0, 5))
5-element Array{Float64,1}:
 0.0
 0.25
 0.5
 0.75
 1.0

julia> y = [:L, :H]
2-element Array{Symbol,1}:
 :L
 :H

julia> z = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> SVector.(x, z')
5×2 Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}:
 [0.0, 1.0]   [0.0, 2.0]
 [0.25, 1.0]  [0.25, 2.0]
 [0.5, 1.0]   [0.5, 2.0]
 [0.75, 1.0]  [0.75, 2.0]
 [1.0, 1.0]   [1.0, 2.0]

julia> SVector.(x, y')
ERROR: MethodError: no method matching transpose(::Symbol)

Basically, this approach of the cartesian products seems to promote discrete dimensions, and doesn't allow symbols?

If I add in a transpose,

julia> Base.transpose(x::Symbol) = x

julia> SVector.(x, y')
5×2 Array{StaticArrays.SArray{Tuple{2},Any,1,2},2}:
 Any[0.0, :L]   Any[0.0, :H]
 Any[0.25, :L]  Any[0.25, :H]
 Any[0.5, :L]   Any[0.5, :H]
 Any[0.75, :L]  Any[0.75, :H]
 Any[1.0, :L]   Any[1.0, :H]

It looks the "promotion" turns these into "Any", which is not what we want.

Isn't the basic problem that SArray are designed to hold a single type?

[dlfivefifty]

Then use tuple. SVector is only needed if you want to add your points.

[jlperla]

OK. So

julia> using IterTools

julia> grids = (1.0:0.1:2.0, [:L, :H])
(1.0:0.1:2.0, Symbol[:L, :H])

julia> [i for i in product(grids...)]
22-element Array{Tuple{Float64,Symbol},1}:
 (1.0, :L)
 (1.1, :L)
 (1.2, :L)
 (1.3, :L)
...

So at least IterTools could be used to get an enumeration on the cartesian problem? I suppose something similar with a Axis generating a NamedTuple would also be possible?

[ChrisRackauckas]

Yes, any of these AbstractVector with scalars or AbstractVectors are a good way to specify the points of a domain. It doesn't matter which you choose, though being able to determine dx will matter for operators.

[jlperla]

Although in the case of discrete states, the dx is not meaningful, since the infinitesimal generator is intensity matrix describing the jump intensity of the continuous-time Markov-chain.

A completely different alternative is that we don't generally let the user define functions of everything, and instead compose up the operators themselves by taking kroneckers of the underlying operators on univariate operators.

[ChrisRackauckas]

Although in the case of discrete states, the dx is not meaningful, since the infinitesimal generator is intensity matrix describing the jump intensity of the continuous-time Markov-chain.

Yes exactly, so it would need an array operator and this space wouldn't be compatible with most of the operators.

A completely different alternative is that we don't generally let the user define functions of everything, and instead compose up the operators themselves by taking kroneckers of the underlying operators on univariate operators.

Yup that's what we're doing.

[jlperla]

@ChrisRackauckas @MSeeker1340 @dlfivefifty I wanted to resurrect this issue since we are making a push on adding differential operators.

I need to have both regular and irregular grids in the very short term, and an irregular grid X a discrete state in the short/medium term, so I think it is worth designing the interface correctly.

My general feeling is that any derivative operator should you should always pass in the "grid", where if it is an AbstractRange you can call the step(r) to get the size of the delta for uniform grids, and if it is AbstractArray then you can assume that it could be irregular. Later on, we could check that the grids === when we compose different operators, or do some fanciness to take the union of them if they don't.

That is, something like

x = linspace(0.0, 1.0, 10) #regular grid since returns an AbstractRange
L = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)

x2 = collect(x) #Or create an irregular grid
L2 = mu * DriftOperator(x2 , :forwards) +sigma_2 * DiffusionOperator(x2 ) #Same basic interface

Alternatively, maybe it makes sense to separate the interface for the grid and the underlying domain/space, etc. but that is a design question above my paygrade for integrating ApproxFun and the DiffEqVerse.

The interface for taking the 2 dimension grid of a discretized continuous state X an enumerable discrete state is the other place that the design of the grids/domains/etc. enters. But lets say I want to have a PDE in a continuous x variable and a discrete variable i in :l and :h where there are different drifts in the stochastic process for the :L vs :H. Then I would love to be able to do something like

x = linspace(0.0, 1.0, 10) 
D = DriftOperator(x, :forwards) 
D2 = DiffusionOperator(x)
sigma = 0.01
mu_L = 0.02
mu_H = 0.05

L_L = mu_L * D  + sigma^2/2 * D2
L_H = mu_H * D + sigma^2/2 * D2
#Stack the independent operators?  Making this up
L = [:L 0; 0 :H] * [L_L; L_H]; %Product of discrete state and continuous state?

Obviously the last state is made-up, but I would how we could have something of similar complexity.

[dlfivefifty]

I don’t remember what the plan was but I’m pretty sure the right design is a software tiers that “lower”. That way ApproxFun+DiffEq can use the appropriate tiers while not forced to integrate everything on the operator level.

So For finite differences we want the following, from highest to lowest:

1. Differential operators and boundary conditions (like you wrote)
2. Lazy banded matrices and lazy constraints

3. Concrete banded matrices (or SymTridiagonal when possible) with concrete constraints.

   The linear algebra happens to incorporate constraints is either at tier 2 or 3.

ApproxFun (without adaptivity) has an extra step:

1. Differential operators and boundary conditions
2. Lazy infinite banded matrices with constraints
3. Lazy finite banded matrices with constraints
4. Concrete almost banded matrices

In the future, there could be another layer on top that represents the differential operators without any discretisation details, that can then be solved via either approach.

The reason for the extra “lazy” tier is that this is going to make distributed memory super easy, as the lazy matrices will be what’s communicated.

[jlperla]

For the creation of the Q boundary extrapolation (if you recall the setup documented in the latex) could the notation be (using https://github.com/JuliaApproximation/ApproxFun.jl/issues/590).

Boundary condition for u'(0) + a * u(0) = 0 and u'(xbar) = 0 be something like

Q = [(𝒟 + a * I )→ 0; 𝒟 → xbar]

Then the composed operator applying the boundary conditions would be

x = linspace(0.0, 1.0, 10) #regular grid since returns an AbstractRange
L1 = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)
L = Q * L1 #Applies boundary extrapolation to L1

Where at that point we can use L as a linear map (or affine, in general)?

ApproxFun.jl, as far as I understand it, does not need to separate boundary conditions from equations, so the other approach is simply

x = linspace(0.0, 1.0, 10) #regular grid since returns an AbstractRange
L1 = mu * DriftOperator(x, :forwards) +sigma_2 * DiffusionOperator(x)
L = [L1;
      (𝒟 + a * I )→ 0;
      𝒟 → xbar];

But I am not sure how that works with the Q approach here. These sort of things are sugar we can wait on, but thinking through the interface might be worthwhile.

[dlfivefifty]

ApproxFun doesn't support time-stepping, though it supports time-evolution linear PDEs like zero Dirichlet Heat equation a la (using the to-be-implemented new syntax for restriction operators, and using u(t,x)):

Dx = Derivative([0,1])
Dt = Derivative([1,0])
[ I → (0..1) × [-1,1];
  I → 0 × (0..1);
  Dt - Dx^2 ]  \ [[0,0], u_0, 0]

In theory, newton and it's auto-diff parser will eventually work for non-linear PDEs, here would be Burgers:

    N=u->[u(:,-1);
                u(:,1);
                u(0,:)-u_0;
                 Dt*u + u*Dx*u - ν*Dx^2*u]
    u = newton(N) # sequence of linear solves

In terms of time-stepping and DiffEq, I'm imagining something like (for Dirichlet heat, copying DAE example):

function heat(out, du, u, _, t)  
    out[1] = 𝒟^2*u - du
   out[2] = ( I → -1) * u
   out[3] = ( I →   1) * u
end
prob = TimeEvolutionPDEProblem(heat, u_0, tspan)
solve(prob, FiniteDifferences(1000)) # Or Ultrapherical() or Fourier() for ultraspherical spectral method a la ApproxFun

and for Burgers

function burgers(out, du, u, ν, t)  
    out[1] = du  + u*(𝒟*u) - ν*𝒟^2*u
   out[2] = ( I → -1) * u
   out[3] = ( I →   1) * u
end
prob = TimeEvolutionPDEProblem(burgers, u_0, tspan)
solve(prob, FiniteDifferences(1000)) # Or Ultraspherical() or Fourier() for ultraspherical spectral method a la ApproxFun

[jlperla]

In terms of time-stepping and DiffEq, I'm imagining something like

I think this is a great thing to think through, especially in the context of https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/260#issuecomment-410939651 s

The thing that stands out to me is the

function heat(out, du, u, _, t)  
    out[1] = 𝒟^2*u - du
   out[2] = ( I → -1) * u
   out[3] = ( I →   1) * u
end

block. When handling about boundary conditions at a high level of generality, there is little distinction between the differential equations and the boundary conditions (except that the latter may be on a manifolds). You just stack `em.

Now with the PDERoadmap math here with the Q in https://github.com/JuliaDiffEq/PDERoadmap/releases/download/v0.2/linear_operators_overview.pdf that is not really the case, since we treat the boundary conditions as something different (interpreted as being a few steps of gaussian elimination, if I recall your interpretation). So shouldn't this seems like the heat equation above could be written as something like,

Q = [( I → -1);  ( I →   1)];
L = 𝒟^2
heat = Q * L
prob = TimeEvolutionPDEProblem(heat, u_0, tspan)

[dlfivefifty]

I don't understand the Q * L part.

I think we need to focus on linear for now, as nonlinear is more complicated (and may require some auto-diff to get working).

Instead of worrying about PDEs first, we should get the lowest level right: which might be some sort of ConstrainedLinearODEProblem that involves only matrices. So maybe something like this would work:

Δ = BandedMatrix{Float64}(undef, (n-2, n), (0,2))
Δ[band(0)] = Δ[band(2)] = 1/h^2
Δ[band(1)] = -2/h^2

B = zeros(2,n)
B[1,1] = B[end,end] = 1

R = eye(n)[2:end-1,:]

prob = ConstrainedLinearODEProblem(B, R, Δ)

solve(prob, GaussReduce(Tsit5(), pivots = (1,n)))

Note the pivots can be any two columns of B that are non-zero, and could also be calculated on the fly.

This is too simple as nothing varies in time. If we have a simple time-varying PDE like u_t = Δ*u + cos(t)*u it could be something like:

prob = ConstrainedLinearODEProblem(t -> B, t -> R, t -> Δ + cos(t)*I)

solve(prob, GaussReduce(Tsit5(), pivots = (1,n)))

Not sure what to do if the pivots themselves need to vary (think cos(t)*Dirichlet + sin(t)*Neumann)...

This is also too simple as there is no support for forcing terms. Something like:

prob = ConstrainedLinearODEProblem(t -> (B,c), t -> R, t -> (Δ + cos(t)*I,f))

solve(prob, GaussReduce(Tsit5(), pivots = (1,n)))

Definitely needs some iteration....

[dlfivefifty]

Following LazyArrays.jl, we could do some broadcast magic (in Julia v1.0) and write things like this:

function heat(du, u, (B, R, Δ, c, f), t)
    Mul(B, u) .= c
    Mul(R, du) .= Mul(Δ + cos(t)*I, u) .+ f
end

Then u and du would be special objects that capture the broadcast.

[ChrisRackauckas]

@dlfivefifty 's syntax up there, at least the first one, has far too much magic and only supports DAEs so I don't like that one. I would like to keep everything as short syntactic operators, or use a DSL. I think it should commit to being automatic or being fully flexible, since anything in between gets complicated but limited.

I don't understand the Q * L part.

That's a generic way to handle linear boundary conditions on any linear operator.

I think we need to focus on linear for now, as nonlinear is more complicated (and may require some auto-diff to get working).

Nonlinear has been working for over a year now and with autodiff.

Instead of worrying about PDEs first, we should get the lowest level right: which might be some sort of ConstrainedLinearODEProblem that involves only matrices.

That's not necessary. A ConstrainedLinearODEProblem is just an ODEProblem with a linear operator and a mass matrix, or it's not fully linear and needs a nonlinear solve anyways on the W matrix.

Anyways, there's no reason to argue at that level of the design because that's very clear anyways. It's the discretization handling for named PDEs that's unclear.

[dlfivefifty]

The whole point is ConstrainedLinearODEProblem is not the same as a mass matrix. One way to solve a ConstrainedLinearODEProblem is to reduce it to a mass matrix formulation by Gaussian elimination, but there exists other methods like IMEX which incorporate the constraints in the implicit solves.

(A special case is when R is a view of the identity matrix, in which case we can reduce further to a standard ODE problem.)

Mixing in the language of finite differences / differential operators / boundary conditions with the actual details of the solve, which only involves linear algebra on a constrained linear ODE, will be limiting.

[ChrisRackauckas]

The whole point is ConstrainedLinearODEProblem is not the same as a mass matrix. One way to solve a ConstrainedLinearODEProblem is to reduce it to a mass matrix formulation by Gaussian elimination, but there exists other methods like IMEX which incorporate the constraints in the implicit solves.

Yes, and that's a SplitODEProblem with a mass matrix. That was exactly my point.

Mixing in the language of finite differences / differential operators / boundary conditions with the actual details of the solve, which only involves linear algebra on a constrained linear ODE, will be limiting.

That's exactly why I don't want to do that.

[dlfivefifty]

It's a SplitODEProblem with a mass matrix (R = view(Eye(n),2:n-1,:) in the case of FD) and also boundary conditions. I don't see how boundary conditions are incorporated in a standard SplitODEProblem without further manipulation.

Also, we should have an idea of what to do with 2D time evolution. Here there's the added difficulty that boundary conditions over-pose the constraints (and work because the RHS is in a special subspace).

[ChrisRackauckas]

Boundary conditions are just handled like any other conditions by the mass matrix.

Also, we should have an idea of what to do with 2D time evolution. Here there's the added difficulty that boundary conditions over-pose the constraints (and work because the RHS is in a special subspace).

That's operator-specific. The finite difference ones can lazy Kron. Since the intermediate level isn't supported, as long as each operator should just do whatever is appropriate for it.

[dlfivefifty]

Boundary conditions are just handled like any other conditions by the mass matrix.

But the user has to do this. I thought the whole point was to make a system that incorporates the boundary conditions automatically, without the user having to set-up their own mass matrices?

[ChrisRackauckas]

No, the whole point is to not have this intermediate level because it's exactly the level that's not possible to do without sacrificing features (and likely performance). At the operator implementation level, you can be given a package for operators and easily build the derivative function. At the high level, you say what kind of discretization you would want on some named equation and get something that solves it with only specifying BCs in an abstract form. At this intermediate level, you'd need every possible operator to satisfy quite a complex interface which won't handle some important things like nonlinearity, so that makes it difficult but not necessarily useful since then using the lower interface is more flexible and the high interface already solves the problem in most cases. So since this is the piece that's the most work with the least payoff I don't want to support it

[dlfivefifty]

At the intermediate level there are no operators, just matrices. It's not a complex interface, it's the AbstractArray interface. And if nonlinearity can be incorporated at all, it can definitely be incorporated at this level.

By doing things on the level of linear algebra, you have something that other discretisation methods can also use.

[ChrisRackauckas]

I'll agree that level is nice in theory and can be useful for composing different pieces, but it's not a level that anyone has really expressed interest in using. The lower level methods are for implementation of discretizations so it's for the developers and those who really just want helpers for building their own solver. The high level methods are what people have expressed a lot of interest in: can I guess say in a high level way that I want you to solve a semilinear diffusion advection equation, and have you come up with how to do it correctly. The middle, expressing the abstract PDE with usable operators that double as code, is really cool and likely where I want to be in the long term but its a lot of work looking for a problem.

[dlfivefifty]

It's the level no one wants because they don't realise everything is just structured linear algebra. When it comes to parallelisation, working in this language will make parallelising a lot easier, because it can be done on the linear algebra level, not on the level of discretisation.

It sounds like you want to go down the same route as FeNICs, Firedrake, Dedalus, etc. which specialise on particular discretisation methods, and jumble up the linear algebra with the details of the discretisation, and with the details of parallelising.

This is fine if you only ever want to solve PDEs, and specific ones at that. But I'm more interested a framework that is flexible enough to handle other problems just as easily (boundary element methods, singular integral equations, Benjamin–Ono equation, fractional differential equations).

[dlfivefifty]

Though maybe our wires are crossed. The really important thing is that SplitODEProblem supports lazy structured matrices (that may depend on time) that are realised as a distributed structured matrix remotedly: that is, the workers have access to the lazy matrix and know how to populate their blocks.

The mass matrix vs boundary conditions conversation is less crucial, and actually you might be right that its hard to incorporate boundary conditions automatically into the mass matrix while retaining structure.