View Issue Details

IDProjectCategoryView StatusLast Update
0003537OpenFOAMBugpublic2020-08-25 17:02
Reportermillszg Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Platformx86-64OSUbuntuOS Version18.04
Summary0003537: faceAgglomerate hangs for certain values of nFacesInCoarsestLevel and featureAngle in parallel
DescriptionfaceAgglomerate will occasionally hang for certain values of nFacesInCoarsestLevel and featureAngle. This only occurs when running in parallel, and the same values that cause it to hang for a certain number of processors will work for a different number of processors. Increasing nFacesInCoarsestLevel will sometimes alleviate it, but oftentimes it needs to be increased so much that very little agglomeration is actually occurring. In runs where the program hangs, it will always print "Agglomerating patch : patch_name" at least one time, but will not always print for each patch being agglomerated. Keeping everything the same, but rearranging the order of patches in the viewfactorsDict file will usually lead to the programming hanging after outputting "Agglomerating patch..." for a different patch. From what I have tested, this occurs in OpenFOAM 6, 7, 8 and dev. I have had this issues with multiple different meshes and on different computers. From what I have found, this is actually caused by a bug in the pairPatchAgglomeration class (see additional information for explanation).
Steps To ReproduceA case which will reproduce the error is attached. It is a single region of a multiregion case, so some of the boundaries are mappedWall. This error was occurring with the entire multi-region mesh was included in the case folder, so the use of mapped walls without the patch being mapped to is not the cause of the issue. The mesh has already been generated and calling faceAgglomerate will cause it to hang. The culprit seems to be innerGas_to_holderBase. When nFacesInCoarsestLevel for the patch is set to 7000, faceAgglomerate hangs. When it is increased to 10000, it will complete successfully. When the innerGas_to_holderBase is the first in the list of patch dictionaries in viewFactorsDict, faceAgglomerate seems to hang on innerGas_to_retort. When it's the fourth dictionary (between innerGas_to_tempDisks and innerGas_to_retort), faceAgglomerate seems to hang on innerGas_to_retortInlet.
Additional InformationFrom debugging it appears the problem is in the pairPatchAgglomeration class. In the agglomerate() function, at least one process hangs while the rest continue, causing the program to be unable to continue. This seems to happen in the do-while loop where the agglomOK variable for processes that hang remains false, while it is set to true for the other processes and they complete the function. The issue seems to be resolved by adding in a call to reduce the agglomOK variable with an orOp operation the reduce operation for nCoarseFaces on line 412 of pairPatchAgglomeration.C (in version 8). I have attached a copy of pairPatchAgglomeration.C that I am using, which is no longer causing the error for me.
TagsNo tags attached.

Activities

millszg

2020-08-21 15:34

reporter  

