/[escript]/branches/domexper/weipa/src/EscriptDataset.cpp
ViewVC logotype

Contents of /branches/domexper/weipa/src/EscriptDataset.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File size: 33552 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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

  ViewVC Help
Powered by ViewVC 1.1.26