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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2810 - (show annotations)
Mon Dec 7 04:13:49 2009 UTC (10 years, 11 months ago) by caltinay
File size: 17548 byte(s)
Reincarnation of the escriptreader as a more flexible escriptexport library:
 - can be initialised with instances of escript::Data and Finley_Mesh
 - is now MPI aware at the EscriptDataset level including Silo writer
 - now uses boost shared pointers
The lib is currently only used by the escriptconvert tool but this is going to
change...

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

  ViewVC Help
Powered by ViewVC 1.1.26