pairPatchAgglomeration.C (15,853 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "pairPatchAgglomeration.H"
#include "meshTools.H"
#include "unitConversion.H"

#include <unistd.h>
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::pairPatchAgglomeration::compactLevels(const label nCreatedLevels)
{
    nFaces_.setSize(nCreatedLevels);
    restrictAddressing_.setSize(nCreatedLevels);
    patchLevels_.setSize(nCreatedLevels);
}


bool Foam::pairPatchAgglomeration::continueAgglomerating
(
    const label nCoarseFaces
)
{
    // Check the need for further agglomeration on all processors
    label localnCoarseFaces = nCoarseFaces;
    bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
    return contAgg;
}


void Foam::pairPatchAgglomeration::setBasedEdgeWeights()
{
    const bPatch& coarsePatch = patchLevels_[0];
    forAll(coarsePatch.edges(), i)
    {
        if (coarsePatch.isInternalEdge(i))
        {
            scalar edgeLength =
                coarsePatch.edges()[i].mag(coarsePatch.localPoints());

            const labelList& eFaces = coarsePatch.edgeFaces()[i];

            if (eFaces.size() == 2)
            {
                scalar cosI =
                    coarsePatch.faceNormals()[eFaces[0]]
                  & coarsePatch.faceNormals()[eFaces[1]];

                const edge edgeCommon = edge(eFaces[0], eFaces[1]);

                if (facePairWeight_.found(edgeCommon))
                {
                    facePairWeight_[edgeCommon] += edgeLength;
                }
                else
                {
                    facePairWeight_.insert(edgeCommon, edgeLength);
                }

                if (cosI < Foam::cos(degToRad(featureAngle_)))
                {
                    facePairWeight_[edgeCommon] = -1.0;
                }
            }
            else
            {
                forAll(eFaces, j)
                {
                    for (label k = j+1; k<eFaces.size(); k++)
                    {
                        facePairWeight_.insert
                        (
                            edge(eFaces[j], eFaces[k]),
                            -1.0
                        );
                    }
                }
            }
        }
    }
}


void Foam::pairPatchAgglomeration::setEdgeWeights
(
    const label fineLevelIndex
)
{

    const bPatch& coarsePatch = patchLevels_[fineLevelIndex];

    const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
    const label nCoarseI =  max(fineToCoarse) + 1;
    labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));

    HashSet<edge, Hash<edge>> fineFeaturedFaces(coarsePatch.nEdges()/10);

    // Map fine faces with featured edge into coarse faces
    forAllConstIter(EdgeMap<scalar>, facePairWeight_, iter)
    {
        if (iter() == -1.0)
        {
            const edge e = iter.key();
            const edge edgeFeatured
            (
                fineToCoarse[e[0]],
                fineToCoarse[e[1]]
            );
            fineFeaturedFaces.insert(edgeFeatured);
        }
    }

    // Clean old weitghs
    facePairWeight_.clear();
    facePairWeight_.resize(coarsePatch.nEdges());

    forAll(coarsePatch.edges(), i)
    {
        if (coarsePatch.isInternalEdge(i))
        {
            scalar edgeLength =
                coarsePatch.edges()[i].mag(coarsePatch.localPoints());

            const labelList& eFaces = coarsePatch.edgeFaces()[i];

            if (eFaces.size() == 2)
            {
                const edge edgeCommon = edge(eFaces[0], eFaces[1]);

                if (facePairWeight_.found(edgeCommon))
                {
                    facePairWeight_[edgeCommon] += edgeLength;
                }
                else
                {
                    facePairWeight_.insert(edgeCommon, edgeLength);
                }

                // If the fine 'pair' faces was featured edge so it is
                // the coarse 'pair'
                if (fineFeaturedFaces.found(edgeCommon))
                {
                    facePairWeight_[edgeCommon] = -1.0;
                }
            }
            else
            {
                // Set edge as barrier by setting weight to -1
                forAll(eFaces, j)
                {
                    for (label k = j+1; k<eFaces.size(); k++)
                    {
                        facePairWeight_.insert
                        (
                            edge(eFaces[j], eFaces[k]),
                            -1.0
                        );
                    }
                }
            }
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::pairPatchAgglomeration::pairPatchAgglomeration
(
    const polyPatch& patch,
    const dictionary& controlDict,
    const bool additionalWeights
)
:
    mergeLevels_
    (
        controlDict.lookupOrDefault<label>("mergeLevels", 2)
    ),
    maxLevels_(50),
    nFacesInCoarsestLevel_
    (
        controlDict.lookup<label>("nFacesInCoarsestLevel")
    ),
    featureAngle_
    (
        controlDict.lookupOrDefault<scalar>("featureAngle", 0)
    ),
    nFaces_(maxLevels_),
    restrictAddressing_(maxLevels_),
    restrictTopBottomAddressing_(identity(patch.size())),
    patchLevels_(maxLevels_),
    facePairWeight_(patch.size())
{
    // Set base fine patch
    patchLevels_.set
    (
        0,
        new bPatch
        (
            patch.localFaces(),
            patch.localPoints()
        )
    );

    // Set number of faces for the base patch
    nFaces_[0] = patch.size();

    // Set edge weights for level 0
    setBasedEdgeWeights();
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::pairPatchAgglomeration::~pairPatchAgglomeration()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

const Foam::bPatch& Foam::pairPatchAgglomeration::patchLevel
(
    const label i
) const
{
    return patchLevels_[i];
}


void Foam::pairPatchAgglomeration::mapBaseToTopAgglom
(
    const label fineLevelIndex
)
{
    const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
    forAll(restrictTopBottomAddressing_, i)
    {
        restrictTopBottomAddressing_[i] =
            fineToCoarse[restrictTopBottomAddressing_[i]];
    }
}


bool Foam::pairPatchAgglomeration::agglomeratePatch
(
    const bPatch& patch,
    const labelList& fineToCoarse,
    const label fineLevelIndex
)
{
    if (min(fineToCoarse) == -1)
    {
        FatalErrorInFunction
            << "min(fineToCoarse) == -1" << exit(FatalError);
    }

    if (fineToCoarse.size() == 0)
    {
        return true;
    }

    if (fineToCoarse.size() != patch.size())
    {
        FatalErrorInFunction
            << "restrict map does not correspond to fine level. " << endl
            << " Sizes: restrictMap: " << fineToCoarse.size()
            << " nEqns: " << patch.size()
            << abort(FatalError);
    }

    const label nCoarseI =  max(fineToCoarse)+1;
    List<face> patchFaces(nCoarseI);


    // Patch faces per agglomeration
    labelListList coarseToFine(invertOneToMany(nCoarseI, fineToCoarse));

    for (label coarseI = 0; coarseI < nCoarseI; coarseI++)
    {
        const labelList& fineFaces = coarseToFine[coarseI];

        // Construct single face
        indirectPrimitivePatch upp
        (
            IndirectList<face>(patch, fineFaces),
            patch.points()
        );

        if (upp.edgeLoops().size() != 1)
        {
            if (fineFaces.size() == 2)
            {
                const edge e(fineFaces[0], fineFaces[1]);
                facePairWeight_[e] = -1.0;
            }
            else if (fineFaces.size() == 3)
            {
                const edge e(fineFaces[0], fineFaces[1]);
                const edge e1(fineFaces[0], fineFaces[2]);
                const edge e2(fineFaces[2], fineFaces[1]);
                facePairWeight_[e] = -1.0;
                facePairWeight_[e1] = -1.0;
                facePairWeight_[e2] = -1.0;
            }
            return false;
        }

        patchFaces[coarseI] = face
        (
            renumber
            (
                upp.meshPoints(),
                upp.edgeLoops()[0]
            )
        );
    }

    patchLevels_.set
    (
        fineLevelIndex,
        new bPatch
        (
            SubList<face>(patchFaces, nCoarseI, 0),
            patch.points()
        )
    );


    return true;

}


void Foam::pairPatchAgglomeration::agglomerate()
{
    label nPairLevels = 0;
    label nCreatedLevels = 1; // 0 level is the base patch
    label nCoarseFaces = 0;
    label nCoarseFacesOld = 0;

    while (nCreatedLevels < maxLevels_)
    {
        const bPatch& patch = patchLevels_[nCreatedLevels - 1];
        tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
        bool agglomOK = false;
        do
        {
            label nCoarseFacesPrev = nCoarseFaces;

            finalAgglomPtr = agglomerateOneLevel
            (
                nCoarseFaces,
                patch
            );
            
            if (nCoarseFaces > 0 && nCoarseFaces != nCoarseFacesPrev)
            {
                if
                (
                    (
                        agglomOK = agglomeratePatch
                        (
                            patch,
                            finalAgglomPtr,
                            nCreatedLevels
                        )
                    )
                )
                {
                    restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
                    mapBaseToTopAgglom(nCreatedLevels);
                    setEdgeWeights(nCreatedLevels);

                    if (nPairLevels % mergeLevels_)
                    {
                        combineLevels(nCreatedLevels);
                    }
                    else
                    {
                        nCreatedLevels++;
                    }

                    nPairLevels++;
                }
            }
            else
            {
                agglomOK = true;
            }

            reduce(nCoarseFaces, sumOp<label>());

            // This line appears to eliminate the hanging issue
            reduce(agglomOK, orOp<bool>());
         
        } while (!agglomOK);

        nFaces_[nCreatedLevels] = nCoarseFaces;

        if
        (
            !continueAgglomerating(nCoarseFaces)
          || (nCoarseFacesOld ==  nCoarseFaces)
        )
        {
            break;
        }
        nCoarseFacesOld = nCoarseFaces;
    }
    compactLevels(nCreatedLevels);
}


Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
(
    label& nCoarseFaces,
    const bPatch& patch
)
{
    const label nFineFaces = patch.size();

    tmp<labelField> tcoarseCellMap(new labelField(nFineFaces, -1));
    labelField& coarseCellMap = tcoarseCellMap.ref();

    const labelListList& faceFaces = patch.faceFaces();

    nCoarseFaces = 0;

    forAll(faceFaces, facei)
    {
        const labelList& fFaces = faceFaces[facei];

        if (coarseCellMap[facei] < 0)
        {
            label matchFaceNo = -1;
            label matchFaceNeibNo = -1;
            scalar maxFaceWeight = -great;

            // Check faces to find ungrouped neighbour with largest face weight
            forAll(fFaces, i)
            {
                label faceNeig = fFaces[i];
                const edge edgeCommon = edge(facei, faceNeig);
                if
                (
                    facePairWeight_[edgeCommon] > maxFaceWeight
                 && coarseCellMap[faceNeig] < 0
                 && facePairWeight_[edgeCommon] != -1.0
                )
                {
                    // Match found. Pick up all the necessary data
                    matchFaceNo = facei;
                    matchFaceNeibNo = faceNeig;
                    maxFaceWeight = facePairWeight_[edgeCommon];
                }
            }

            if (matchFaceNo >= 0)
            {
                // Make a new group
                coarseCellMap[matchFaceNo] = nCoarseFaces;
                coarseCellMap[matchFaceNeibNo] = nCoarseFaces;
                nCoarseFaces++;
            }
            else
            {
                // No match. Find the best neighbouring cluster and
                // put the cell there
                label clusterMatchFaceNo = -1;
                scalar clusterMaxFaceCoeff = -great;

                forAll(fFaces, i)
                {
                    label faceNeig = fFaces[i];
                    const edge edgeCommon = edge(facei, faceNeig);
                    if
                    (
                        facePairWeight_[edgeCommon] > clusterMaxFaceCoeff
                        && facePairWeight_[edgeCommon] != -1.0
                        && coarseCellMap[faceNeig] > 0
                    )
                    {
                        clusterMatchFaceNo = faceNeig;
                        clusterMaxFaceCoeff = facePairWeight_[edgeCommon];
                    }
                }

                if (clusterMatchFaceNo >= 0)
                {
                    // Add the cell to the best cluster
                    coarseCellMap[facei] = coarseCellMap[clusterMatchFaceNo];
                }
                else
                {
                    // if not create single-cell "clusters" for each
                    coarseCellMap[facei] = nCoarseFaces;
                    nCoarseFaces++;
                }
            }
        }
    }

    // Check that all faces are part of clusters,

    for (label facei=0; facei<nFineFaces; facei++)
    {
        if (coarseCellMap[facei] < 0)
        {
            FatalErrorInFunction
              << " face " << facei
            << " is not part of a cluster"
            << exit(FatalError);
        }
    }

    return tcoarseCellMap;
}


void Foam::pairPatchAgglomeration::combineLevels(const label curLevel)
{
    label prevLevel = curLevel - 1;

    // Set the previous level nCells to the current
    nFaces_[prevLevel] = nFaces_[curLevel];

    // Map the restrictAddressing from the coarser level into the previous
    // finer level

    const labelList& curResAddr = restrictAddressing_[curLevel];
    labelList& prevResAddr = restrictAddressing_[prevLevel];

    forAll(prevResAddr, i)
    {
        prevResAddr[i] = curResAddr[prevResAddr[i]];
    }

    // Delete the restrictAddressing for the coarser level
    restrictAddressing_.set(curLevel, nullptr);

    patchLevels_.set(prevLevel, patchLevels_.set(curLevel, nullptr));
}


// ************************************************************************* //
pairPatchAgglomeration.C (15,853 bytes)   

henry

2020-08-25 17:02

manager   ~0011478

Thanks for the report and correction, I agree that with the current rather convoluted logic an additional reduction is needed

Resolved in OpenFOAM-8 by commit dd46ea14cf84c7d19a07f2f8d633d7c74025d5dd
Resolved in OpenFOAM-dev by commit 2ad3fa242aa2c0c4de9cab5fbb275df093e3ac29

Issue History

Date Modified Username Field Change
2020-08-21 15:34 millszg New Issue
2020-08-21 15:34 millszg File Added: pairPatchAgglomeration.C
2020-08-25 17:02 henry Assigned To => henry
2020-08-25 17:02 henry Status new => resolved
2020-08-25 17:02 henry Resolution open => fixed
2020-08-25 17:02 henry Note Added: 0011478