FSI: Implicit subcycling incorrect displacement

[vidulejs]

While working on issue #58, I discovered the real reason for instability in FSI simulations with subcycling. Initially it ticked me off, that the mesh was changing during subcycling steps. Maybe this was intended but I was not sure.

However, I realized that at each subcycling step the displacement is read again. This is the total displacement and not multiplied by subcycling dt. This means that for N subcycling steps, the total displacement will be amplified N times.

EDIT: What I thought before was wrong and not the cause of the issue. The core of the problem is that the implicit subcycling simulations converge to different solutions than the non-subcycling ones. See my new comment below.

[davidscn]

However, I realized that at each subcycling step the displacement is read again

For FSI and subcycling, this would be the expected behavior.

This is the total displacement and not multiplied by subcycling dt.

That's the expected behavior either. preCICE takes care of the time interpolation between the beginning and the end of the time window, so each intermediate state will just be in between those two displacement states. The read displacement is total, yes, but the state of the fluid grid within a simulation is already deformed, so there is only an incremental portion left to reach this absolute state. I think this might go wrong when using relative displacements. But from your description, I understand you are dealing with absolute displacements, right?

This is the total displacement and not multiplied by subcycling dt

Why should the displacement be multiplied by dt? What would be the unit of this quantity?

This means that for N subcycling steps, the total displacement will be amplified N times.

I dont understand this and disagree, could explain this more in detail?

[vidulejs]

Excuse the pasted text, I'm in a slight rush. I was in the process of doing a convergence study on the perpendicular flap with subcycling when I wrote most of this text.

The idea was finding a timestep, time-window and subcycling timestep that converges to the same solution. The natural assumption is that as the timestep goes down, the solutions will converge. What I found instead was that there are two solutions that the simulations converge to; (1) the non-implicit-subcycling cases and (2) the implicit-subcycling cases. See the plots below.

Note that explicit subcycling converges to the first group IIRC, so the problem is strictly with implicit subcycling. However I cannot conclude whether the problem is on the OpenFOAM adapter side or somewhere else (the solid participant is not the problem though, same issue with Calculix and Deal II solid).

Timestep convergence analysis

First I perform the convergence study on the timestep. Reference case is parallel-implict fluid-openfoam solid-calculix with time-window value="0.01" and waveform-degree="1". Both participants have equal timestep same as time-window, no subcycling.

<?xml version="1.0" encoding="UTF-8" ?>
<precice-configuration>
    
  <!--  -->

  <data:vector name="Force" waveform-degree="1"/>
  <data:vector name="Displacement" waveform-degree="1"/>
  
  <!--  -->

  <coupling-scheme:parallel-implicit>
    <time-window-size value="0.01" />
    <max-time value="5" />
    <participants first="Fluid" second="Solid" />
    <exchange data="Force" mesh="Solid-Mesh" from="Fluid" to="Solid" substeps="true"/>
    <exchange data="Displacement" mesh="Solid-Mesh" from="Solid" to="Fluid" substeps="true"/>
    <max-iterations value="50" />
    <relative-convergence-measure limit="5e-3" data="Displacement" mesh="Solid-Mesh" />
    <relative-convergence-measure limit="5e-3" data="Force" mesh="Solid-Mesh" />
    <acceleration:IQN-ILS>
      <data name="Displacement" mesh="Solid-Mesh" />
      <data name="Force" mesh="Solid-Mesh" />
      <preconditioner type="residual-sum" update-on-threshold="false" />
      <filter type="QR2" limit="1e-2" />
      <initial-relaxation value="0.5" />
      <max-used-iterations value="100" />
      <time-windows-reused value="15" />
    </acceleration:IQN-ILS>
  </coupling-scheme:parallel-implicit>
</precice-configuration>

[...]

Subcycling convergence analysis

Displacements

Subcycling does have an issue after all. All the subcycling cases converge to a different solution than the non-subcycling ones - independent of timestep!

The subcycling cases are much smoother. That is a clue to something. I've previously seen that the boundary conditions are not applied properly during substeps (but that may have to do with openfoam write time and paraview visualization).

Forces

smooth the forces proportional to 1/dt number of points (smoothing the curves by higher frequency)

Deal II solid solver

