Mesh points & volumes checkpointing

[vidulejs]

Overview

• Removed comments and vectors to mesh fields which don't exist
• Fixed meshOldPoints checkpointing. Previously mesh would move from incorrect interim state to meshPoints, instead of meshOldPoints (leads to incorrect mesh flux).
• Mesh volume checkpointing (needed for higher order time derivative schemes)

Detailed Changes

Cleanup and refactor

This relates to the bulk of line removals in this branch.

Just like checkpointed fields of physical variables, there were different types of checkpointed mesh fields, however they were unnecessary and empty. There were no meshVolVectorFields_ or meshSurfaceVectorFields_ (none that needed checkpointing) . The only checkpointed mesh field is face motion flux mesh_.phi() of type surfaceScalarField. As pointed to by in a comment in setupMeshCheckpointing():

// The other mesh Fields:
// C
// Cf
// Sf
// magSf
// delta
// are updated by the function fvMesh::movePoints. Only the meshPhi needs checkpointing.

Mesh Points Checkpointing

• Store meshPoints_ and meshOldPoints_ as pointers instead of objects for consistency with the rest of the codebase. These are actually pointers to the copies managed by the adapter.
• The main logic of restoring mesh checkpoints is like this:

> // const_cast<pointField&>(mesh_.points()) = oldMeshPoints_;
> 
> // Reload mesh points
> const_cast(mesh_).movePoints(*meshPoints_)

Currently the first line is commented out and probably forgotten by the previous developer Derek (and I do not understand the meaning of this comment).
> // Switch oldpoints on for pure physics. (is this required?). Switch off for better mesh deformation capabilities?
> // const_cast<pointField&>(mesh_.points()) = oldMeshPoints_;

I had independently come to the conclusion that this needs to be enabled, when looking at the fvMesh.movePoints() method. Perhaps we need a longer discussion on this.

EDIT*

• The main logic of restoring mesh checkpoints is like this:

    // Reload mesh points
    const_cast<Foam::fvMesh&>(mesh_).movePoints(*meshPoints_);

    // polyMesh.movePoints will only update oldPoints
    // if (curMotionTimeIndex_ != time().timeIndex())
    // I.e., first implicit iteration
    const_cast<pointField&>(mesh_.oldPoints()) = *meshOldPoints_;

mesh_.oldPoints() will not get updated inside of fvMesh.movePoints() -> polyMesh.movePoints() in implicit iterations. See the discussion in the comments. Related issue - FSI: Subcycling with implicit coupling #58.

Mesh Volumes Checkpointing

For solvers that use higher-order time derivatives (e.g., Euler ddtScheme), checkpoint cell volumes (V, V0, V00).

• Implement new methods setupMeshVolCheckpointing, readMeshVolCheckpoint, writeMeshVolCheckpoint
• A counter volumeCheckpointCounter is used to get V0, V00, when they become available. Otherwise mesh_.V00() access will throw an error.
• Only store copies of the volume fields (meshVolFieldCopies_). The references given by mesh_.V0() may be invalidated by OpenFOAM, so we cannot store pointers to them. Instead directly call mesh_.V0() and so on.
• readMeshVolCheckpoint() will const_cast the const reference and overwrite the memory, e.g.,

const_cast<volScalarField::Internal&>(mesh_.V()) = *(meshVolFieldCopies_.at(0));

TODO list:

• I updated the documentation in docs/
• I added a changelog entry in changelog-entries/ (create directory if missing)

[MakisH]

Looks good overall, thank you for all the cleanup effort! This is an important step.

Please also:

• Add a changelog entry
• Update the PR description with the motivation and linking to related issues/PRs.

I am triggering the system tests mainly as a routine check, but we will (indeed) need to update the reference results after merging.

[MakisH]

Please remind me, it's been a while: What is the reason that here (and in addMeshCheckpointField()) we only need to store the copies, but for surfaceScalarFields we needed two objects (main and copy)?

Can we get rid of one of the two in the other fields as well?

[vidulejs]

We do not need to explicitly store pointers to the fvMesh.points() and fvMesh.oldPoints() because we can just call these methods when we need to access them. The logic for these fields is handled specifically as opposed to generically. Same goes for the mesh volumes.

