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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3397 - (show annotations)
Sun Dec 5 22:46:25 2010 UTC (8 years, 8 months ago) by caltinay
File size: 34320 byte(s)
Fixed an issue in saveVTK/Silo.

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 <weipa/EscriptDataset.h>
15 #include <weipa/DataVar.h>
16 #include <weipa/ElementData.h>
17 #include <weipa/FileWriter.h>
18 #include <weipa/FinleyDomain.h>
19 #include <weipa/NodeData.h>
20
21 #ifndef VISIT_PLUGIN
22 #include <escript/Data.h>
23 #include <dudley/CppAdapter/MeshAdapter.h>
24 #include <finley/CppAdapter/MeshAdapter.h>
25 #endif
26
27 #include <cstring>
28 #include <iostream>
29 #include <numeric> // for std::accumulate
30 #include <sstream> // for std::ostringstream
31
32 #if USE_SILO
33 #include <silo.h>
34
35 #if HAVE_MPI
36 #include <pmpio.h>
37 #endif
38
39 #endif // USE_SILO
40
41 using namespace std;
42
43 namespace weipa {
44
45 const char* MESH_VARS = "mesh_vars/";
46 const int NUM_SILO_FILES = 1;
47
48 //
49 // Default constructor
50 //
51 EscriptDataset::EscriptDataset() :
52 cycle(0),
53 time(0.),
54 externalDomain(false),
55 mpiRank(0),
56 mpiSize(1)
57 {
58 }
59
60 //
61 // Constructor with communicator
62 //
63 #if HAVE_MPI
64 EscriptDataset::EscriptDataset(MPI_Comm comm) :
65 cycle(0),
66 time(0.),
67 externalDomain(false),
68 mpiComm(comm)
69 {
70 MPI_Comm_rank(mpiComm, &mpiRank);
71 MPI_Comm_size(mpiComm, &mpiSize);
72 }
73 #endif
74
75 //
76 // Destructor
77 //
78 EscriptDataset::~EscriptDataset()
79 {
80 }
81
82 //
83 // Sets the domain using an escript domain instance.
84 //
85 bool EscriptDataset::setDomain(const escript::AbstractDomain* domain)
86 {
87 #ifndef VISIT_PLUGIN
88 int myError = 0, gError;
89
90 // fail if the domain has already been set
91 if (domainChunks.size() > 0) {
92 cerr << "Domain has already been set!" << endl;
93 myError = 1;
94 } else if (!domain) {
95 cerr << "Domain is NULL!" << endl;
96 myError = 1;
97 } else {
98 #if HAVE_MPI
99 mpiComm = domain->getMPIComm();
100 mpiRank = domain->getMPIRank();
101 mpiSize = domain->getMPISize();
102 #endif
103 // FinleyDomain can handle both finley and dudley
104 if (dynamic_cast<const finley::MeshAdapter*>(domain)
105 || dynamic_cast<const dudley::MeshAdapter*>(domain)) {
106 DomainChunk_ptr dom(new FinleyDomain());
107 if (dom->initFromEscript(domain)) {
108 if (mpiSize > 1)
109 dom->reorderGhostZones(mpiRank);
110 domainChunks.push_back(dom);
111 } else {
112 cerr << "Error initializing domain!" << endl;
113 myError = 2;
114 }
115 } else {
116 cerr << "Unsupported domain type!" << endl;
117 myError = 2;
118 }
119 }
120
121 if (mpiSize > 1) {
122 #if HAVE_MPI
123 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
124 #else
125 gError = myError;
126 #endif
127 } else {
128 gError = myError;
129 }
130
131 if (gError>1) {
132 domainChunks.clear();
133 } else if (gError==0) {
134 // Convert mesh data to variables
135 convertMeshVariables();
136 }
137 return (gError==0);
138
139 #else // VISIT_PLUGIN
140 return false;
141 #endif
142 }
143
144 //
145 //
146 //
147 bool EscriptDataset::addData(escript::Data& data, const string name,
148 const string units)
149 {
150 #ifndef VISIT_PLUGIN
151 bool success = true;
152
153 // fail if no domain has been set
154 if (domainChunks.size() == 0) {
155 success = false;
156 } else {
157 // initialize variable
158 VarInfo vi;
159 vi.varName = name;
160 vi.units = units;
161
162 DataVar_ptr var(new DataVar(vi.varName));
163 if (var->initFromEscript(data, domainChunks[0])) {
164 vi.dataChunks.push_back(var);
165 updateSampleDistribution(vi);
166 vi.valid = true;
167 } else {
168 var.reset();
169 vi.valid = false;
170 }
171 variables.push_back(vi);
172 }
173 return success;
174
175 #else // VISIT_PLUGIN
176 return false;
177 #endif
178 }
179
180 //
181 //
182 //
183 bool EscriptDataset::loadNetCDF(const string domainFilePattern,
184 const StringVec& varFiles,
185 const StringVec& varNames, int nChunks)
186 {
187 // sanity check
188 if (varFiles.size() != varNames.size()) {
189 return false;
190 }
191
192 // load the domain files
193 if (!loadDomain(domainFilePattern, nChunks)) {
194 return false;
195 }
196
197 // load the variables
198 StringVec::const_iterator fileIt = varFiles.begin();
199 StringVec::const_iterator nameIt = varNames.begin();
200 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
201 loadData(*fileIt, *nameIt, "");
202 }
203
204 return true;
205 }
206
207 //
208 // Load only variables using provided domain
209 //
210 bool EscriptDataset::loadNetCDF(const DomainChunks& domain,
211 const StringVec& varFiles,
212 const StringVec& varNames)
213 {
214 // sanity check
215 if (varFiles.size() != varNames.size()) {
216 return false;
217 }
218
219 // set the domain
220 if (!setExternalDomain(domain)) {
221 return false;
222 }
223
224 // load the variables
225 StringVec::const_iterator fileIt = varFiles.begin();
226 StringVec::const_iterator nameIt = varNames.begin();
227 for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
228 loadData(*fileIt, *nameIt, "");
229 }
230
231 return true;
232 }
233
234 //
235 //
236 //
237 bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)
238 {
239 #if USE_SILO
240 if (domainChunks.size() == 0)
241 return false;
242
243 const char* blockDirFmt = "/block%04d";
244 string siloPath;
245 DBfile* dbfile = NULL;
246 int driver = DB_HDF5; // prefer HDF5 if available
247 //FIXME: Silo's HDF5 driver and NetCDF 4 are currently incompatible because
248 //NetCDF calls H5close() when all its files are closed.
249 //Unidata has been contacted, Ticket ID: YTC-894489.
250 //When this issue is resolved, remove the following line.
251 driver = DB_PDB;
252 #if HAVE_MPI
253 PMPIO_baton_t* baton = NULL;
254 #endif
255
256 if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
257 fileName+=".silo";
258 }
259
260 if (mpiSize > 1) {
261 #if HAVE_MPI
262 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
263 mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
264 PMPIO_DefaultClose, (void*)&driver);
265 // try the fallback driver in case of error
266 if (!baton && driver != DB_PDB) {
267 driver = DB_PDB;
268 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE,
269 mpiComm, 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
270 PMPIO_DefaultClose, (void*)&driver);
271 }
272 if (baton) {
273 char str[64];
274 snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, mpiRank));
275 siloPath = str;
276 dbfile = (DBfile*) PMPIO_WaitForBaton(
277 baton, fileName.c_str(), siloPath.c_str());
278 }
279 #endif
280 } else {
281 dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
282 "escriptData", driver);
283 // try the fallback driver in case of error
284 if (!dbfile && driver != DB_PDB) {
285 driver = DB_PDB;
286 dbfile = DBCreate(fileName.c_str(), DB_CLOBBER, DB_LOCAL,
287 "escriptData", driver);
288 }
289 }
290
291 if (!dbfile) {
292 cerr << "Could not create Silo file." << endl;
293 if (mpiSize > 1) {
294 #if HAVE_MPI
295 PMPIO_HandOffBaton(baton, dbfile);
296 PMPIO_Finish(baton);
297 #endif
298 }
299 return false;
300 }
301
302 if (driver==DB_HDF5) {
303 // gzip level 1 already provides good compression with minimal
304 // performance penalty. Some tests showed that gzip levels >3 performed
305 // rather badly on escript data both in terms of time and space
306 DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
307 }
308
309 DomainChunks::iterator domIt;
310 VarVector::iterator viIt;
311 int idx = 0;
312 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
313 if (mpiSize == 1) {
314 char str[64];
315 snprintf(str, 64, blockDirFmt, idx);
316 siloPath = str;
317 DBMkdir(dbfile, siloPath.c_str());
318 }
319 // write block of the mesh if we don't use an external domain
320 if (!externalDomain) {
321 if (! (*domIt)->writeToSilo(
322 dbfile, siloPath, meshLabels, meshUnits)) {
323 cerr << "Error writing block " << idx
324 << " of mesh to Silo file!" << endl;
325 break;
326 }
327 }
328
329 // write variables for current domain chunk
330 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
331 // do not attempt to write this variable if previous steps failed
332 if (!viIt->valid) continue;
333 DataVar_ptr var = viIt->dataChunks[idx];
334 if (!var->writeToSilo(dbfile, siloPath, viIt->units)) {
335 cerr << "Error writing block " << idx << " of '"
336 << var->getName() << "' to Silo file!" << endl;
337 viIt->valid = false;
338 }
339 }
340 }
341
342 // rank 0 writes additional data that describe how the parts fit together
343 if (mpiRank == 0) {
344 if (useMultiMesh) {
345 const StringVec& meshNames = domainChunks[0]->getMeshNames();
346 StringVec::const_iterator it;
347 for (it = meshNames.begin(); it != meshNames.end(); it++)
348 putSiloMultiMesh(dbfile, *it);
349
350 DBMkdir(dbfile, MESH_VARS);
351 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++)
352 putSiloMultiVar(dbfile, *viIt, true);
353
354 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
355 if (!viIt->valid) continue;
356 DataVar_ptr var = viIt->dataChunks[0];
357 if (var->getRank() < 2)
358 putSiloMultiVar(dbfile, *viIt);
359 else
360 putSiloMultiTensor(dbfile, *viIt);
361 }
362 }
363
364 vector<char*> tensorNames;
365 vector<string> tensorDefStrings;
366 vector<char*> tensorDefs;
367
368 // collect tensors for their Silo definitions
369 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
370 if (!viIt->valid) continue;
371 DataVar_ptr var = viIt->dataChunks[0];
372 if (var->getRank() == 2) {
373 tensorDefStrings.push_back(var->getTensorDef());
374 tensorDefs.push_back((char*)tensorDefStrings.back().c_str());
375 tensorNames.push_back((char*)var->getName().c_str());
376 }
377 }
378
379 if (tensorDefs.size()) {
380 DBSetDir(dbfile, "/");
381 DBoptlist* optList = DBMakeOptlist(2);
382 DBAddOption(optList, DBOPT_CYCLE, &cycle);
383 DBAddOption(optList, DBOPT_DTIME, &time);
384 vector<DBoptlist*> defOpts(tensorDefs.size(), optList);
385 vector<int> defTypes(tensorDefs.size(), DB_VARTYPE_TENSOR);
386 DBPutDefvars(dbfile, "tensors", tensorDefs.size(), &tensorNames[0],
387 &defTypes[0], &tensorDefs[0], &defOpts[0]);
388 DBFreeOptlist(optList);
389 }
390 }
391
392 if (mpiSize > 1) {
393 #if HAVE_MPI
394 PMPIO_HandOffBaton(baton, dbfile);
395 PMPIO_Finish(baton);
396 #endif
397 } else {
398 DBClose(dbfile);
399 }
400
401 return true;
402
403 #else // !USE_SILO
404 return false;
405 #endif
406 }
407
408 //
409 //
410 //
411 bool EscriptDataset::saveVTK(string fileName)
412 {
413 if (domainChunks.size() == 0)
414 return false;
415
416 map<string,VarVector> varsPerMesh;
417
418 // determine meshes and variables to write
419 VarVector::iterator viIt;
420 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
421 // skip empty variable
422 int numSamples = accumulate(viIt->sampleDistribution.begin(),
423 viIt->sampleDistribution.end(), 0);
424 if (numSamples == 0 || !viIt->valid) {
425 continue;
426 }
427 string meshName = viIt->dataChunks[0]->getMeshName();
428
429 // ranks with 0 samples can't know the correct mesh name.
430 // We assume rank 0 always has samples, if this turns out to be a
431 // wrong assumption then a bit more work is needed to get the correct
432 // mesh name to all ranks.
433 #if HAVE_MPI
434 if (mpiSize > 1) {
435 char name[100];
436 if (mpiRank == 0) {
437 strncpy(&name[0], meshName.c_str(), 100);
438 }
439
440 MPI_Bcast(name, 100, MPI_CHAR, 0, mpiComm);
441 meshName = name;
442 }
443 #endif
444
445 map<string,VarVector>::iterator it = varsPerMesh.find(meshName);
446 if (it != varsPerMesh.end()) {
447 it->second.push_back(*viIt);
448 } else {
449 VarVector v;
450 v.push_back(*viIt);
451 varsPerMesh[meshName] = v;
452 }
453 }
454
455 if (fileName.length() < 5 || fileName.compare(fileName.length()-4, 4, ".vtu") != 0) {
456 fileName+=".vtu";
457 }
458
459 bool ret = true;
460 if (varsPerMesh.empty()) {
461 // no valid variables so just write default mesh
462 ret = saveVTKsingle(fileName, "Elements", VarVector());
463 } else {
464 // write one file per required mesh
465 string newName(fileName);
466 string filePrefix(fileName.substr(0, fileName.length()-4));
467 bool prependMeshName=(varsPerMesh.size()>1);
468 map<string,VarVector>::const_iterator vpmIt;
469 for (vpmIt=varsPerMesh.begin(); vpmIt!=varsPerMesh.end(); vpmIt++) {
470 if (prependMeshName) {
471 newName=filePrefix+"_"+vpmIt->first+".vtu";
472 }
473 // attempt to write all files even if one fails
474 if (!saveVTKsingle(newName, vpmIt->first, vpmIt->second)) {
475 ret = false;
476 }
477 }
478 }
479 return ret;
480 }
481
482 //
483 //
484 //
485 void EscriptDataset::setMeshLabels(const string x, const string y, const string z)
486 {
487 meshLabels.clear();
488 meshLabels.push_back(x);
489 meshLabels.push_back(y);
490 if (z.length()>0)
491 meshLabels.push_back(z);
492 }
493
494 //
495 //
496 //
497 void EscriptDataset::setMeshUnits(const string x, const string y, const string z)
498 {
499 meshUnits.clear();
500 meshUnits.push_back(x);
501 meshUnits.push_back(y);
502 if (z.length()>0)
503 meshUnits.push_back(z);
504 }
505
506 //
507 //
508 //
509 bool EscriptDataset::saveVTKsingle(const string& fileName,
510 const string& meshName,
511 const VarVector& vars)
512 {
513 VarVector nodalVars, cellVars;
514 VarVector::const_iterator viIt;
515 for (viIt = vars.begin(); viIt != vars.end(); viIt++) {
516 const DataChunks& varChunks = viIt->dataChunks;
517 if (varChunks[0]->isNodeCentered()) {
518 nodalVars.push_back(*viIt);
519 } else {
520 cellVars.push_back(*viIt);
521 }
522 }
523
524 // add mesh variables
525 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
526 DataVar_ptr var = viIt->dataChunks[0];
527 if (meshName == var->getMeshName()) {
528 VarInfo vi = *viIt;
529 vi.varName = string(MESH_VARS)+vi.varName;
530 if (var->isNodeCentered()) {
531 nodalVars.push_back(vi);
532 } else {
533 cellVars.push_back(vi);
534 }
535 }
536 }
537
538 DomainChunks::iterator domIt;
539 int gNumPoints;
540 int gNumCells = 0;
541 int gCellSizeAndType[2] = { 0, 0 };
542
543 FileWriter* fw = NULL;
544
545 if (mpiSize > 1) {
546 #if HAVE_MPI
547 fw = new FileWriter(mpiComm);
548 domainChunks[0]->removeGhostZones(mpiRank);
549 ElementData_ptr elements = domainChunks[0]->getElementsByName(meshName);
550 int myNumCells = 0;
551 if (elements) {
552 myNumCells = elements->getNumElements();
553 }
554 MPI_Reduce(&myNumCells, &gNumCells, 1, MPI_INT, MPI_SUM, 0, mpiComm);
555
556 int myCellSizeAndType[2] = { 0, 0 };
557 if (elements) {
558 myCellSizeAndType[0] = elements->getNodesPerElement();
559 myCellSizeAndType[1] = elements->getType();
560 }
561
562 // rank 0 needs to know element type and size but it's possible that
563 // this information is only available on other ranks (value=0) so
564 // retrieve it
565 MPI_Reduce(&myCellSizeAndType, &gCellSizeAndType, 2, MPI_INT, MPI_MAX,
566 0, mpiComm);
567 #endif
568 } else {
569 fw = new FileWriter();
570 int idx = 0;
571 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
572 if (domainChunks.size() > 1)
573 (*domIt)->removeGhostZones(idx);
574 ElementData_ptr elements = (*domIt)->getElementsByName(meshName);
575 gNumCells += elements->getNumElements();
576 if (gCellSizeAndType[0] == 0)
577 gCellSizeAndType[0] = elements->getNodesPerElement();
578 if (gCellSizeAndType[1] == 0)
579 gCellSizeAndType[1] = elements->getType();
580 }
581 }
582
583 gNumPoints = domainChunks[0]->getNodes()->getGlobalNumNodes();
584
585 ostringstream oss;
586 oss.setf(ios_base::scientific, ios_base::floatfield);
587
588 if (!fw->openFile(fileName)) {
589 return false;
590 }
591
592 if (mpiRank == 0) {
593 #ifdef _DEBUG
594 cout << meshName << ": pts=" << gNumPoints << ", cells=" << gNumCells
595 << ", cellsize=" << gCellSizeAndType[0] << ", type="
596 << gCellSizeAndType[1] << ", numNodalVars=" << nodalVars.size()
597 << ", numCellVars=" << cellVars.size()
598 << endl;
599 #endif
600
601 oss << "<?xml version=\"1.0\"?>" << endl;
602 oss << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
603 if (mdSchema.length()>0) {
604 oss << " " << mdSchema;
605 }
606 oss << ">" << endl;
607 if (mdString.length()>0) {
608 oss << "<MetaData>" << endl << mdString << endl
609 << "</MetaData>" << endl;
610 }
611 oss << "<UnstructuredGrid>" << endl;
612
613 // write time and cycle values
614 oss << "<FieldData>" << endl;
615 oss << "<DataArray Name=\"TIME\" type=\"Float64\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
616 oss << time << endl;
617 oss << "</DataArray>" << endl << "</FieldData>" << endl;
618 oss << "<FieldData>" << endl;
619 oss << "<DataArray Name=\"CYCLE\" type=\"Int32\" format=\"ascii\" NumberOfTuples=\"1\">" << endl;
620 oss << cycle << endl;
621 oss << "</DataArray>" << endl << "</FieldData>" << endl;
622
623 // write mesh header
624 oss << "<Piece NumberOfPoints=\"" << gNumPoints
625 << "\" NumberOfCells=\"" << gNumCells << "\">" << endl;
626 oss << "<Points>" << endl;
627 oss << "<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">" << endl;
628 }
629
630 //
631 // coordinates (note that we are using the original nodes)
632 //
633 int blockNum = (mpiSize>1 ? mpiRank : 0);
634 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, blockNum++) {
635 (*domIt)->getNodes()->writeCoordinatesVTK(oss, blockNum);
636 }
637
638 if (!fw->writeOrdered(oss)) {
639 cerr << "Warning, ignoring file write error!" << endl;
640 }
641
642 if (mpiRank == 0) {
643 oss << "</DataArray>" << endl << "</Points>" << endl;
644 oss << "<Cells>" << endl;
645 oss << "<DataArray Name=\"connectivity\" type=\"Int32\" format=\"ascii\">" << endl;
646 }
647
648 //
649 // connectivity
650 //
651 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++) {
652 ElementData_ptr el = (*domIt)->getElementsByName(meshName);
653 if (el) {
654 el->writeConnectivityVTK(oss);
655 }
656 }
657
658 if (!fw->writeOrdered(oss)) {
659 cerr << "Warning, ignoring file write error!" << endl;
660 }
661
662 //
663 // offsets & types
664 //
665 if (mpiRank == 0) {
666 oss << "</DataArray>" << endl;
667 oss << "<DataArray Name=\"offsets\" type=\"Int32\" format=\"ascii\">" << endl;
668 for (int i=1; i < gNumCells+1; i++) {
669 oss << i*gCellSizeAndType[0] << endl;
670 }
671 oss << "</DataArray>" << endl;
672 oss << "<DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">" << endl;
673 for (int i=1; i < gNumCells+1; i++) {
674 oss << gCellSizeAndType[1] << endl;
675 }
676 oss << "</DataArray>" << endl << "</Cells>" << endl;
677 if (!fw->writeShared(oss)) {
678 cerr << "Warning, ignoring file write error!" << endl;
679 }
680 }
681
682 // now write all variables - first the nodal data, then cell data
683
684 // write nodal data if any
685 if (!nodalVars.empty()) {
686 if (mpiRank == 0)
687 oss << "<PointData>" << endl;
688 for (viIt = nodalVars.begin(); viIt != nodalVars.end(); viIt++) {
689 writeVarToVTK(*viIt, oss);
690 if (!fw->writeOrdered(oss)) {
691 cerr << "Warning, ignoring file write error!" << endl;
692 }
693 if (mpiRank == 0)
694 oss << "</DataArray>" << endl;
695 }
696 if (mpiRank == 0)
697 oss << "</PointData>" << endl;
698 }
699
700 // write cell data if any
701 if (!cellVars.empty()) {
702 if (mpiRank == 0)
703 oss << "<CellData>" << endl;
704 for (viIt = cellVars.begin(); viIt != cellVars.end(); viIt++) {
705 writeVarToVTK(*viIt, oss);
706 if (!fw->writeOrdered(oss)) {
707 cerr << "Warning, ignoring file write error!" << endl;
708 }
709 if (mpiRank == 0)
710 oss << "</DataArray>" << endl;
711 }
712 if (mpiRank == 0)
713 oss << "</CellData>" << endl;
714 }
715
716 if (mpiRank == 0) {
717 oss << "</Piece>" << endl << "</UnstructuredGrid>" << endl
718 << "</VTKFile>" << endl;
719 if (!fw->writeShared(oss)) {
720 cerr << "Warning, ignoring file write error!" << endl;
721 }
722 }
723
724 fw->close();
725 delete fw;
726 return true;
727 }
728
729 //
730 //
731 //
732 void EscriptDataset::writeVarToVTK(const VarInfo& varInfo, ostream& os)
733 {
734 const DataChunks& varChunks = varInfo.dataChunks;
735 int rank = varChunks[0]->getRank();
736 int numComps = 1;
737 if (rank > 0)
738 numComps *= 3;
739 if (rank > 1)
740 numComps *= 3;
741
742 if (mpiRank == 0) {
743 os << "<DataArray Name=\"" << varInfo.varName
744 << "\" type=\"Float64\" NumberOfComponents=\"" << numComps
745 << "\" format=\"ascii\">" << endl;
746 }
747
748 // this is required in case we read a dataset with more than one chunk on
749 // one rank
750 int blockNum = (mpiSize>1 ? mpiRank : 0);
751 DataChunks::const_iterator blockIt;
752 for (blockIt = varChunks.begin(); blockIt != varChunks.end(); blockIt++, blockNum++) {
753 (*blockIt)->writeToVTK(os, blockNum);
754 }
755 }
756
757 //
758 // Sets the domain from dump files.
759 //
760 bool EscriptDataset::loadDomain(const string filePattern, int nChunks)
761 {
762 int myError = 0, gError;
763
764 if (mpiSize > 1 && nChunks != mpiSize) {
765 cerr << "Cannot load " << nChunks << " chunks on " << mpiSize
766 << " MPI ranks!" << endl;
767 myError = 1;
768
769 } else if (domainChunks.size() > 0) {
770 cerr << "Domain has already been set!" << endl;
771 myError = 1;
772
773 } else {
774 char* str = new char[filePattern.length()+10];
775 // FIXME: This assumes a finley domain!
776 if (mpiSize > 1) {
777 DomainChunk_ptr chunk(new FinleyDomain());
778 sprintf(str, filePattern.c_str(), mpiRank);
779 string domainFile = str;
780 if (chunk->initFromFile(domainFile)) {
781 chunk->reorderGhostZones(mpiRank);
782 domainChunks.push_back(chunk);
783 } else {
784 cerr << "Error initializing domain!" << endl;
785 myError = 2;
786 }
787 } else {
788 for (int idx=0; idx < nChunks; idx++) {
789 DomainChunk_ptr chunk(new FinleyDomain());
790 sprintf(str, filePattern.c_str(), idx);
791 string domainFile = str;
792 if (chunk->initFromFile(domainFile)) {
793 if (nChunks > 1)
794 chunk->reorderGhostZones(idx);
795 domainChunks.push_back(chunk);
796 } else {
797 cerr << "Error initializing domain block " << idx << endl;
798 myError = 2;
799 break;
800 }
801 }
802 }
803 delete[] str;
804 }
805
806 if (mpiSize > 1) {
807 #if HAVE_MPI
808 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
809 #else
810 gError = myError;
811 #endif
812 } else {
813 gError = myError;
814 }
815
816 if (gError>1) {
817 domainChunks.clear();
818 } else if (gError==0) {
819 // Convert mesh data to variables
820 convertMeshVariables();
821 }
822 return (gError==0);
823 }
824
825 //
826 // Sets an already converted domain.
827 //
828 bool EscriptDataset::setExternalDomain(const DomainChunks& domain)
829 {
830 int myError = 0, gError;
831
832 if (mpiSize > 1 && domain.size() > 1) {
833 cerr << "Can only add one domain block per rank when using MPI!"
834 << endl;
835 myError = 1;
836
837 } else if (domainChunks.size() > 0) {
838 cerr << "Domain has already been set!" << endl;
839 myError = 1;
840
841 }
842
843 if (mpiSize > 1) {
844 #if HAVE_MPI
845 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
846 #else
847 gError = myError;
848 #endif
849 } else {
850 gError = myError;
851 }
852
853 if (!gError) {
854 externalDomain = true;
855 domainChunks = domain;
856 }
857
858 return !gError;
859 }
860
861 //
862 //
863 //
864 bool EscriptDataset::loadData(const string filePattern, const string name,
865 const string units)
866 {
867 int myError = 0, gError;
868
869 // fail if no domain has been set
870 if (domainChunks.size() == 0) {
871 gError = 1;
872
873 } else {
874 // initialize variable
875 VarInfo vi;
876 vi.varName = name;
877 vi.units = units;
878 vi.valid = true;
879 char* str = new char[filePattern.length()+10];
880
881 // read all parts of the variable
882 DomainChunks::iterator domIt;
883 int idx = (mpiSize > 1) ? mpiRank : 0;
884 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++, idx++) {
885 sprintf(str, filePattern.c_str(), idx);
886 string dfile = str;
887 DataVar_ptr var(new DataVar(name));
888 if (var->initFromFile(dfile, *domIt))
889 vi.dataChunks.push_back(var);
890 else {
891 cerr << "Error reading " << dfile << endl;
892 myError = 1;
893 break;
894 }
895 }
896 delete[] str;
897
898 if (mpiSize > 1) {
899 #if HAVE_MPI
900 MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
901 #else
902 gError = myError;
903 #endif
904 } else {
905 gError = myError;
906 }
907
908 if (!gError) {
909 // only add variable if all chunks have been read without error
910 updateSampleDistribution(vi);
911 variables.push_back(vi);
912 }
913 }
914
915 return !gError;
916 }
917
918 //
919 //
920 //
921 void EscriptDataset::convertMeshVariables()
922 {
923 const StringVec& varNames = domainChunks[0]->getVarNames();
924 StringVec::const_iterator it;
925 for (it = varNames.begin(); it != varNames.end(); it++) {
926 VarInfo vi;
927 vi.varName = *it;
928 vi.valid = true;
929 // get all parts of current variable
930 DomainChunks::iterator domIt;
931 for (domIt = domainChunks.begin(); domIt != domainChunks.end(); domIt++) {
932 DataVar_ptr var = (*domIt)->getDataVarByName(*it);
933 if (var != NULL) {
934 vi.dataChunks.push_back(var);
935 } else {
936 cerr << "Error converting mesh variable " << *it << endl;
937 vi.valid = false;
938 break;
939 }
940 }
941 updateSampleDistribution(vi);
942 meshVariables.push_back(vi);
943 }
944 }
945
946 // retrieves the number of samples at each block - used to determine which
947 // blocks contribute to given variable if any.
948 void EscriptDataset::updateSampleDistribution(VarInfo& vi)
949 {
950 IntVec sampleDist;
951 const DataChunks& varChunks = vi.dataChunks;
952
953 if (mpiSize > 1) {
954 #if HAVE_MPI
955 int myNumSamples = varChunks[0]->getNumberOfSamples();
956 sampleDist.insert(sampleDist.end(), mpiSize, 0);
957 MPI_Allgather(
958 &myNumSamples, 1, MPI_INT, &sampleDist[0], 1, MPI_INT, mpiComm);
959 #endif
960 } else {
961 DataChunks::const_iterator it;
962 for (it = varChunks.begin(); it != varChunks.end(); it++) {
963 sampleDist.push_back((*it)->getNumberOfSamples());
964 }
965 }
966 vi.sampleDistribution = sampleDist;
967 }
968
969 //
970 //
971 //
972 void EscriptDataset::putSiloMultiMesh(DBfile* dbfile, const string& meshName)
973 {
974 #if USE_SILO
975 vector<int> meshtypes;
976 vector<string> tempstrings;
977 vector<char*> meshnames;
978 string pathPrefix;
979
980 int ppIndex = domainChunks[0]->getSiloPath().find(':');
981 if (ppIndex != string::npos) {
982 pathPrefix = domainChunks[0]->getSiloPath().substr(0, ppIndex+1);
983 }
984
985 // find a variable that is defined on this mesh to get the sample
986 // distribution (which tells us which ranks contribute to this mesh).
987 // Try mesh variables first, then regular ones.
988 VarVector::const_iterator viIt;
989 for (viIt = meshVariables.begin(); viIt != meshVariables.end(); viIt++) {
990 if (meshName == viIt->dataChunks[0]->getMeshName())
991 break;
992 }
993
994 if (viIt == meshVariables.end()) {
995 for (viIt = variables.begin(); viIt != variables.end(); viIt++) {
996 if (meshName == viIt->dataChunks[0]->getMeshName())
997 break;
998 }
999 }
1000 // this probably means that the mesh is empty
1001 if (viIt == variables.end()) {
1002 return;
1003 }
1004
1005 for (size_t idx = 0; idx < viIt->sampleDistribution.size(); idx++) {
1006 if (viIt->sampleDistribution[idx] > 0) {
1007 stringstream siloPath;
1008 siloPath << pathPrefix << "/block";
1009 int prevWidth = siloPath.width(4);
1010 char prevFill = siloPath.fill('0');
1011 siloPath << right << idx;
1012 siloPath.width(prevWidth);
1013 siloPath.fill(prevFill);
1014 siloPath << "/";
1015 siloPath << meshName;
1016 tempstrings.push_back(siloPath.str());
1017 meshnames.push_back((char*)tempstrings.back().c_str());
1018 meshtypes.push_back(DB_UCDMESH);
1019 }
1020 }
1021
1022 // ignore empty mesh
1023 if (!meshnames.empty()) {
1024 DBSetDir(dbfile, "/");
1025 DBoptlist* optList = DBMakeOptlist(2);
1026 DBAddOption(optList, DBOPT_CYCLE, &cycle);
1027 DBAddOption(optList, DBOPT_DTIME, &time);
1028 DBPutMultimesh(dbfile, meshName.c_str(), meshnames.size(),
1029 &meshnames[0], &meshtypes[0], optList);
1030 DBFreeOptlist(optList);
1031 }
1032 #endif
1033 }
1034
1035 //
1036 //
1037 //
1038 void EscriptDataset::putSiloMultiVar(DBfile* dbfile, const VarInfo& vi,
1039 bool useMeshFile)
1040 {
1041 #if USE_SILO
1042 vector<int> vartypes;
1043 vector<string> tempstrings;
1044 vector<char*> varnames;
1045 string pathPrefix;
1046 if (useMeshFile) {
1047 int ppIndex = domainChunks[0]->getSiloPath().find(':');
1048 if (ppIndex != string::npos) {
1049 pathPrefix = domainChunks[0]->getSiloPath().substr(0, ppIndex+1);
1050 }
1051 }
1052
1053 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
1054 if (vi.sampleDistribution[idx] > 0) {
1055 stringstream siloPath;
1056 siloPath << pathPrefix << "/block";
1057 int prevWidth = siloPath.width(4);
1058 char prevFill = siloPath.fill('0');
1059 siloPath << right << idx;
1060 siloPath.width(prevWidth);
1061 siloPath.fill(prevFill);
1062 siloPath << "/";
1063 siloPath << vi.varName;
1064 tempstrings.push_back(siloPath.str());
1065 varnames.push_back((char*)tempstrings.back().c_str());
1066 vartypes.push_back(DB_UCDVAR);
1067 }
1068 }
1069
1070 // ignore empty variables
1071 if (!varnames.empty()) {
1072 DBSetDir(dbfile, "/");
1073 DBoptlist* optList = DBMakeOptlist(2);
1074 DBAddOption(optList, DBOPT_CYCLE, &cycle);
1075 DBAddOption(optList, DBOPT_DTIME, &time);
1076 if (useMeshFile) {
1077 string vpath = string(MESH_VARS)+vi.varName;
1078 DBPutMultivar(dbfile, vpath.c_str(), varnames.size(),
1079 &varnames[0], &vartypes[0], optList);
1080 } else {
1081 DBPutMultivar(dbfile, vi.varName.c_str(), varnames.size(),
1082 &varnames[0], &vartypes[0], optList);
1083 }
1084 DBFreeOptlist(optList);
1085 }
1086 #endif
1087 }
1088
1089 //
1090 //
1091 //
1092 void EscriptDataset::putSiloMultiTensor(DBfile* dbfile, const VarInfo& vi)
1093 {
1094 #if USE_SILO
1095 string tensorDir = vi.varName+string("_comps/");
1096 DBSetDir(dbfile, "/");
1097 DBMkdir(dbfile, tensorDir.c_str());
1098 int one = 1;
1099 DBoptlist* optList = DBMakeOptlist(3);
1100 DBAddOption(optList, DBOPT_CYCLE, &cycle);
1101 DBAddOption(optList, DBOPT_DTIME, &time);
1102 DBAddOption(optList, DBOPT_HIDE_FROM_GUI, &one);
1103 const IntVec& shape = vi.dataChunks[0]->getShape();
1104
1105 for (int i=0; i<shape[1]; i++) {
1106 for (int j=0; j<shape[0]; j++) {
1107 vector<string> tempstrings;
1108 vector<char*> varnames;
1109 vector<int> vartypes;
1110 stringstream comp;
1111 comp << vi.varName << "_comps/a_";
1112 comp << i;
1113 comp << j;
1114 for (size_t idx = 0; idx < vi.sampleDistribution.size(); idx++) {
1115 if (vi.sampleDistribution[idx] > 0) {
1116 stringstream siloPath;
1117 siloPath << "/block";
1118 int prevWidth = siloPath.width(4);
1119 char prevFill = siloPath.fill('0');
1120 siloPath << right << idx;
1121 siloPath.width(prevWidth);
1122 siloPath.fill(prevFill);
1123 siloPath << "/" << comp.str();
1124 tempstrings.push_back(siloPath.str());
1125 varnames.push_back((char*)tempstrings.back().c_str());
1126 vartypes.push_back(DB_UCDVAR);
1127 }
1128 }
1129 if (!varnames.empty()) {
1130 DBPutMultivar(dbfile, comp.str().c_str(), varnames.size(),
1131 &varnames[0], &vartypes[0], optList);
1132 }
1133 }
1134 }
1135 DBFreeOptlist(optList);
1136 #endif
1137 }
1138
1139 } // namespace weipa
1140

  ViewVC Help
Powered by ViewVC 1.1.26