Please don't pay attention to the waveform-degree-0 cases. I understand that (constant time interpolation) waveform-degree-0 is not meant to be used for real simulations.

Displacements

Waveform degree 0 cases show the same behaviour as seen with Calculix. dealii-dt-0.01-subcycling-dt-0.005-waveform-0 diverge exponentially just as before, with a similar curve. dealii-dt-0.01-subcycling-dt-0.001-waveform-0 doesn't diverge, just as before. What is special about this combination?

Difference from solid-Calculix is that the waveform degree 1 cases don't converge to something. Maybe requires even smaller timestep.

Forces

---

For reproducibility I created a tutorials and adapter fork with visualization scripts and force and displacement debug logging.

Tutorials fork with the modified solver and python scripts

The modified pimpleFoam solver writes probedLocations.dat which contains the coordinate of the flap upper right corner point.

Modified OpenFOAM adapter which writes forces and displacements from preCICE to files

There's a lot of things that I've tried to diagnose the cause of this problem from the adapter side, but that's a more detailed discussion for which I need more time to write it down.

[MakisH]

Thanks for documenting! Looking forward to more details.

So, as you write: the clear information at the moment is that, for the perpendicular-flap tutorial (FSI), subcycling itself (on or off, independent of the step) seems to affect the values. This does suggest a technical issue, in my opinion, and I still believe it should be in the way we do checkpointing in the adapter.

Assuming that the non-subcycling coupling group is the "correct" one, and trying to brainstorm further directions:

1. Have we shown that there is no similar issue in conjugate heat transfer cases? Let's take the flow-over-heated-plate tutorial and run a similar study, to have evidence (once you have time again).
   A) If same behavior: then, the issue is not specific to the displacements or forces (computation or checkpointing).
   B) If different behavior: it might be.
2. If 1B: Does maybe the number of nCorrectors, nNonOrthogonalCorrectors, or nOuterCorrectors (fvSchemes) affect the subcycling plots? Since the latter two are already 1, I would first try modifying the nCorrectors 4.

@uekerman @davidscn do you see something else interesting in these plots? Do you have any further ideas on what we could try to better pinpoint the cause?

