Implement adjoint sensitivity analysis

[schmoelder]

Currently, only first order derivatives are provided using Forward Mode Algorithmic Differentiation.
If implemented, the efficiency for first order would be much higher and second order would be possible with reasonable effort.

[sleweke]

There are several ways to do it:

• Fully discrete adjoints by differentiation of the time integrator as done by the group of U. Naumann at RWTH Aachen
• (Semi-)Discrete adjoints as in Hahn et al. 2014
• Continuous adjoints as in Osterroth et al. 2020

Personally, I'm in favor of continuous adjoints, but it is not clear how to do rapid-equilibrium binding.
Discrete adjoints are doable, especially if we have removed domain decomposition and switched to DG #22. For the backwards integration, we need to solve with the transposed Jacobian.