View Issue Details

IDProjectCategoryView StatusLast Update
0001203OpenFOAM[All Projects] Bugpublic2015-02-06 00:42
Reporteruser503Assigned Tohenry 
PrioritynormalSeveritycrashReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSUbuntuOS Version12.04
Product Version 
Fixed in Version 
Summary0001203: Several bugs in mesh distribution prevent run-time parallel load balancing
DescriptionI have constructed a runtime parallel load balancing library for OpenFOAM derived from dynamicRefineFvMesh and in the process identified several bugs in the distribution of fields during the mesh rebalancing operations. These issues are present in 2.1.x to 2.3.x, and without fixing them there are a variety of fatal crashes when using the fvMeshDistribute::distribute function.

There are 6 bugs in this report, described in detail below. The line numbers are based on 2.1.x and may vary slightly from the original source.
Steps To ReproduceFor all bugs except those related to DimensionedFields, attempt to use load balancing on an interDyMFoam case. For the DimensionedField bug, use load balancing on a solver with a dynamic mesh and a chemistry model.
Additional Information*******************************************************************************
Issue 1
*******************************************************************************
Add a guard in src/dynamicMesh/polyTopoChange/refinementHistory.C in
refinementHistory::distribute at line 927 to catch
when fully un-refined cells are transferred
     
    if( newVisibleCells[i] >= 0 )
    {
        visibleCells_[constructMap[i]] = newVisibleCells[i] + offset;
    }
    
    
*******************************************************************************
Issue 2
*******************************************************************************
Add `if (debug)` guard around print statements in
`refinementHistory::countProc` in the same file as Issue 1 to prevent excessive
printing to stdout during mesh balancing operations


*******************************************************************************
Issue 3
*******************************************************************************
Add the following to src/finiteVolume/fvMesh/wallDist/nearWallDist.C
at `nearWallDist::correct()` inside the
`mesh_.changing()` condition but before the patch sizes are set. If the
total number of patches on a processor increases, there is a crash unless
we increase the size of nearWallDist first. There may be a more elegant way
of doing this too.
      
    // If the number of patches on this processor increased, we need to
    // increase the number of patches in nearWallDist to match
    if( mesh_.boundary().size() != size() )
    {
        //Resize nearWallDist to match mesh patches
        resize(mesh_.boundary().size());

        //Set any newly created (unset) patches
        forAll(*this, patchI)
        {
            if( !set(patchI) )
            {
                set
                (
                    patchI,
                    fvPatchField<scalar>::New
                    (
                        calculatedFvPatchScalarField::typeName,
                        mesh_.boundary()[patchI],
                        mesh_.V()
                    )
                );
            }
        }
    }

*******************************************************************************
Issue 4
*******************************************************************************
Mapping of patches had an error for mixed patches. In
src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C in function
`fvMeshAdder::MapVolField` around line 263 there is a section where imported patch values
are mapped to existing patches. This correctly maps the patch value, but
omits mapping of the refValue, refGradient, and other stored fields in
mixed patches. The refValue is then left
set to a section of allocated but unset memory if all the original faces
of that patch get moved to another processor. A similar edit may be required in
the same file near line 578 for the MapSurfaceField function, although errors from this have not shown up in my tests so far.

Remove the manual loop near the end of the function(s) and use the rmap
function with a reversed map by replacing
 
    forAll(newFld, i)
    {
        label oldFaceI = newToAdded[i];
    
        if (oldFaceI >= 0 && oldFaceI < addedFld.size())
        {
            newFld[i] = addedFld[oldFaceI];
        }
    }
    
with

    labelList addedToNew(addedFld.size(),-1);
    forAll(newFld, i)
    {
        label oldFaceI = newToAdded[i];

        if (oldFaceI >= 0 && oldFaceI < addedFld.size())
        {
            addedToNew[oldFaceI] = i;
        }
    }
    
    newFld.rmap(addedFld, addedToNew);

*******************************************************************************
Issue 5
*******************************************************************************
DimensionedFields are not distributed. See the attached file for the necessary fixes (this is the most involved of the bugs). This becomes more of a problem in 2.2.x and newer since more of the chemistry fields were switched to DimensionedFields.

