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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2822 - (show annotations)
Mon Dec 14 03:10:16 2009 UTC (9 years, 9 months ago) by caltinay
File size: 18676 byte(s)
dataexporter: Ensure that NetCDF files created using N ranks are only loaded
when mpi size equals N or 1.

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

  ViewVC Help
Powered by ViewVC 1.1.26