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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File size: 33572 byte(s)
Merging dudley and scons updates from branches

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

  ViewVC Help
Powered by ViewVC 1.1.26