Implement stage-level coupling for split integration

[efaulhaber]

This PR implements stage-level coupling as opposed to the previously used step-level coupling.
In step-level coupling, the fluid is advanced a full step before the structure is advanced a full step. This reduces the stability of the main time integration, reducing the maximum stable time step by a factor of 2 (in my tests with Carpenter-Kennedy).
Stage-level coupling saves the state of the structure before a time step, and then advances it from there to the intermediate stage times of the time integrator for every stage. This means we throw away the structure solution after every stage and start integrating again from the previous time step, which requires ~3x more sub-steps than step-level coupling (in my tests with Carpenter-Kennedy). The advantage is that the stability properties of the time integrator are preserved.

Because of these downsides and the fact that this is probably only required for large time step ratios between fluid and structure, this option is turned off by default.

Δt	Step-level coupling	Stage-level coupling
9e-4
1e-3
1.8e-3

For 1.9e-3, I get "instability detected" with stage-level coupling. The maximum stable time step of 1.8e-3 with stage-level coupling is the same as when making the ball a moving solid wall boundary (non-elastic).
Note that we cannot do this with a CFL because the StepsizeCallback has an upper limit of 1.2e-3 due to #1048.

Note that something is still not right, as the larger time step makes the ball fall deeper with stage-level coupling.

[efaulhaber]

I now implemented a slightly different benchmark simulation. A TLSPH square with E=1e8 and rho=1200 is fully submerged in a fluid with rho=1000. It starts with zero velocity and is then slowly sinking. I measure the stability in the form of the largest stable time step and the accuracy in the form of deviation from a simulation with a very small time step.

The version I showed above was stage-coupling with "restart" (the first row in the table below), which means I reset the sub-integrator in every stage to the state of the previous time step. Without "restart", the sub-integrator is integrated to the first stage time with the fluid state prediction of the first stage, and then continues from there in the second stage, etc.
"Predict" means I apply an explicit Euler step with the structure velocity to predict the structure position at the stage time as u += v * (t_new - t_previous). Then I use this prediction to compute $F_\text{fluid}$.
"Deviation" is the (normalized) difference in the y-coordinates of the square at t=1.5 between the split integration and the non-split integration (fluid and structure integrated together at Δt=1.44e-4), which I consider the reference for this simulation (time integration error almost zero). A negative deviation means the square sank too fast, a positive deviation means it didn't sink fast enough.

Stage-coupling	Restart	Predict	max Δt	Deviation @ Δt=8e-4	Deviation @ Δt=1.6e-3	#Sub-steps @ Δt=8e-4	#Sub-steps @ Δt=1.6e-3
✅	✅	❌	1.6e-3	-8.47e-2	-1.91e-1	38k	35.2k
✅	✅	✅	1.6e-3	-2.50e-4	1.47e-3	38k	35.2k
✅	❌	❌	1.6e-3	-2.01e-2	-3.86e-2	15k	12.3k
✅	❌	✅	1.6e-3	-2.91e-4	-3.95e-4	15k	12.3k
❌	—	❌	8.0e-4	-4.39e-2	—	11.2k	—
❌	—	✅	8.0e-4	3.70e-2	—	11.2k	—

Interpreting these results, we can see that the methods without position prediction all have a significant negative deviation (obviously even more pronounced with the larger time step, for the stage-level coupling that is stable at this time step), indicating an underestimation of $F_\text{fluid}$. This is reasonable because $F_\text{fluid}$ is computed at the previous time step (or stage for the "restart" methods), at which the square is higher up, so the forces from the fluid are underestimated. This effect is smaller for non-restart methods because they compute $F_\text{fluid}$ at the previous stage instead of going further back to the previous step. Prediction significantly increases the accuracy.

The "restart" methods compute the final structure state by integrating from $t_n$ to $t_{n+1}$ with a constant $F_\text{fluid}$ computed at $t_n$. The non-restart methods instead compute $F_\text{fluid}$ based on the predicted fluid state for each stage (which can be considered less accurate since it's a prediction), but re-compute $F_\text{fluid}$ for each state (which, in turn, is more accurate). The resulting method has similar accuracy in this benchmark (higher even for the larger time step) and the same stability, but the number of sub-steps is close to the step-coupling method, whereas the restart methods require 2-3x more stub-steps.
Note that a larger ratio between fluid and structure time step further reduces the overhead of the non-restart methods compared to the step-coupling method.

In summary, the non-restart stage-coupling method with prediction is more stable (allows for a 2x larger fluid time step), more accurate, and not more expensive (does not require significantly more structure sub-steps) than the previous implementation of step-level coupling. For testing purposes, I added kwargs, so all methods can still be tested. I am not sure if the non-restart will work as well with every time integrator as it did with Carpenter-Kennedy.