Mass matrix change

[MSeeker1340]

These are the proposed changes to how mass matrix should be handled in DiffEq. Basically, the mass_matrix field is no longer stored as a field of *DEProblem but instead *DEFunction.

The main reason for these changes is to make sparse/lazy W operator handling in OrdinaryDiffEq simpler SciML/OrdinaryDiffEq.jl#416 (comment). Apart from this, there are also be cases where knowing what the mass matrix is solely by querying the *DEFunction might help. One case I can think of is converting a problem with non-unity mass matrix to a regular problem, i.e. Mu' = f(u) -> u' = M^(-1)f(u) = g(u). Since a lot of integrators do not use the mass matrix explicitly (e.g. the ExpRK integrators) this can come handy. By storing mass_matrix in ODEFunction we can do something like:

function convert_canonical(fun::ODEFunction)
  f = (u,p,t) -> M \ fun.f(u,p,t)
  ODEFunction(f; mass_matrix=I, analytic=fun.analytic, ...)
end

An additional benefit is when a single *DEProblem type support multiple *DEFucntion. For example, DynamicalProblem, SecondOrderODEProblem and SplitODEProblem are all ODEProblem in disguise, each with their own *DEFunction type. Whereas previously ODEProblem.mass_matrix need to handle all these cases, now the mass matrices are handled by the functions and can thus allow specialization (e.g. in calc_W!).

These changes are only the first step though. @ChrisRackauckas can you create a new branch for DiffEqBase with these changes, so that I can update OrdinaryDiffEq and other packages using it as a base?

[ChrisRackauckas]

This needs a tuple mass matrix (mm,I)

[ChrisRackauckas]

Or formalize it as two mass matrix fields: mm1 and mm2 in DynamicalODEFunction. This might be more clear.

[MSeeker1340]

DynamicalODEFunction should always have a tuple (mm1, mm2) as its mass matrix? (I can probably make this a requirement in the type declaration).

If I understand it correctly, the functions here are intended to convert a DynamicalODEProblem to a SecondOrderODEProblem right? If this is the case, then I believe we should inherit the mass matrix from the original problem/function.

[ChrisRackauckas]

A Dynamical ODE is a 2-partitioned ODE, so it could have two mass matrices:

http://docs.juliadiffeq.org/latest/types/dynamical_types.html#Mathematical-Specification-of-a-Dynamical-ODE-Problem-1

[MSeeker1340]

A Dynamical ODE is a 2-partitioned ODE, so it could have two mass matrices

I'm slightly confused. If what you're referring to is separate mass matrix for the u and v parts, then the current implementation already uses a 2-tuple mass matrix and I did the same with DynamicalODEFunction.

I think keep a single mass_matrix field makes the interface clearer. After all, having f.mm1 and f.mm2 is functionally the same as having f.mass_matrix = (mm1, mm2). Of course, if we use a getter method instead, then none of these are needed as we can just get_mm(f::DynamicalODEFunction) = (get_mm(f.f1), get_mm(f.f2)).

[ChrisRackauckas]

mass_matrix=f.mass_matrix that's passing on a tuple?

[ChrisRackauckas]

oh misread it, my bad.

[ChrisRackauckas]

Nah let's take this out.

[ChrisRackauckas]

Let's take this out.

[MSeeker1340]

One thing I forgot to mention is how the default mass matrix for SplitFunction and DynamicalODEFunction can be automatically set in the constructor: f2.mass_matrix for SplitFunction (instead of I) and (f1.mass_matrix, f2.mass_matrix) for DynamicalODEFunction (instead of (I, I)). The constructor for these two are currently a bit messy though, so we can probably do this in a future refactor of diffeqfunction.jl (we need this anyways to drop the outdated trait interface).

Of course there's also the alternative option of dropping f.mass_matrix in favor of a gettor method, e.g.

get_mm(f::ODEFunction) = f.mass_matrix
get_mm(f::SplitFunction) = get_mm(f.f2)`
get_mm(f::DynamicalODEFunction) = (get_mm(f.f1), get_mm(f.f2))

This is a bit restrictive since we cannot override the mass matrix for SplitFunction and DynamicalODEFunction. Not sure if this is needed though.

[ChrisRackauckas]

SplitFunction should really only have one mass matrix though?

[MSeeker1340]

SplitFunction should really only have one mass matrix though?

Seeing that SplitFunction, or in general an N-split function, is just a sum of component functions, I would say yes. This is the case before the changes as well, for example, in https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/126/files#diff-86615b6bd63cae214249239a6d0e855a.

I think that, just like seeing a regular ODEFunction with mass matrix $M$ as $M^{-1}f$, a N-split function can be seen as $M^{-1}(f_1 + f_2 + \cdots + f_m)$. We can always generalize and allow different mass matrix for different components, i.e. $M_1^{-1}f_1 + \cdots + M_m^{-1}f_m$, but then this cannot be converted to the original definition of mass matrix $Mu' = f(u)$.

[ChrisRackauckas]

I think that, just like seeing a regular ODEFunction with mass matrix $M$ as $M^{-1}f$, a N-split function can be seen as $M^{-1}(f_1 + f_2 + \cdots + f_m)$. We can always generalize and allow different mass matrix for different components, i.e. $M_1^{-1}f_1 + \cdots + M_m^{-1}f_m$, but then this cannot be converted to the original definition of mass matrix $Mu' = f(u)$.

I'm just not really sure what that generalization would even mean physically.

[ChrisRackauckas]

For split functions, it's IMEX so f1 should be the mass matrix of interest?

[coveralls]

Coverage decreased (-0.02%) to 91.908% when pulling 7338265 on MSeeker1340:mass-matrix into a5a94e0 on JuliaDiffEq:master.

[ChrisRackauckas]

LGTM. I don't see a mm -> f.f2 assumption anywhere so it sounds like we're agreeing.