FreeFOAM The Cross-Platform CFD Toolkit
decomposePar.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  decomposePar
26 
27 Description
28  Automatically decomposes a mesh and fields of a case for parallel
29  execution of OpenFOAM.
30 
31 Usage
32 
33  - decomposePar [OPTION]
34 
35  @param -cellDist \n
36  Write the cell distribution as a labelList for use with 'manual'
37  decomposition method and as a volScalarField for post-processing.
38 
39  @param -region regionName \n
40  Decompose named region. Does not check for existence of processor*.
41 
42  @param -copyUniform \n
43  Copy any @a uniform directories too.
44 
45  @param -fields \n
46  Use existing geometry decomposition and convert fields only.
47 
48  @param -filterPatches \n
49  Remove empty patches when decomposing the geometry.
50 
51  @param -force \n
52  Remove any existing @a processor subdirectories before decomposing the
53  geometry.
54 
55  @param -ifRequired \n
56  Only decompose the geometry if the number of domains has changed from a
57  previous decomposition. No @a processor subdirectories will be removed
58  unless the @a -force option is also specified. This option can be used
59  to avoid redundant geometry decomposition (eg, in scripts), but should
60  be used with caution when the underlying (serial) geometry or the
61  decomposition method etc. have been changed between decompositions.
62 
63  @param -case <dir> \n
64  Case directory.
65 
66  @param -help \n
67  Display help message.
68 
69  @param -doc \n
70  Display Doxygen API documentation page for this application.
71 
72  @param -srcDoc \n
73  Display Doxygen source documentation page for this application.
74 
75 \*---------------------------------------------------------------------------*/
76 
77 #include <OpenFOAM/OSspecific.H>
78 #include <finiteVolume/fvCFD.H>
79 #include <OpenFOAM/IOobjectList.H>
81 #include "domainDecomposition.H"
82 #include <OpenFOAM/labelIOField.H>
83 #include <OpenFOAM/scalarIOField.H>
84 #include <OpenFOAM/vectorIOField.H>
87 #include <OpenFOAM/tensorIOField.H>
88 #include <OpenFOAM/pointFields.H>
89 
90 #include "readFields.H"
91 #include "fvFieldDecomposer.H"
92 #include "pointFieldDecomposer.H"
94 
95 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
96 
97 int main(int argc, char *argv[])
98 {
99  // enable -constant
100  timeSelector::addOptions(true, false);
101  argList::noParallel();
102 # include <OpenFOAM/addRegionOption.H>
103  argList::validOptions.insert("cellDist", "");
104  argList::validOptions.insert("copyUniform", "");
105  argList::validOptions.insert("fields", "");
106  argList::validOptions.insert("filterPatches", "");
107  argList::validOptions.insert("force", "");
108  argList::validOptions.insert("ifRequired", "");
109 
110 # include <OpenFOAM/setRootCase.H>
111 
112  word regionName = fvMesh::defaultRegion;
113  word regionDir = word::null;
114 
115  if (args.optionFound("region"))
116  {
117  regionName = args.option("region");
118  regionDir = regionName;
119  Info<< "Decomposing mesh " << regionName << nl << endl;
120  }
121 
122 
123  bool writeCellDist = args.optionFound("cellDist");
124  bool copyUniform = args.optionFound("copyUniform");
125  bool decomposeFieldsOnly = args.optionFound("fields");
126  bool filterPatches = args.optionFound("filterPatches");
127  bool forceOverwrite = args.optionFound("force");
128  bool ifRequiredDecomposition = args.optionFound("ifRequired");
129 
130 # include <OpenFOAM/createTime.H>
131 
132  // Allow -constant to override controlDict settings.
133  if (args.optionFound("constant"))
134  {
135  instantList timeDirs = timeSelector::select0(runTime, args);
136  if (runTime.timeName() != runTime.constant())
137  {
139  << "No '" << runTime.constant() << "' time present." << endl
140  << "Valid times are " << runTime.times()
141  << exit(FatalError);
142  }
143  }
144 
145 
146  Info<< "Time = " << runTime.timeName() << endl;
147 
148  // determine the existing processor count directly
149  label nProcs = 0;
150  while
151  (
152  isDir
153  (
154  runTime.path()
155  /(word("processor") + name(nProcs))
156  /runTime.constant()
157  /regionDir
158  /polyMesh::meshSubDir
159  )
160  )
161  {
162  ++nProcs;
163  }
164 
165  // get requested numberOfSubdomains
166  label nDomains = 0;
167  {
168  IOdictionary decompDict
169  (
170  IOobject
171  (
172  "decomposeParDict",
173  runTime.time().system(),
174  regionDir, // use region if non-standard
175  runTime,
176  IOobject::MUST_READ,
177  IOobject::NO_WRITE,
178  false
179  )
180  );
181 
182  decompDict.lookup("numberOfSubdomains") >> nDomains;
183  }
184 
185  if (decomposeFieldsOnly)
186  {
187  // Sanity check on previously decomposed case
188  if (nProcs != nDomains)
189  {
191  << "Specified -fields, but the case was decomposed with "
192  << nProcs << " domains"
193  << nl
194  << "instead of " << nDomains
195  << " domains as specified in decomposeParDict"
196  << nl
197  << exit(FatalError);
198  }
199  }
200  else if (nProcs)
201  {
202  bool procDirsProblem = true;
203 
204  if (ifRequiredDecomposition && nProcs == nDomains)
205  {
206  // we can reuse the decomposition
207  decomposeFieldsOnly = true;
208  procDirsProblem = false;
209  forceOverwrite = false;
210 
211  Info<< "Using existing processor directories" << nl;
212  }
213 
214  if (forceOverwrite)
215  {
216  Info<< "Removing " << nProcs
217  << " existing processor directories" << endl;
218 
219  // remove existing processor dirs
220  // reverse order to avoid gaps if someone interrupts the process
221  for (label procI = nProcs-1; procI >= 0; --procI)
222  {
223  fileName procDir
224  (
225  runTime.path()/(word("processor") + name(procI))
226  );
227 
228  rmDir(procDir);
229  }
230 
231  procDirsProblem = false;
232  }
233 
234  if (procDirsProblem)
235  {
237  << "Case is already decomposed with " << nProcs
238  << " domains, use the -force option or manually" << nl
239  << "remove processor directories before decomposing. e.g.,"
240  << nl
241  << " rm -rf " << runTime.path().c_str() << "/processor*"
242  << nl
243  << exit(FatalError);
244  }
245  }
246 
247  Info<< "Create mesh" << endl;
249  (
250  IOobject
251  (
252  regionName,
253  runTime.timeName(),
254  runTime
255  )
256  );
257 
258  // Decompose the mesh
259  if (!decomposeFieldsOnly)
260  {
261  mesh.decomposeMesh(filterPatches);
262 
263  mesh.writeDecomposition();
264 
265  if (writeCellDist)
266  {
267  const labelList& procIds = mesh.cellToProc();
268 
269  // Write the decomposition as labelList for use with 'manual'
270  // decomposition method.
271  labelIOList cellDecomposition
272  (
273  IOobject
274  (
275  "cellDecomposition",
277  mesh,
278  IOobject::NO_READ,
279  IOobject::NO_WRITE,
280  false
281  ),
282  procIds
283  );
284  cellDecomposition.write();
285 
286  Info<< nl << "Wrote decomposition to "
287  << cellDecomposition.objectPath()
288  << " for use in manual decomposition." << endl;
289 
290  // Write as volScalarField for postprocessing.
291  volScalarField cellDist
292  (
293  IOobject
294  (
295  "cellDist",
296  runTime.timeName(),
297  mesh,
298  IOobject::NO_READ,
299  IOobject::AUTO_WRITE
300  ),
301  mesh,
302  dimensionedScalar("cellDist", dimless, 0),
303  zeroGradientFvPatchScalarField::typeName
304  );
305 
306  forAll(procIds, celli)
307  {
308  cellDist[celli] = procIds[celli];
309  }
310 
311  cellDist.write();
312 
313  Info<< nl << "Wrote decomposition as volScalarField to "
314  << cellDist.name() << " for use in postprocessing."
315  << endl;
316  }
317  }
318 
319 
320  // Search for list of objects for this time
321  IOobjectList objects(mesh, runTime.timeName());
322 
323  // Construct the vol fields
324  // ~~~~~~~~~~~~~~~~~~~~~~~~
325  PtrList<volScalarField> volScalarFields;
326  readFields(mesh, objects, volScalarFields);
327 
328  PtrList<volVectorField> volVectorFields;
329  readFields(mesh, objects, volVectorFields);
330 
331  PtrList<volSphericalTensorField> volSphericalTensorFields;
332  readFields(mesh, objects, volSphericalTensorFields);
333 
334  PtrList<volSymmTensorField> volSymmTensorFields;
335  readFields(mesh, objects, volSymmTensorFields);
336 
337  PtrList<volTensorField> volTensorFields;
338  readFields(mesh, objects, volTensorFields);
339 
340 
341  // Construct the surface fields
342  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
343  PtrList<surfaceScalarField> surfaceScalarFields;
344  readFields(mesh, objects, surfaceScalarFields);
345  PtrList<surfaceVectorField> surfaceVectorFields;
346  readFields(mesh, objects, surfaceVectorFields);
347  PtrList<surfaceSphericalTensorField> surfaceSphericalTensorFields;
348  readFields(mesh, objects, surfaceSphericalTensorFields);
349  PtrList<surfaceSymmTensorField> surfaceSymmTensorFields;
350  readFields(mesh, objects, surfaceSymmTensorFields);
351  PtrList<surfaceTensorField> surfaceTensorFields;
352  readFields(mesh, objects, surfaceTensorFields);
353 
354 
355  // Construct the point fields
356  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
357  pointMesh pMesh(mesh);
358 
359  PtrList<pointScalarField> pointScalarFields;
360  readFields(pMesh, objects, pointScalarFields);
361 
362  PtrList<pointVectorField> pointVectorFields;
363  readFields(pMesh, objects, pointVectorFields);
364 
365  PtrList<pointSphericalTensorField> pointSphericalTensorFields;
366  readFields(pMesh, objects, pointSphericalTensorFields);
367 
368  PtrList<pointSymmTensorField> pointSymmTensorFields;
369  readFields(pMesh, objects, pointSymmTensorFields);
370 
371  PtrList<pointTensorField> pointTensorFields;
372  readFields(pMesh, objects, pointTensorFields);
373 
374 
375  // Construct the Lagrangian fields
376  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
377 
378  fileNameList cloudDirs
379  (
380  readDir(runTime.timePath()/cloud::prefix, fileName::DIRECTORY)
381  );
382 
383  // Particles
384  PtrList<Cloud<indexedParticle> > lagrangianPositions(cloudDirs.size());
385  // Particles per cell
386  PtrList< List<SLList<indexedParticle*>*> > cellParticles(cloudDirs.size());
387 
388  PtrList<PtrList<labelIOField> > lagrangianLabelFields(cloudDirs.size());
389  PtrList<PtrList<scalarIOField> > lagrangianScalarFields(cloudDirs.size());
390  PtrList<PtrList<vectorIOField> > lagrangianVectorFields(cloudDirs.size());
391  PtrList<PtrList<sphericalTensorIOField> > lagrangianSphericalTensorFields
392  (
393  cloudDirs.size()
394  );
395  PtrList<PtrList<symmTensorIOField> > lagrangianSymmTensorFields
396  (
397  cloudDirs.size()
398  );
399  PtrList<PtrList<tensorIOField> > lagrangianTensorFields
400  (
401  cloudDirs.size()
402  );
403 
404  label cloudI = 0;
405 
406  forAll(cloudDirs, i)
407  {
408  IOobjectList sprayObjs
409  (
410  mesh,
411  runTime.timeName(),
412  cloud::prefix/cloudDirs[i]
413  );
414 
415  IOobject* positionsPtr = sprayObjs.lookup("positions");
416 
417  if (positionsPtr)
418  {
419  // Read lagrangian particles
420  // ~~~~~~~~~~~~~~~~~~~~~~~~~
421 
422  Info<< "Identified lagrangian data set: " << cloudDirs[i] << endl;
423 
424  lagrangianPositions.set
425  (
426  cloudI,
428  (
429  mesh,
430  cloudDirs[i],
431  false
432  )
433  );
434 
435 
436  // Sort particles per cell
437  // ~~~~~~~~~~~~~~~~~~~~~~~
438 
439  cellParticles.set
440  (
441  cloudI,
443  (
444  mesh.nCells(),
445  static_cast<SLList<indexedParticle*>*>(NULL)
446  )
447  );
448 
449  label i = 0;
450 
451  forAllIter
452  (
454  lagrangianPositions[cloudI],
455  iter
456  )
457  {
458  iter().index() = i++;
459 
460  label celli = iter().cell();
461 
462  // Check
463  if (celli < 0 || celli >= mesh.nCells())
464  {
466  << "Illegal cell number " << celli
467  << " for particle with index " << iter().index()
468  << " at position " << iter().position() << nl
469  << "Cell number should be between 0 and "
470  << mesh.nCells()-1 << nl
471  << "On this mesh the particle should be in cell "
472  << mesh.findCell(iter().position())
473  << exit(FatalError);
474  }
475 
476  if (!cellParticles[cloudI][celli])
477  {
478  cellParticles[cloudI][celli] = new SLList<indexedParticle*>
479  ();
480  }
481 
482  cellParticles[cloudI][celli]->append(&iter());
483  }
484 
485  // Read fields
486  // ~~~~~~~~~~~
487 
488  IOobjectList lagrangianObjects
489  (
490  mesh,
491  runTime.timeName(),
492  cloud::prefix/cloudDirs[cloudI]
493  );
494 
496  (
497  cloudI,
498  lagrangianObjects,
499  lagrangianLabelFields
500  );
501 
503  (
504  cloudI,
505  lagrangianObjects,
506  lagrangianScalarFields
507  );
508 
510  (
511  cloudI,
512  lagrangianObjects,
513  lagrangianVectorFields
514  );
515 
517  (
518  cloudI,
519  lagrangianObjects,
520  lagrangianSphericalTensorFields
521  );
522 
524  (
525  cloudI,
526  lagrangianObjects,
527  lagrangianSymmTensorFields
528  );
529 
531  (
532  cloudI,
533  lagrangianObjects,
534  lagrangianTensorFields
535  );
536 
537  cloudI++;
538  }
539  }
540 
541  lagrangianPositions.setSize(cloudI);
542  cellParticles.setSize(cloudI);
543  lagrangianLabelFields.setSize(cloudI);
544  lagrangianScalarFields.setSize(cloudI);
545  lagrangianVectorFields.setSize(cloudI);
546  lagrangianSphericalTensorFields.setSize(cloudI);
547  lagrangianSymmTensorFields.setSize(cloudI);
548  lagrangianTensorFields.setSize(cloudI);
549 
550 
551  // Any uniform data to copy/link?
552  fileName uniformDir("uniform");
553 
554  if (isDir(runTime.timePath()/uniformDir))
555  {
556  Info<< "Detected additional non-decomposed files in "
557  << runTime.timePath()/uniformDir
558  << endl;
559  }
560  else
561  {
562  uniformDir.clear();
563  }
564 
565  Info<< endl;
566 
567  // split the fields over processors
568  for (label procI = 0; procI < mesh.nProcs(); procI++)
569  {
570  Info<< "Processor " << procI << ": field transfer" << endl;
571 
572  // open the database
573  Time processorDb
574  (
575  Time::controlDictName,
576  args.rootPath(),
577  args.caseName()/fileName(word("processor") + name(procI))
578  );
579 
580  processorDb.setTime(runTime);
581 
582  // remove files remnants that can cause horrible problems
583  // - mut and nut are used to mark the new turbulence models,
584  // their existence prevents old models from being upgraded
585  {
586  fileName timeDir(processorDb.path()/processorDb.timeName());
587 
588  rm(timeDir/"mut");
589  rm(timeDir/"nut");
590  }
591 
592  // read the mesh
593  fvMesh procMesh
594  (
595  IOobject
596  (
597  regionName,
598  processorDb.timeName(),
599  processorDb
600  )
601  );
602 
603  labelIOList cellProcAddressing
604  (
605  IOobject
606  (
607  "cellProcAddressing",
608  procMesh.facesInstance(),
609  procMesh.meshSubDir,
610  procMesh,
611  IOobject::MUST_READ,
612  IOobject::NO_WRITE
613  )
614  );
615 
616  labelIOList boundaryProcAddressing
617  (
618  IOobject
619  (
620  "boundaryProcAddressing",
621  procMesh.facesInstance(),
622  procMesh.meshSubDir,
623  procMesh,
624  IOobject::MUST_READ,
625  IOobject::NO_WRITE
626  )
627  );
628 
629  // FV fields
630  if
631  (
632  volScalarFields.size()
633  || volVectorFields.size()
634  || volSphericalTensorFields.size()
635  || volSymmTensorFields.size()
636  || volTensorFields.size()
637  || surfaceScalarFields.size()
638  || surfaceVectorFields.size()
639  || surfaceSphericalTensorFields.size()
640  || surfaceSymmTensorFields.size()
641  || surfaceTensorFields.size()
642  )
643  {
644  labelIOList faceProcAddressing
645  (
646  IOobject
647  (
648  "faceProcAddressing",
649  procMesh.facesInstance(),
650  procMesh.meshSubDir,
651  procMesh,
652  IOobject::MUST_READ,
653  IOobject::NO_WRITE
654  )
655  );
656 
657  fvFieldDecomposer fieldDecomposer
658  (
659  mesh,
660  procMesh,
661  faceProcAddressing,
662  cellProcAddressing,
663  boundaryProcAddressing
664  );
665 
666  fieldDecomposer.decomposeFields(volScalarFields);
667  fieldDecomposer.decomposeFields(volVectorFields);
668  fieldDecomposer.decomposeFields(volSphericalTensorFields);
669  fieldDecomposer.decomposeFields(volSymmTensorFields);
670  fieldDecomposer.decomposeFields(volTensorFields);
671 
672  fieldDecomposer.decomposeFields(surfaceScalarFields);
673  fieldDecomposer.decomposeFields(surfaceVectorFields);
674  fieldDecomposer.decomposeFields(surfaceSphericalTensorFields);
675  fieldDecomposer.decomposeFields(surfaceSymmTensorFields);
676  fieldDecomposer.decomposeFields(surfaceTensorFields);
677  }
678 
679 
680  // Point fields
681  if
682  (
683  pointScalarFields.size()
684  || pointVectorFields.size()
685  || pointSphericalTensorFields.size()
686  || pointSymmTensorFields.size()
687  || pointTensorFields.size()
688  )
689  {
690  labelIOList pointProcAddressing
691  (
692  IOobject
693  (
694  "pointProcAddressing",
695  procMesh.facesInstance(),
696  procMesh.meshSubDir,
697  procMesh,
698  IOobject::MUST_READ,
699  IOobject::NO_WRITE
700  )
701  );
702 
703  pointMesh procPMesh(procMesh, true);
704 
705  pointFieldDecomposer fieldDecomposer
706  (
707  pMesh,
708  procPMesh,
709  pointProcAddressing,
710  boundaryProcAddressing
711  );
712 
713  fieldDecomposer.decomposeFields(pointScalarFields);
714  fieldDecomposer.decomposeFields(pointVectorFields);
715  fieldDecomposer.decomposeFields(pointSphericalTensorFields);
716  fieldDecomposer.decomposeFields(pointSymmTensorFields);
717  fieldDecomposer.decomposeFields(pointTensorFields);
718  }
719 
720 
721  // If there is lagrangian data write it out
722  forAll(lagrangianPositions, cloudI)
723  {
724  if (lagrangianPositions[cloudI].size())
725  {
726  lagrangianFieldDecomposer fieldDecomposer
727  (
728  mesh,
729  procMesh,
730  cellProcAddressing,
731  cloudDirs[cloudI],
732  lagrangianPositions[cloudI],
733  cellParticles[cloudI]
734  );
735 
736  // Lagrangian fields
737  if
738  (
739  lagrangianLabelFields[cloudI].size()
740  || lagrangianScalarFields[cloudI].size()
741  || lagrangianVectorFields[cloudI].size()
742  || lagrangianSphericalTensorFields[cloudI].size()
743  || lagrangianSymmTensorFields[cloudI].size()
744  || lagrangianTensorFields[cloudI].size()
745  )
746  {
747  fieldDecomposer.decomposeFields
748  (
749  cloudDirs[cloudI],
750  lagrangianLabelFields[cloudI]
751  );
752  fieldDecomposer.decomposeFields
753  (
754  cloudDirs[cloudI],
755  lagrangianScalarFields[cloudI]
756  );
757  fieldDecomposer.decomposeFields
758  (
759  cloudDirs[cloudI],
760  lagrangianVectorFields[cloudI]
761  );
762  fieldDecomposer.decomposeFields
763  (
764  cloudDirs[cloudI],
765  lagrangianSphericalTensorFields[cloudI]
766  );
767  fieldDecomposer.decomposeFields
768  (
769  cloudDirs[cloudI],
770  lagrangianSymmTensorFields[cloudI]
771  );
772  fieldDecomposer.decomposeFields
773  (
774  cloudDirs[cloudI],
775  lagrangianTensorFields[cloudI]
776  );
777  }
778  }
779  }
780 
781 
782  // Any non-decomposed data to copy?
783  if (uniformDir.size())
784  {
785  const fileName timePath = processorDb.timePath();
786 
787  if (copyUniform || mesh.distributed())
788  {
789  cp
790  (
791  runTime.timePath()/uniformDir,
792  timePath/uniformDir
793  );
794  }
795  else
796  {
797  // link with relative paths
798  const string parentPath = string("..")/"..";
799 
800  fileName currentDir(cwd());
801  chDir(timePath);
802  ln
803  (
804  parentPath/runTime.timeName()/uniformDir,
805  uniformDir
806  );
807  chDir(currentDir);
808  }
809  }
810  }
811 
812 
813  Info<< "\nEnd.\n" << endl;
814 
815  return 0;
816 }
817 
818 
819 // ************************ vim: set sw=4 sts=4 et: ************************ //