*******************************************************************************
Issue 6 (non-crash)
*******************************************************************************
In the current implementation of dynamicRefineFvMesh, a refined cell can be
coarsened if the minimum refinementField value in any of its child cells is
less than the threshhold in the dictionary. This is leads to oscillatory
refinement when the refinement field is sharply defined. Consider the following cell, showing the value of the field on which refinement is based in each cell
  
     +------+------+
     | | |
     | 0 | 9.9 |
     | | |
     +------+------+
     | | |
     | 9.9 | 9.9 |
     | | |
     +------+------+
     
where refinement is triggered at a value of 10 and unrefinement at a value of, say, 0.5. This cell, rather than being left alone, will be coarsened. It should use the maximum value rather than the minimum value. To change this, in `dynamicRefineFvMesh.H` on line 112 change `minCellField`
to `maxCellField`. In `dynamicRefineFvMesh.C`, make the following 4 changes. Line numbers here may not match your line numbers exactly.
Changes 1-3 will all be in the existing `minCellField` function and change 4 is the only location in the file that referenced the `minCellField` function.
 
 1. Line 646: change `minCellField` to `maxCellField`
 2. Line 648: change `GREAT` to `-GREAT`
 3. Line 656: change `min` to `max`
 4. Line 1236: change `minCellField` to `maxCellField`
TagsNo tags attached.

Activities

user503

2014-03-02 18:05

 

Dimensioned Field Fixes (11,420 bytes)
Steps to Fix Distribution of DimensionedFields

===============================================================================
Step 1 - Add a constructor for DimensionedField which is consistent with the one
         used in fvMeshSubset for geometric fields, with inputs of
         IOobject, mesh, and dictionary
         
         It's a good idea to re-compile src after this step.
===============================================================================

*******************************************************************************
src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H

line 151ish, add:

    //- Construct from dictionary
    DimensionedField
    (
        const IOobject&,
        const Mesh& mesh,
        const dictionary& fieldDict,
        const word& fieldDictEntry="value"
    );
    
*******************************************************************************
src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C

line 82ish, add:

    template<class Type, class GeoMesh>
    Foam::DimensionedField<Type, GeoMesh>::DimensionedField
    (
        const IOobject& io,
        const Mesh& mesh,
        const dictionary& fieldDict,
        const word& fieldDictEntry
    )
    :
        regIOobject(io),
        Field<Type>(0),
        mesh_(mesh),
        dimensions_(dimless)
    {
        readField(fieldDict, fieldDictEntry);
    }

===============================================================================
Step 2 - Add an interpolate template to fvMeshSubset that can handle
         a DimensionedField<type,volMesh> input
===============================================================================

*******************************************************************************
src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H

line 345ish, add:

    //- Map dimensioned fields
    template<class Type>
    static tmp<DimensionedField<Type, volMesh> >
    interpolate
    (
        const DimensionedField<Type, volMesh>&,
        const fvMesh& sMesh,
        const labelList& cellMap
    );
    
    template<class Type>
    tmp<DimensionedField<Type, volMesh> >
    interpolate
    (
        const DimensionedField<Type, volMesh>&
    ) const;
    
*******************************************************************************
src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetInterpolate.C