(Related to or duplicate of #58)

[uekerman]

To make things clearer: Some notation we typically use:

• dt_F, dt_S: the time step sizes
• Dt: the time window size

I agree that the results look buggy. What I would do next is to only let one side subcycle to maybe get a grasp on whether reading or writing data is an issue. Meaning

• dt_S = Dt fixed, dt_F = Dt / N, N = 1,2,4,8
• dt_F = Dt fixed, dt_S = Dt / N, N = 1,2,4,8
• As a reference: 8 times smaller Dt and dt_F=dt_S=Dt

All use waveform-degree="1" and substeps="true".

If 1B: Does maybe the number of nCorrectors, nNonOrthogonalCorrectors, or nOuterCorrectors (fvSchemes) affect the subcycling plots? Since the latter two are already 1, I would first try modifying the nCorrectors 4.

All this should make no difference if the coupling converges. Please always check that.

[MakisH]

* dt_F = Dt fixed, dt_S = Dt / N, N = 1,2,4,8

Unfortunately, we do have a clear issue with CalculiX subcycling, and we already confirmed that with plots (not shown here): precice/calculix-adapter#9

We could still use OpenFOAM for both sides.

[davidscn]

One could also emulate one (the solid) participant using ASTE's replay mode. Works only for explicit coupling schemes at the moment, but this way, you can also check reading and writing data separately.

Have we shown that there is no similar issue in conjugate heat transfer cases? Let's take the flow-over-heated-plate tutorial and run a similar study, to have evidence (once you have time again).

Maybe the coupled Laplace OpenFOAM solver would be better suited: there, we even have an analytical solution and an error computation included.

[vidulejs]

Hello all!

Thank you for the suggestions! I am back on working on this issue.

@uekerman from now on I will be using the notation you suggested. Since solid-CalculiX does not support subcycling, all the cases I've tested have been subcycling only on the fluid-OpenFOAM side. I think we can conclude that the problem is on the OpenFOAM side when reading displacement. I've redrawn the displacement plot with dashed lines for the dt_F < Dt (fluid subcycling) cases.

solid subcycling

However, since I've only tested the perpendicular-flap case, I cannot say if this issue is with OpenFOAM subcycling in general or just with FSI. CalculiX adapter does not support subcycling and the Deal II linear elasticity solver does not either:

An error occurred in line <666> of file </home/dan/Desktop/preCICE/dealii/source/linear_elasticity/linear_elasticity.cc> in function
    void Linear_Elasticity::ElastoDynamics<dim>::run() [with int dim = 2]
The violated condition was: 
    std::abs(time.get_delta_t() - adapter.precice.getMaxTimeStepSize()) < 1e-10
Additional information: 
    This solver supports only constant time-step sizes.Configured time
    step size in deal.II parameter file: 0.005000. Time-window size from
    preCICE: 0.010000.

Which other solid solver would you suggest for subcycling dt_S?

• I could try solid-OpenFOAM, which also uses our custom-written solver, but then I can't isolate if the issue is in the OpenFOAM-adapter.
• I could try to implement subcycling in the Deal II linear_elasticity.cc solver.

---

I want to make two notes here, which I think might be relevant. I may be mistaken, but I believe the displacement read from preCICE may not necessarily be one-to-one with the OpenFOAM displacement. Admittedly, the residual of the OpenFOAM displacement solver is very low, but this could be something to consider.

smoothSolver:  Solving for cellDisplacementx, Initial residual = 5.96942e-17, Final residual = 1.30367e-17, No Iterations 2
smoothSolver:  Solving for cellDisplacementy, Initial residual = 3.2827e-17, Final residual = 5.20547e-18, No Iterations 2

Secondly, I agree that nCorrectors of pimpleFoam should not make a difference as long as it converges. But what if the boundary conditions change inside the pimpleFoam loop? Displacements are not affected either way and force is written only after the OpenFOAM solver steps; but the boundary conditions on the interface should be set according to received preCICE data at each solver iteration. I recently encountered this problem with an implicit solver, so I'm just thinking out loud here. However, this also does not explain why explicit coupling with subcycling does not exhibit this problem in comparison to the implicit coupling with subcycling.

• I need to draw the explicit coupling with subcycling cases, but IIRC they looked fine.
• Would you suggest I try a different fluid-OpenFOAM solver, e.g., icoFoam?

pimpleFoam log

Time = 4.949

PIMPLE: iteration 1
---[pimpleFoam] PIMPLE loop iteration: 1
---[pimpleFoam] (pimple.firstIter() || moveMeshOuterCorrectors) is true, entering mesh update
smoothSolver: Solving for cellDisplacementx, Initial residual = 5.96942e-17, Final residual = 1.30367e-17, No Iterations 2
smoothSolver: Solving for cellDisplacementy, Initial residual = 3.2827e-17, Final residual = 5.20547e-18, No Iterations 2
---[pimpleFoam] mesh.controlledUpdate()
---[pimpleFoam] (correctPhi) is true, calculate absolute flux: phi = mesh.Sf() & Uf();
DICPCG: Solving for pcorr, Initial residual = 1, Final residual = 0.000856733, No Iterations 65
DICPCG: Solving for pcorr, Initial residual = 0.00468318, Final residual = 9.79052e-09, No Iterations 82
time step continuity errors : sum local = 2.28403e-15, global = 2.30905e-17, cumulative = -1.14727e-08
--- [pimpleFoam] CorrectPhi
---[pimpleFoam] fvc::makeRelative(phi, U) Make flux relative to mesh motion
---[pimpleFoam] (pimple.momentumPredictor()) is true, solving momentum equation
smoothSolver: Solving for Ux, Initial residual = 3.98434e-05, Final residual = 1.80632e-07, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.000129526, Final residual = 8.26559e-07, No Iterations 2
---[pimpleFoam] fvOptions.correct(U)
---[pimpleFoam] (pimple.consistent()) is true, adjusting phiHbyA
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 0.000125338, Final residual = 1.01508e-07, No Iterations 56
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 1.77351e-06, Final residual = 9.56868e-09, No Iterations 54
---[pimpleFoam] (pimple.finalNonOrthogonalIter()) is true, update phi
time step continuity errors : sum local = 3.74345e-12, global = 3.80413e-14, cumulative = -1.14727e-08
---[pimpleFoam] Correcting Uf fvc::correctUf(Uf, U, phi)
---[pimpleFoam] Make phi relative to mesh motion fvc::makeRelative(phi, U)
---[pimpleFoam] (pimple.consistent()) is true, adjusting phiHbyA
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 2.26121e-05, Final residual = 2.16286e-08, No Iterations 57
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 2.07847e-07, Final residual = 9.9894e-09, No Iterations 17
---[pimpleFoam] (pimple.finalNonOrthogonalIter()) is true, update phi
time step continuity errors : sum local = 3.90797e-12, global = -1.65963e-13, cumulative = -1.14728e-08
---[pimpleFoam] Correcting Uf fvc::correctUf(Uf, U, phi)
---[pimpleFoam] Make phi relative to mesh motion fvc::makeRelative(phi, U)
---[pimpleFoam] (pimple.consistent()) is true, adjusting phiHbyA
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 1.07586e-05, Final residual = 1.07388e-08, No Iterations 56
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 9.1636e-08, Final residual = 9.76858e-09, No Iterations 10
---[pimpleFoam] (pimple.finalNonOrthogonalIter()) is true, update phi
time step continuity errors : sum local = 3.82154e-12, global = 3.17849e-14, cumulative = -1.14728e-08
---[pimpleFoam] Correcting Uf fvc::correctUf(Uf, U, phi)
---[pimpleFoam] Make phi relative to mesh motion fvc::makeRelative(phi, U)
---[pimpleFoam] (pimple.consistent()) is true, adjusting phiHbyA
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 5.77247e-06, Final residual = 9.37586e-09, No Iterations 55
---[pimpleFoam] In non-orthogonal pressure corrector loop
DICPCG: Solving for p, Initial residual = 5.20589e-08, Final residual = 9.09224e-09, No Iterations 3
---[pimpleFoam] (pimple.finalNonOrthogonalIter()) is true, update phi
time step continuity errors : sum local = 3.55694e-12, global = 1.13276e-13, cumulative = -1.14727e-08
---[pimpleFoam] Correcting Uf fvc::correctUf(Uf, U, phi)
---[pimpleFoam] Make phi relative to mesh motion fvc::makeRelative(phi, U)
ExecutionTime = 341.35 s ClockTime = 346 s

other tutorial cases

Else I can try a different tutorial case like CHT like @MakisH suggested. Then we can say that the problem is either with FSI or moving mesh in the OpenFOAM-adapter. Are you referring to a tutorial case with the Laplacian OpenFOAM solver, @davidscn ?

[davidscn]

I think the deal.II solver does in fact support subcycling. It does not support variable time-step sizes, so you would need to pre-calculate the subcycling factor.
The OpenFOAM solver is in this case https://github.com/precice/tutorials/tree/develop/partitioned-heat-conduction and if that's not enough you can also find it under https://journal.openfoam.com/index.php/ofj/article/view/88 (see the attached dataset).

[vidulejs]

I've conducted three additional experiments to try and pinpoint the cause of the issue. But I am not certain what the conclusions are. Please see the results and see if you have any insights.

1. Checked if perpendicular-flap solid subcycling also exhibits the behavior.
2. Tried partitioned-heat-conduction to see if OpenFOAM implicit subcycling is at fault.
3. Isolated the checkpointing mechanism in OpenFOAM-Adapter, more info below.

Firstly, I implemented a hard coded solid subcycling dt_S in the linear_elasticity.cc solver in Deal II, which doesn't seem to differ from the non-subcycling case as drastically. However, I read only the constant data from the start of the window, i.e., no waveform interpolation. That was an oversight on my part, and I can re-do the simulations if anyone thinks it will make a difference. Anyways, here is the displacement plot with the solid subcycling cases (dotted lines), which overlap and hide behind the reference non-subcycling pink line.

---

Secondly, I tried the implicit subcycling experiment on the partitioned-heat-conduction case. The Dirichlet left side is Nutils and the Neumann right side is OpenFOAM with our custom heatTransfer solver. Excuse the Paraview plots, I didn't bother to improve the readability much.

If we zoom in closer we can see that the Nutils coarse grained solution differs from the analytical reference. Also the lines are disconnected at the interface, but I assume this is normal and due to the separate discretizations. The OpenFOAM subcycling and non-subcycling cases seem to be identical at first glance.

However, if we zoom even further, there is a miniscule difference in the subcycling and non-subcycling cases. They are grouped separately just like in the perpendicular flap case! What does this tell us?

---

Checkpointing Sanity Check

Thirdly, Makis suggested we try to see if the OpenFOAM-Adapter checkpointing implementation is missing something. For example, if we do not correctly checkpoint all fields, then the reload would lead to a different state than before and we would accumulate errors over the iterations. In order to test if the checkpointing works correctly, we can pass the same data in each subsequent iteration and we should get the same results if acceleration is disabled. PreCICE throws an error message for "constant acceleration" with relaxation coefficient 0, but it does allow an arbitrarily small value like 1e-16. That configuration results in the first plot below.

  <coupling-scheme:parallel-implicit>
    <!-- ...  -->
    <max-iterations value="5" />
    <relative-convergence-measure limit="1e-12" data="Displacement" mesh="Solid-Mesh" />
    <relative-convergence-measure limit="1e-12" data="Force" mesh="Solid-Mesh" />
    <acceleration:constant>
      <relaxation value="1e-16"/>
    </acceleration:constant>
  </coupling-scheme:parallel-implicit>

These plots might look strange to you, but I had an easier time reading them when time is on the vertical axis, sort of like a spacetime diagram. In this plot you can see the displacement over the iterations and subcycling steps (here 10 subcycling steps), where the step number is color coded. You can see the checkpoint reloading. Here maximum iterations are 5 per time window. We receive the same data at each iteration, and since at the first iteration of each parallel implicit coupling window the data is not yet received and set to 0, the curves for a time window are constant (straight vertical lines).

The first plot is not interesting and the data exchange is not really valid. I implemented a mechanism in my fork of the OpenFOAM adapter, where the second iteration coupling data is cached and repeated for all subsequent iterations (modified readCouplingData function and additionally track iterations and substeps). Therefore, we repeat the data at the second iteration and see if there's any deviations between the iterations.

I plotted displacement over the iterations and also checked the log file and saw that there are no differences. I will post it here for completeness, but I didn't see any numerical differences in the iterations. That's good, it means that the checkpointing is working. But what else could be the issue then?

Here's the plot and log file for the full simulation with Dt=0.1 and dt_F=0.01.
perpendicular-flap_watchpoint.txt

For reference, a case with the standard OF-adapter and acceleration 1, where the points clearly differ between iterations.

[MakisH]

Notes from out meeting:

• Checkpointing for displacements looks indeed fine, but what about forces?
• deal.II subcycling (implemented for the experiment) works fine, so the issue must be in the OpenFOAM adapter side.
• In the partitioned heat conduction, there is a small (much smaller) difference between the subcycling and non-subcycling cases. Might be the same issue, but there might also be two different issues, both affecting the perpendicular-flap. Let's focus on the FSI case for now.

As a next step, let's get the same plots for forces or fields that are used to compute the force.

[vidulejs]

Good news and Bad news

• Good news
  The checkpointing is working.
• Bad news
  I have no idea why the implicit subcycling is converging to a different solution than non-subcycling.

Just like before, I performed the checkpointing sanity check (whether restoring checkpoints works as expected), but this time comparing the forces. The read data from preCICE (displacements) on the second iteration is repeated for the next iterations. I modified the OpenFOAM-adapter to cache and repeat the second iteration. If OpenFOAM-adapter checkpointing is working correctly, then the initial state at iterations 2, 3, 4 is identical and thus we should get identical outputs at all those iterations. That is the case, and I believe we can conclude that the checkpointing is working correctly.

  <coupling-scheme:parallel-implicit>
    <time-window-size value="0.1" />
    <max-time value="5" />

    <participants first="Fluid" second="Solid" />
    <exchange data="Force" mesh="Solid-Mesh" from="Fluid" to="Solid" substeps="true"/>
    <exchange data="Displacement" mesh="Solid-Mesh" from="Solid" to="Fluid" substeps="true"/>
    <max-iterations value="5" />
    <relative-convergence-measure limit="1e-12" data="Displacement" mesh="Solid-Mesh" />
    <relative-convergence-measure limit="1e-12" data="Force" mesh="Solid-Mesh" />

    <acceleration:constant>
      <relaxation value="1"/>
    </acceleration:constant>
  </coupling-scheme:parallel-implicit>
</precice-configuration>

Dt=0.1, dt_F=0.01

Full log: forces_comparison.log

Look on the right column, where the 2nd iteration is repeated for all subsequent iterations. The forces written by OpenFOAM are identical, as expected.

Time         Force_x         Repeat 2nd it  
---------- Time window 0.10 ----------
0.01         2183.930        | 2183.930       
0.02         598.770         | 598.770        
0.03         371.720         | 371.720        
0.04         325.926         | 325.926        
0.05         305.977         | 305.977        
0.06         293.806         | 293.806        
0.07         285.540         | 285.540        
0.08         279.506         | 279.506        
0.09         274.895         | 274.895        
0.10         271.263         | 271.263        
------------ Iteration 2 -------------
0.01         2183.930        | 2183.930       
0.02         598.770         | 598.770        
0.03         371.720         | 371.720        
0.04         325.926         | 325.926        
0.05         305.977         | 305.977        
0.06         293.806         | 293.806        
0.07         285.540         | 285.540        
0.08         279.506         | 279.506        
0.09         274.895         | 274.895        
0.10         271.263         | 271.263        
------------ Iteration 3 -------------
0.01         2180.250        | 2183.930       
0.02         597.161         | 598.770        
0.03         370.505         | 371.720        
0.04         324.735         | 325.926        
0.05         304.777         | 305.977        
0.06         292.594         | 293.806        
0.07         284.316         | 285.540        
0.08         278.272         | 279.506        
0.09         273.653         | 274.895        
0.10         270.014         | 271.263        
------------ Iteration 4 -------------
...

...
---------- Time window 5.00 ----------
4.91         243.847         | 244.169        
4.92         245.699         | 245.530        
4.93         245.566         | 245.539        
4.94         245.283         | 245.273        
4.95         245.026         | 245.020        
4.96         244.817         | 244.812        
4.97         244.653         | 244.650        
4.98         244.530         | 244.528        
4.99         244.441         | 244.440        
5.00         244.379         | 244.378        
------------ Iteration 2 -------------
4.91         239.751         | 240.058        
4.92         244.321         | 244.139        
4.93         244.976         | 244.937        
4.94         245.157         | 245.135        
4.95         245.214         | 245.196        
4.96         245.199         | 245.183        
4.97         245.130         | 245.115        
4.98         245.018         | 245.005        
4.99         244.870         | 244.857        
5.00         244.694         | 244.682        
------------ Iteration 3 -------------
4.91         239.790         | 240.058        
4.92         244.337         | 244.139        
4.93         244.989         | 244.937        
4.94         245.169         | 245.135        
4.95         245.225         | 245.196        
4.96         245.210         | 245.183        
4.97         245.141         | 245.115        
4.98         245.030         | 245.005        
4.99         244.882         | 244.857        
5.00         244.706         | 244.682        
------------ Iteration 4 -------------
4.91         239.790         | 240.058        
4.92         244.338         | 244.139        
4.93         244.989         | 244.937        
4.94         245.169         | 245.135        
4.95         245.226         | 245.196        
4.96         245.211         | 245.183        
4.97         245.141         | 245.115        
4.98         245.030         | 245.005        
4.99         244.882         | 244.857        
5.00         244.706         | 244.682

[MakisH]

The most interesting aspect in the above plots remains that even the mere presence of subcycling (two steps instead of one) drastically changes the results, while the number of steps does not seem to play as strong a role.

I ran a similar study with SU2 (with and without subcycling), and the values that SU2 converges to are essentially the same:

This still points out that the reason is technical.

Good to know that the checkpointing of displacements and forces is correct. This does not guarantee that the complete checkpointing is correct, but it should be enough.

You also wrote earlier that:

Note that explicit subcycling converges to the first group IIRC, so the problem is strictly with implicit subcycling.

Apparently, repeating a step causes side effects to OpenFOAM. As a next step, I think we need to go step by step through the code and think of potential side-effects. There are some areas marked with "if implicit coupling" or "if subcycling". I would start from there.