[WIP] Compute reactions only on selected fields

[anne-glerum]

No need to comment at the moment, I want to see what the tester says first.

This PR adds an argument to the compute_reactions() function so that reactions can be computed on only a subset of fields. In this PR, this is used to eliminate fields tracked by particles from the operator splitting. For example in the case of tracking viscoelastic stresses on particles, the operator split takes place directly on the particles, and should not happen on the fields constructed from the particle properties.

This PR is also in preparation for #6909, where we want to compute reactions for some fields before the nonlinear solver loop, and for some fields after the solver loop (e.g. porosity vs stress updates).

I have added a test, but as noted in the prm file header, I am not sure whether the result is correct - the porosity in the beam is higher than the test I based it on says it should be. However, moving the beam to be the last compositional field does not change the results, so it seems not related to the new indexing in compute_reactions(). Also, both reaction solvers produce the higher porosity.

For all pull requests:

• I have followed the instructions for indenting my code.
• I followed the AI Policy:
  • Significant parts of this PR were written by AI (please add a sentence below).
  • AI has not been used in significant ways.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.