Subcycling leads to time drift leading to desync

[fsimonis]

Describe your setup

Operating system (e.g. Linux distribution and version): Ubuntu 00.00
preCICE Version: develop at 05541eb

Describe the problem

Using subcycling leads to a time drift which can result in a desync of participants and a crash once only one participant decides that isCouplingOngoing() == false.

Problem by @BenjaminRodenberg:

max_time   = 1.0
precice_dt = 0.2
solver_dt  = precice_dt / 640 # = 0.0003125

t = 0
NUMERICAL_ZERO_DIFFERENCE = 1e-14

while abs(t - max_time) > NUMERICAL_ZERO_DIFFERENCE:
    i = 0
    t_loc = 0
    while abs(t_loc - precice_dt) > 10e-14:
        i+=1
        t_loc += solver_dt
    t += t_loc
    print(f"Reached {t} after {i} time steps")
    if(t > max_time):
        break

We are using the computed time window to advance timeWindowStart in the BaseCouplingScheme. This should always be moved by _timeWindowSize if available.
This leads to a drift which builds up over time.

Step To Reproduce

1. max-time = 1.0, time-window-size=0.2, P1 does 1 timestep, P2 does 640 timesteps of size 0.2/640.
2. Run the simulation
3. P1 finalizes while P2 crashes as it is waiting for data, but the communication was closed.

Expected behaviour

No drift, no crash.

Additional context

This is related to but not caused by #1788, which advances the time window size correctly.
This used to simply hang with the participant wait in finalize, which is now disabled by default #1600.

There are 3 relevant cases at the end of a time window:

