FSI: Subcycling with implicit coupling

[MakisH]

Additional fields should be checkpointed, for stability. Meanwhile, these fields should not be checkpointed when we don't subcycle.

See the method storeMeshPoints() in Adapter.C for more. @derekrisseeuw could you maybe provide a few more details?

Relevant to #7.

[derekrisseeuw]

In the openfoam simulation 3 parts need updating. A distinction is being made between fluid fields and mesh fields, which both need checkpointing. The Time is the third element which needs updating. Fluid fields are all the Fields, such as the velocity which is a volVectorField. These can include old fields.

For the mesh only the points and oldpoints (location of the mesh) and the mesh flux area read. This is

The fields are stored and read when required by preCICE by evaluating the bool isWriteCheckpointRequired and isReadCheckpointRequired. If so, the respective functions writeCheckpoint or readCheckpoint are called which follow similar procedures but one is for writing and the other for reading.

In writeCheckpoint, first the time is checkpointed, next the mesh is written. This includes writing the checkpoint, and the meshPhi. This is valid for both implicit simulations with the same timestep as simulations with subcycling. The catch is found in the checkpointing of the mesh (old)volumes: V, V0 (and V00 dependent on the time scheme). Without subcycling these are not to be touched because openFOAM computes them correctly for the single timestep. However, for subcycling they should be set back to the original checkpointed value, see the comments in the Adapter.C for this.

[MakisH]

In other words (correct me if I am wrong):

• When there is no coupling at all, the solver updates the old values at every time step, doing something like V00 = V0, V0 = V, V = compute().

• Without subcycling, the solver only executes one step before having to return to the previous coupling time and reload the checkpoint, so it has executed the same update as above (once). When reloading the mesh, we should also return V to its value before compute(). The important question is: when does OpenFOAM update the V0 and V00? If in the meantime they have not been updated yet at the time we reload the checkpoint, then we should leave them as is. But why not store them and reload them eitherways? (apart from any performance reasons)

• When subcycling, the solver executes multiple small time steps until it has to reload a checkpoint, so it has executed multiple updates. If we do not restore the V0 and V00 as well, there is a problem.

• In explicit coupling, it doesn't matter (there is no checkpointing).

By the way, this is related to #55.

[MakisH]

After discussing this with @derekrisseeuw and trying to understand the problem a bit better (what I am writing above for "without subcycling" is not exactly the case), a few notes:

1. The old volumes are actually updated inside the fvMesh::movePoints(), which calls the fvMesh::storeOldVol(), only if the time has been incremented:

   Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
   {
    // Grab old time volumes if the time has been incremented
    // This will update V0, V00
    if (curTimeIndex_ < time().timeIndex())
    {
        storeOldVol(V());
    }
    ...
   }

2. The fvMesh::storeOldVol() itself then updates the old volumes to point to a range of time from

   void Foam::fvMesh::storeOldVol(const scalarField& V)
   {
       if (curTimeIndex_ < time().timeIndex())
       {
           if (debug)
           {
               InfoInFunction
                   << " Storing old time volumes since from time " << curTimeIndex_
                   << " and time now " << time().timeIndex()
                   << " V:" << V.size()
                   << endl;
           }
   
           if (V00Ptr_ && V0Ptr_)
           {
               // Copy V0 into V00 storage
               *V00Ptr_ = *V0Ptr_;
           }
   
           if (V0Ptr_)
           {
               // Copy V into V0 storage
               V0Ptr_->scalarField::operator=(V);
           }
           else
           {
               ...
           }
   
           curTimeIndex_ = time().timeIndex();
   
           ...
       }
   }

3. When we reload a checkpoint, we reload the time, and then we call movePoints(). This is the sequence of events:

   void preciceAdapter::Adapter::readCheckpoint()
   {
       // Reload the runTime
       reloadCheckpointTime();
       // In there: const_cast<Time&>(runTime_).setTime(couplingIterationTimeValue_, couplingIterationTimeIndex_);
   
       // Reload the meshPoints (if FSI is enabled)
       if (FSIenabled_)
       {
           reloadMeshPoints();
           // In there: const_cast<fvMesh&>(mesh_).movePoints(meshPoints_);
           // and: readMeshCheckpoint();
           // which restores the old values
       }
   
       ... // Restore fields
   }

The problem (or one of the problems) is that, when we reset the time value and index, we do not reload the curTimeIndex_, which is a private member of fvMesh, and set to a value larger or equal to the (global/normal) time index. Under these conditions, the storeOldVol() does nothing and therefore does not update the old volumes.

One workaround could be that we implement our own storeOldVol(), which does not depend on the curTimeIndex_.

See also the method preciceAdapter::Adapter::setupMeshVolCheckpointing().

I will add more details eventually. For now, I would assign a lower priority to this, as its effect is mainly performance-related (the fact that we cannot use subcycling).

[derekrisseeuw]

This is an accurate description @MakisH, thanks for taking the time to write it all out.

For additional information: The old volumes are used in the time scheme. Therefore, using the incorrect (old) volumes would induce errors in the time integration.

As a third option we also discussed creating an additional method of the fvMesh class, which would set back the curTimeIndex_. Then we could create this method somewhere in the adapter and call it whenever the global time checkpoint is read. It would reset curTimeIndex_ to timeIndex - 1, to assure that the storeOldVols would be called by movePoints.

Of course the mesh volume fields must then be checkpointed too, as you mention by the function setupMeshVolCheckpointing()

[precice-bot]

This issue has been mentioned on preCICE Forum on Discourse. There might be relevant details there:

https://precice.discourse.group/t/openfoam-calculix-fsi-diverging-sim-never-ends/645/21

[precice-bot]

This issue has been mentioned on preCICE Forum on Discourse. There might be relevant details there:

https://precice.discourse.group/t/strange-phenomenon-during-multiphase-fsi/802/20