OpenFOAM throwing mesh errors from packs

[acatwithacomputer]

I'm using the .STL of a sphere pack with snappyHexMesh to generate a mesh for OpenFOAM. I originally had snapping on, and verified with checkMesh that there were no disconnected regions and that the mesh was fine. However, as soon as I try to run the case, I get this error regarding the mesh:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 10
Exec   : simpleFoamMod -parallel
Date   : Feb 26 2025
Time   : 14:38:45
Host   : "brittany-XPS-13-9370"
PID    : 71013
I/O    : uncollated
Case   : /home/brittany/sedTransp/seq10/monoDense
nProcs : 4
Slaves : 
3
(
"brittany-XPS-13-9370.71014"
"brittany-XPS-13-9370.71015"
"brittany-XPS-13-9370.71016"
)

Pstream initialised with:
    floatTransfer      : 0
    nProcsSimpleSum    : 0
    commsType          : nonBlocking
    polling iterations : 0
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

[2] processorPolyPatch::calcGeometry : Writing my 1443 faces to "/home/brittany/sedTransp/seq10/monoDense/processor2/procBoundary2to0throughyL_faces.obj"
[3] processorPolyPatch::calcGeometry : Writing my 1098 faces to "/home/brittany/sedTransp/seq10/monoDense/processor3/procBoundary3to1throughyL_faces.obj"
[0] processorPolyPatch::calcGeometry : Writing my 1813 faces to "/home/brittany/sedTransp/seq10/monoDense/processor0/procBoundary0to1throughx0_faces.obj"
[1] processorPolyPatch::calcGeometry : Writing my 1813 faces to "/home/brittany/sedTransp/seq10/monoDense/processor1/procBoundary1to0throughxL_faces.obj"
[3] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary3to1throughyL_faceCentresConnections.obj
[2] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary2to0throughyL_faceCentresConnections.obj
[1] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary1to0throughxL_faceCentresConnections.obj
[0] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary0to1throughx0_faceCentresConnections.obj
[3] 
[3] 
[3] --> FOAM FATAL ERROR: 
[3] face 57 area does not match neighbour by 0.0124883% -- possible face ordering problem.
patch:procBoundary3to1throughyL my area:0.000368015 neighbour area:0.000368061 matching tolerance:2.87081e-08
Mesh face:231314 vertices:4((1.36457 2.12138 0.715095) (1.33914 2.12385 0.699719) (1.33857 2.1272 0.725713) (1.35883 2.12728 0.71873))
If you are certain your matching is correct you can increase the 'matchTolerance' setting in the patch dictionary in the boundary file.
Rerun with processor debug flag set for more information.
[3] 
[3]     From function virtual void Foam::processorPolyPatch::calcGeometry(Foam::PstreamBuffers&)
[3]     in file meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C at line 269.
[3] 
FOAM parallel run exiting
[3] 
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 3 in communicator MPI COMMUNICATOR 3 SPLIT FROM 0
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
[brittany-XPS-13-9370:71009] PMIX ERROR: UNREACHABLE in file ../../../src/server/pmix_server.c at line 2193
[2] 
[2] 
[2] --> FOAM FATAL ERROR: 
[2] face 82 area does not match neighbour by 0.00551079% -- possible face ordering problem.
patch:procBoundary2to0throughyL my area:0.000330345 neighbour area:0.000330363 matching tolerance:1.78836e-08
Mesh face:235236 vertices:4((0.493972 2.12737 1.93407) (0.509532 2.12517 1.93394) (0.509855 2.12525 1.91275) (0.494651 2.12761 1.91277))
If you are certain your matching is correct you can increase the 'matchTolerance' setting in the patch dictionary in the boundary file.
Rerun with processor debug flag set for more information.
[2] 
[2]     From function virtual void Foam::processorPolyPatch::calcGeometry(Foam::PstreamBuffers&)
[2]     in file meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C at line 269.
[2] 
FOAM parallel run exiting
[2] 
[brittany-XPS-13-9370:71009] 1 more process has sent help message help-mpi-api.txt / mpi-abort
[brittany-XPS-13-9370:71009] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

In order to rule out anything wrong with the case, I've successfully run the case with:

1. A single sphere in the middle of the domain (rule out anything other than the mesh)
2. A single sphere, but broken into octants at each of the corners of the domain; stl was generated with this function:

def makeCornerTestSpheres(length,r):
    lx = [0, length[0]]
    ly = [0, length[1]]
    lz = [0, length[2]]
    with pygmsh.geo.Geometry() as geom:
        for i in range(0,2):
            for j in range(0,2):
                for k in range(0,2):
                    pos = [lx[k],ly[j],lz[i]]
                    geom.add_ball(pos,r,mesh_size=r * 0.1)
        mesh = geom.generate_mesh()
    mesh.write("./corners.stl")

I then moved a 100 sphere pack to sycamore (too large to run locally) with the case only to find it fail. I also toggled off the snapping and am still getting the exact same issue. I then generated a smaller pack with 11 spheres (small enough to run locally). It runs on a single processor (snapping on), but then gives the error above when I try to run it in parallel. This is also the same type of error that the 100 sphere pack gives on sycamore (face x area does not match neighbour by 0.0????% -- possible face ordering problem.). I tried googling the error, and it always seems to be about something else. I'm not sure if this is a purely openFoam problem or something about how the meshes from this code is interacting with OpenFoam.

[acatwithacomputer]