1. the next time step is too small, so we move to the real end Change valid-digits to min-timestep #1788
2. the next time step size is numerically 0 (this doesn't move to the next time window correctly)
3. there is no time window size (nothing to correct)

Possible solution

Always "snap" the provided time of the last timestep to the end of the time window. This also requires correcting the time stamps of the data samples.

precice/src/cplscheme/BaseCouplingScheme.cpp

Lines 298 to 325 in 05541eb

	void BaseCouplingScheme::firstExchange()
	{
	PRECICE_TRACE(_timeWindows, getTime());
	checkCompletenessRequiredActions();
	PRECICE_ASSERT(_isInitialized, "Before calling advance() coupling scheme has to be initialized via initialize().");
	_hasDataBeenReceived = false;
	_isTimeWindowComplete = false;
	PRECICE_ASSERT(_couplingMode != Undefined);
	if (reachedEndOfTimeWindow()) {
	_timeWindows += 1; // increment window counter. If not converged, will be decremented again later.
	//If preCICE has stopped before the end of the time window we have to duplicate the last available sample and put it at the end of the time window.
	// We have to exclude the case where coupling scheme does not have a time window size, since this will cause problem with the interpolation later on
	if (getNextTimeStepMaxSize() > math::NUMERICAL_ZERO_DIFFERENCE && hasTimeWindowSize()) {
	addTimeStepAtWindowEnd();
	// Update the _computedTimeWindowPart in order to keep the time within preCICE synchronised
	// Has to be done before the second exchange, since the serial coupling scheme moves to the new time window before updating _timeWindowStartTime
	_computedTimeWindowPart = _timeWindowSize;
	}
	exchangeFirstData();
	}
	}

if (reachedEndOfTimeWindow()) {

    _timeWindows += 1; // increment window counter. If not converged, will be decremented again later.

    if (hasTimeWindowSize()) {

      //If preCICE has stopped before the end of the time window we have to duplicate the last available sample and put it at the end of the time window.
      // We have to exclude the case where coupling scheme does not have a time window size, since this will cause problem with the interpolation later on
      if (getNextTimeStepMaxSize() > math::NUMERICAL_ZERO_DIFFERENCE) {

        addTimeStepAtWindowEnd();

        // Update the _computedTimeWindowPart in order to keep the time within preCICE synchronised
        // Has to be done before the second exchange, since the serial coupling scheme moves to the new time window before updating _timeWindowStartTime
        _computedTimeWindowPart = _timeWindowSize;
      } else {
        // snap final data samples to the end
        snapDataToTimeWindowEnd();
        _computedTimeWindowPart = _timeWindowSize;
      }
    }

    exchangeFirstData();
  }

[BenjaminRodenberg]

The python script above basically summarizes what's happening with the CPP solverdummy, if you choose to use a maximum time of 1.0, a time window size of 0.2 and one participant subcycles with 640 steps, while the other uses a single step. Note that you also need to comment out reading and writing of data, since otherwise an assertion is triggered before reaching the end of the simulation.

[fsimonis]

Easier than updating the timestamps of write-data Stamples is to update the time before generating the Stamples.

This is something I already do in the upcoming fix for the compositional coupling scheme #1705

[BenjaminRodenberg]

With #1897 it seems like the problem here has shifted a bit. As we are tracking our time differently now, we don't run into any observable trouble for the case described above anymore. The python script below represents the updated situation:

max_time   = 1.0
window_dt = 0.2
solver_dt  = window_dt / 640

t_start = t = 0
NUMERICAL_ZERO_DIFFERENCE = 1e-14

while abs(t - max_time) > NUMERICAL_ZERO_DIFFERENCE:
    i = 0
    precice_dt = window_dt
    while abs(t - t_start - window_dt) > NUMERICAL_ZERO_DIFFERENCE:  # BaseCouplingScheme::reachedEndOfTimeWindow()
        i += 1
        dt = min(solver_dt, precice_dt)
        t += dt
        precice_dt = t_start + window_dt - t
    t_start += window_dt
    t = t_start
    print(f"Reached {t} after {i} time steps")
    if(t > max_time):
        break

This works and we get 640 steps for every window. So the new way for keeping track of time seems to be more robust.

However, if we use 6400 steps for every window, we can also see a window, where 6401 steps are performed:

Reached 0.2 after 6400 time steps
Reached 0.4 after 6401 time steps
Reached 0.6000000000000001 after 6400 time steps
Reached 0.8 after 6400 time steps
Reached 1.0 after 6400 time steps

Additionally, I'm currently looking for a good strategy that might help us to return an error to the user in such a situation.

Note that a user can still quite easily implement a fix in the adapter by following this idea: If there is a danger for precice_dt getting really small in the next step, allow to use a dt that is slightly larger than the solver_dt we would like to use:

max_time   = 1.0
window_dt = 0.2
solver_dt  = window_dt / 6400

t_start = t = 0
NUMERICAL_ZERO_DIFFERENCE = 1e-14

while abs(t - max_time) > NUMERICAL_ZERO_DIFFERENCE:
    i = 0
    precice_dt = window_dt
    while abs(t - t_start - window_dt) > NUMERICAL_ZERO_DIFFERENCE:  # BaseCouplingScheme::reachedEndOfTimeWindow()
        i += 1
        if abs(solver_dt - precice_dt) < 10 * NUMERICAL_ZERO_DIFFERENCE:
            dt = precice_dt
        else:
            dt = min(solver_dt, precice_dt)
        t += dt
        precice_dt = t_start + window_dt - t
    t_start += window_dt
    t = t_start
    print(f"Reached {t} after {i} time steps")
    if(t > max_time):
        break

Important question w.r.t #1867 do we need to be able to detect this situation? If not, we can close this issue, because the actual problem is the adapter code and not preCICE.

[uekerman]

Is the remaining problem really worth investing more time? It is extremely niche and even if someone hits it, it does not give real troubles (crash in finalize).

[BenjaminRodenberg]

@uekerman I agree. Let's close this for now. The crash in finalize has as far as I've tested it so far even disappeared with #1897. There is only the (very niche) possibility of performing +-few substeps, if you are doing a large number of substeps due to round-off errors. But this is less critical as a crash and we also still have to learn a lot about how users will actually use subcycling in the end (I assume doing exactly 6400 substeps is not the main use-case).