Suppose I have some time dependent PDE, which can be written in the strong form as
$$ \frac{\partial u}{\partial t} + \mathcal{L}(u) = f $$
Where $\mathcal{L}$ is some differential operator. If I write the weak form of the problem without discretizing the time term, I end up with the standard weak form problem of finding $u_h \in X_h$ such that
$$ m\left(\frac{\partial u_h}{\partial t},\, v\right) + a(u_h,\,v) = f(v), \qquad \forall v \in X_h $$
At this point, I usually discretize in time in the following way (using backward-euler here) as
$$ m\left(\frac{u_h^{k+1} - u_h^k}{\Delta t},\, v\right) + a(u_h^{k+1 },\,v) = f(v), \qquad \forall v \in X_h $$
at which point, I know how to solve the problem by putting the $m(u_h^k, v)$ term to the right hand side, choosing a test, trial space, discretizing and solving.
But what if I wanted to start from the second equation and structure my code to use a black box integrator? Is there a standard way of doing this? I know that the right hand side terms can become much more complicated with higher-order time discretizations.
Edit:
It seems that the answer is to first discretize in space, and express the discrete system as a matrix system $B\dot{U_h} + AU_h = F$. However, if the black box time integrator library does not support DAE systems (allowing user to input $B$) matrix, you're out of luck (this seems to be the case with scipy.integrate.ode). If $B$ can be inverted cheaply, then the black box integrator can handle the equivalent problem $\dot{U_h} = -B^{-1}(AU_h - F)$. Alternatively, mass lumping can be used, but this seems ridiculous if you were interested in implicit time-marching schemes to begin with.