We store pointers to the OpenFOAM field objects and our copies for convenience. In setupCheckpointing() we loop over the OpenFOAM object registry to get all the fields, then in addMeshCheckpointField() we store the pointers in *Fields_ and *FieldCopies_ vectors for each type of field.

The mesh face motion flux meshPhi is not indexed in the object registry, therefore it is added manually by calling it's method.

void preciceAdapter::Adapter::setupMeshCheckpointing()
{
    DEBUG(adapterInfo("Creating a list of the mesh checkpointed fields..."));
    // Add meshPhi (Face motion flux)
    addMeshCheckpointField(const_cast<surfaceScalarField&>(mesh_.phi()));

    DEBUG(adapterInfo("Added " + mesh_.phi().name() + " to the list of checkpointed fields."));
}

Afterwards it is accessed from meshSurfaceScalarFields_ and meshSurfaceScalarFieldCopies_, however these vectors contain only meshPhi and no other surfaceScalarFields.

Perhaps we can refactor these parts for consistency with the rest of the codebase.

[MakisH]

Suggested change

	// First timestep oldPoints doesn't exist yet
	// First timestep, oldPoints doesn't exist yet

But I don't understand: if it is an issue that it does not yet exist, how to copy that?

[MakisH]

@vidulejs you probably also want to change this one to some kind of more meaningful flag. Unless you think that a counter makes more sense here.

How does the code behave for >=3?

[vidulejs]

Thanks for the discussion @MakisH! I have now had a look with fresh eyes, it's been a while since I made the original commits.

There's something I need to bring attention to. Although this PR is not about implicit subcycling (WIP PR #370), the two are closely related. Bear with me for a moment as I try to explain.

The old issue on FSI: Subcycling with implicit coupling #58 mentions the private member variable curTimeIndex_ of fvMesh, which will prevent the storeOldVol() method from executing in subsequent implicit iterations.

 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());
  }
  // (...)
}

When we restore a checkpoint, we restore the Time object, which is in the past and less than curTImeIndex_. The fvMesh class is responsible for the cell volumes V, V0, V00.
Similarly, the polyMesh class is responsible for points and oldPoints, and it also performs such a check:

void Foam::polyMesh::movePoints(const pointField& newPoints)
{
    // (...)

    // Pick up old points
    if (curMotionTimeIndex_ != time().timeIndex())
    {

        // Mesh motion in the new time step
        oldPointsPtr_.reset(nullptr);
        oldPointsPtr_.reset(new pointField(points_));
        curMotionTimeIndex_ = time().timeIndex();
    }

    points_ = newPoints;

    // (...)
}

Now, what does this mean? In implicit coupling w/o subcycling this is inconsequential (!) We perform a solver step spanning the whole window. In the first coupling iteration curTimeIndex_ is updated and we store the old volumes and old points corresponding to the start of the window. The new points are updated irregardless. In subsequent implicit coupling iterations the state of old volumes and old points is correct, i.e., corresponds to the start of the window. Also mesh face flux meshPhi is correct.

However, for implicit coupling w/ subcycling, polyMesh::movePoints() will store the old points from the last (unconverged/interim) subcycling step .

We perform multiple subcycling steps within the window, the old points will always update since curMotionTimeIndex_ != time().timeIndex().

On the next implicit coupling iteration, the old points will be from the last subcycling step - not from the start of the window! The mesh face motion flux meshPhi will be incorrect and introduce spurious motion!

Solution

Therefore, instead of relying on fvMesh.movePoints() to store the old points (sometimes correctly), we have to make sure they are always updated correctly like this.

void preciceAdapter::Adapter::reloadMeshPoints()
{
    if (!mesh_.moving())
    {
        DEBUG(adapterInfo("Mesh points not moved as the mesh is not moving"));
        return;
    }

-    // fvMesh.movePoints overwrites the pointer to the oldPoints
-    const_cast<pointField&>(mesh_.points()) = *meshOldPoints_;

     // Reload mesh points
     const_cast<Foam::fvMesh&>(mesh_).movePoints(*meshPoints_);

+    // polyMesh.movePoints will only update oldPoints
+    // if (curMotionTimeIndex_ != time().timeIndex())
+    const_cast<pointField&>(mesh_.oldPoints()) = *meshOldPoints_;

     readMeshCheckpoint();
     
    DEBUG(adapterInfo("Moved mesh points to their previous locations."));

    readMeshVolCheckpoint();
}