line 38ish, add:

    template<class Type>
    tmp<DimensionedField<Type, volMesh> > fvMeshSubset::interpolate
    (
        const DimensionedField<Type, volMesh>& df,
        const fvMesh& sMesh,
        const labelList& cellMap
    )
    {
        // Create and map the internal-field values
        Field<Type> internalField(df, cellMap);

        // Create the complete field from the pieces
        tmp<DimensionedField<Type, volMesh> > tresF
        (
            new DimensionedField<Type, volMesh>
            (
                IOobject
                (
                    "subset"+df.name(),
                    sMesh.time().timeName(),
                    sMesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                sMesh,
                df.dimensions(),
                internalField
            )
        );

        return tresF;
    }


    template<class Type>
    tmp<DimensionedField<Type, volMesh> > fvMeshSubset::interpolate
    (
        const DimensionedField<Type, volMesh>& df
    ) const
    {
        return interpolate
        (
            df,
            subMesh(),
            cellMap()
        );
    }
    
===============================================================================
Step 3 - Update the field mapping in fvMeshAdder to include the templates for
         DimensionedFields and add the calls to map those fields to the "add"
         function
===============================================================================

*******************************************************************************
src/dynamicMesh/fvMeshAdder/fvMeshAdder.H

line 148ish, add:

    //- Update single dimensioned Field.
    template<class Type>
    static void MapDimField
    (
        const mapAddedPolyMesh& meshMap,

        DimensionedField<Type, volMesh>& fld,
        const DimensionedField<Type, volMesh>& fldToAdd
    );

line 190ish, add:

    //- Map all dimensionedFields of Type
    template<class Type>
    static void MapDimFields
    (
        const mapAddedPolyMesh&,
        const fvMesh& mesh,
        const fvMesh& meshToAdd
    );

*******************************************************************************
src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C

line 353ish, add:

    template<class Type>
    void Foam::fvMeshAdder::MapDimField
    (
        const mapAddedPolyMesh& meshMap,

        DimensionedField<Type, volMesh>& fld,
        const DimensionedField<Type, volMesh>& fldToAdd
    )
    {
        const fvMesh& mesh = fld.mesh();

        // Internal field
        // ~~~~~~~~~~~~~~

        {
            // Store old internal field
            Field<Type> oldInternalField(fld);

            // Modify internal field
            Field<Type>& intFld = fld;

            intFld.setSize(mesh.nCells());

            map(oldInternalField, meshMap.oldCellMap(), intFld);
            map(fldToAdd, meshMap.addedCellMap(), intFld);
        }
    }


    template<class Type>
    void Foam::fvMeshAdder::MapDimFields
    (
        const mapAddedPolyMesh& meshMap,
        const fvMesh& mesh,
        const fvMesh& meshToAdd
    )
    {
        HashTable<const DimensionedField<Type, volMesh>*> fields
        (
            mesh.objectRegistry::lookupClass
            <DimensionedField<Type, volMesh> >
            ()
        );

        HashTable<const DimensionedField<Type, volMesh>*> fieldsToAdd
        (
            meshToAdd.objectRegistry::lookupClass
            <DimensionedField<Type, volMesh> >
            ()
        );


        for
        (
            typename HashTable<const DimensionedField<Type, volMesh>*>::
                iterator fieldIter = fields.begin();
            fieldIter != fields.end();
            ++fieldIter
        )
        {
            DimensionedField<Type, volMesh>& fld =
                const_cast<DimensionedField<Type, volMesh>&>
                (
                    *fieldIter()
                );

            if (fieldsToAdd.found(fld.name()))
            {
                const DimensionedField<Type, volMesh>& fldToAdd =
                    *fieldsToAdd[fld.name()];

                MapDimField<Type>(meshMap, fld, fldToAdd);
            }
            else
            {
                WarningIn("fvMeshAdder::MapDimFields(..)")
                    << "Not mapping field " << fld.name()
                    << " since not present on mesh to add"
                    << endl;
            }
        }
    }

*******************************************************************************
src/dynamicMesh/fvMeshAdder/fvMeshAdder.C

line 104ish, add:

    fvMeshAdder::MapDimFields<scalar>(mapPtr, mesh0, mesh1);
    fvMeshAdder::MapDimFields<vector>(mapPtr, mesh0, mesh1);
    fvMeshAdder::MapDimFields<sphericalTensor>(mapPtr, mesh0, mesh1);
    fvMeshAdder::MapDimFields<symmTensor>(mapPtr, mesh0, mesh1);
    fvMeshAdder::MapDimFields<tensor>(mapPtr, mesh0, mesh1);
    


===============================================================================
Step 4 - Add in the actual distribution if the different types of
         DimensionedFields in fvMeshDistribute.
         
         Some typedefs of DimensionedFields like there are of geometric fields
         would be nice. Something like:
         
           typedef DimensionedField<scalar,volmesh> dimensionedScalarField;
           typedef DimensionedField<vector,volmesh> dimensionedVectorField;
           
           and soforth
           
===============================================================================

*******************************************************************************
src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
        
line 1735ish, add:
    
    const wordList dimScalarFields
    (
        mesh_.names(DimensionedField<scalar, volMesh>::typeName)
    );
    checkEqualWordList("dimScalarFields", dimScalarFields);
    
    const wordList dimVectorFields
    (
        mesh_.names(DimensionedField<vector, volMesh>::typeName)
    );
    checkEqualWordList("dimVectorFields", dimVectorFields);
    
    const wordList dimSphericalTensorFields
    (
        mesh_.names(DimensionedField<sphericalTensor, volMesh>::typeName)
    );
    checkEqualWordList("dimSphericalTensorFields", dimSphericalTensorFields);
    
    const wordList dimSymmTensorFields
    (
        mesh_.names(DimensionedField<symmTensor, volMesh>::typeName)
    );
    checkEqualWordList("dimSymmTensorFields", dimSymmTensorFields);
    
    const wordList dimTensorFields
    (
        mesh_.names(DimensionedField<tensor, volMesh>::typeName)
    );
    checkEqualWordList("dimTensorFields", dimTensorFields);
    
    
    
line 1997ish (immediately after the sendFields<surfaceTensorField>, add:
    
    sendFields<DimensionedField<scalar,volMesh> >
    (
        recvProc,
        dimScalarFields,
        subsetter,
        str
    );
    sendFields<DimensionedField<vector,volMesh> >
    (
        recvProc,
        dimVectorFields,
        subsetter,
        str
    );
    sendFields<DimensionedField<sphericalTensor,volMesh> >
    (
        recvProc,
        dimSphericalTensorFields,
        subsetter,
        str
    );
    sendFields<DimensionedField<symmTensor,volMesh> >
    (
        recvProc,
        dimSymmTensorFields,
        subsetter,
        str
    );
    sendFields<DimensionedField<tensor,volMesh> >
    (
        recvProc,
        dimTensorFields,
        subsetter,
        str
    );
    
line 2180ish (after the PtrList<surfaceTensorField> line), add:

    PtrList<DimensionedField<scalar,volMesh> > isf;
    PtrList<DimensionedField<vector,volMesh> > ivf;
    PtrList<DimensionedField<sphericalTensor,volMesh> > istf;
    PtrList<DimensionedField<symmTensor,volMesh> > isytf;
    PtrList<DimensionedField<tensor,volMesh> > itf;
    
line 2290ish (after the receiveFields<surfaceTensorField> call), add:

    receiveFields<DimensionedField<scalar,volMesh> >
    (
        sendProc,
        dimScalarFields,
        domainMesh,
        isf,
        fieldDicts.subDict(DimensionedField<scalar,volMesh>::typeName)
    );
    receiveFields<DimensionedField<vector,volMesh> >
    (
        sendProc,
        dimVectorFields,
        domainMesh,
        ivf,
        fieldDicts.subDict(DimensionedField<vector,volMesh>::typeName)
    );
    receiveFields<DimensionedField<sphericalTensor,volMesh> >
    (
        sendProc,
        dimSphericalTensorFields,
        domainMesh,
        istf,
        fieldDicts.subDict(DimensionedField<sphericalTensor,volMesh>::typeName)
    );
    receiveFields<DimensionedField<symmTensor,volMesh> >
    (
        sendProc,
        dimSymmTensorFields,
        domainMesh,
        isytf,
        fieldDicts.subDict(DimensionedField<symmTensor,volMesh>::typeName)
    );
    receiveFields<DimensionedField<tensor,volMesh> >
    (
        sendProc,
        dimTensorFields,
        domainMesh,
        itf,
        fieldDicts.subDict(DimensionedField<tensor,volMesh>::typeName)
    );
    

Dimensioned Field Fixes (11,420 bytes)

user503

2014-03-03 02:34

  ~0002921

Small update, for 2.3.x the MapDimField function should use the newer rmap functions instead of 'map', to be consistent with MapVolField

intFld.rmap(oldInternalField, meshMap.oldCellMap());
intFld.rmap(fldToAdd, meshMap.addedCellMap());

user4

2014-05-23 16:11

  ~0003082

Dear tgvosk,

thanks for the bugfixes

1,2 fixed in 23x
3 fixed in the next version
4 fixed in 23x
5 fixed in the next version (and also pre,postprocessing of dimensioned fields)
6 currently testing

user4

2014-05-28 15:26

  ~0003091

6. resolved in 50304f561dfd61a94abe99d8a02639acb167431b in 23x

user4

2014-05-30 16:47

  ~0003099

3. resolving in 14edbacd706552be1f7f65d987d8a71f57026001.

Thanks for reporting.

Issue History

Date Modified Username Field Change
2014-03-02 18:05 user503 New Issue
2014-03-02 18:05 user503 File Added: Dimensioned Field Fixes
2014-03-03 02:34 user503 Note Added: 0002921
2014-05-23 16:11 user4 Note Added: 0003082
2014-05-28 15:26 user4 Note Added: 0003091
2014-05-30 16:47 user4 Note Added: 0003099
2015-02-05 09:15 henry Status new => resolved
2015-02-05 09:15 henry Resolution open => fixed
2015-02-05 09:15 henry Assigned To => henry