[Kellis]

Good afternoon,

I am trying to write a piece of code that will track the center of a patch as it moves within the dynamic mesh. The center location is needed for calculations in the solver during runtime. Using the following code, I get the expected results when running the case serially:

        // Get current patch
        scalar patchNum = mesh.boundaryMesh().findPatchID(fricPatchList[patchi]);
        const polyPatch& currPatch = mesh.boundaryMesh()[patchNum];

        // Initialize center vector and area
        vector patchCenter = Zero;
        scalar patchArea = VSMALL;

        // Create pointers to face surface area and centers (used by both loops)
        const surfaceScalarField& magSf = mesh.magSf();
        const surfaceVectorField& Cf = mesh.Cf();

        // Calculate area-weighted patch center
        forAll(currPatch, facei)
        {
            // Add area-weighted face center and face area to totals
            patchCenter += Cf.boundaryField()[patchNum][facei]*magSf.boundaryField()[patchNum][facei];
            patchArea += magSf.boundaryField()[patchNum][facei];
        }

        // Calculate patch centroid
        patchCenter = patchCenter / patchArea;

        Info << "Center for patch " << fricPatchList[patchi] << " found to be at " << patchCenter << endl;

On my test case, a circle patch with origin at (0 0 0.0015), this accurately calculates the center to within rounding error. However, issues begin after decomposing the case and trying to run it in parallel. The calculated center is immediately wrong, and by adding a counter to the loop I can see that it's only looping over the faces of the patch owned by processor0, with the centroid ending up in the middle of this group of faces.

So, is there a way to access the full patch data while running in parallel? Or, will I have to come up with a workaround, such as keeping the entire patch on one processor? This would be less than ideal, as the patch is directly next to a finely discretized region in my mesh that should be split up into several parts.

Any help would be appreciated. Thanks.

Kellis

[Kellis]

Bump. Has anyone been able to overcome this issue before?

[olesen]

Originally Posted by Kellis

Bump. Has anyone been able to overcome this issue before?

Lots of places in OpenFOAM have that type of thing.
For example,

vector patchCentre = Zero;

scalar patchArea = 0;


// Do loops here ...


reduce(patchCentre, sumOp<vector>());
reduce(patchArea, sumOp<area>());


patchCentre /= (patchArea + VSMALL);

Not sure why this doesn't work for you.

[Kellis]

Olesen,


Thank you for chiming in. I'm not quite clear on the purpose of the "reduce" function. Does it somehow sum the patchCentre calculations for all processors? In that case, would it be placed after the forAll loops in my sample code? I might be missing something obvious, but I would appreciate a bit more help with the implementation.


Thanks,
Kellis

[Kellis]

Nevermind, I plugged the code you mentioned into my solver and it seems to be working perfectly now - thank you for pointing me in the right direction. If anyone needs to replicate this in the future, the working patch center calculation looks like this:

        // Get current patch
        scalar patchNum = mesh.boundaryMesh().findPatchID(fricPatchList[patchi]);
        const polyPatch& currPatch = mesh.boundaryMesh()[patchNum];

        // Initialize center vector and area
        vector patchCenter = Zero;
        scalar patchArea = VSMALL;

        // Create pointers to face surface area and centers (used by both loops)
        const surfaceScalarField& magSf = mesh.magSf();
        const surfaceVectorField& Cf = mesh.Cf();

        // Calculate area-weighted patch center
        forAll(currPatch, facei)
        {
            // Add area-weighted face center and face area to totals
            patchCenter += Cf.boundaryField()[patchNum][facei]*magSf.boundaryField()[patchNum][facei];
            patchArea += magSf.boundaryField()[patchNum][facei];
        }

        reduce(patchCenter, sumOp<vector>());
        reduce(patchArea, sumOp<scalar>());

        // Calculate patch centroid
        patchCenter = patchCenter / patchArea;

         Info << "Center for patch " << fricPatchList[patchi] << " found to be at " << patchCenter << endl;

Thanks again olesen!