How to describe the 'state' of a simulation.

[bangerth]

@luca-heltai and I were having discussions about time steppers that advance multiphysics codes, like the ones in the SUNDIALS multirate splitting methods at https://sundials.readthedocs.io/en/latest/arkode/Mathematics_link.html#arkode-mathematics-splittingstep . That also made me think about our roll-back scheme for when a time step fails. Let me outline what I see as the issues, and then perhaps propose how we should encode this.

• When advancing one physics in a multiphysics code, integrators such as those in https://sundials.readthedocs.io/en/latest/arkode/Mathematics_link.html#arkode-mathematics-splittingstep see the time evolution as an ODE of the form $\dot y = f_1(t,y) + f_2(t,y) + \cdots + f_N(t,y)$ where each of the terms $f_i$ encodes one physical aspect. We can think of $f_1$ as the flow fields, for example, then $f_2,\ldots$ as the compositional fields, then $f_M$ as the particles, then the free surface, and then some more components of the ASPECT universe. Programmatically, for each of these physics, we need access to a "state vector" that encodes the current state of the simulation. In our current scheme, whenever we advance one part of the universe, we update the current state variable solution (or sometimes current_linearization_point) or the particle solution vector, or the free surface, or ... That gives rise to a specific operator splitting scheme, which we should think of as the Lie scheme or a variation. But really, this is conceptually difficult to understand because it isn't ever quite clear at which time point each part of the "solution" actually is. And it isn't possible to specify at around which solution vector a particular function that advances a physics function should linearize things. In other words, in the formula above, we write things as $\dot y = f_1(t) + f_2(t) + \cdots + f_N(t)$ and each of the $f_i$ functions has some sort of implicit notion what it takes as $y$, and then obtains it through the environment (i.e., this->).
• @tjhei started an effort to allow for "roll-back" when a time step fails in Timestepping: allow repeating of time steps #3842. This is made difficult because we don't ever collect what the "state" actually consists of, and there are undoubtedly bugs in the code base because we forget to roll back one part of the state -- implement timestep repeating with mesh deformation #3923 is one example; Unify restoration of particle backup between solver schemes #6858 addresses another place where code was duplicated but correct. The core issue here, as above, is that if $y$ denotes the "state" of the simulation, then that state is scattered across many variables in multiple objects, and we don't have a way to say

  saved_state = state;
  ...do updates...
  if (repeat_time_step)
  {
    state = saved_state;   // roll-back changes
    goto repeat;
  }

It is not easy for me to think of a scheme that would make all of that easy. But perhaps one way to start doing that would be to introduce a class SimulatorState in which we collect the current state -- or whatever parts of the state we can reasonably store in one place. Thoughts about this issue?

[bangerth]

I find this sort of confusion of what state is and where it is stored also exemplified in this comment and the ones below: #3932 (comment)