precice / openfoam-adapter

OpenFOAM-preCICE adapter
https://precice.org/adapter-openfoam-overview.html
GNU General Public License v3.0
135 stars 83 forks source link

Checkpointing the old fields and mesh #43

Closed derekrisseeuw closed 5 years ago

derekrisseeuw commented 6 years ago

Field checkpointing OpenFOAM uses various timeschemes which can use the old field and old old field. <field>.oldTime() and <field>.oldTime().oldTime() The correct value of these old fields is required for correct timestepping.

When writing the checkpoint, pointers to the old fields are actually saved, but when reading the checkpoint these are not automatically updated. Therefore, this must be performed explicitly in the adapter. In the Adapter.C I added the following lines for all fields. (here for the volScalarFields, line 1260)

        int nOldTimes(volScalarFields_.at(i)->nOldTimes());
        if (nOldTimes >= 1)
        {
            volScalarFields_.at(i)->oldTime() == volScalarFieldCopies_.at(i)->oldTime();        
        }
        if (nOldTimes == 2)
        {
            volScalarFields_.at(i)->oldTime().oldTime() == volScalarFieldCopies_.at(i)->oldTime().oldTime();
        }

In the openfoam sourcecode there is one mention of the oldTime().oldTime().oldTime(), but this is used for a structural solver. Therefore I only include two 'layers'.

Mesh checkpointing A similar approach can be taken for the mesh update. When the function movepoint is called on line 764, the pointer to the mesh.oldPoints() is reset to the current mesh.

const_cast<fvMesh&>(mesh_).movePoints(meshPoints_); This leads to a reversion of the points of the mesh in time, which is not physical. therefore I added the pointer to the old mesh in a checkpointed field, which is read later in the reloadMeshPoints function:

void preciceAdapter::Adapter::storeMeshPoints()
{
    DEBUG(adapterInfo("Storing mesh points..."));
    // TODO: In foam-extend, we would need "allPoints()". Check if this gives the same data.
    meshPoints_ = mesh_.points();
    oldMeshPoints_ = mesh_.oldPoints();
    DEBUG(adapterInfo("Stored mesh points."));
}

void preciceAdapter::Adapter::reloadMeshPoints()
{
    // In Foam::polyMesh::movePoints.
    // TODO: This function overwrites the pointer to the old mesh. Therefore, if you revert the mesh, the oldpointer will be set to the points, which are the new values. 
    DEBUG(adapterInfo("Moving mesh points to their previous locations..."));
    const_cast<fvMesh&>(mesh_).movePoints(meshPoints_);
    const_cast<pointField&>(mesh_.oldPoints()) = oldMeshPoints_;
    DEBUG(adapterInfo("Moved mesh points to their previous locations."));
}

Please leave any feedback on the validity of this approach. I am curious to other users experience with this issue.

MakisH commented 6 years ago

This is a very interesting issue, both because it is pretty much non-obvious and because it helped @derekrisseeuw to get a stable Turek&Hron FSI3 benchmark simulation.

This is relevant to #7 (FSI module), but also to CHT. Moreover, all of these are changes required in the Adapter.C and not in the modules.

derekrisseeuw commented 6 years ago

One comment on the mesh part of this issue.

The mesh deformation method only stores the oldPoints. Therefore checkpointing it is not required as I see it now: Every sub-iteration the fields must be returned to their checkpointed fields. For the mesh this includes that the .points() reference and .oldpoints() are updated. However, for the next subiteration the mesh is moved again to a new estimate, which means that the checkpointed .points() is set to the .oldpoints() field. So the checkpointed .oldpoints() is not actually used.

I would like feedback if my reasoning is correct.

The reason that it might be convenient to neglect checkpointing the mesh .oldpoints() (other than computations effort) is that the native mesh deformations in openfoam are not very good at handling large displacements, and poor at handling rotation. These cause problems in typical FSI problems. RBF mesh deformation would be a superior option.

The performance of the openfoam mesh deformation is actually improved if the mesh is not updated, because then the laplacian equation is solved every sub-iteration on the new mesh which in term diffuses more of the accumulated motion.

uekerman commented 6 years ago

I might misunderstand some of the technical details, but your reasoning sounds good. If .oldpoints() is never used, there is no reason to checkpoint it. But why does this now allow for RBF mesh deformation and checkpointing .oldpoints() does not?

And how is this connected to:

The performance of the openfoam mesh deformation is actually improved if the mesh is not updated, because then the laplacian equation is solved every sub-iteration on the new mesh which in term diffuses more of the accumulated motion.

?

derekrisseeuw commented 6 years ago

Hi Benjamin,

Good to hear that you don't think this field is necessary.

On the part of the RBF deformation you misunderstood me. I mean to say that the mesh deformation method in openfoam are quite poor, especially at handling rotation. Therefore using RBF deformation would be better.

I need to check, but I thought that Not updating the oldpoints actually improved the capability of the laplacian mesh deformation somewhat. So this might be a reason to leave them out of the updating.

MakisH commented 5 years ago

@derekrisseeuw do you consider this issue closed by now? I understand that there is still an issue with checkpointing the mesh volumes, but maybe we could now isolate it in a separate issue.

derekrisseeuw commented 5 years ago

@MakisH we can close the issue.

The 'issue' regarding the mesh volumes is not relevant for coupled simulations with the same timestep, but only when subcycling is considered. Even more so, the volume fields should NOT be altered when NOT subcycling.