FreeFOAM The Cross-Platform CFD Toolkit
redistributeMeshPar.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Application
25  redistributeMeshPar
26 
27 Description
28  Redistributes existing decomposed mesh and fields according to the current
29  settings in the decomposeParDict file.
30 
31  Must be run on maximum number of source and destination processors.
32  Balances mesh and writes new mesh to new time directory.
33 
34  Can also work like decomposePar:
35 
36  # Create empty processors
37  mkdir processor0
38  ..
39  mkdir processorN
40 
41  # Copy undecomposed polyMesh
42  cp -r constant processor0
43 
44  # Distribute
45  mpirun -np ddd redistributeMeshPar -parallel
46 
47 Note
48  You might want to unset FOAM_SIGFPE and FOAM_SETNAN since
49  patchfields that hold additional data might not be initialised
50  (since mapped from 0 faces)
51 
52 Usage
53  - redistributeMeshPar [OPTIONS]
54 
55  @param -mergeTol <number>\n
56  Relative merge distance.
57 
58  @param -case <dir>\n
59  Case directory.
60 
61  @param -parallel \n
62  Run in parallel.
63 
64  @param -help \n
65  Display help message.
66 
67  @param -doc \n
68  Display Doxygen API documentation page for this application.
69 
70  @param -srcDoc \n
71  Display Doxygen source documentation page for this application.
72 
73 \*---------------------------------------------------------------------------*/
74 
75 #include <finiteVolume/fvMesh.H>
78 #include <finiteVolume/fvCFD.H>
81 #include <OpenFOAM/IOobjectList.H>
82 
83 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
84 
85 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
86 // usually meshes get written with limited precision (6 digits)
87 static const scalar defaultMergeTol = 1E-6;
88 
89 
90 // Read mesh if available. Otherwise create empty mesh with same non-proc
91 // patches as proc0 mesh. Requires all processors to have all patches
92 // (and in same order).
93 autoPtr<fvMesh> createMesh
94 (
95  const Time& runTime,
96  const word& regionName,
97  const fileName& instDir,
98  const bool haveMesh
99 )
100 {
101  Pout<< "Create mesh for time = "
102  << runTime.timeName() << nl << endl;
103 
104  IOobject io
105  (
106  regionName,
107  instDir,
108  runTime,
109  IOobject::MUST_READ
110  );
111 
112  if (!haveMesh)
113  {
114  // Create dummy mesh. Only used on procs that don't have mesh.
115  fvMesh dummyMesh
116  (
117  io,
118  xferCopy(pointField()),
119  xferCopy(faceList()),
120  xferCopy(labelList()),
121  xferCopy(labelList()),
122  false
123  );
124  Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
125  << endl;
126  dummyMesh.write();
127  }
128 
129  Pout<< "Reading mesh from " << io.objectPath() << endl;
130  autoPtr<fvMesh> meshPtr(new fvMesh(io));
131  fvMesh& mesh = meshPtr();
132 
133 
134  // Determine patches.
135  if (Pstream::master())
136  {
137  // Send patches
138  for
139  (
140  int slave=Pstream::firstSlave();
141  slave<=Pstream::lastSlave();
142  slave++
143  )
144  {
145  OPstream toSlave(Pstream::blocking, slave);
146  toSlave << mesh.boundaryMesh();
147  }
148  }
149  else
150  {
151  // Receive patches
152  IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
153  PtrList<entry> patchEntries(fromMaster);
154 
155  if (haveMesh)
156  {
157  // Check master names against mine
158 
159  const polyBoundaryMesh& patches = mesh.boundaryMesh();
160 
161  forAll(patchEntries, patchI)
162  {
163  const entry& e = patchEntries[patchI];
164  const word type(e.dict().lookup("type"));
165  const word& name = e.keyword();
166 
167  if (type == processorPolyPatch::typeName)
168  {
169  break;
170  }
171 
172  if (patchI >= patches.size())
173  {
175  (
176  "createMesh(const Time&, const fileName&, const bool)"
177  ) << "Non-processor patches not synchronised."
178  << endl
179  << "Processor " << Pstream::myProcNo()
180  << " has only " << patches.size()
181  << " patches, master has "
182  << patchI
183  << exit(FatalError);
184  }
185 
186  if
187  (
188  type != patches[patchI].type()
189  || name != patches[patchI].name()
190  )
191  {
193  (
194  "createMesh(const Time&, const fileName&, const bool)"
195  ) << "Non-processor patches not synchronised."
196  << endl
197  << "Master patch " << patchI
198  << " name:" << type
199  << " type:" << type << endl
200  << "Processor " << Pstream::myProcNo()
201  << " patch " << patchI
202  << " has name:" << patches[patchI].name()
203  << " type:" << patches[patchI].type()
204  << exit(FatalError);
205  }
206  }
207  }
208  else
209  {
210  // Add patch
211  List<polyPatch*> patches(patchEntries.size());
212  label nPatches = 0;
213 
214  forAll(patchEntries, patchI)
215  {
216  const entry& e = patchEntries[patchI];
217  const word type(e.dict().lookup("type"));
218  const word& name = e.keyword();
219 
220  if (type == processorPolyPatch::typeName)
221  {
222  break;
223  }
224 
225  Pout<< "Adding patch:" << nPatches
226  << " name:" << name
227  << " type:" << type << endl;
228 
229  dictionary patchDict(e.dict());
230  patchDict.remove("nFaces");
231  patchDict.add("nFaces", 0);
232  patchDict.remove("startFace");
233  patchDict.add("startFace", 0);
234 
235  patches[patchI] = polyPatch::New
236  (
237  name,
238  patchDict,
239  nPatches++,
240  mesh.boundaryMesh()
241  ).ptr();
242  }
243  patches.setSize(nPatches);
244  mesh.addFvPatches(patches, false); // no parallel comms
245 
247  //meshPtr().write();
248  }
249  }
250 
251  if (!haveMesh)
252  {
253  // We created a dummy mesh file above. Delete it.
254  Pout<< "Removing dummy mesh " << io.objectPath()
255  << endl;
256  rmDir(io.objectPath());
257  }
258 
259  // Force recreation of globalMeshData.
260  mesh.clearOut();
261  mesh.globalData();
262 
263  return meshPtr;
264 }
265 
266 
267 // Get merging distance when matching face centres
268 scalar getMergeDistance
269 (
270  const argList& args,
271  const Time& runTime,
272  const boundBox& bb
273 )
274 {
275  scalar mergeTol = defaultMergeTol;
276  args.optionReadIfPresent("mergeTol", mergeTol);
277 
278  scalar writeTol =
279  Foam::pow(scalar(10.0), -scalar(IOstream::defaultPrecision()));
280 
281  Info<< "Merge tolerance : " << mergeTol << nl
282  << "Write tolerance : " << writeTol << endl;
283 
284  if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
285  {
286  FatalErrorIn("getMergeDistance")
287  << "Your current settings specify ASCII writing with "
288  << IOstream::defaultPrecision() << " digits precision." << endl
289  << "Your merging tolerance (" << mergeTol << ") is finer than this."
290  << endl
291  << "Please change your writeFormat to binary"
292  << " or increase the writePrecision" << endl
293  << "or adjust the merge tolerance (-mergeTol)."
294  << exit(FatalError);
295  }
296 
297  scalar mergeDist = mergeTol * bb.mag();
298 
299  Info<< "Overall meshes bounding box : " << bb << nl
300  << "Relative tolerance : " << mergeTol << nl
301  << "Absolute matching distance : " << mergeDist << nl
302  << endl;
303 
304  return mergeDist;
305 }
306 
307 
308 void printMeshData(Ostream& os, const polyMesh& mesh)
309 {
310  os << "Number of points: " << mesh.points().size() << nl
311  << " faces: " << mesh.faces().size() << nl
312  << " internal faces: " << mesh.faceNeighbour().size() << nl
313  << " cells: " << mesh.cells().size() << nl
314  << " boundary patches: " << mesh.boundaryMesh().size() << nl
315  << " point zones: " << mesh.pointZones().size() << nl
316  << " face zones: " << mesh.faceZones().size() << nl
317  << " cell zones: " << mesh.cellZones().size() << nl;
318 }
319 
320 
321 // Debugging: write volScalarField with decomposition for post processing.
322 void writeDecomposition
323 (
324  const word& name,
325  const fvMesh& mesh,
326  const labelList& decomp
327 )
328 {
329  Info<< "Writing wanted cell distribution to volScalarField " << name
330  << " for postprocessing purposes." << nl << endl;
331 
332  volScalarField procCells
333  (
334  IOobject
335  (
336  name,
337  mesh.time().timeName(),
338  mesh,
339  IOobject::NO_READ,
340  IOobject::AUTO_WRITE,
341  false // do not register
342  ),
343  mesh,
344  dimensionedScalar(name, dimless, -1),
345  zeroGradientFvPatchScalarField::typeName
346  );
347 
348  forAll(procCells, cI)
349  {
350  procCells[cI] = decomp[cI];
351  }
352  procCells.write();
353 }
354 
355 
356 // Read vol or surface fields
357 //template<class T, class Mesh>
358 template<class GeoField>
359 void readFields
360 (
361  const boolList& haveMesh,
362  const fvMesh& mesh,
363  const autoPtr<fvMeshSubset>& subsetterPtr,
364  IOobjectList& allObjects,
366 )
367 {
368  //typedef GeometricField<T, fvPatchField, Mesh> fldType;
369 
370  // Get my objects of type
371  IOobjectList objects(allObjects.lookupClass(GeoField::typeName));
372  // Check that we all have all objects
373  wordList objectNames = objects.toc();
374  // Get master names
375  wordList masterNames(objectNames);
376  Pstream::scatter(masterNames);
377 
378  if (haveMesh[Pstream::myProcNo()] && objectNames != masterNames)
379  {
380  FatalErrorIn("readFields()")
381  << "differing fields of type " << GeoField::typeName
382  << " on processors." << endl
383  << "Master has:" << masterNames << endl
384  << Pstream::myProcNo() << " has:" << objectNames
385  << abort(FatalError);
386  }
387 
388  fields.setSize(masterNames.size());
389 
390  // Have master send all fields to processors that don't have a mesh
391  if (Pstream::master())
392  {
393  forAll(masterNames, i)
394  {
395  const word& name = masterNames[i];
396  IOobject& io = *objects[name];
397  io.writeOpt() = IOobject::AUTO_WRITE;
398 
399  // Load field
400  fields.set(i, new GeoField(io, mesh));
401 
402  // Create zero sized field and send
403  if (subsetterPtr.valid())
404  {
405  tmp<GeoField> tsubfld = subsetterPtr().interpolate(fields[i]);
406 
407  // Send to all processors that don't have a mesh
408  for (label procI = 1; procI < Pstream::nProcs(); procI++)
409  {
410  if (!haveMesh[procI])
411  {
412  OPstream toProc(Pstream::blocking, procI);
413  toProc<< tsubfld();
414  }
415  }
416  }
417  }
418  }
419  else if (!haveMesh[Pstream::myProcNo()])
420  {
421  // Don't have mesh (nor fields). Receive empty field from master.
422 
423  forAll(masterNames, i)
424  {
425  const word& name = masterNames[i];
426 
427  // Receive field
428  IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
429 
430  fields.set
431  (
432  i,
433  new GeoField
434  (
435  IOobject
436  (
437  name,
438  mesh.time().timeName(),
439  mesh,
440  IOobject::NO_READ,
441  IOobject::AUTO_WRITE
442  ),
443  mesh,
444  fromMaster
445  )
446  );
447 
449  //fields[i].write();
450  }
451  }
452  else
453  {
454  // Have mesh so just try to load
455  forAll(masterNames, i)
456  {
457  const word& name = masterNames[i];
458  IOobject& io = *objects[name];
459  io.writeOpt() = IOobject::AUTO_WRITE;
460 
461  // Load field
462  fields.set(i, new GeoField(io, mesh));
463  }
464  }
465 }
466 
467 
468 // Debugging: compare two fields.
469 void compareFields
470 (
471  const scalar tolDim,
472  const volVectorField& a,
473  const volVectorField& b
474 )
475 {
476  forAll(a, cellI)
477  {
478  if (mag(b[cellI] - a[cellI]) > tolDim)
479  {
481  (
482  "compareFields"
483  "(const scalar, const volVectorField&, const volVectorField&)"
484  ) << "Did not map volVectorField correctly:" << nl
485  << "cell:" << cellI
486  << " transfer b:" << b[cellI]
487  << " real cc:" << a[cellI]
488  << abort(FatalError);
489  }
490  }
491  forAll(a.boundaryField(), patchI)
492  {
493  // We have real mesh cellcentre and
494  // mapped original cell centre.
495 
496  const fvPatchVectorField& aBoundary =
497  a.boundaryField()[patchI];
498 
499  const fvPatchVectorField& bBoundary =
500  b.boundaryField()[patchI];
501 
502  if (!bBoundary.coupled())
503  {
504  forAll(aBoundary, i)
505  {
506  if (mag(aBoundary[i] - bBoundary[i]) > tolDim)
507  {
509  (
510  "compareFields"
511  "(const scalar, const volVectorField&"
512  ", const volVectorField&)"
513  ) << "Did not map volVectorField correctly:"
514  << endl
515  << "patch:" << patchI << " patchFace:" << i
516  << " cc:" << endl
517  << " real :" << aBoundary[i] << endl
518  << " mapped :" << bBoundary[i] << endl
519  << abort(FatalError);
520  }
521  }
522  }
523  }
524 }
525 
526 
527 // Main program:
528 
529 int main(int argc, char *argv[])
530 {
531 # include <OpenFOAM/addRegionOption.H>
532  argList::validOptions.insert("mergeTol", "relative merge distance");
533 
534 # include <OpenFOAM/setRootCase.H>
535 
536  // Create processor directory if non-existing
537  if (!Pstream::master() && !isDir(args.path()))
538  {
539  Pout<< "Creating case directory " << args.path() << endl;
540  mkDir(args.path());
541  }
542 
543 # include <OpenFOAM/createTime.H>
544 
545  word regionName = polyMesh::defaultRegion;
546  fileName meshSubDir;
547 
548  if (args.optionReadIfPresent("region", regionName))
549  {
550  meshSubDir = regionName/polyMesh::meshSubDir;
551  }
552  else
553  {
554  meshSubDir = polyMesh::meshSubDir;
555  }
556  Info<< "Using mesh subdirectory " << meshSubDir << nl << endl;
557 
558 
559  // Get time instance directory. Since not all processors have meshes
560  // just use the master one everywhere.
561 
562  fileName masterInstDir;
563  if (Pstream::master())
564  {
565  masterInstDir = runTime.findInstance(meshSubDir, "points");
566  }
567  Pstream::scatter(masterInstDir);
568 
569  // Check who has a mesh
570  const fileName meshPath = runTime.path()/masterInstDir/meshSubDir;
571 
572  Info<< "Found points in " << meshPath << nl << endl;
573 
574 
575  boolList haveMesh(Pstream::nProcs(), false);
576  haveMesh[Pstream::myProcNo()] = isDir(meshPath);
577  Pstream::gatherList(haveMesh);
578  Pstream::scatterList(haveMesh);
579  Info<< "Per processor mesh availability : " << haveMesh << endl;
580  const bool allHaveMesh = (findIndex(haveMesh, false) == -1);
581 
582  // Create mesh
583  autoPtr<fvMesh> meshPtr = createMesh
584  (
585  runTime,
586  regionName,
587  masterInstDir,
588  haveMesh[Pstream::myProcNo()]
589  );
590  fvMesh& mesh = meshPtr();
591 
592  Pout<< "Read mesh:" << endl;
593  printMeshData(Pout, mesh);
594  Pout<< endl;
595 
596  IOdictionary decompositionDict
597  (
598  IOobject
599  (
600  "decomposeParDict",
601  runTime.system(),
602  mesh,
603  IOobject::MUST_READ,
604  IOobject::NO_WRITE
605  )
606  );
607 
608  labelList finalDecomp;
609 
610  // Create decompositionMethod and new decomposition
611  {
613  (
614  decompositionMethod::New
615  (
616  decompositionDict,
617  mesh
618  )
619  );
620 
621  if (!decomposer().parallelAware())
622  {
623  WarningIn(args.executable())
624  << "You have selected decomposition method "
625  << decomposer().typeName
626  << " which does" << endl
627  << "not synchronise the decomposition across"
628  << " processor patches." << endl
629  << " You might want to select a decomposition method which"
630  << " is aware of this. Continuing."
631  << endl;
632  }
633 
634  finalDecomp = decomposer().decompose(mesh.cellCentres());
635  }
636 
637  // Dump decomposition to volScalarField
638  writeDecomposition("decomposition", mesh, finalDecomp);
639 
640 
641  // Create 0 sized mesh to do all the generation of zero sized
642  // fields on processors that have zero sized meshes. Note that this is
643  // only nessecary on master but since polyMesh construction with
644  // Pstream::parRun does parallel comms we have to do it on all
645  // processors
646  autoPtr<fvMeshSubset> subsetterPtr;
647 
648  if (!allHaveMesh)
649  {
650  // Find last non-processor patch.
651  const polyBoundaryMesh& patches = mesh.boundaryMesh();
652 
653  label nonProcI = -1;
654 
655  forAll(patches, patchI)
656  {
657  if (isA<processorPolyPatch>(patches[patchI]))
658  {
659  break;
660  }
661  nonProcI++;
662  }
663 
664  if (nonProcI == -1)
665  {
666  FatalErrorIn(args.executable())
667  << "Cannot find non-processor patch on processor "
668  << Pstream::myProcNo() << endl
669  << " Current patches:" << patches.names() << abort(FatalError);
670  }
671 
672  // Subset 0 cells, no parallel comms. This is used to create zero-sized
673  // fields.
674  subsetterPtr.reset(new fvMeshSubset(mesh));
675  subsetterPtr().setLargeCellSubset(labelHashSet(0), nonProcI, false);
676  }
677 
678 
679  // Get original objects (before incrementing time!)
680  IOobjectList objects(mesh, runTime.timeName());
681  // We don't want to map the decomposition (mapping already tested when
682  // mapping the cell centre field)
683  IOobjectList::iterator iter = objects.find("decomposition");
684  if (iter != objects.end())
685  {
686  objects.erase(iter);
687  }
688 
689  PtrList<volScalarField> volScalarFields;
690  readFields
691  (
692  haveMesh,
693  mesh,
694  subsetterPtr,
695  objects,
696  volScalarFields
697  );
698 
699  PtrList<volVectorField> volVectorFields;
700  readFields
701  (
702  haveMesh,
703  mesh,
704  subsetterPtr,
705  objects,
706  volVectorFields
707  );
708 
709  PtrList<volSphericalTensorField> volSphereTensorFields;
710  readFields
711  (
712  haveMesh,
713  mesh,
714  subsetterPtr,
715  objects,
716  volSphereTensorFields
717  );
718 
719  PtrList<volSymmTensorField> volSymmTensorFields;
720  readFields
721  (
722  haveMesh,
723  mesh,
724  subsetterPtr,
725  objects,
726  volSymmTensorFields
727  );
728 
729  PtrList<volTensorField> volTensorFields;
730  readFields
731  (
732  haveMesh,
733  mesh,
734  subsetterPtr,
735  objects,
736  volTensorFields
737  );
738 
739  PtrList<surfaceScalarField> surfScalarFields;
740  readFields
741  (
742  haveMesh,
743  mesh,
744  subsetterPtr,
745  objects,
746  surfScalarFields
747  );
748 
749  PtrList<surfaceVectorField> surfVectorFields;
750  readFields
751  (
752  haveMesh,
753  mesh,
754  subsetterPtr,
755  objects,
756  surfVectorFields
757  );
758 
759  PtrList<surfaceSphericalTensorField> surfSphereTensorFields;
760  readFields
761  (
762  haveMesh,
763  mesh,
764  subsetterPtr,
765  objects,
766  surfSphereTensorFields
767  );
768 
769  PtrList<surfaceSymmTensorField> surfSymmTensorFields;
770  readFields
771  (
772  haveMesh,
773  mesh,
774  subsetterPtr,
775  objects,
776  surfSymmTensorFields
777  );
778 
779  PtrList<surfaceTensorField> surfTensorFields;
780  readFields
781  (
782  haveMesh,
783  mesh,
784  subsetterPtr,
785  objects,
786  surfTensorFields
787  );
788 
789 
790  // Debugging: Create additional volField that will be mapped.
791  // Used to test correctness of mapping
792  volVectorField mapCc("mapCc", 1*mesh.C());
793 
794  // Global matching tolerance
795  const scalar tolDim = getMergeDistance
796  (
797  args,
798  runTime,
799  mesh.globalData().bb()
800  );
801 
802  // Mesh distribution engine
803  fvMeshDistribute distributor(mesh, tolDim);
804 
805  Pout<< "Wanted distribution:"
806  << distributor.countCells(finalDecomp) << nl << endl;
807 
808  // Do actual sending/receiving of mesh
809  autoPtr<mapDistributePolyMesh> map = distributor.distribute(finalDecomp);
810 
812  //map().distributeFaceData(faceCc);
813 
814 
815  // Print a bit
816  Pout<< "After distribution mesh:" << endl;
817  printMeshData(Pout, mesh);
818  Pout<< endl;
819 
820  runTime++;
821  Pout<< "Writing redistributed mesh to " << runTime.timeName() << nl << endl;
822  mesh.write();
823 
824 
825  // Debugging: test mapped cellcentre field.
826  compareFields(tolDim, mesh.C(), mapCc);
827 
828  // Print nice message
829  // ~~~~~~~~~~~~~~~~~~
830 
831  // Determine which processors remain so we can print nice final message.
832  labelList nFaces(Pstream::nProcs());
833  nFaces[Pstream::myProcNo()] = mesh.nFaces();
834  Pstream::gatherList(nFaces);
835  Pstream::scatterList(nFaces);
836 
837  Info<< nl
838  << "You can pick up the redecomposed mesh from the polyMesh directory"
839  << " in " << runTime.timeName() << "." << nl
840  << "If you redecomposed the mesh to less processors you can delete"
841  << nl
842  << "the processor directories with 0 sized meshes in them." << nl
843  << "Below is a sample set of commands to do this."
844  << " Take care when issuing these" << nl
845  << "commands." << nl << endl;
846 
847  forAll(nFaces, procI)
848  {
849  fileName procDir = "processor" + name(procI);
850 
851  if (nFaces[procI] == 0)
852  {
853  Info<< " rm -r " << procDir.c_str() << nl;
854  }
855  else
856  {
857  fileName timeDir = procDir/runTime.timeName()/meshSubDir;
858  fileName constDir = procDir/runTime.constant()/meshSubDir;
859 
860  Info<< " rm -r " << constDir.c_str() << nl
861  << " mv " << timeDir.c_str()
862  << ' ' << constDir.c_str() << nl;
863  }
864  }
865  Info<< endl;
866 
867 
868  Pout<< "End\n" << endl;
869 
870  return 0;
871 }
872 
873 
874 // ************************ vim: set sw=4 sts=4 et: ************************ //