[WIP] Enable operator split reactions at end of the timestep as well

[anne-glerum]

For discussion with @gassmoeller

This PR enables the user to choose to apply the operator split at the end of a timestep instead of at the beginning. The PR also requires this order when elasticity is included.

Since the time interpolation of the current linearization point is no longer called after the reaction solve, the compute_reactions function requires no changes I think.

I ran the simple stress relaxation benchmark, and this shows the same results. In this benchmark, the stresses only decay with a factor \eta_ve / \eta_e every timestep applied through the reactions. As the timestep is fixed, it doesn't matter whether the reaction solve is applied at the beginning on the previous stress, or at the end of the timestep on the still unchanged stress.

Question 1: When we discussed making this change, we decided to put the compute_reactions after the new timestep is computed and to pass this timestep to the reaction solve. Looking at the code now, I wonder why this would be necessary. Wouldn't we want to use the same timestep that is used for the advection and rotation of the stresses?

Question 2: In main, atm no reaction solves occur during timestep number 0. But in timestep number 1, a reaction solve does occur. Since this reaction actually belongs to t0, should it even happen?

TODO:

• Check whole of elasticity.cc to update the comments
• Check that the reaction rates are computed at the right timesteps
• Also reorder the update for the particles, so far only for fields. For particles, the update happens every time after they are restored to their state of the beginning of the time step. But as the calculation of the reaction rates doesn't have to be changed, they are still correct. The signal the particle update connects to could be changed, so that the update is also applied at the end of the time step.
• Run more benchmarks

For all pull requests:

• I have followed the instructions for indenting my code.

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.

[gassmoeller]

I like it! I have a suggestion for how to make this more future proof (see below).

To your questions:

When we discussed making this change, we decided to put the compute_reactions after the new timestep is computed and to pass this timestep to the reaction solve. Looking at the code now, I wonder why this would be necessary. Wouldn't we want to use the same timestep that is used for the advection and rotation of the stresses?

I think the old way of computing reactions for elastic stresses was abusing the operator splitting in a way it was not intended for. It applied the update of timestep i in timestep i+1, while ordinarily the reactions in timestep i are intended to be applied in timestep i. So yes, I think you are right, we should be using the current time step length.

In main, atm no reaction solves occur during timestep number 0. But in timestep number 1, a reaction solve does occur. Since this reaction actually belongs to t0, should it even happen?

That is for the same reason I mentioned above. The reactions in time step number 1 are intended to be the reactions applying to time step 1, not the ones applying to time step 0. For all reactions except for the elastic stress fields that was always the case, we just used it differently for the elastic stress fields.

[gassmoeller]

This inconsistency seems concerning and prone to bugs in the future. Would it be possible to change the particles the same way?

[anne-glerum]

Yes, I agree, this is still on the TODO list at the top of the PR. I didn't manage to include it last week.

[anne-glerum]

Hi Rene, I like the flexibility with the struct, I'll update the changes accordingly.

In main, atm no reaction solves occur during timestep number 0. But in timestep number 1, a reaction solve does occur. Since this reaction actually belongs to t0, should it even happen?

That is for the same reason I mentioned above. The reactions in time step number 1 are intended to be the reactions applying to time step 1, not the ones applying to time step 0. For all reactions except for the elastic stress fields that was always the case, we just used it differently for the elastic stress fields.

To clarify, my question only pertains to the use of reactions for elastic stresses. In line with the "misuse" of the reactions to update the stress of the previous timestep, I don't think the reactions should be computed in t1, as they would belong to t0. However, which is my doing, the function fill_reaction_rates in elasticity.cc only checks whether the time step number is not zero, not whether it is larger than 1. The compute_reactions function itself also only checks for the time step number not being zero - which is of course correct when using the reactions as intended.
Therefore I think the current reactions for updating the stresses when elasticity is included, in main, should only be computed at t2. But such a change would be superseded by this PR.

[anne-glerum]

I have included the Strategy struct instead of the boolean. I called it ReactionStrategy etc to match the other parameters that use 'reaction' rather than 'operator splitting', but this can be changed.

I also looked at a way to update the particles at the same time as the fields, instead after each restore_particles call. There is a signal post_nonlinear_solver called at the end of each nonlinear solver scheme, but this communicates Solver information. I don't think we want it to include a particle_manager. Therefore I have added a new signal in the same location as the after_nonlinear_solver location for the operator split. The previous post_restore_particles signal could be removed, as the elastic_stress particle property was the only thing using it.

The annoying thing is, the stresses on the particles are updated after the Stokes solve, which means their properties are not interpolated onto the fields again before postprocessing. The compositional fields therefore hold different values from the particles.

I'll think about it some more, but any thoughts are appreciated.

[anne-glerum]

Options/todos discussed off-GitHub:

• Use the existing post_linear_solver, and get the particle manager through SimulatorAccess
• Keep the post_restore_particles signal
• Run 3D simple shear + rotation benchmark

[gassmoeller]

@anne-glerum see my PR to your branch here anne-glerum#19