Force calculation depends on number of ranks

We currently face an issue with the force functionObject: depending on the total number of ranks in parallel runs the obtained forces differ.

We use pimpleFoam in combination with dynamic meshes to perform our simulations. The issue appeared originally in a rather complex application case and the best guess for the rank dependency there was the FV mesh motion solver itself. The behavior can be reproduced using the tutorial case movingCone of the latest release v2006. The tutorial includes a uniform mesh motion in x-direction, which is solved by the velocityComponentLaplacian solver of the dynamicMotionSolverFvMesh dynamic mesh. The force is evaluated at the moving interface with the libforces.so library, as usual. Since wedge type patches cannot be partitioned as expected with the latest OpenFOAM version, the out-of-plane patches have been switched from type wedge to the type patch.

We performed the simulation on a various number of ranks and looked at the resulting force at the moving interface. Comparing a serial run with a simulation on 5 x 5 = 25 ranks shows on the first time level a deviation around O(10^{-5}). At the final time level 3e-3s, a deviation of around O(10^{-6}) can be observed. Both measures are absolute deviations, so that the accuracy is clearly below machine accuracy. The force diff for all time steps is also attached diff_serial_25_procs.log. diff_serial_25_procs.log

In a second step, the mesh movement has been removed. Since the non-zero fluid velocity is entirely triggered by the moving obstacle, a uniform inlet flow velocity of strength 1 is placed on the left side and an outlet is placed on the right side of the domain. Again, the force for a serial run and a 5 x 5 = 25 rank parallel run is considered. The deviation of the calculated force values is initially around O(10^{-6}) and at the final time level around O(10^{-9}). Thus, the differences do not disappear in case the mesh motion is deactivated. The complete log is also attached diff_serial_25_procs_no_motion.log. diff_serial_25_procs_no_motion.log

I'm a bit puzzled about the reason. Although the attached log files contain only the force.dat file comparison, I also checked the console output for more than the default write precision in order to exclude round-off errors while creating the force.dat file. Also, selecting very (too) strict residual controls doesn't affect the observed differences.

Are there any hints for these differences?

Background: We are performing black-box coupled FSI simulations with OpenFOAM using the openfoam-adapter and these differences can lead to severe differences for coupled simulations so that the entire coupled system converges to different states. In contrast to the tutorial setup, as described above, our original problem setup applies the displacement based FV mesh motion solver with time varying boundary conditions.