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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 4 months ago) by jfenwick
File size: 19574 byte(s)
Don't panic.
Updating copyright stamps

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 DBSetDir(dbfile, "/");
302 DBoptlist* optList = DBMakeOptlist(2);
303 DBAddOption(optList, DBOPT_CYCLE, &cycle);
304 DBAddOption(optList, DBOPT_DTIME, &time);
305 vector<DBoptlist*> defOpts(tensorDefs.size(), optList);
306 vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
307 DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
308 &defTypes[0], &tensorDefs[0], &defOpts[0]);
309 DBFreeOptlist(optList);
310 }
311 }
312
313 if (mpiSize > 1) {
314 #if HAVE_MPI
315 PMPIO_HandOffBaton(baton, dbfile);
316 PMPIO_Finish(baton);
317 #endif
318 } else {
319 DBClose(dbfile);
320 }
321
322 return true;
323
324 #else // !USE_SILO
325 return false;
326 #endif
327 }
328
329 //
330 //
331 //
332 bool EscriptDataset::saveVTK(string fileName)
333 {
334 #if 0
335 if (numParts == 0)
336 return false;
337
338 int globalNumPointsAndCells[2] = { 0, 0 };
339
340 #if HAVE_MPI
341 int myNumPointsAndCells[2];
342 meshBlocks[0]->removeGhostZones(mpiRank);
343 ElementData_ptr elements = meshBlocks[0]->getElements();
344 myNumPointsAndCells[0] = elements->getNodeMesh()->getNumNodes();
345 myNumPointsAndCells[1] = elements->getNumElements();
346 MPI_Reduce(&myNumPointsAndCells[0], &globalNumPointsAndCells[0], 2,
347 MPI_INT, MPI_SUM, 0, mpiComm);
348 #else
349 MeshBlocks::iterator meshIt;
350 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {
351 ElementData_ptr elements = (*meshIt)->getElements();
352 globalNumPointsAndCells[0] += elements->getNodeMesh()->getNumNodes();
353 globalNumPointsAndCells[1] += elements->getNumElements();
354 }
355 #endif
356 #endif
357 return false;
358 }
359
360 //
361 //
362 //
363 bool EscriptDataset::loadMeshFromNetCDF()
364 {
365 bool ok = true;
366 char* str = new char[meshFmt.length()+10];
367 if (mpiSize > 1) {
368 FinleyMesh_ptr meshPart(new FinleyMesh());
369 sprintf(str, meshFmt.c_str(), mpiRank);
370 string meshfile = str;
371 if (meshPart->initFromNetCDF(meshfile)) {
372 meshPart->reorderGhostZones(mpiRank);
373 meshBlocks.push_back(meshPart);
374 } else {
375 meshPart.reset();
376 ok = false;
377 }
378 } else {
379 for (int idx=0; idx < numParts; idx++) {
380 FinleyMesh_ptr meshPart(new FinleyMesh());
381 sprintf(str, meshFmt.c_str(), idx);
382 string meshfile = str;
383 if (meshPart->initFromNetCDF(meshfile)) {
384 if (numParts > 1)
385 meshPart->reorderGhostZones(idx);
386 meshBlocks.push_back(meshPart);
387 } else {
388 meshPart.reset();
389 ok = false;
390 break;
391 }
392 }
393 }
394 delete[] str;
395 return ok;
396 }
397
398 //
399 //
400 //
401 void EscriptDataset::convertMeshVariables()
402 {
403 const StringVec& varNames = meshBlocks[0]->getVarNames();
404 StringVec::const_iterator it;
405 for (it = varNames.begin(); it != varNames.end(); it++) {
406 VarInfo vi;
407 vi.varName = *it;
408 vi.valid = true;
409 // get all parts of current variable
410 MeshBlocks::iterator mIt;
411 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++) {
412 DataVar_ptr var(new DataVar(*it));
413 if (var->initFromMesh(*mIt)) {
414 vi.dataBlocks.push_back(var);
415 } else {
416 cerr << "Error converting mesh variable " << *it << endl;
417 vi.valid = false;
418 break;
419 }
420 }
421 updateSampleDistribution(vi);
422 meshVariables.push_back(vi);
423 }
424 }
425
426 //
427 //
428 //
429 bool EscriptDataset::loadVarFromNetCDF(const string& fileName,
430 const string& varName)
431 {
432 int myError = false;
433 char* str = new char[fileName.length()+10];
434 VarInfo vi;
435
436 vi.varName = varName;
437 vi.valid = true;
438
439 // read all parts of the variable
440 MeshBlocks::iterator mIt;
441 int idx = (mpiSize > 1) ? mpiRank : 0;
442 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++, idx++) {
443 sprintf(str, fileName.c_str(), idx);
444 string dfile = str;
445 DataVar_ptr var(new DataVar(varName));
446 if (var->initFromNetCDF(dfile, *mIt))
447 vi.dataBlocks.push_back(var);
448 else {
449 cerr << "Error reading " << dfile << endl;
450 var.reset();
451 myError = true;
452 break;
453 }
454 }
455 delete[] str;
456
457 int gError;
458 #if HAVE_MPI
459 MPI_Allreduce(&myError, &gError, 1, MPI_LOGICAL, MPI_LOR, mpiComm);
460 #else
461 gError = myError;
462 #endif
463
464 if (gError) {
465 // at least one chunk was not read correctly
466 vi.dataBlocks.clear();
467 vi.valid = false;
468 }
469
470 updateSampleDistribution(vi);
471 variables.push_back(vi);
472 return vi.valid;
473 }
474
475 // returns the number of samples at each block - used to determine which
476 // blocks contribute to given variable if any.
477 void EscriptDataset::updateSampleDistribution(VarInfo& vi)
478 {
479 IntVec sampleDist;
480 const DataBlocks& varBlocks = vi.dataBlocks;
481
482 if (mpiSize > 1) {
483 #if HAVE_MPI
484 int myNumSamples = varBlocks[0]->getNumberOfSamples();
485 sampleDist.insert(sampleDist.end(), mpiSize, 0);
486 MPI_Allgather(
487 &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, mpiComm);
488 #endif
489 } else {
490 DataBlocks::const_iterator it;
491 for (it = varBlocks.begin(); it != varBlocks.end(); it++) {
492 sampleDist.push_back((*it)->getNumberOfSamples());
493 }
494 }
495 vi.sampleDistribution = sampleDist;
496 }
497
498 //
499 //
500 //
501 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
502 {
503 #if USE_SILO
504 vector<int> meshtypes;
505 vector<string> tempstrings;
506 vector<char*> meshnames;
507 string pathPrefix;
508
509 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
510 if (ppIndex != string::npos) {
511 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
512 }
513
514 // find a variable belonging to this mesh to get the sample
515 // distribution (which tells us which ranks contribute to this mesh).
516 // Try mesh variables first, then regular ones.
517 VarVector::const_iterator viIt;
518 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
519 if (meshName == viIt->dataBlocks[0]->getMeshName())
520 break;
521 }
522
523 if (viIt == meshVariables.end()) {
524 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
525 if (meshName == viIt->dataBlocks[0]->getMeshName())
526 break;
527 }
528 }
529 // this probably means that the mesh is empty
530 if (viIt == variables.end()) {
531 return;
532 }
533
534 for (size_t idx = 0; idx < viIt->sampleDistribution.size(); idx++) {
535 if (viIt->sampleDistribution[idx] > 0) {
536 stringstream siloPath;
537 siloPath << pathPrefix << "/block";
538 int prevWidth = siloPath.width(4);
539 char prevFill = siloPath.fill('0');
540 siloPath << right << idx;
541 siloPath.width(prevWidth);
542 siloPath.fill(prevFill);
543 siloPath << "/";
544 siloPath << meshName;
545 tempstrings.push_back(siloPath.str());
546 meshnames.push_back((char*)tempstrings.back().c_str());
547 meshtypes.push_back(DB_UCDMESH);
548 }
549 }
550
551 // ignore empty mesh
552 if (meshnames.size() > 0) {
553 DBSetDir(dbfile, "/");
554 DBoptlist* optList = DBMakeOptlist(2);
555 DBAddOption(optList, DBOPT_CYCLE, &cycle);
556 DBAddOption(optList, DBOPT_DTIME, &time);
557 DBPutMultimesh(dbfile, meshName.c_str(), meshnames.size(),
558 &meshnames[0], &meshtypes[0], optList);
559 DBFreeOptlist(optList);
560 }
561 #endif
562 }
563
564 //
565 //
566 //
567 void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
568 bool useMeshFile)
569 {
570 #if USE_SILO
571 vector<int> vartypes;
572 vector<string> tempstrings;
573 vector<char*> varnames;
574 string pathPrefix;
575 if (useMeshFile) {
576 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
577 if (ppIndex != string::npos) {
578 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
579 }
580 }
581
582 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
583 if (vi.sampleDistribution[idx] > 0) {
584 stringstream siloPath;
585 siloPath << pathPrefix << "/block";
586 int prevWidth = siloPath.width(4);
587 char prevFill = siloPath.fill('0');
588 siloPath << right << idx;
589 siloPath.width(prevWidth);
590 siloPath.fill(prevFill);
591 siloPath << "/";
592 siloPath << vi.varName;
593 tempstrings.push_back(siloPath.str());
594 varnames.push_back((char*)tempstrings.back().c_str());
595 vartypes.push_back(DB_UCDVAR);
596 }
597 }
598
599 // ignore empty variables
600 if (varnames.size() > 0) {
601 DBSetDir(dbfile, "/");
602 DBoptlist* optList = DBMakeOptlist(2);
603 DBAddOption(optList, DBOPT_CYCLE, &cycle);
604 DBAddOption(optList, DBOPT_DTIME, &time);
605 if (useMeshFile) {
606 string vpath = string(MESH_VARS)+vi.varName;
607 DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
608 &varnames[0], &vartypes[0], optList);
609 } else {
610 DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
611 &varnames[0], &vartypes[0], optList);
612 }
613 DBFreeOptlist(optList);
614 }
615 #endif
616 }
617
618 //
619 //
620 //
621 void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
622 {
623 #if USE_SILO
624 string tensorDir = vi.varName+string("_comps/");
625 DBSetDir(dbfile, "/");
626 DBMkdir(dbfile, tensorDir.c_str());
627 int one = 1;
628 DBoptlist* optList = DBMakeOptlist(3);
629 DBAddOption(optList, DBOPT_CYCLE, &cycle);
630 DBAddOption(optList, DBOPT_DTIME, &time);
631 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
632 const IntVec& shape = vi.dataBlocks[0]->getShape();
633
634 for (int i=0; i<shape[1]; i++) {
635 for (int j=0; j<shape[0]; j++) {
636 vector<string> tempstrings;
637 vector<char*> varnames;
638 vector<int> vartypes;
639 stringstream comp;
640 comp << vi.varName << "_comps/a_";
641 comp << i;
642 comp << j;
643 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
644 if (vi.sampleDistribution[idx] > 0) {
645 stringstream siloPath;
646 siloPath << "/block";
647 int prevWidth = siloPath.width(4);
648 char prevFill = siloPath.fill('0');
649 siloPath << right << idx;
650 siloPath.width(prevWidth);
651 siloPath.fill(prevFill);
652 siloPath << "/" << comp.str();
653 tempstrings.push_back(siloPath.str());
654 varnames.push_back((char*)tempstrings.back().c_str());
655 vartypes.push_back(DB_UCDVAR);
656 }
657 }
658 if (varnames.size() > 0) {
659 DBPutMultivar(dbfile, comp.str().c_str(), varnames.size(),
660 &varnames[0], &vartypes[0], optList);
661 }
662 }
663 }
664 DBFreeOptlist(optList);
665 #endif
666 }
667
668 } // namespace escriptexport
669

  ViewVC Help
Powered by ViewVC 1.1.26