Fix different dimensions error / Checkpointing issue #158

[vidulejs]

In issue #158 it was discovered that using the kOmegaSST turbulence model together with implicit coupling causes a dimension mismatch error in OpenFOAM, when starting the second timestep / reading the checkpointed fields. The perpetrator was found to be the field 'grad(U)', which was created in the first timestep during the execution of the kOmegaSST turbulence model. The field is created only temporarily and afterwards is not stored in the OpenFOAM registry / fields database, thus causing an error when trying to read it in the next timesteps. A failing case was posted by a user:

Hi makis, I have made a case related to this problem, you may run this case to catch the error.
3000_40Mpa_orig.zip

Originally posted by @c7888026 in #158 (comment)

However the problem is more general since other functions can create temporary checkpointable fields, such as 'yPsi' used for 'wallDist' calculation as posted by user:

... The log is attached above , I found that the variable yPSi that is calculated in the fvSolution and required by the method of the wall distance in the fvSchemes is the problem, so it behaves as the gradU as you described and not passed correctly through the check point . I changed the wallDist method to meshWave as @c7888026 did in his case , and it worked fine even with K-omega sst and in parallel setting using mpirun , so it is matter of variables that are not passing well ( may be ) through the checkpoint step.

Originally posted by @KariimAhmed in #158 (comment)

I could not yet replicate this issue, because for some reason OpenFOAM is still selecting 'meshWave' method for 'wallDist' instead of 'Poisson', which I specified in fvSchemes. Also there are no preCICE tutorials without 'meshWave'.

The general solution I found is to check which fields are stored in the OpenFOAM database and which ones are checkpointed, then remove those checkpointed fields which are no longer used by OpenFOAM. This new function 'pruneCheckpointedFields()' is executed before the readCheckpoint() function, where the copy was loaded *field = *fieldCopy;. I am not sure, but I assume the field pointer was set to another field by OpenFOAM, otherwise it would be a memory error and not a dimensions error.

void preciceAdapter::Adapter::pruneCheckpointedFields()
{
    // Check if checkpointed fields exist in OpenFOAM database
    // If not, remove them from the checkpointed fields list
	// [...]
}

void preciceAdapter::Adapter::readCheckpoint()
{
	// [...]
	
    // Reload all the fields of type volVectorField
    for (const auto& kv : volVectorFields_)
    {
        const std::string& key = kv.first;
        volVectorField* field = kv.second;
        volVectorField* fieldCopy = volVectorFieldCopies_.at(key);

        // Load the volume field
        *field = *fieldCopy;

        int nOldTimes = field->nOldTimes();
        if (nOldTimes >= 1)
        {
            field->oldTime() = fieldCopy->oldTime();
        }
        if (nOldTimes == 2)
        {
            field->oldTime().oldTime() = fieldCopy->oldTime().oldTime();
        }
    }
	// [...]
}

Originally I thought it was an issue with storing the field pointers in a vector and using the wrong index to access this vector, so I tried to change the storage to an unordered_map and use the field name as the key:

    //- Checkpointed volScalarField fields
    std::unordered_map<std::string, Foam::volScalarField*> volScalarFields_;

However, the problem remained unchanged. Thus I can revert the storage back to vectors for simplicity and access the name of the field from the field object when accessing it from the vector.

Also, the new function 'pruneCheckpointedFields()' uses a macro and is not the prettiest, so I am open for a refactor. I included the '' header in 'Adapter.H' to search the fields database, but it can be done without it. When a field is removed from checkpointing, a new debug message is logged:

[2] ---[preciceAdapter] [DEBUG] Removed volTensorField : grad(U) from the checkpointed fields list.

TODO list:

• Revert the fields storage from vectors to unordered_map ?
• Can be done without using macro (it's fine as it is - @MakisH)
• Can be done without including header
• Replicate and test issue with a different wallDist method than meshWave (e.g. Poisson)
• I updated the documentation in docs/ -> N/A
• I added a changelog entry in changelog-entries/ (create directory if missing)

[KariimAhmed]

Thank you Daniels for this work. Do you suggest trying your modified Adapter.C and .H for my Poisson wall distance issue now ? or there still other modifications to do ?
since i will run it on a cluster so i have to be sure before changing the original Precice installed there.

[vidulejs]

Hello, Karim!
As I wrote before I was not able to test the bugfix yet with the 'yPsi' Poisson wall distance issue. I tried to modify the fvSchemes wallDist method to Poisson, but for some reason OpenFOAM still selects meshWave. I did not yet create a seperate test case for this specifically. I might try to do this tonight.

If you wish, you may try recompiling with this code and try to run your case. But there may be changes in the final pull request.

[MakisH]

Already good that it solves the original issue. If you cannot replicate the second one, I would not make that a road-blocker. We will figure out eventually if that is a separate issue.

[MakisH]

The code looks clean, and the provided test case seems to work, at least for 10 time windows. In that case, the grad(U) is removed.

I ran the flow-over-heated-plate and the perpendicular-flap tutorials and they also seem to work.

@davidscn do you want to check anything else?

[vidulejs]

Hello again and sorry for the delay! I briefly had an issue with my system and couldn't run preCICE, but now I resolved it and could test the code. Thanks for your comments David @davidscn ! I added the STL algorithm header in Adapter.h and made the changes you suggested.

[davidscn]

Looks good, I haven't tested it though.

[MakisH]

After some iterations of improvements, I think this is now ready to merge. I tested a previous state of this and it worked, and I ran the system tests as well (which do not trigger the actual issue). @vidulejs tested the latest state as well.

I am triggering the system tests once more for completeness, and merging.

[MakisH]

I forgot that the system tests cannot currently run in forks (see precice/tutorials#580). I still think that this is fine to merge. I will trigger the system tests in develop.