Small update: The processor configuration that got the error above was 2x2x1. I had noticed that the points it prints out don't seem to be internal boundaries, so I reran it with a processor configuration 4x1x1 to obtain this error:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 10
Exec   : simpleFoamMod -parallel
Date   : Feb 26 2025
Time   : 15:16:53
Host   : "brittany-XPS-13-9370"
PID    : 74558
I/O    : uncollated
Case   : /home/brittany/sedTransp/seq10/monoDense
nProcs : 4
Slaves : 
3
(
"brittany-XPS-13-9370.74559"
"brittany-XPS-13-9370.74560"
"brittany-XPS-13-9370.74561"
)

Pstream initialised with:
    floatTransfer      : 0
    nProcsSimpleSum    : 0
    commsType          : nonBlocking
    polling iterations : 0
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

[3] processorPolyPatch::calcGeometry : Writing my 3176 faces to "/home/brittany/sedTransp/seq10/monoDense/processor3/procBoundary3to0throughxL_faces.obj"
[0] processorPolyPatch::calcGeometry : Writing my 3176 faces to "/home/brittany/sedTransp/seq10/monoDense/processor0/procBoundary0to3throughx0_faces.obj"
[3] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary3to0throughxL_faceCentresConnections.obj
[0] processorPolyPatch::calcGeometry : Dumping lines between corresponding face centres to procBoundary0to3throughx0_faceCentresConnections.obj
[3] 
[3] 
[3] --> FOAM FATAL ERROR: 
[3] face 207 area does not match neighbour by 0.0208105% -- possible face ordering problem.
patch:procBoundary3to0throughxL my area:0.000203829 neighbour area:0.000203787 matching tolerance:3.24415e-08
Mesh face:231603 vertices:4((2.13247 1.98177 0.265764) (2.12521 1.95518 0.277314) (2.12493 1.95453 0.26394) (2.1257 1.97587 0.265102))
If you are certain your matching is correct you can increase the 'matchTolerance' setting in the patch dictionary in the boundary file.
Rerun with processor debug flag set for more information.
[3] 
[3]     From function virtual void Foam::processorPolyPatch::calcGeometry(Foam::PstreamBuffers&)
[3]     in file meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C at line 269.
[3] 
FOAM parallel run exiting
[3] 
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 3 in communicator MPI COMMUNICATOR 3 SPLIT FROM 0
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

My guess is that it may have something to do with the periodic boundaries going from one processor to the next.

ps: sorry, I got confused about how the button interacts with the comment button (turns out it doesn't..?) and accidentally closed it.

[tmweigand]

Yes this is due to how snappyHexMesh discretizes in parallel. It is not an issue with the stl file.

In this work, I used this sphere pack code, had spheres crossing boundaries and did not have this issue. I was using OpenFoam v2406. Here is how I was generating the domains.

I suggest switching to that version and seeing what is different in your process vs what I used previously.

Edit: @acatwithacomputer I would first look at my workflow and parameters and then switch to different version if that did not solve the issue.

[acatwithacomputer]

Saw the edit. Ok, I tried this:

echo "Sourcing OF-10 bashrc"
source ~/OpenFOAM/OpenFOAM-10/etc/bashrc

blockMesh > log.blockMesh
snappyHexMesh -overwrite > log.snappyHexMesh
renumberMesh -overwrite > log.renumberMesh
checkMesh > log.checkMesh

before running, but still get same error with slightly different numbers. I did see on v2406 that they said they made some snappyHexMesh and renumberMesh improvements, so I'm installing v2406 now.

[tmweigand]

Comparing to these files may also be good to do: https://github.com/tmweigand/tcat-1p/tree/main/applications/sim-steady-compressible/system

[acatwithacomputer]

Yes, it was helpful. It's running on v2406 now, and that seems to have resolved the issue. No idea what's up with OpenFOAM 10.... Thanks!

[acatwithacomputer]

The error is back again when I try to move up from 50 spheres and use finer discretizations.

Edit to add: I've also verified that the case that fails in parallel does seem to be running in serial with no issue.

[acatwithacomputer]

Since it runs in serial, but not in parallel, I don't disagree with the notion that this is fundamentally an issue with OpenFOAM's decompositon algorithm. However, a couple observations that I wanted to add is that:

1. I've seen across both versions of OpenFOAM is that when I get this error, I've noticed that the absolute value of the differences as consistently around 1e-8.
2. Sometimes python adds on these little machine epsilons. If there's a square root anywhere, that could be enough to throw this off.
   I'm suspecting that these small truncation errors are enough to trigger whatever is going on with OpenFOAM.

[acatwithacomputer]

Brief update:
First, I found the mergeTolerance the error was referring to, and it's found in the boundary file of the mesh (polyMesh/boundary). Moreover, the mergeTolerance setting between boundaries to be adjusted is automatically generated from decomposePar, so I don't think it's practical to change it after the fact.

I have a possible solution that I'm testing at the moment and will return with an update once that is done.

[acatwithacomputer]

I believe I found a solution.

First, I think the "correct" way to go about it is to set the match tolerance in blockMeshDict and then add the preservePatch constraint in decomposeParDict to keep the matching cyclic patches on the same processor. Unfortunately, decomposePar seems to hang on sycamore, so I wrote a python script to just replace the matchTolerance in all of the processor*/constant/polyMesh/boundary files and it seems to be running consistently now.

edit: of course, as soon as I write that, one doesn't seem to be working