View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0001203 | OpenFOAM | Bug | public | 2014-03-02 18:05 | 2015-02-06 00:42 |
Reporter | Assigned To | henry | |||
Priority | normal | Severity | crash | Reproducibility | always |
Status | resolved | Resolution | fixed | ||
Platform | Linux | OS | Ubuntu | OS Version | 12.04 |
Summary | 0001203: Several bugs in mesh distribution prevent run-time parallel load balancing | ||||
Description | I 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 Reproduce | For 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` | ||||
Tags | No tags attached. | ||||
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) ); |
|
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()); |
|
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 |
|
6. resolved in 50304f561dfd61a94abe99d8a02639acb167431b in 23x |
|
3. resolving in 14edbacd706552be1f7f65d987d8a71f57026001. Thanks for reporting. |
Date Modified | Username | Field | Change |
---|---|---|---|
2014-03-02 18:05 |
|
New Issue | |
2014-03-02 18:05 |
|
File Added: Dimensioned Field Fixes | |
2014-03-03 02:34 |
|
Note Added: 0002921 | |
2014-05-23 16:11 |
|
Note Added: 0003082 | |
2014-05-28 15:26 |
|
Note Added: 0003091 | |
2014-05-30 16:47 |
|
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 |