/[escript]/trunk/dataexporter/src/EscriptDataset.cpp
ViewVC logotype

Contents of /trunk/dataexporter/src/EscriptDataset.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2847 - (show annotations)
Fri Jan 15 03:35:05 2010 UTC (9 years, 6 months ago) by caltinay
File size: 19335 byte(s)
Improved error handling when using MPI and fixed Silo output when not all
ranks contribute to a data variable.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include <escriptexport/EscriptDataset.h>
15 #include <escriptexport/DataVar.h>
16 #include <escriptexport/ElementData.h>
17 #include <escriptexport/FinleyMesh.h>
18 #include <escriptexport/NodeData.h>
19 #include <escript/Data.h>
20
21 #include <iostream>
22
23 #if USE_SILO
24 #include <silo.h>
25
26 #if HAVE_MPI
27 #include <pmpio.h>
28 #endif
29
30 #endif // USE_SILO
31
32 using namespace std;
33
34 namespace escriptexport {
35
36 const char* MESH_VARS = "mesh_vars/";
37 const int NUM_SILO_FILES = 1;
38
39 //
40 // Default constructor
41 //
42 EscriptDataset::EscriptDataset() :
43 numParts(0),
44 cycle(0),
45 time(0.),
46 keepMesh(false),
47 externalMesh(false),
48 mpiRank(0),
49 mpiSize(1)
50 {
51 }
52
53 //
54 // Constructor with communicator
55 //
56 #if HAVE_MPI
57 EscriptDataset::EscriptDataset(MPI_Comm comm) :
58 numParts(0),
59 cycle(0),
60 time(0.),
61 keepMesh(false),
62 externalMesh(false),
63 mpiComm(comm)
64 {
65 MPI_Comm_rank(mpiComm, &mpiRank);
66 MPI_Comm_size(mpiComm, &mpiSize);
67 }
68 #endif
69
70 //
71 // Destructor
72 //
73 EscriptDataset::~EscriptDataset()
74 {
75 }
76
77 //
78 //
79 //
80 bool EscriptDataset::initFromEscript(escript::const_Domain_ptr escriptDomain,
81 DataVec& escriptVars,
82 const StringVec& varNames)
83 {
84 int myError;
85 numParts = 1;
86 FinleyMesh_ptr mesh(new FinleyMesh());
87 if (mesh->initFromEscript(escriptDomain)) {
88 if (mpiSize > 1)
89 mesh->reorderGhostZones(mpiRank);
90 meshBlocks.push_back(mesh);
91 myError = false;
92 } else {
93 mesh.reset();
94 myError = true;
95 }
96
97 int gError;
98 #if HAVE_MPI
99 MPI_Allreduce(&myError, &gError, 1, MPI_LOGICAL, MPI_LOR, mpiComm);
100 #else
101 gError = myError;
102 #endif
103
104 if (!gError) {
105 // initialize variables
106 DataVec::iterator varIt = escriptVars.begin();
107 StringVec::const_iterator nameIt = varNames.begin();
108 for (; varIt != escriptVars.end(); varIt++, nameIt++) {
109 VarInfo vi;
110 vi.varName = *nameIt;
111
112 DataVar_ptr var(new DataVar(vi.varName));
113 if (var->initFromEscript(*varIt, mesh)) {
114 vi.dataBlocks.push_back(var);
115 updateSampleDistribution(vi);
116 vi.valid = true;
117 } else {
118 var.reset();
119 vi.valid = false;
120 }
121 variables.push_back(vi);
122 }
123
124 // Convert mesh data to variables
125 convertMeshVariables();
126 }
127
128 return !gError;
129 }
130
131 //
132 //
133 //
134 bool EscriptDataset::loadNetCDF(const string meshFilePattern,
135 const StringVec& varFiles,
136 const StringVec& varNames, int nBlocks)
137 {
138 if (mpiSize > 1 && nBlocks != mpiSize) {
139 cerr << "Cannot load " << nBlocks << " chunks on " << mpiSize
140 << " MPI ranks!" << endl;
141 return false;
142 }
143
144 numParts = nBlocks;
145 meshFmt = meshFilePattern;
146
147 // Load the mesh files
148 if (!loadMeshFromNetCDF()) {
149 cerr << "Reading the domain failed." << endl;
150 return false;
151 }
152
153 // Convert mesh data to variables
154 convertMeshVariables();
155
156 // Load the variables
157 StringVec::const_iterator fileIt = varFiles.begin();
158 StringVec::const_iterator nameIt = varNames.begin();
159 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
160 loadVarFromNetCDF(*fileIt, *nameIt);
161 }
162
163 return true;
164 }
165
166 //
167 // Load only variables using provided mesh
168 //
169 bool EscriptDataset::loadNetCDF(const MeshBlocks& mesh,
170 const StringVec& varFiles,
171 const StringVec& varNames)
172 {
173 if (mpiSize > 1 && mesh.size() > 1) {
174 cerr << "Can only read one mesh block per rank when using MPI!" << endl;
175 return false;
176 }
177
178 externalMesh = true;
179 meshBlocks = mesh;
180 numParts = meshBlocks.size();
181 keepMesh = true;
182
183 // Load the variables
184 StringVec::const_iterator fileIt = varFiles.begin();
185 StringVec::const_iterator nameIt = varNames.begin();
186 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
187 loadVarFromNetCDF(*fileIt, *nameIt);
188 }
189
190 return true;
191 }
192
193 //
194 //
195 //
196 bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)
197 {
198 #if USE_SILO
199 if (numParts == 0)
200 return false;
201
202 const char* blockDirFmt = "/block%04d";
203 string siloPath;
204 DBfile* dbfile = NULL;
205 #if HAVE_MPI
206 PMPIO_baton_t* baton = NULL;
207 #endif
208
209 if (mpiSize > 1) {
210 #if HAVE_MPI
211 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
212 mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
213 PMPIO_DefaultClose, NULL);
214 if (baton) {
215 char str[64];
216 snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, mpiRank));
217 siloPath = str;
218 dbfile = (DBfile*) PMPIO_WaitForBaton(
219 baton, fileName.c_str(), siloPath.c_str());
220 }
221 #endif
222 } else {
223 dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
224 "escriptData", DB_PDB);
225 }
226
227 if (!dbfile) {
228 cerr << "Could not create Silo file." << endl;
229 return false;
230 }
231
232 MeshBlocks::iterator meshIt;
233 VarVector::iterator viIt;
234 int idx = 0;
235 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
236 if (mpiSize == 1) {
237 char str[64];
238 snprintf(str, 64, blockDirFmt, idx);
239 siloPath = str;
240 DBMkdir(dbfile, siloPath.c_str());
241 }
242 // write block of the mesh if we don't use an external mesh
243 if (!externalMesh) {
244 if (! (*meshIt)->writeToSilo(dbfile, siloPath)) {
245 cerr << "Error writing block " << idx << " of mesh to Silo file!\n";
246 break;
247 }
248 }
249
250 // write variables for current mesh block
251 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
252 // do not attempt to write this variable if previous steps failed
253 if (!(*viIt).valid) continue;
254 DataVar_ptr var = (*viIt).dataBlocks[idx];
255 if (!var->writeToSilo(dbfile, siloPath)) {
256 cerr << "Error writing block " << idx << " of '"
257 << var->getName() << "' to Silo file!" << endl;
258 (*viIt).valid = false;
259 }
260 }
261 }
262
263 // rank 0 writes additional data that describe how the parts fit together
264 if (mpiRank == 0) {
265 if (useMultiMesh) {
266 const StringVec& meshNames = meshBlocks[0]->getMeshNames();
267 StringVec::const_iterator it;
268 for (it = meshNames.begin(); it != meshNames.end(); it++)
269 putSiloMultiMesh(dbfile, *it);
270
271 DBMkdir(dbfile, MESH_VARS);
272 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
273 putSiloMultiVar(dbfile, *viIt, true);
274
275 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
276 if (!(*viIt).valid) continue;
277 DataVar_ptr var = (*viIt).dataBlocks[0];
278 if (var->getRank() < 2)
279 putSiloMultiVar(dbfile, *viIt);
280 else
281 putSiloMultiTensor(dbfile, *viIt);
282 }
283 }
284
285 vector<char*> tensorNames;
286 vector<string> tensorDefStrings;
287 vector<char*> tensorDefs;
288
289 // collect tensors for their Silo definitions
290 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
291 if (!(*viIt).valid) continue;
292 DataVar_ptr var = (*viIt).dataBlocks[0];
293 if (var->getRank() == 2) {
294 tensorDefStrings.push_back(var->getTensorDef());
295 tensorDefs.push_back((char*)tensorDefStrings.back().c_str());
296 tensorNames.push_back((char*)var->getName().c_str());
297 }
298 }
299
300 if (tensorDefs.size()) {
301 DBoptlist* optList = DBMakeOptlist(2);
302 DBAddOption(optList, DBOPT_CYCLE, &cycle);
303 DBAddOption(optList, DBOPT_DTIME, &time);
304 vector<DBoptlist*> defOpts(tensorDefs.size(), optList);
305 vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
306 DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
307 &defTypes[0], &tensorDefs[0], &defOpts[0]);
308 DBFreeOptlist(optList);
309 }
310 }
311
312 if (mpiSize > 1) {
313 #if HAVE_MPI
314 PMPIO_HandOffBaton(baton, dbfile);
315 PMPIO_Finish(baton);
316 #endif
317 } else {
318 DBClose(dbfile);
319 }
320
321 return true;
322
323 #else // !USE_SILO
324 return false;
325 #endif
326 }
327
328 //
329 //
330 //
331 bool EscriptDataset::saveVTK(string fileName)
332 {
333 #if 0
334 if (numParts == 0)
335 return false;
336
337 int globalNumPointsAndCells[2] = { 0, 0 };
338
339 #if HAVE_MPI
340 int myNumPointsAndCells[2];
341 meshBlocks[0]->removeGhostZones(mpiRank);
342 ElementData_ptr elements = meshBlocks[0]->getElements();
343 myNumPointsAndCells[0] = elements->getNodeMesh()->getNumNodes();
344 myNumPointsAndCells[1] = elements->getNumElements();
345 MPI_Reduce(&myNumPointsAndCells[0], &globalNumPointsAndCells[0], 2,
346 MPI_INT, MPI_SUM, 0, mpiComm);
347 #else
348 MeshBlocks::iterator meshIt;
349 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {
350 ElementData_ptr elements = (*meshIt)->getElements();
351 globalNumPointsAndCells[0] += elements->getNodeMesh()->getNumNodes();
352 globalNumPointsAndCells[1] += elements->getNumElements();
353 }
354 #endif
355 #endif
356 return false;
357 }
358
359 //
360 //
361 //
362 bool EscriptDataset::loadMeshFromNetCDF()
363 {
364 bool ok = true;
365 char* str = new char[meshFmt.length()+10];
366 if (mpiSize > 1) {
367 FinleyMesh_ptr meshPart(new FinleyMesh());
368 sprintf(str, meshFmt.c_str(), mpiRank);
369 string meshfile = str;
370 if (meshPart->initFromNetCDF(meshfile)) {
371 meshPart->reorderGhostZones(mpiRank);
372 meshBlocks.push_back(meshPart);
373 } else {
374 meshPart.reset();
375 ok = false;
376 }
377 } else {
378 for (int idx=0; idx < numParts; idx++) {
379 FinleyMesh_ptr meshPart(new FinleyMesh());
380 sprintf(str, meshFmt.c_str(), idx);
381 string meshfile = str;
382 if (meshPart->initFromNetCDF(meshfile)) {
383 if (numParts > 1)
384 meshPart->reorderGhostZones(idx);
385 meshBlocks.push_back(meshPart);
386 } else {
387 meshPart.reset();
388 ok = false;
389 break;
390 }
391 }
392 }
393 delete[] str;
394 return ok;
395 }
396
397 //
398 //
399 //
400 void EscriptDataset::convertMeshVariables()
401 {
402 const StringVec& varNames = meshBlocks[0]->getVarNames();
403 StringVec::const_iterator it;
404 for (it = varNames.begin(); it != varNames.end(); it++) {
405 VarInfo vi;
406 vi.varName = *it;
407 vi.valid = true;
408 // get all parts of current variable
409 MeshBlocks::iterator mIt;
410 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++) {
411 DataVar_ptr var(new DataVar(*it));
412 if (var->initFromMesh(*mIt)) {
413 vi.dataBlocks.push_back(var);
414 } else {
415 cerr << "Error converting mesh variable " << *it << endl;
416 vi.valid = false;
417 break;
418 }
419 }
420 updateSampleDistribution(vi);
421 meshVariables.push_back(vi);
422 }
423 }
424
425 //
426 //
427 //
428 bool EscriptDataset::loadVarFromNetCDF(const string& fileName,
429 const string& varName)
430 {
431 int myError = false;
432 char* str = new char[fileName.length()+10];
433 VarInfo vi;
434
435 vi.varName = varName;
436 vi.valid = true;
437
438 // read all parts of the variable
439 MeshBlocks::iterator mIt;
440 int idx = (mpiSize > 1) ? mpiRank : 0;
441 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++, idx++) {
442 sprintf(str, fileName.c_str(), idx);
443 string dfile = str;
444 DataVar_ptr var(new DataVar(varName));
445 if (var->initFromNetCDF(dfile, *mIt))
446 vi.dataBlocks.push_back(var);
447 else {
448 cerr << "Error reading " << dfile << endl;
449 var.reset();
450 myError = true;
451 break;
452 }
453 }
454 delete[] str;
455
456 int gError;
457 #if HAVE_MPI
458 MPI_Allreduce(&myError, &gError, 1, MPI_LOGICAL, MPI_LOR, mpiComm);
459 #else
460 gError = myError;
461 #endif
462
463 if (gError) {
464 // at least one chunk was not read correctly
465 vi.dataBlocks.clear();
466 vi.valid = false;
467 }
468
469 updateSampleDistribution(vi);
470 variables.push_back(vi);
471 return vi.valid;
472 }
473
474 // returns the number of samples at each block - used to determine which
475 // blocks contribute to given variable if any.
476 void EscriptDataset::updateSampleDistribution(VarInfo& vi)
477 {
478 IntVec sampleDist;
479 int numBlocks = 0;
480 const DataBlocks& varBlocks = vi.dataBlocks;
481
482 if (mpiSize > 1) {
483 #if HAVE_MPI
484 int myNumSamples = varBlocks[0]->getNumberOfSamples();
485 if (mpiRank == 0) {
486 numBlocks = mpiSize;
487 sampleDist.insert(sampleDist.end(), numBlocks, 0);
488 }
489 MPI_Gather(
490 &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, 0, mpiComm);
491 #endif
492 } else {
493 numBlocks = varBlocks.size();
494 DataBlocks::const_iterator it;
495 for (it = varBlocks.begin(); it != varBlocks.end(); it++) {
496 sampleDist.push_back((*it)->getNumberOfSamples());
497 }
498 }
499 vi.sampleDistribution = sampleDist;
500 }
501
502 //
503 //
504 //
505 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
506 {
507 #if USE_SILO
508 int numBlocks = 0;
509
510 if (mpiSize > 1) {
511 // FIXME: empty ranks are not accounted for
512 numBlocks = mpiSize;
513 } else {
514 MeshBlocks::iterator meshIt;
515 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {
516 const StringVec& meshNames = (*meshIt)->getMeshNames();
517 if (find(meshNames.begin(), meshNames.end(), meshName) != meshNames.end()) {
518 numBlocks++;
519 }
520 }
521 }
522 vector<int> meshtypes(numBlocks, DB_UCDMESH);
523 vector<string> tempstrings;
524 vector<char*> meshnames;
525 string pathPrefix;
526 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
527 if (ppIndex != string::npos) {
528 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
529 }
530
531 for (size_t idx = 0; idx < numBlocks; idx++) {
532 stringstream siloPath;
533 siloPath << pathPrefix << "/block";
534 int prevWidth = siloPath.width(4);
535 char prevFill = siloPath.fill('0');
536 siloPath << right << idx;
537 siloPath.width(prevWidth);
538 siloPath.fill(prevFill);
539 siloPath << "/";
540 siloPath << meshName;
541 tempstrings.push_back(siloPath.str());
542 meshnames.push_back((char*)tempstrings.back().c_str());
543 }
544
545 // ignore empty mesh
546 if (meshnames.size() > 0) {
547 DBoptlist* optList = DBMakeOptlist(2);
548 DBAddOption(optList, DBOPT_CYCLE, &cycle);
549 DBAddOption(optList, DBOPT_DTIME, &time);
550 DBPutMultimesh(dbfile, meshName.c_str(), numBlocks, &meshnames[0],
551 &meshtypes[0], optList);
552 DBFreeOptlist(optList);
553 }
554 #endif
555 }
556
557 //
558 //
559 //
560 void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
561 bool useMeshFile)
562 {
563 #if USE_SILO
564 vector<int> vartypes;
565 vector<string> tempstrings;
566 vector<char*> varnames;
567 string pathPrefix;
568 if (useMeshFile) {
569 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
570 if (ppIndex != string::npos) {
571 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
572 }
573 }
574
575 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
576 if (vi.sampleDistribution[idx] > 0) {
577 stringstream siloPath;
578 siloPath << pathPrefix << "/block";
579 int prevWidth = siloPath.width(4);
580 char prevFill = siloPath.fill('0');
581 siloPath << right << idx;
582 siloPath.width(prevWidth);
583 siloPath.fill(prevFill);
584 siloPath << "/";
585 siloPath << vi.varName;
586 tempstrings.push_back(siloPath.str());
587 varnames.push_back((char*)tempstrings.back().c_str());
588 vartypes.push_back(DB_UCDVAR);
589 }
590 }
591
592 // ignore empty variables
593 if (varnames.size() > 0) {
594 DBoptlist* optList = DBMakeOptlist(2);
595 DBAddOption(optList, DBOPT_CYCLE, &cycle);
596 DBAddOption(optList, DBOPT_DTIME, &time);
597 if (useMeshFile) {
598 string vpath = string(MESH_VARS)+vi.varName;
599 DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
600 &varnames[0], &vartypes[0], optList);
601 } else {
602 DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
603 &varnames[0], &vartypes[0], optList);
604 }
605 DBFreeOptlist(optList);
606 }
607 #endif
608 }
609
610 //
611 //
612 //
613 void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
614 {
615 #if USE_SILO
616 int numBlocks = 0;
617 if (mpiSize > 1) {
618 // FIXME: empty ranks are not accounted for
619 numBlocks = mpiSize;
620 } else {
621 DataBlocks::const_iterator it;
622 for (it = vi.dataBlocks.begin(); it != vi.dataBlocks.end(); it++) {
623 if ((*it)->getNumberOfSamples() > 0)
624 numBlocks++;
625 }
626 }
627 string tensorDir = vi.varName+string("_comps/");
628 DBSetDir(dbfile, "/");
629 DBMkdir(dbfile, tensorDir.c_str());
630 int one = 1;
631 DBoptlist* optList = DBMakeOptlist(3);
632 DBAddOption(optList, DBOPT_CYCLE, &cycle);
633 DBAddOption(optList, DBOPT_DTIME, &time);
634 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
635 vector<int> vartypes(numBlocks, DB_UCDVAR);
636 const IntVec& shape = vi.dataBlocks[0]->getShape();
637 for (int i=0; i<shape[1]; i++) {
638 for (int j=0; j<shape[0]; j++) {
639 vector<string> tempstrings;
640 vector<char*> varnames;
641 stringstream comp;
642 comp << vi.varName << "_comps/a_";
643 comp << i;
644 comp << j;
645 for (size_t idx = 0; idx < numBlocks; idx++) {
646 stringstream siloPath;
647 siloPath << "/block";
648 int prevWidth = siloPath.width(4);
649 char prevFill = siloPath.fill('0');
650 siloPath << right << idx;
651 siloPath.width(prevWidth);
652 siloPath.fill(prevFill);
653 siloPath << "/" << comp.str();
654 tempstrings.push_back(siloPath.str());
655 varnames.push_back((char*)tempstrings.back().c_str());
656 }
657 DBPutMultivar(dbfile, comp.str().c_str(), numBlocks, &varnames[0],
658 &vartypes[0], optList);
659 }
660 }
661 DBFreeOptlist(optList);
662 #endif
663 }
664
665 } // namespace escriptexport
666

  ViewVC Help
Powered by ViewVC 1.1.26