Handling conservation laws via solver constraints

[dweindl]

The CVODES constraint handling (see also #2327) would also allow us to enforce conservation laws, and probably make our code quite a bit simpler.

To be discussed whether we'd want to keep the current conservation law handling anyways in order to support other solver at some point, but it seems rather unlikely that this will happen.

This would also solve the currently missing/difficult handling of non-negativity constraints in combination with conservation laws.

This will, however, first require upgrading SUNDIALS (#1565).

[1] Note: relevant memory leaks fixed in v6.4.0