diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev index 416d297f77b6b6afb6323e47dc85ac7c0e2c55bc..6af144a309419dd0cbbae5abfb5391b36dcb18df 100644 --- a/ReleaseNotes-dev +++ b/ReleaseNotes-dev @@ -104,6 +104,8 @@ taking film into account + Parallel aware *** *New* ptscotch decomposition method +*** *Updated* decomposePar maps polyPatches instead of recreating them so + polyPatches holding data can map the data. *** *Updated* particle tracking algorithm + uses non-blocking parallel transfers + does 'minimum-tet' decomposition of face to work with warped faces (snappyHexMesh!) @@ -150,7 +152,12 @@ OpenFOAM. + *New* wall functions: + kappatJayatillekeWallFunction: incompressible RAS thermal wall function - + + directMappedFixedValue: + + takes interpolationScheme so can interpolate instead of always getting + cell value + + takes optional fieldName to sample + + directMapped patch added 'normal' method to calculate sample points + to e.g. sample fields just above wall (e.g. for streaklines) * Utilities There have been some utilities added and updated in this release. *** *New* utilities @@ -172,8 +179,10 @@ *** Updated utilities + =setFields=: optionally use faceSets to set patch values (see e.g. hotRoom tutorial). + =blockMesh=: specify patches via dictionary instead of type only. This - makes rereading the boundary superfluous. see + makes rereading the boundary file superfluous. see e.g. pitzDailyDirectMapped tutorial. + + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument + to execute the commands for multiple times. * Post-processing + =foamToEnsight=: parallel continuous data. new =-nodeValues= option to generate and output nodal field data. @@ -190,6 +199,12 @@ specified time has been reached, e.g. to automagically change fvSchemes and fvSolution during a calculation + =streamLine=: generate streamlines; ouputs both trajectory and field data + + =surfaceInterpolateFields=: constructs face interpolate of registered + volFields for any future functionObjects that need surfaceFields. + + =readFields=: reads field if not yet registered. Makes functionObjects + useable through standalone execFlowFunctionObjects. + + =faceSource=: can now calculate on a sampledSurface (e.g. flow through a + triSurfaceMesh) * New tutorials There is a large number of new tutorials for existing and new solvers in the diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict index 81b087456d2567e87467d07ce36c1e1e9f8c84c6..afa8e1e21dec194410e9235ed2a07362effc5dcd 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposeParDict +++ b/applications/utilities/parallelProcessing/decomposePar/decomposeParDict @@ -19,15 +19,17 @@ FoamFile numberOfSubdomains 8; - //- Keep owner and neighbour on same processor for faces in zones: // preserveFaceZones (heater solid1 solid3); - //- Keep owner and neighbour on same processor for faces in patches: // (makes sense only for cyclic patches) //preservePatches (cyclic_half0 cyclic_half1); +//- Use the volScalarField named here as a weight for each cell in the +// decomposition. For example, use a particle population field to decompose +// for a balanced number of particles in a lagrangian simulation. +// weightField dsmcRhoNMean; method scotch; // method hierarchical; @@ -59,11 +61,8 @@ multiLevelCoeffs } } - // Desired output - - simpleCoeffs { n (2 1 1); diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C index 785cf8e975b57095b65dc3ee526fe492a5b34fc1..d69c2c4dfd0b456571bd28dbde3e676ab1d63d02 100644 --- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C +++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C @@ -24,7 +24,6 @@ License \*---------------------------------------------------------------------------*/ #include "domainDecomposition.H" -#include "Time.H" #include "dictionary.H" #include "labelIOList.H" #include "processorPolyPatch.H" @@ -341,10 +340,10 @@ bool Foam::domainDecomposition::writeDecomposition() const labelList& curProcessorPatchStarts = procProcessorPatchStartIndex_[procI]; - const labelListList& curSubPatchIDs = + const labelListList& curSubPatchIDs = procProcessorPatchSubPatchIDs_[procI]; - const labelListList& curSubStarts = + const labelListList& curSubStarts = procProcessorPatchSubPatchStarts_[procI]; const polyPatchList& meshPatches = boundaryMesh(); diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H index d916d11ab1ec21f8171632896be7552cf982daa0..74ffcbfaf755888de0856eb5ac53ce5d6914d542 100644 --- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H +++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.H @@ -41,6 +41,9 @@ SourceFiles #include "SLList.H" #include "PtrList.H" #include "point.H" +#include "Time.H" +#include "volFields.H" + namespace Foam { @@ -80,7 +83,7 @@ class domainDecomposition // original face. In order to do this properly, all face // indices will be incremented by 1 and the decremented as // necessary to avoid the problem of face number zero having no - // sign. + // sign. List<DynamicList<label> > procFaceAddressing_; //- Labels of cells for each processor diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C index 6544ab769f39765aedf83c306fbbc3fad7a85ef3..2c0067df2fa7511f3747a648ed2fb9254c7acf86 100644 --- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C +++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C @@ -114,7 +114,35 @@ void Foam::domainDecomposition::distributeCells() if (sameProcFaces.empty()) { - cellToProc_ = decomposePtr().decompose(*this, cellCentres()); + if (decompositionDict_.found("weightField")) + { + word weightName = decompositionDict_.lookup("weightField"); + + volScalarField weights + ( + IOobject + ( + weightName, + time().timeName(), + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + *this + ); + + cellToProc_ = decomposePtr().decompose + ( + *this, + cellCentres(), + weights.internalField() + ); + } + else + { + cellToProc_ = decomposePtr().decompose(*this, cellCentres()); + } + } else { @@ -173,12 +201,49 @@ void Foam::domainDecomposition::distributeCells() // Do decomposition on agglomeration // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - cellToProc_ = decomposePtr().decompose - ( - *this, - globalRegion, - regionCentres - ); + if (decompositionDict_.found("weightField")) + { + scalarField regionWeights(globalRegion.nRegions(), 0); + + word weightName = decompositionDict_.lookup("weightField"); + + volScalarField weights + ( + IOobject + ( + weightName, + time().timeName(), + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + *this + ); + + forAll(globalRegion, cellI) + { + label regionI = globalRegion[cellI]; + + regionWeights[regionI] += weights.internalField()[cellI]; + } + + cellToProc_ = decomposePtr().decompose + ( + *this, + globalRegion, + regionCentres, + regionWeights + ); + } + else + { + cellToProc_ = decomposePtr().decompose + ( + *this, + globalRegion, + regionCentres + ); + } } Info<< "\nFinished decomposition in " diff --git a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C index 4bb427665e81a2a3b937a8a2681b8950f166639f..985349019f1a3f1043c94364c282b74266fe38f5 100644 --- a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C +++ b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C @@ -86,9 +86,35 @@ autoPtr<fvMesh> createMesh if (!haveMesh) { // Create dummy mesh. Only used on procs that don't have mesh. + // WIP: how to avoid parallel comms when loading IOdictionaries? + // For now just give error message. + if + ( + regIOobject::fileModificationChecking + == regIOobject::timeStampMaster + || regIOobject::fileModificationChecking + == regIOobject::inotifyMaster + ) + { + FatalErrorIn("createMesh(..)") + << "Cannot use 'fileModificationChecking' mode " + << regIOobject::fileCheckTypesNames + [ + regIOobject::fileModificationChecking + ] + << " since this uses parallel communication." + << exit(FatalError); + } + fvMesh dummyMesh ( - io, + IOobject + ( + regionName, + instDir, + runTime, + IOobject::NO_READ + ), xferCopy(pointField()), xferCopy(faceList()), xferCopy(labelList()), @@ -510,16 +536,14 @@ int main(int argc, char *argv[]) "specify the merge distance relative to the bounding box size " "(default 1E-6)" ); - // Create argList. This will check for non-existing processor dirs. # include "setRootCase.H" - //- Not useful anymore. See above. - //// Create processor directory if non-existing - //if (!Pstream::master() && !isDir(args.path())) - //{ - // Pout<< "Creating case directory " << args.path() << endl; - // mkDir(args.path()); - //} + // Create processor directory if non-existing + if (!Pstream::master() && !isDir(args.path())) + { + Pout<< "Creating case directory " << args.path() << endl; + mkDir(args.path()); + } # include "createTime.H" diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C index 069ef6de9344c71e4c47fc6081bf37add74ec669..329e52ccc0a68ed9515126fba3f7d12ca611de1a 100644 --- a/src/OSspecific/POSIX/fileMonitor.C +++ b/src/OSspecific/POSIX/fileMonitor.C @@ -319,7 +319,7 @@ void Foam::fileMonitor::checkFiles() const if (ready < 0) { - FatalErrorIn("fileMonitor::updateStates()") + FatalErrorIn("fileMonitor::checkFiles()") << "Problem in issuing select." << abort(FatalError); } @@ -335,7 +335,7 @@ void Foam::fileMonitor::checkFiles() const if (nBytes < 0) { - FatalErrorIn("fileMonitor::updateStates(const fileName&)") + FatalErrorIn("fileMonitor::checkFiles()") << "read of " << watcher_->inotifyFd_ << " failed with " << label(nBytes) << abort(FatalError); @@ -374,7 +374,7 @@ void Foam::fileMonitor::checkFiles() const ) { // Correct directory and name - state_[i] = MODIFIED; + localState_[i] = MODIFIED; } } } @@ -403,18 +403,17 @@ void Foam::fileMonitor::checkFiles() const if (newTime == 0) { - state_[watchFd] = DELETED; + localState_[watchFd] = DELETED; } else { if (newTime > (oldTime + regIOobject::fileModificationSkew)) { - watcher_->lastMod_[watchFd] = newTime; - state_[watchFd] = MODIFIED; + localState_[watchFd] = MODIFIED; } else { - state_[watchFd] = UNMODIFIED; + localState_[watchFd] = UNMODIFIED; } } } @@ -422,12 +421,14 @@ void Foam::fileMonitor::checkFiles() const } } + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::fileMonitor::fileMonitor(const bool useInotify) : useInotify_(useInotify), + localState_(20), state_(20), watchFile_(20), freeWatchFds_(2), @@ -476,6 +477,7 @@ Foam::label Foam::fileMonitor::addWatch(const fileName& fName) } else { + localState_(watchFd) = UNMODIFIED; state_(watchFd) = UNMODIFIED; watchFile_(watchFd) = fName; } @@ -517,30 +519,26 @@ void Foam::fileMonitor::updateStates { if (Pstream::master() || !masterOnly) { + // Update the localState_ checkFiles(); } if (syncPar) { - // Pack current state (might be on master only) + // Pack local state (might be on master only) PackedList<2> stats(state_.size(), MODIFIED); if (Pstream::master() || !masterOnly) { forAll(state_, watchFd) { - stats[watchFd] = static_cast<unsigned int>(state_[watchFd]); + stats[watchFd] = static_cast<unsigned int> + ( + localState_[watchFd] + ); } } - // Save local state for warning message below - PackedList<2> thisProcStats; - if (!masterOnly) - { - thisProcStats = stats; - } - - // Scatter or reduce to synchronise state if (masterOnly) { @@ -573,43 +571,49 @@ void Foam::fileMonitor::updateStates } - // Update local state + // Update synchronised state forAll(state_, watchFd) { - if (masterOnly) - { - // No need to check for inconsistent state. Just assign. - unsigned int stat = stats[watchFd]; - state_[watchFd] = fileState(stat); - } - else + // Assign synchronised state + unsigned int stat = stats[watchFd]; + state_[watchFd] = fileState(stat); + + if (!masterOnly) { - // Check for inconsistent state before assigning. - if (thisProcStats[watchFd] != UNMODIFIED) + // Give warning for inconsistent state + if (state_[watchFd] != localState_[watchFd]) { - if (stats[watchFd] == UNMODIFIED) + if (debug) { - WarningIn("fileMonitor::updateStates(const bool) const") - << "Delaying reading " << watchFile_[watchFd] + Pout<< "fileMonitor : Delaying reading " + << watchFile_[watchFd] << " due to inconsistent " "file time-stamps between processors" << endl; } - else - { - unsigned int stat = stats[watchFd]; - state_[watchFd] = fileState(stat); - } + + WarningIn + ( + "fileMonitor::updateStates" + "(const bool, const bool) const" + ) << "Delaying reading " << watchFile_[watchFd] + << " due to inconsistent " + "file time-stamps between processors" << endl; } } } } + else + { + state_ = localState_; + } } void Foam::fileMonitor::setUnmodified(const label watchFd) { state_[watchFd] = UNMODIFIED; + localState_[watchFd] = UNMODIFIED; if (!useInotify_) { diff --git a/src/OSspecific/POSIX/fileMonitor.H b/src/OSspecific/POSIX/fileMonitor.H index d19c3b0e90fb35c4aa6f26949a55fc4ac7898799..bdb0dd2bb2f75a82473ef12ffb87c04f844317d8 100644 --- a/src/OSspecific/POSIX/fileMonitor.H +++ b/src/OSspecific/POSIX/fileMonitor.H @@ -71,7 +71,7 @@ public: { UNMODIFIED = 0, MODIFIED = 1, - DELETED = 2, + DELETED = 2 }; static const NamedEnum<fileState, 3> fileStateNames_; @@ -82,7 +82,10 @@ private: //- Whether to use inotify (requires -DFOAM_USE_INOTIFY, see above) const bool useInotify_; - //- State for all watchFds + //- State for all watchFds based on local files + mutable DynamicList<fileState> localState_; + + //- State for all watchFds - synchronised mutable DynamicList<fileState> state_; //- Filename for all watchFds @@ -97,7 +100,7 @@ private: // Private Member Functions - //- Update state_ from any events. + //- Update localState_ from any events. void checkFiles() const; //- Disallow default bitwise copy construct diff --git a/src/OpenFOAM/db/regIOobject/regIOobject.C b/src/OpenFOAM/db/regIOobject/regIOobject.C index 2367d79c60dbc9c5bdcc3066b647920b385d72a1..f43cb4f44ef22261fccfc03393e8797f6aa6c068 100644 --- a/src/OpenFOAM/db/regIOobject/regIOobject.C +++ b/src/OpenFOAM/db/regIOobject/regIOobject.C @@ -60,7 +60,6 @@ Foam::regIOobject::fileCheckTypes Foam::regIOobject::fileModificationChecking debug::optimisationSwitches().lookup ( "fileModificationChecking" - //Foam::regIOobject::timeStamp ) ) ); diff --git a/src/OpenFOAM/global/argList/argList.C b/src/OpenFOAM/global/argList/argList.C index d1ee0a63f75f4516402baa1d50de965f92ff3c63..ee623a30bbeed8732506b3796364c07b18251060 100644 --- a/src/OpenFOAM/global/argList/argList.C +++ b/src/OpenFOAM/global/argList/argList.C @@ -31,6 +31,7 @@ License #include "IOobject.H" #include "JobInfo.H" #include "labelList.H" +#include "regIOobject.H" #include <cctype> @@ -767,6 +768,16 @@ Foam::argList::argList sigQuit_.set(bannerEnabled); sigSegv_.set(bannerEnabled); + if (bannerEnabled) + { + Info<< "Monitoring run-time modified files using " + << regIOobject::fileCheckTypesNames + [ + regIOobject::fileModificationChecking + ] + << endl; + } + if (Pstream::master() && bannerEnabled) { Info<< endl; diff --git a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C index 9e59239fd1f7000c6d7c9e01e144f9bcc75fdfd4..d4c2802d588b7e0135d5b2b7dbaa600f52d0c9f2 100644 --- a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C +++ b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C @@ -58,13 +58,35 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) {} +Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +( + const List<int>& initxadj, + const List<int>& initadjncy, + const scalarField& initcWeights, + + List<int>& finalDecomp +) const +{ + FatalErrorIn + ( + "label ptscotchDecomp::decompose" + "(" + "const List<int>&, " + "const List<int>&, " + "const scalarField&, " + "List<int>&" + ")" + ) << notImplementedMessage << exit(FatalError); + + return -1; +} Foam::label Foam::ptscotchDecomp::decompose ( - List<int>& adjncy, - List<int>& xadj, + const List<int>& adjncy, + const List<int>& xadj, const scalarField& cWeights, List<int>& finalDecomp -) +) const { FatalErrorIn ( diff --git a/src/lagrangian/basic/Particle/Particle.C b/src/lagrangian/basic/Particle/Particle.C index 1fb312298ec6425fa0b8506ef6f37d8a2df8f455..b0bc7a8d0ceba789c6d43621aefc2f3469f430a9 100644 --- a/src/lagrangian/basic/Particle/Particle.C +++ b/src/lagrangian/basic/Particle/Particle.C @@ -62,37 +62,29 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer if (!ppp.parallel()) { - if (ppp.forwardT().size() == 1) - { - const tensor& T = ppp.forwardT()[0]; - transformPosition(T); - static_cast<ParticleType&>(*this).transformProperties(T); - } - else - { - const tensor& T = ppp.forwardT()[faceI_]; - transformPosition(T); - static_cast<ParticleType&>(*this).transformProperties(T); - } + const tensor& T = + ( + ppp.forwardT().size() == 1 + ? ppp.forwardT()[0] + : ppp.forwardT()[faceI_] + ); + + transformPosition(T); + static_cast<ParticleType&>(*this).transformProperties(T); } else if (ppp.separated()) { - if (ppp.separation().size() == 1) - { - position_ -= ppp.separation()[0]; - static_cast<ParticleType&>(*this).transformProperties - ( - -ppp.separation()[0] - ); - } - else - { - position_ -= ppp.separation()[faceI_]; - static_cast<ParticleType&>(*this).transformProperties - ( - -ppp.separation()[faceI_] - ); - } + const vector& s = + ( + (ppp.separation().size() == 1) + ? ppp.separation()[0] + : ppp.separation()[faceI_] + ); + position_ -= s; + static_cast<ParticleType&>(*this).transformProperties + ( + -s + ); } tetFaceI_ = faceI_ + ppp.start(); @@ -815,21 +807,34 @@ void Foam::Particle<ParticleType>::hitCyclicPatch // See note in correctAfterParallelTransfer for tetPtI_ addressing. tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_; + const cyclicPolyPatch& receiveCpp = cpp.neighbPatch(); + // Now the particle is on the receiving side - if (!cpp.parallel()) + if (!receiveCpp.parallel()) { - const tensor& T = cpp.reverseT()[0]; + const tensor& T = + ( + receiveCpp.forwardT().size() == 1 + ? receiveCpp.forwardT()[0] + : receiveCpp.forwardT()[receiveCpp.whichFace(faceI_)] + ); transformPosition(T); static_cast<ParticleType&>(*this).transformProperties(T); } - else if (cpp.separated()) + else if (receiveCpp.separated()) { - position_ += cpp.separation()[0]; + const vector& s = + ( + (receiveCpp.separation().size() == 1) + ? receiveCpp.separation()[0] + : receiveCpp.separation()[receiveCpp.whichFace(faceI_)] + ); + position_ -= s; static_cast<ParticleType&>(*this).transformProperties ( - cpp.separation()[0] + -s ); } } diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C index ff262c5645d28388a71529f20f6c0a2dfaceed91..f08ab6536369c1286e9433b841372a1d30e75ea8 100644 --- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C @@ -108,6 +108,7 @@ License #include "addToRunTimeSelectionTable.H" #include "Time.H" #include "OFstream.H" +#include "globalIndex.H" extern "C" { @@ -157,16 +158,204 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) } +//- Does prevention of 0 cell domains and calls ptscotch. +Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +( + const List<int>& initadjncy, + const List<int>& initxadj, + const scalarField& initcWeights, + + List<int>& finalDecomp +) const +{ + globalIndex globalCells(initxadj.size()-1); + + bool hasZeroDomain = false; + for (label procI = 0; procI < Pstream::nProcs(); procI++) + { + if (globalCells.localSize(procI) == 0) + { + hasZeroDomain = true; + break; + } + } + + if (!hasZeroDomain) + { + return decompose + ( + initadjncy, + initxadj, + initcWeights, + finalDecomp + ); + } + + + if (debug) + { + Info<< "ptscotchDecomp : have graphs with locally 0 cells." + << " trickling down." << endl; + } + + // Make sure every domain has at least one cell + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // (scotch does not like zero sized domains) + // Trickle cells from processors that have them up to those that + // don't. + + + // Number of cells to send to the next processor + // (is same as number of cells next processor has to receive) + List<int> nSendCells(Pstream::nProcs(), 0); + + for (label procI = nSendCells.size()-1; procI >=1; procI--) + { + label nLocalCells = globalCells.localSize(procI); + if (nLocalCells-nSendCells[procI] < 1) + { + nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1; + } + } + + // First receive (so increasing the sizes of all arrays) + + Field<int> xadj(initxadj); + Field<int> adjncy(initadjncy); + scalarField cWeights(initcWeights); + + if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) + { + // Receive cells from previous processor + IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1); + + Field<int> prevXadj(fromPrevProc); + Field<int> prevAdjncy(fromPrevProc); + scalarField prevCellWeights(fromPrevProc); + + if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) + { + FatalErrorIn("ptscotchDecomp::decompose(..)") + << "Expected from processor " << Pstream::myProcNo()-1 + << " connectivity for " << nSendCells[Pstream::myProcNo()-1] + << " nCells but only received " << prevXadj.size() + << abort(FatalError); + } + + // Insert adjncy + prepend(prevAdjncy, adjncy); + // Adapt offsets and prepend xadj + xadj += prevAdjncy.size(); + prepend(prevXadj, xadj); + // Weights + prepend(prevCellWeights, cWeights); + } + + + // Send to my next processor + + if (nSendCells[Pstream::myProcNo()] > 0) + { + // Send cells to next processor + OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1); + + int nCells = nSendCells[Pstream::myProcNo()]; + int startCell = xadj.size()-1 - nCells; + int startFace = xadj[startCell]; + int nFaces = adjncy.size()-startFace; + + // Send for all cell data: last nCells elements + // Send for all face data: last nFaces elements + toNextProc + << Field<int>::subField(xadj, nCells, startCell)-startFace + << Field<int>::subField(adjncy, nFaces, startFace) + << + ( + cWeights.size() + ? scalarField::subField(cWeights, nCells, startCell) + : scalarField(0) + ); + + // Remove data that has been sent + if (cWeights.size()) + { + cWeights.setSize(cWeights.size()-nCells); + } + adjncy.setSize(adjncy.size()-nFaces); + xadj.setSize(xadj.size() - nCells); + } + + + // Do decomposition as normal. Sets finalDecomp. + label result = decompose(adjncy, xadj, cWeights, finalDecomp); + + + if (debug) + { + Info<< "ptscotchDecomp : have graphs with locally 0 cells." + << " trickling up." << endl; + } + + + // If we sent cells across make sure we undo it + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // Receive back from next processor if I sent something + if (nSendCells[Pstream::myProcNo()] > 0) + { + IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1); + + List<int> nextFinalDecomp(fromNextProc); + + if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) + { + FatalErrorIn("parMetisDecomp::decompose(..)") + << "Expected from processor " << Pstream::myProcNo()+1 + << " decomposition for " << nSendCells[Pstream::myProcNo()] + << " nCells but only received " << nextFinalDecomp.size() + << abort(FatalError); + } + + append(nextFinalDecomp, finalDecomp); + } + + // Send back to previous processor. + if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) + { + OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1); + + int nToPrevious = nSendCells[Pstream::myProcNo()-1]; + + toPrevProc << + SubList<int> + ( + finalDecomp, + nToPrevious, + finalDecomp.size()-nToPrevious + ); + + // Remove locally what has been sent + finalDecomp.setSize(finalDecomp.size()-nToPrevious); + } + return result; +} + + // Call scotch with options from dictionary. Foam::label Foam::ptscotchDecomp::decompose ( - List<int>& adjncy, - List<int>& xadj, + const List<int>& adjncy, + const List<int>& xadj, const scalarField& cWeights, List<int>& finalDecomp -) +) const { + if (debug) + { + Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl; + } + // // Dump graph // if (decompositionDict_.found("ptscotchCoeffs")) // { @@ -262,8 +451,7 @@ Foam::label Foam::ptscotchDecomp::decompose { WarningIn ( - "ptscotchDecomp::decompose" - "(const pointField&, const scalarField&)" + "ptscotchDecomp::decompose(..)" ) << "Illegal minimum weight " << minWeights << endl; } @@ -272,8 +460,7 @@ Foam::label Foam::ptscotchDecomp::decompose { FatalErrorIn ( - "ptscotchDecomp::decompose" - "(const pointField&, const scalarField&)" + "ptscotchDecomp::decompose(..)" ) << "Number of cell weights " << cWeights.size() << " does not equal number of cells " << xadj.size()-1 << exit(FatalError); @@ -289,8 +476,25 @@ Foam::label Foam::ptscotchDecomp::decompose + if (debug) + { + Pout<< "SCOTCH_dgraphInit" << endl; + } SCOTCH_Dgraph grafdat; check(SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD), "SCOTCH_dgraphInit"); + + + if (debug) + { + Pout<< "SCOTCH_dgraphBuild with:" << nl + << "xadj.size() : " << xadj.size()-1 << nl + << "xadj : " << long(xadj.begin()) << nl + << "velotab : " << long(velotab.begin()) << nl + << "adjncy.size() : " << adjncy.size() << nl + << "adjncy : " << long(adjncy.begin()) << nl + << endl; + } + check ( SCOTCH_dgraphBuild @@ -302,19 +506,25 @@ Foam::label Foam::ptscotchDecomp::decompose const_cast<SCOTCH_Num*>(xadj.begin()), // vertloctab, start index per cell into // adjncy - &xadj[1], // vendloctab, end index ,, + const_cast<SCOTCH_Num*>(&xadj[1]),// vendloctab, end index ,, - velotab.begin(), // veloloctab, vertex weights + const_cast<SCOTCH_Num*>(velotab.begin()),// veloloctab, vtx weights NULL, // vlblloctab adjncy.size(), // edgelocnbr, number of arcs adjncy.size(), // edgelocsiz - adjncy.begin(), // edgeloctab + const_cast<SCOTCH_Num*>(adjncy.begin()), // edgeloctab NULL, // edgegsttab NULL // edlotab, edge weights ), "SCOTCH_dgraphBuild" ); + + + if (debug) + { + Pout<< "SCOTCH_dgraphCheck" << endl; + } check(SCOTCH_dgraphCheck(&grafdat), "SCOTCH_dgraphCheck"); @@ -322,6 +532,10 @@ Foam::label Foam::ptscotchDecomp::decompose // ~~~~~~~~~~~~ // (fully connected network topology since using switch) + if (debug) + { + Pout<< "SCOTCH_archInit" << endl; + } SCOTCH_Arch archdat; check(SCOTCH_archInit(&archdat), "SCOTCH_archInit"); @@ -349,6 +563,10 @@ Foam::label Foam::ptscotchDecomp::decompose } else { + if (debug) + { + Pout<< "SCOTCH_archCmplt" << endl; + } check ( SCOTCH_archCmplt(&archdat, nProcessors_), @@ -373,6 +591,10 @@ Foam::label Foam::ptscotchDecomp::decompose ); # endif + if (debug) + { + Pout<< "SCOTCH_dgraphMap" << endl; + } finalDecomp.setSize(xadj.size()-1); finalDecomp = 0; check @@ -406,6 +628,10 @@ Foam::label Foam::ptscotchDecomp::decompose // "SCOTCH_graphPart" //); + if (debug) + { + Pout<< "SCOTCH_dgraphExit" << endl; + } // Release storage for graph SCOTCH_dgraphExit(&grafdat); // Release storage for strategy @@ -465,7 +691,13 @@ Foam::labelList Foam::ptscotchDecomp::decompose // Decompose using default weights List<int> finalDecomp; - decompose(cellCells.m(), cellCells.offsets(), pointWeights, finalDecomp); + decomposeZeroDomains + ( + cellCells.m(), + cellCells.offsets(), + pointWeights, + finalDecomp + ); // Copy back to labelList labelList decomp(finalDecomp.size()); @@ -510,7 +742,13 @@ Foam::labelList Foam::ptscotchDecomp::decompose // Decompose using weights List<int> finalDecomp; - decompose(cellCells.m(), cellCells.offsets(), pointWeights, finalDecomp); + decomposeZeroDomains + ( + cellCells.m(), + cellCells.offsets(), + pointWeights, + finalDecomp + ); // Rework back into decomposition for original mesh labelList fineDistribution(agglom.size()); @@ -557,7 +795,13 @@ Foam::labelList Foam::ptscotchDecomp::decompose // Decompose using weights List<int> finalDecomp; - decompose(cellCells.m(), cellCells.offsets(), cWeights, finalDecomp); + decomposeZeroDomains + ( + cellCells.m(), + cellCells.offsets(), + cWeights, + finalDecomp + ); // Copy back to labelList labelList decomp(finalDecomp.size()); diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H index 7471f5ddd79540d112e902eeea99602790751df7..6a6937ac79a0d50a2ce66a5a8d90ba29af22390b 100644 --- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H @@ -50,16 +50,33 @@ class ptscotchDecomp { // Private Member Functions + //- Insert list in front of list. + template<class Type> + static void prepend(const UList<Type>&, List<Type>&); + //- Insert list at end of list. + template<class Type> + static void append(const UList<Type>&, List<Type>&); + //- Check and print error message static void check(const int, const char*); + //- Decompose with optionally zero sized domains + label decomposeZeroDomains + ( + const List<int>& initadjncy, + const List<int>& initxadj, + const scalarField& initcWeights, + List<int>& finalDecomp + ) const; + + //- Decompose label decompose ( - List<int>& adjncy, - List<int>& xadj, + const List<int>& adjncy, + const List<int>& xadj, const scalarField& cWeights, List<int>& finalDecomp - ); + ) const; //- Disallow default bitwise copy construct and assignment void operator=(const ptscotchDecomp&); @@ -135,6 +152,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "ptscotchDecompTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecompTemplates.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecompTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..e133dc9ceb60ada703affd55735f6962b4fa65fe --- /dev/null +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecompTemplates.C @@ -0,0 +1,76 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\/ 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 "ptscotchDecomp.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +// Insert at front of list +template<class Type> +void Foam::ptscotchDecomp::prepend +( + const UList<Type>& extraLst, + List<Type>& lst +) +{ + label nExtra = extraLst.size(); + + // Make space for initial elements + lst.setSize(lst.size() + nExtra); + for (label i = lst.size()-1; i >= nExtra; i--) + { + lst[i] = lst[i-nExtra]; + } + + // Insert at front + forAll(extraLst, i) + { + lst[i] = extraLst[i]; + } +} + + +// Insert at back of list +template<class Type> +void Foam::ptscotchDecomp::append +( + const UList<Type>& extraLst, + List<Type>& lst +) +{ + label sz = lst.size(); + + // Make space for initial elements + lst.setSize(sz + extraLst.size()); + + // Insert at back + forAll(extraLst, i) + { + lst[sz++] = extraLst[i]; + } +} + + +// ************************************************************************* //