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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2937 - (show annotations)
Tue Feb 16 00:01:52 2010 UTC (10 years, 8 months ago) by caltinay
File size: 28777 byte(s)
Corrected an MPI datatype

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_INT, 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 = myError;
725
726 if (mpiSize > 1) {
727 #if HAVE_MPI
728 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_LOR, mpiComm);
729 #endif
730 }
731
732 if (gError) {
733 // at least one chunk was not read correctly
734 vi.dataBlocks.clear();
735 vi.valid = false;
736 }
737
738 updateSampleDistribution(vi);
739 variables.push_back(vi);
740 return vi.valid;
741 }
742
743 //
744 //
745 //
746 void EscriptDataset::convertMeshVariables()
747 {
748 const StringVec& varNames = meshBlocks[0]->getVarNames();
749 StringVec::const_iterator it;
750 for (it = varNames.begin(); it != varNames.end(); it++) {
751 VarInfo vi;
752 vi.varName = *it;
753 vi.valid = true;
754 // get all parts of current variable
755 MeshBlocks::iterator mIt;
756 for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++) {
757 DataVar_ptr var(new DataVar(*it));
758 if (var->initFromMesh(*mIt)) {
759 vi.dataBlocks.push_back(var);
760 } else {
761 cerr << "Error converting mesh variable " << *it << endl;
762 vi.valid = false;
763 break;
764 }
765 }
766 updateSampleDistribution(vi);
767 meshVariables.push_back(vi);
768 }
769 }
770
771 // retrieves the number of samples at each block - used to determine which
772 // blocks contribute to given variable if any.
773 void EscriptDataset::updateSampleDistribution(VarInfo& vi)
774 {
775 IntVec sampleDist;
776 const DataBlocks& varBlocks = vi.dataBlocks;
777
778 if (mpiSize > 1) {
779 #if HAVE_MPI
780 int myNumSamples = varBlocks[0]->getNumberOfSamples();
781 sampleDist.insert(sampleDist.end(), mpiSize, 0);
782 MPI_Allgather(
783 &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, mpiComm);
784 #endif
785 } else {
786 DataBlocks::const_iterator it;
787 for (it = varBlocks.begin(); it != varBlocks.end(); it++) {
788 sampleDist.push_back((*it)->getNumberOfSamples());
789 }
790 }
791 vi.sampleDistribution = sampleDist;
792 }
793
794 //
795 //
796 //
797 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
798 {
799 #if USE_SILO
800 vector<int> meshtypes;
801 vector<string> tempstrings;
802 vector<char*> meshnames;
803 string pathPrefix;
804
805 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
806 if (ppIndex != string::npos) {
807 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
808 }
809
810 // find a variable belonging to this mesh to get the sample
811 // distribution (which tells us which ranks contribute to this mesh).
812 // Try mesh variables first, then regular ones.
813 VarVector::const_iterator viIt;
814 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
815 if (meshName == viIt->dataBlocks[0]->getMeshName())
816 break;
817 }
818
819 if (viIt == meshVariables.end()) {
820 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
821 if (meshName == viIt->dataBlocks[0]->getMeshName())
822 break;
823 }
824 }
825 // this probably means that the mesh is empty
826 if (viIt == variables.end()) {
827 return;
828 }
829
830 for (size_t idx = 0; idx < viIt->sampleDistribution.size(); idx++) {
831 if (viIt->sampleDistribution[idx] > 0) {
832 stringstream siloPath;
833 siloPath << pathPrefix << "/block";
834 int prevWidth = siloPath.width(4);
835 char prevFill = siloPath.fill('0');
836 siloPath << right << idx;
837 siloPath.width(prevWidth);
838 siloPath.fill(prevFill);
839 siloPath << "/";
840 siloPath << meshName;
841 tempstrings.push_back(siloPath.str());
842 meshnames.push_back((char*)tempstrings.back().c_str());
843 meshtypes.push_back(DB_UCDMESH);
844 }
845 }
846
847 // ignore empty mesh
848 if (meshnames.size() > 0) {
849 DBSetDir(dbfile, "/");
850 DBoptlist* optList = DBMakeOptlist(2);
851 DBAddOption(optList, DBOPT_CYCLE, &cycle);
852 DBAddOption(optList, DBOPT_DTIME, &time);
853 DBPutMultimesh(dbfile, meshName.c_str(), meshnames.size(),
854 &meshnames[0], &meshtypes[0], optList);
855 DBFreeOptlist(optList);
856 }
857 #endif
858 }
859
860 //
861 //
862 //
863 void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
864 bool useMeshFile)
865 {
866 #if USE_SILO
867 vector<int> vartypes;
868 vector<string> tempstrings;
869 vector<char*> varnames;
870 string pathPrefix;
871 if (useMeshFile) {
872 int ppIndex = meshBlocks[0]->getSiloPath().find(':');
873 if (ppIndex != string::npos) {
874 pathPrefix = meshBlocks[0]->getSiloPath().substr(0, ppIndex+1);
875 }
876 }
877
878 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
879 if (vi.sampleDistribution[idx] > 0) {
880 stringstream siloPath;
881 siloPath << pathPrefix << "/block";
882 int prevWidth = siloPath.width(4);
883 char prevFill = siloPath.fill('0');
884 siloPath << right << idx;
885 siloPath.width(prevWidth);
886 siloPath.fill(prevFill);
887 siloPath << "/";
888 siloPath << vi.varName;
889 tempstrings.push_back(siloPath.str());
890 varnames.push_back((char*)tempstrings.back().c_str());
891 vartypes.push_back(DB_UCDVAR);
892 }
893 }
894
895 // ignore empty variables
896 if (varnames.size() > 0) {
897 DBSetDir(dbfile, "/");
898 DBoptlist* optList = DBMakeOptlist(2);
899 DBAddOption(optList, DBOPT_CYCLE, &cycle);
900 DBAddOption(optList, DBOPT_DTIME, &time);
901 if (useMeshFile) {
902 string vpath = string(MESH_VARS)+vi.varName;
903 DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
904 &varnames[0], &vartypes[0], optList);
905 } else {
906 DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
907 &varnames[0], &vartypes[0], optList);
908 }
909 DBFreeOptlist(optList);
910 }
911 #endif
912 }
913
914 //
915 //
916 //
917 void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
918 {
919 #if USE_SILO
920 string tensorDir = vi.varName+string("_comps/");
921 DBSetDir(dbfile, "/");
922 DBMkdir(dbfile, tensorDir.c_str());
923 int one = 1;
924 DBoptlist* optList = DBMakeOptlist(3);
925 DBAddOption(optList, DBOPT_CYCLE, &cycle);
926 DBAddOption(optList, DBOPT_DTIME, &time);
927 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
928 const IntVec& shape = vi.dataBlocks[0]->getShape();
929
930 for (int i=0; i<shape[1]; i++) {
931 for (int j=0; j<shape[0]; j++) {
932 vector<string> tempstrings;
933 vector<char*> varnames;
934 vector<int> vartypes;
935 stringstream comp;
936 comp << vi.varName << "_comps/a_";
937 comp << i;
938 comp << j;
939 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
940 if (vi.sampleDistribution[idx] > 0) {
941 stringstream siloPath;
942 siloPath << "/block";
943 int prevWidth = siloPath.width(4);
944 char prevFill = siloPath.fill('0');
945 siloPath << right << idx;
946 siloPath.width(prevWidth);
947 siloPath.fill(prevFill);
948 siloPath << "/" << comp.str();
949 tempstrings.push_back(siloPath.str());
950 varnames.push_back((char*)tempstrings.back().c_str());
951 vartypes.push_back(DB_UCDVAR);
952 }
953 }
954 if (varnames.size() > 0) {
955 DBPutMultivar(dbfile, comp.str().c_str(), varnames.size(),
956 &varnames[0], &vartypes[0], optList);
957 }
958 }
959 }
960 DBFreeOptlist(optList);
961 #endif
962 }
963
964 } // namespace escriptexport
965

  ViewVC Help
Powered by ViewVC 1.1.26