Refactoring of the ODE module

[AlexandreMagueresse]

As the ODE module of Gridap is expanding, its structure is becoming a little messy and redundant. @santiagobadia and I propose the following changes

ODETools

• [x] Replace the current OperatorType subtypes by the more generic types MassLinear and ConstantMass. These traits are not fully taken advantage of for now.
• [x] Reindex the jacobians (from zero) to match the time derivative order.
• [x] Unify the signature of allocate_residual and allocate_jacobian to accept AbstractVectors and Tuple{Vararg{AbstractVector}}.
• [x] Enable allocate_cache to accept supplementary arguments (in the main signature).
• [x] Bring all the butcher tableaus under an AbstractButcherTableau{T} type, where T is a trait used to tell whether the tableau corresponds to an explicit / (diagonally) implicit / implicit-explicit scheme.
• [ ] For now, Runge-Kutta schemes will only be defined for mass-linear operators, and will throw a @notimplemented for general transient FE operators.

TransientFESpaces

• [ ] Replace Transient(...)RKFEOperatorFromWeakForm by a more general TransientMassLinearFEOperatorFromWeakForm to represent ODE operators that are linear in $\partial_{t} u$. This type was originally created for Runge-Kutta schemes, but it represents a much more general setting that can benefit many other schemes.
• [ ] Mass-linear operators are currently represented by lhs and rhs, which are somewhat confusing names. These are currently defined through writing the ODE in the form $\mathrm{LHS}(t, u) \partial_{t} u = \mathrm{RHS}(t, u)$, so in fact LHS corresponds to the existing jacobian_t (the mass matrix). The new representation for these operators is $$\mathrm{mass}(\partial_{t} u, v) + \mathrm{res}(t, u, v) = 0,$$ where $\mathrm{mass}$ is linear with respect to $\partial_{t} u$ and $\mathrm{res}$ does not depend on $\partial_{t} u$.
• [ ] IMEX Runge-Kutta schemes need a special TransientFEOperator as we need to split the residual in two or three parts:
• General case: $\mathrm{res_implicit}(t, u, v) + \mathrm{res_explicit}(t, u, v) = 0$.
• Semilinear case: $\mathrm{mass}(\partial_{t} u, v) + \mathrm{res_implicit}(t, u, v) + \mathrm{res_explicit}(t, u, v) = 0$, where here again the mass bilinear form is linear with respect to $\partial_{t} u$, and neither residual can depend on $\partial_{t} u$.