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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2910 - (show annotations)
Wed Feb 3 03:22:31 2010 UTC (9 years, 6 months ago) by caltinay
File size: 28813 byte(s)
Added a libescriptreader.a target which does not depend on escript or finley
at runtime. This will be used for the VisIt plugin so is not built by default.

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
20 #ifndef VISIT_PLUGIN
21 #include <escript/Data.h>
22 #endif
23
24 #include <fstream>
25 #include <iostream>
26 #include <numeric> // for std::accumulate
27 #include <sstream>
28
29 #if USE_SILO
30 #include <silo.h>
31
32 #if HAVE_MPI
33 #include <pmpio.h>
34 #endif
35
36 #endif // USE_SILO
37
38 using namespace std;
39
40 namespace escriptexport {
41
42 const char* MESH_VARS = "mesh_vars/";
43 const int NUM_SILO_FILES = 1;
44
45 //
46 // Default constructor
47 //
48 EscriptDataset::EscriptDataset() :
49 numParts(0),
50 cycle(0),
51 time(0.),
52 keepMesh(false),
53 externalMesh(false),
54 mpiRank(0),
55 mpiSize(1)
56 {
57 }
58
59 //
60 // Constructor with communicator
61 //
62 #if HAVE_MPI
63 EscriptDataset::EscriptDataset(MPI_Comm comm) :
64 numParts(0),
65 cycle(0),
66 time(0.),
67 keepMesh(false),
68 externalMesh(false),
69 mpiComm(comm)
70 {
71 MPI_Comm_rank(mpiComm, &mpiRank);
72 MPI_Comm_size(mpiComm, &mpiSize);
73 }
74 #endif
75
76 //
77 // Destructor
78 //
79 EscriptDataset::~EscriptDataset()
80 {
81 }
82
83 //
84 //
85 //
86 bool EscriptDataset::initFromEscript(
87 const escript::AbstractDomain* escriptDomain,
88 DataVec& escriptVars,
89 const StringVec& varNames)
90 {
91 #ifndef VISIT_PLUGIN
92 int myError;
93 numParts = 1;
94 FinleyMesh_ptr mesh(new FinleyMesh());
95 if (mesh->initFromEscript(escriptDomain)) {
96 if (mpiSize > 1)
97 mesh->reorderGhostZones(mpiRank);
98 meshBlocks.push_back(mesh);
99 myError = false;
100 } else {
101 mesh.reset();
102 myError = true;
103 }
104
105 int gError;
106 #if HAVE_MPI
107 MPI_Allreduce(&myError, &gError, 1, MPI_LOGICAL, MPI_LOR, mpiComm);
108 #else
109 gError = myError;
110 #endif
111
112 if (!gError) {
113 // initialize variables
114 DataVec::iterator varIt = escriptVars.begin();
115 StringVec::const_iterator nameIt = varNames.begin();
116 for (; varIt != escriptVars.end(); varIt++, nameIt++) {
117 VarInfo vi;
118 vi.varName = *nameIt;
119
120 DataVar_ptr var(new DataVar(vi.varName));
121 if (var->initFromEscript(*varIt, mesh)) {
122 vi.dataBlocks.push_back(var);
123 updateSampleDistribution(vi);
124 vi.valid = true;
125 } else {
126 var.reset();
127 vi.valid = false;
128 }
129 variables.push_back(vi);
130 }
131
132 // Convert mesh data to variables
133 convertMeshVariables();
134 }
135
136 return !gError;
137
138 #else // VISIT_PLUGIN
139 return false;
140 #endif
141 }
142
143 //
144 //
145 //
146 bool EscriptDataset::loadNetCDF(const string meshFilePattern,
147 const StringVec& varFiles,
148 const StringVec& varNames, int nBlocks)
149 {
150 if (mpiSize > 1 && nBlocks != mpiSize) {
151 cerr << "Cannot load " << nBlocks << " chunks on " << mpiSize
152 << " MPI ranks!" << endl;
153 return false;
154 }
155
156 numParts = nBlocks;
157 meshFmt = meshFilePattern;
158
159 // Load the mesh files
160 if (!loadMeshFromNetCDF()) {
161 cerr << "Reading the domain failed." << endl;
162 return false;
163 }
164
165 // Convert mesh data to variables
166 convertMeshVariables();
167
168 // Load the variables
169 StringVec::const_iterator fileIt = varFiles.begin();
170 StringVec::const_iterator nameIt = varNames.begin();
171 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
172 loadVarFromNetCDF(*fileIt, *nameIt);
173 }
174
175 return true;
176 }
177
178 //
179 // Load only variables using provided mesh
180 //
181 bool EscriptDataset::loadNetCDF(const MeshBlocks& mesh,
182 const StringVec& varFiles,
183 const StringVec& varNames)
184 {
185 if (mpiSize > 1 && mesh.size() > 1) {
186 cerr << "Can only read one mesh block per rank when using MPI!" << endl;
187 return false;
188 }
189
190 externalMesh = true;
191 meshBlocks = mesh;
192 numParts = meshBlocks.size();
193 keepMesh = true;
194
195 // Load the variables
196 StringVec::const_iterator fileIt = varFiles.begin();
197 StringVec::const_iterator nameIt = varNames.begin();
198 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
199 loadVarFromNetCDF(*fileIt, *nameIt);
200 }
201
202 return true;
203 }
204
205 //
206 //
207 //
208 bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)
209 {
210 #if USE_SILO
211 if (numParts == 0)
212 return false;
213
214 const char* blockDirFmt = "/block%04d";
215 string siloPath;
216 DBfile* dbfile = NULL;
217 #if HAVE_MPI
218 PMPIO_baton_t* baton = NULL;
219 #endif
220
221 if (mpiSize > 1) {
222 #if HAVE_MPI
223 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
224 mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
225 PMPIO_DefaultClose, NULL);
226 if (baton) {
227 char str[64];
228 snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, mpiRank));
229 siloPath = str;
230 dbfile = (DBfile*) PMPIO_WaitForBaton(
231 baton, fileName.c_str(), siloPath.c_str());
232 }
233 #endif
234 } else {
235 dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
236 "escriptData", DB_PDB);
237 }
238
239 if (!dbfile) {
240 cerr << "Could not create Silo file." << endl;
241 return false;
242 }
243
244 MeshBlocks::iterator meshIt;
245 VarVector::iterator viIt;
246 int idx = 0;
247 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
248 if (mpiSize == 1) {
249 char str[64];
250 snprintf(str, 64, blockDirFmt, idx);
251 siloPath = str;
252 DBMkdir(dbfile, siloPath.c_str());
253 }
254 // write block of the mesh if we don't use an external mesh
255 if (!externalMesh) {
256 if (! (*meshIt)->writeToSilo(dbfile, siloPath)) {
257 cerr << "Error writing block " << idx << " of mesh to Silo file!\n";
258 break;
259 }
260 }
261
262 // write variables for current mesh block
263 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
264 // do not attempt to write this variable if previous steps failed
265 if (!(*viIt).valid) continue;
266 DataVar_ptr var = (*viIt).dataBlocks[idx];
267 if (!var->writeToSilo(dbfile, siloPath)) {
268 cerr << "Error writing block " << idx << " of '"
269 << var->getName() << "' to Silo file!" << endl;
270 (*viIt).valid = false;
271 }
272 }
273 }
274
275 // rank 0 writes additional data that describe how the parts fit together
276 if (mpiRank == 0) {
277 if (useMultiMesh) {
278 const StringVec& meshNames = meshBlocks[0]->getMeshNames();
279 StringVec::const_iterator it;
280 for (it = meshNames.begin(); it != meshNames.end(); it++)
281 putSiloMultiMesh(dbfile, *it);
282
283 DBMkdir(dbfile, MESH_VARS);
284 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
285 putSiloMultiVar(dbfile, *viIt, true);
286
287 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
288 if (!(*viIt).valid) continue;
289 DataVar_ptr var = (*viIt).dataBlocks[0];
290 if (var->getRank() < 2)
291 putSiloMultiVar(dbfile, *viIt);
292 else
293 putSiloMultiTensor(dbfile, *viIt);
294 }
295 }
296
297 vector<char*> tensorNames;
298 vector<string> tensorDefStrings;
299 vector<char*> tensorDefs;
300
301 // collect tensors for their Silo definitions
302 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
303 if (!(*viIt).valid) continue;
304 DataVar_ptr var = (*viIt).dataBlocks[0];
305 if (var->getRank() == 2) {
306 tensorDefStrings.push_back(var->getTensorDef());
307 tensorDefs.push_back((char*)tensorDefStrings.back().c_str());
308 tensorNames.push_back((char*)var->getName().c_str());
309 }
310 }
311
312 if (tensorDefs.size()) {
313 DBSetDir(dbfile, "/");
314 DBoptlist* optList = DBMakeOptlist(2);
315 DBAddOption(optList, DBOPT_CYCLE, &cycle);
316 DBAddOption(optList, DBOPT_DTIME, &time);
317 vector<DBoptlist*> defOpts(tensorDefs.size(), optList);
318 vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
319 DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
320 &defTypes[0], &tensorDefs[0], &defOpts[0]);
321 DBFreeOptlist(optList);
322 }
323 }
324
325 if (mpiSize > 1) {
326 #if HAVE_MPI
327 PMPIO_HandOffBaton(baton, dbfile);
328 PMPIO_Finish(baton);
329 #endif
330 } else {
331 DBClose(dbfile);
332 }
333
334 return true;
335
336 #else // !USE_SILO
337 return false;
338 #endif
339 }
340
341 //
342 //
343 //
344 bool EscriptDataset::saveVTK(string fileName)
345 {
346 if (numParts == 0)
347 return false;
348
349 string meshName;
350
351 // determine mesh type and variables to write
352 VarVector nodalVars, cellVars;
353 VarVector::iterator viIt;
354 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
355 const DataBlocks& varBlocks = viIt->dataBlocks;
356 // skip empty variable
357 int numSamples = accumulate(viIt->sampleDistribution.begin(),
358 viIt->sampleDistribution.end(), 0);
359 if (numSamples == 0 || !viIt->valid) {
360 continue;
361 }
362 string tmpName = varBlocks[0]->getMeshName();
363 if (meshName != "") {
364 if (meshName != tmpName) {
365 cerr << "VTK supports only one mesh! Skipping variable "
366 << varBlocks[0]->getName() << " on " << tmpName << endl;
367 continue;
368 }
369 } else {
370 meshName = tmpName;
371 }
372
373 if (varBlocks[0]->isNodeCentered()) {
374 nodalVars.push_back(*viIt);
375 } else {
376 cellVars.push_back(*viIt);
377 }
378 }
379
380 // no valid variables so use default mesh
381 if (meshName == "")
382 meshName = "Elements";
383
384 // add mesh variables
385 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
386 DataVar_ptr var = viIt->dataBlocks[0];
387 if (meshName == var->getMeshName()) {
388 VarInfo vi = *viIt;
389 vi.varName = string(MESH_VARS)+vi.varName;
390 if (var->isNodeCentered()) {
391 nodalVars.push_back(vi);
392 } else {
393 cellVars.push_back(vi);
394 }
395 }
396 }
397
398 MeshBlocks::iterator meshIt;
399 int gNumPoints;
400 int gNumCells = 0;
401 int gCellSizeAndType[2] = { 0, 0 };
402
403 #if HAVE_MPI
404 // remove file first if it exists
405 int error = 0;
406 int mpiErr;
407 if (mpiRank == 0) {
408 ifstream f(fileName.c_str());
409 if (f.is_open()) {
410 f.close();
411 if (remove(fileName.c_str())) {
412 error=1;
413 }
414 }
415 }
416 MPI_Allreduce(&error, &mpiErr, 1, MPI_INT, MPI_MAX, mpiComm);
417 if (mpiErr != 0) {
418 cerr << "Error removing " << fileName << "!" << endl;
419 return false;
420 }
421
422 MPI_File mpiFileHandle;
423 MPI_Info mpiInfo = MPI_INFO_NULL;
424 int amode = MPI_MODE_CREATE|MPI_MODE_WRONLY|MPI_MODE_UNIQUE_OPEN;
425 mpiErr = MPI_File_open(mpiComm, const_cast<char*>(fileName.c_str()),
426 amode, mpiInfo, &mpiFileHandle);
427 if (mpiErr == MPI_SUCCESS) {
428 mpiErr = MPI_File_set_view(mpiFileHandle, 0, MPI_CHAR, MPI_CHAR,
429 const_cast<char*>("native"), mpiInfo);
430 }
431 if (mpiErr != MPI_SUCCESS) {
432 cerr << "Error opening " << fileName << " for parallel writing!" << endl;
433 return false;
434 }
435
436 int myNumCells;
437 if (mpiSize > 1)
438 meshBlocks[0]->removeGhostZones(mpiRank);
439 ElementData_ptr elements = meshBlocks[0]->getElementsByName(meshName);
440 myNumCells = elements->getNumElements();
441 MPI_Reduce(&myNumCells, &gNumCells, 1, MPI_INT, MPI_SUM, 0, mpiComm);
442
443 int myCellSizeAndType[2];
444 myCellSizeAndType[0] = elements->getNodesPerElement();
445 myCellSizeAndType[1] = elements->getType();
446
447 // rank 0 needs to know element type and size but it's possible that this
448 // information is only available on other ranks (value=0) so retrieve it
449 MPI_Reduce(&myCellSizeAndType, &gCellSizeAndType, 2, MPI_INT, MPI_MAX, 0,
450 mpiComm);
451
452 // these definitions make the code more readable and life easier -
453 // WRITE_ORDERED causes all ranks to write the contents of the stream while
454 // WRITE_SHARED should only be called for rank 0.
455 #define WRITE_ORDERED(_stream) \
456 do {\
457 MPI_Status mpiStatus;\
458 string contents = _stream.str();\
459 mpiErr = MPI_File_write_ordered(\
460 mpiFileHandle, const_cast<char*>(contents.c_str()),\
461 contents.length(), MPI_CHAR, &mpiStatus);\
462 _stream.str("");\
463 } while(0)
464
465 #define WRITE_SHARED(_stream) \
466 do {\
467 MPI_Status mpiStatus;\
468 string contents = _stream.str();\
469 mpiErr = MPI_File_write_shared(\
470 mpiFileHandle, const_cast<char*>(contents.c_str()),\
471 contents.length(), MPI_CHAR, &mpiStatus);\
472 _stream.str("");\
473 } while(0)
474
475 #else /*********** non-MPI version follows ************/
476
477 ofstream ofs(fileName.c_str());
478
479 int idx = 0;
480 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
481 if (numParts > 1)
482 (*meshIt)->removeGhostZones(idx);
483 ElementData_ptr elements = (*meshIt)->getElementsByName(meshName);
484 gNumCells += elements->getNumElements();
485 if (gCellSizeAndType[0] == 0)
486 gCellSizeAndType[0] = elements->getNodesPerElement();
487 if (gCellSizeAndType[1] == 0)
488 gCellSizeAndType[1] = elements->getType();
489 }
490
491 // the non-MPI versions of WRITE_ORDERED and WRITE_SHARED are
492 // identical.
493 #define WRITE_ORDERED(_stream) \
494 do {\
495 ofs << _stream.str();\
496 _stream.str("");\
497 } while(0)
498
499 #define WRITE_SHARED(_stream) \
500 do {\
501 ofs << _stream.str();\
502 _stream.str("");\
503 } while(0)
504
505 #endif
506
507 gNumPoints = meshBlocks[0]->getNodes()->getGlobalNumNodes();
508
509 ostringstream oss;
510 oss.setf(ios_base::scientific, ios_base::floatfield);
511
512 if (mpiRank == 0) {
513 #ifdef _DEBUG
514 cout << meshName << ": pts=" << gNumPoints << ", cells=" << gNumCells
515 << ", cellsize=" << gCellSizeAndType[0] << ", type="
516 << gCellSizeAndType[1] << ", numNodalVars=" << nodalVars.size()
517 << ", numCellVars=" << cellVars.size()
518 << endl;
519 #endif
520
521 oss << "<?xml version=\"1.0\"?>" << endl;
522 oss << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">" << endl;
523 oss << "<UnstructuredGrid>" << endl;
524
525 // write time and cycle values
526 oss << "<FieldData>" << endl;
527 oss << "<DataArray Name=\"TIME\" type=\"Float64\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
528 oss << time << endl;
529 oss << "</DataArray>" << endl << "</FieldData>" << endl;
530 oss << "<FieldData>" << endl;
531 oss << "<DataArray Name=\"CYCLE\" type=\"Int32\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
532 oss << cycle << endl;
533 oss << "</DataArray>" << endl << "</FieldData>" << endl;
534
535 // write mesh header
536 oss << "<Piece NumberOfPoints=\"" << gNumPoints
537 << "\" NumberOfCells=\"" << gNumCells << "\">" << endl;
538 oss << "<Points>" << endl;
539 oss << "<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">" << endl;
540 }
541
542 //
543 // coordinates (note that we are using the original nodes)
544 //
545 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {
546 (*meshIt)->getNodes()->writeCoordinatesVTK(oss, mpiRank);
547 }
548
549 WRITE_ORDERED(oss);
550
551 if (mpiRank == 0) {
552 oss << "</DataArray>" << endl << "</Points>" << endl;
553 oss << "<Cells>" << endl;
554 oss << "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">" << endl;
555 }
556
557 //
558 // connectivity
559 //
560 for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {
561 ElementData_ptr el = (*meshIt)->getElementsByName(meshName);
562 el->writeConnectivityVTK(oss);
563 }
564
565 WRITE_ORDERED(oss);
566
567 //
568 // offsets & types
569 //
570 if (mpiRank == 0) {
571 oss << "</DataArray>" << endl;
572 oss << "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">" << endl;
573 for (int i=1; i < gNumCells+1; i++) {
574 oss << i*gCellSizeAndType[0] << endl;
575 }
576 oss << "</DataArray>" << endl;
577 oss << "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">" << endl;
578 for (int i=1; i < gNumCells+1; i++) {
579 oss << gCellSizeAndType[1] << endl;
580 }
581 oss << "</DataArray>" << endl << "</Cells>" << endl;
582 WRITE_SHARED(oss);
583 }
584
585 // now write all variables - first the nodal data, then cell data
586
587 // write nodal data if any
588 if (nodalVars.size() > 0) {
589 if (mpiRank == 0)
590 oss << "<PointData>" << endl;
591 for (viIt = nodalVars.begin(); viIt != nodalVars.end(); viIt++) {
592 writeVarToVTK(*viIt, oss);
593 WRITE_ORDERED(oss);
594 if (mpiRank == 0)
595 oss << "</DataArray>" << endl;
596 }
597 if (mpiRank == 0)
598 oss << "</PointData>" << endl;
599 }
600
601 // write cell data if any
602 if (cellVars.size() > 0) {
603 if (mpiRank == 0)
604 oss << "<CellData>" << endl;
605 for (viIt = cellVars.begin(); viIt != cellVars.end(); viIt++) {
606 writeVarToVTK(*viIt, oss);
607 WRITE_ORDERED(oss);
608 if (mpiRank == 0)
609 oss << "</DataArray>" << endl;
610 }
611 if (mpiRank == 0)
612 oss << "</CellData>" << endl;
613 }
614
615 if (mpiRank == 0) {
616 oss << "</Piece>" << endl << "</UnstructuredGrid>" << endl
617 << "</VTKFile>" << endl;
618 WRITE_SHARED(oss);
619 }
620
621 #if HAVE_MPI
622 mpiErr = MPI_File_close(&mpiFileHandle);
623 #else
624 ofs.close();
625 #endif
626
627 return true;
628 }
629
630 //
631 //
632 //
633 void EscriptDataset::writeVarToVTK(const VarInfo& varInfo, ostream& os)
634 {
635 const DataBlocks& varBlocks = varInfo.dataBlocks;
636 int rank = varBlocks[0]->getRank();
637 int numComps = 1;
638 if (rank > 0)
639 numComps *= 3;
640 if (rank > 1)
641 numComps *= 3;
642
643 if (mpiRank == 0) {
644 os << "<DataArray Name=\"" << varInfo.varName
645 << "\" type=\"Float64\" NumberOfComponents=\"" << numComps
646 << "\" format=\"ascii\">" << endl;
647 }
648
649 DataBlocks::const_iterator blockIt;
650 for (blockIt = varBlocks.begin(); blockIt != varBlocks.end(); blockIt++) {
651 (*blockIt)->writeToVTK(os, mpiRank);
652 }
653 }
654
655 //
656 //
657 //
658 bool EscriptDataset::loadMeshFromNetCDF()
659 {
660 bool ok = true;
661 char* str = new char[meshFmt.length()+10];
662 if (mpiSize > 1) {
663 FinleyMesh_ptr meshPart(new FinleyMesh());
664 sprintf(str, meshFmt.c_str(), mpiRank);
665 string meshfile = str;
666 if (meshPart->initFromNetCDF(meshfile)) {
667 meshPart->reorderGhostZones(mpiRank);
668 meshBlocks.push_back(meshPart);
669 } else {
670 meshPart.reset();
671 ok = false;
672 }
673 } else {
674 for (int idx=0; idx < numParts; idx++) {
675 FinleyMesh_ptr meshPart(new FinleyMesh());
676 sprintf(str, meshFmt.c_str(), idx);
677 string meshfile = str;
678 if (meshPart->initFromNetCDF(meshfile)) {
679 if (numParts > 1)
680 meshPart->reorderGhostZones(idx);
681 meshBlocks.push_back(meshPart);
682 } else {
683 meshPart.reset();
684 ok = false;
685 break;
686 }
687 }
688 }
689 delete[] str;
690 return ok;
691 }
692
693 //
694 //
695 //
696 bool EscriptDataset::loadVarFromNetCDF(const string& fileName,
697 const string& varName)
698 {
699 int myError = false;
700 char* str = new char[fileName.length()+10];
701 VarInfo vi;
702
703 vi.varName = varName;
704 vi.valid = true;
705
706 // read all parts of the variable
707 MeshBlocks::iterator mIt;
708 int idx = (mpiSize > 1) ? mpiRank : 0;
709 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++, idx++) {
710 sprintf(str, fileName.c_str(), idx);
711 string dfile = str;
712 DataVar_ptr var(new DataVar(varName));
713 if (var->initFromNetCDF(dfile, *mIt))
714 vi.dataBlocks.push_back(var);
715 else {
716 cerr << "Error reading " << dfile << endl;
717 var.reset();
718 myError = true;
719 break;
720 }
721 }
722 delete[] str;
723
724 int gError;
725 if (mpiSize > 1) {
726 #if HAVE_MPI
727 MPI_Allreduce(&myError, &gError, 1, MPI_LOGICAL, MPI_LOR, mpiComm);
728 #endif
729 } else {
730 gError = myError;
731 }
732
733 if (gError) {
734 // at least one chunk was not read correctly
735 vi.dataBlocks.clear();
736 vi.valid = false;
737 }
738
739 updateSampleDistribution(vi);
740 variables.push_back(vi);
741 return vi.valid;
742 }
743
744 //
745 //
746 //
747 void EscriptDataset::convertMeshVariables()
748 {
749 const StringVec& varNames = meshBlocks[0]->getVarNames();
750 StringVec::const_iterator it;
751 for (it = varNames.begin(); it != varNames.end(); it++) {
752 VarInfo vi;
753 vi.varName = *it;
754 vi.valid = true;
755 // get all parts of current variable
756 MeshBlocks::iterator mIt;
757 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++) {
758 DataVar_ptr var(new DataVar(*it));
759 if (var->initFromMesh(*mIt)) {
760 vi.dataBlocks.push_back(var);
761 } else {
762 cerr << "Error converting mesh variable " << *it << endl;
763 vi.valid = false;
764 break;
765 }
766 }
767 updateSampleDistribution(vi);
768 meshVariables.push_back(vi);
769 }
770 }
771
772 // retrieves the number of samples at each block - used to determine which
773 // blocks contribute to given variable if any.
774 void EscriptDataset::updateSampleDistribution(VarInfo& vi)
775 {
776 IntVec sampleDist;
777 const DataBlocks& varBlocks = vi.dataBlocks;
778
779 if (mpiSize > 1) {
780 #if HAVE_MPI
781 int myNumSamples = varBlocks[0]->getNumberOfSamples();
782 sampleDist.insert(sampleDist.end(), mpiSize, 0);
783 MPI_Allgather(
784 &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, mpiComm);
785 #endif
786 } else {
787 DataBlocks::const_iterator it;
788 for (it = varBlocks.begin(); it != varBlocks.end(); it++) {
789 sampleDist.push_back((*it)->getNumberOfSamples());
790 }
791 }
792 vi.sampleDistribution = sampleDist;
793 }
794
795 //
796 //
797 //
798 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
799 {
800 #if USE_SILO
801 vector<int> meshtypes;
802 vector<string> tempstrings;
803 vector<char*> meshnames;
804 string pathPrefix;
805
806 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
807 if (ppIndex != string::npos) {
808 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
809 }
810
811 // find a variable belonging to this mesh to get the sample
812 // distribution (which tells us which ranks contribute to this mesh).
813 // Try mesh variables first, then regular ones.
814 VarVector::const_iterator viIt;
815 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
816 if (meshName == viIt->dataBlocks[0]->getMeshName())
817 break;
818 }
819
820 if (viIt == meshVariables.end()) {
821 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
822 if (meshName == viIt->dataBlocks[0]->getMeshName())
823 break;
824 }
825 }
826 // this probably means that the mesh is empty
827 if (viIt == variables.end()) {
828 return;
829 }
830
831 for (size_t idx = 0; idx < viIt->sampleDistribution.size(); idx++) {
832 if (viIt->sampleDistribution[idx] > 0) {
833 stringstream siloPath;
834 siloPath << pathPrefix << "/block";
835 int prevWidth = siloPath.width(4);
836 char prevFill = siloPath.fill('0');
837 siloPath << right << idx;
838 siloPath.width(prevWidth);
839 siloPath.fill(prevFill);
840 siloPath << "/";
841 siloPath << meshName;
842 tempstrings.push_back(siloPath.str());
843 meshnames.push_back((char*)tempstrings.back().c_str());
844 meshtypes.push_back(DB_UCDMESH);
845 }
846 }
847
848 // ignore empty mesh
849 if (meshnames.size() > 0) {
850 DBSetDir(dbfile, "/");
851 DBoptlist* optList = DBMakeOptlist(2);
852 DBAddOption(optList, DBOPT_CYCLE, &cycle);
853 DBAddOption(optList, DBOPT_DTIME, &time);
854 DBPutMultimesh(dbfile, meshName.c_str(), meshnames.size(),
855 &meshnames[0], &meshtypes[0], optList);
856 DBFreeOptlist(optList);
857 }
858 #endif
859 }
860
861 //
862 //
863 //
864 void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
865 bool useMeshFile)
866 {
867 #if USE_SILO
868 vector<int> vartypes;
869 vector<string> tempstrings;
870 vector<char*> varnames;
871 string pathPrefix;
872 if (useMeshFile) {
873 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
874 if (ppIndex != string::npos) {
875 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
876 }
877 }
878
879 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
880 if (vi.sampleDistribution[idx] > 0) {
881 stringstream siloPath;
882 siloPath << pathPrefix << "/block";
883 int prevWidth = siloPath.width(4);
884 char prevFill = siloPath.fill('0');
885 siloPath << right << idx;
886 siloPath.width(prevWidth);
887 siloPath.fill(prevFill);
888 siloPath << "/";
889 siloPath << vi.varName;
890 tempstrings.push_back(siloPath.str());
891 varnames.push_back((char*)tempstrings.back().c_str());
892 vartypes.push_back(DB_UCDVAR);
893 }
894 }
895
896 // ignore empty variables
897 if (varnames.size() > 0) {
898 DBSetDir(dbfile, "/");
899 DBoptlist* optList = DBMakeOptlist(2);
900 DBAddOption(optList, DBOPT_CYCLE, &cycle);
901 DBAddOption(optList, DBOPT_DTIME, &time);
902 if (useMeshFile) {
903 string vpath = string(MESH_VARS)+vi.varName;
904 DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
905 &varnames[0], &vartypes[0], optList);
906 } else {
907 DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
908 &varnames[0], &vartypes[0], optList);
909 }
910 DBFreeOptlist(optList);
911 }
912 #endif
913 }
914
915 //
916 //
917 //
918 void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
919 {
920 #if USE_SILO
921 string tensorDir = vi.varName+string("_comps/");
922 DBSetDir(dbfile, "/");
923 DBMkdir(dbfile, tensorDir.c_str());
924 int one = 1;
925 DBoptlist* optList = DBMakeOptlist(3);
926 DBAddOption(optList, DBOPT_CYCLE, &cycle);
927 DBAddOption(optList, DBOPT_DTIME, &time);
928 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
929 const IntVec& shape = vi.dataBlocks[0]->getShape();
930
931 for (int i=0; i<shape[1]; i++) {
932 for (int j=0; j<shape[0]; j++) {
933 vector<string> tempstrings;
934 vector<char*> varnames;
935 vector<int> vartypes;
936 stringstream comp;
937 comp << vi.varName << "_comps/a_";
938 comp << i;
939 comp << j;
940 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
941 if (vi.sampleDistribution[idx] > 0) {
942 stringstream siloPath;
943 siloPath << "/block";
944 int prevWidth = siloPath.width(4);
945 char prevFill = siloPath.fill('0');
946 siloPath << right << idx;
947 siloPath.width(prevWidth);
948 siloPath.fill(prevFill);
949 siloPath << "/" << comp.str();
950 tempstrings.push_back(siloPath.str());
951 varnames.push_back((char*)tempstrings.back().c_str());
952 vartypes.push_back(DB_UCDVAR);
953 }
954 }
955 if (varnames.size() > 0) {
956 DBPutMultivar(dbfile, comp.str().c_str(), varnames.size(),
957 &varnames[0], &vartypes[0], optList);
958 }
959 }
960 }
961 DBFreeOptlist(optList);
962 #endif
963 }
964
965 } // namespace escriptexport
966

  ViewVC Help
Powered by ViewVC 1.1.26