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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3090 by caltinay, Mon Jun 14 00:59:30 2010 UTC revision 3091 by caltinay, Fri Aug 13 07:47:57 2010 UTC
# Line 46  const int NUM_SILO_FILES = 1; Line 46  const int NUM_SILO_FILES = 1;
46  // Default constructor  // Default constructor
47  //  //
48  EscriptDataset::EscriptDataset() :  EscriptDataset::EscriptDataset() :
     numParts(0),  
49      cycle(0),      cycle(0),
50      time(0.),      time(0.),
     keepMesh(false),  
51      externalMesh(false),      externalMesh(false),
52      mpiRank(0),      mpiRank(0),
53      mpiSize(1)      mpiSize(1)
# Line 61  EscriptDataset::EscriptDataset() : Line 59  EscriptDataset::EscriptDataset() :
59  //  //
60  #if HAVE_MPI  #if HAVE_MPI
61  EscriptDataset::EscriptDataset(MPI_Comm comm) :  EscriptDataset::EscriptDataset(MPI_Comm comm) :
     numParts(0),  
62      cycle(0),      cycle(0),
63      time(0.),      time(0.),
     keepMesh(false),  
64      externalMesh(false),      externalMesh(false),
65      mpiComm(comm)      mpiComm(comm)
66  {  {
# Line 88  bool EscriptDataset::initFromEscript( Line 84  bool EscriptDataset::initFromEscript(
84          DataVec& escriptVars,          DataVec& escriptVars,
85          const StringVec& varNames)          const StringVec& varNames)
86  {  {
87  #ifndef VISIT_PLUGIN      // Set the domain
88      int myError;      if (!setDomain(escriptDomain)) {
89      numParts = 1;          return false;
     FinleyMesh_ptr mesh(new FinleyMesh());  
     if (mesh->initFromEscript(escriptDomain)) {  
         if (mpiSize > 1)  
             mesh->reorderGhostZones(mpiRank);  
         meshBlocks.push_back(mesh);  
         myError = false;  
     } else {  
         mesh.reset();  
         myError = true;  
90      }      }
91    
92      int gError;      // initialize variables
93  #if HAVE_MPI      DataVec::iterator varIt = escriptVars.begin();
94      MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_LOR, mpiComm);      StringVec::const_iterator nameIt = varNames.begin();
95  #else      for (; varIt != escriptVars.end(); varIt++, nameIt++) {
96      gError = myError;          addData(*varIt, *nameIt);
 #endif  
   
     if (!gError) {  
         // initialize variables  
         DataVec::iterator varIt = escriptVars.begin();  
         StringVec::const_iterator nameIt = varNames.begin();  
         for (; varIt != escriptVars.end(); varIt++, nameIt++) {  
             VarInfo vi;  
             vi.varName = *nameIt;  
   
             DataVar_ptr var(new DataVar(vi.varName));  
             if (var->initFromEscript(*varIt, mesh)) {  
                 vi.dataBlocks.push_back(var);  
                 updateSampleDistribution(vi);  
                 vi.valid = true;  
             } else {  
                 var.reset();  
                 vi.valid = false;  
             }  
             variables.push_back(vi);  
         }  
   
         // Convert mesh data to variables  
         convertMeshVariables();  
97      }      }
98    
99      return !gError;      return true;
   
 #else // VISIT_PLUGIN  
     return false;  
 #endif  
100  }  }
101    
102  //  //
# Line 147  bool EscriptDataset::loadNetCDF(const st Line 106  bool EscriptDataset::loadNetCDF(const st
106                                  const StringVec& varFiles,                                  const StringVec& varFiles,
107                                  const StringVec& varNames, int nBlocks)                                  const StringVec& varNames, int nBlocks)
108  {  {
109      if (mpiSize > 1 && nBlocks != mpiSize) {      // Load the domain files
110          cerr << "Cannot load " << nBlocks << " chunks on " << mpiSize      if (!setDomain(meshFilePattern, nBlocks)) {
             << " MPI ranks!" << endl;  
111          return false;          return false;
112      }      }
113    
     numParts = nBlocks;  
     meshFmt = meshFilePattern;  
       
     // Load the mesh files  
     if (!loadMeshFromNetCDF()) {  
         cerr << "Reading the domain failed." << endl;  
         return false;  
     }  
   
     // Convert mesh data to variables  
     convertMeshVariables();  
   
114      // Load the variables      // Load the variables
115      StringVec::const_iterator fileIt = varFiles.begin();      StringVec::const_iterator fileIt = varFiles.begin();
116      StringVec::const_iterator nameIt = varNames.begin();      StringVec::const_iterator nameIt = varNames.begin();
117      for (; fileIt != varFiles.end(); fileIt++, nameIt++) {      for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
118          loadVarFromNetCDF(*fileIt, *nameIt);          addData(*fileIt, *nameIt);
119      }      }
120    
121      return true;      return true;
122  }  }
123    
124  //  //
125  // Load only variables using provided mesh  // Load only variables using provided domain
126  //  //
127  bool EscriptDataset::loadNetCDF(const MeshBlocks& mesh,  bool EscriptDataset::loadNetCDF(const MeshBlocks& mesh,
128                                  const StringVec& varFiles,                                  const StringVec& varFiles,
129                                  const StringVec& varNames)                                  const StringVec& varNames)
130  {  {
131      if (mpiSize > 1 && mesh.size() > 1) {      // Set the domain
132          cerr << "Can only read one mesh block per rank when using MPI!" << endl;      if (!setDomain(mesh)) {
133          return false;          return false;
134      }      }
135    
     externalMesh = true;  
     meshBlocks = mesh;  
     numParts = meshBlocks.size();  
     keepMesh = true;  
   
136      // Load the variables      // Load the variables
137      StringVec::const_iterator fileIt = varFiles.begin();      StringVec::const_iterator fileIt = varFiles.begin();
138      StringVec::const_iterator nameIt = varNames.begin();      StringVec::const_iterator nameIt = varNames.begin();
139      for (; fileIt != varFiles.end(); fileIt++, nameIt++) {      for (; fileIt != varFiles.end(); fileIt++, nameIt++) {
140          loadVarFromNetCDF(*fileIt, *nameIt);          addData(*fileIt, *nameIt);
141      }      }
142    
143      return true;      return true;
# Line 208  bool EscriptDataset::loadNetCDF(const Me Line 149  bool EscriptDataset::loadNetCDF(const Me
149  bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)  bool EscriptDataset::saveSilo(string fileName, bool useMultiMesh)
150  {  {
151  #if USE_SILO  #if USE_SILO
152      if (numParts == 0)      if (meshBlocks.size() == 0)
153          return false;          return false;
154    
155      const char* blockDirFmt = "/block%04d";      const char* blockDirFmt = "/block%04d";
# Line 343  bool EscriptDataset::saveSilo(string fil Line 284  bool EscriptDataset::saveSilo(string fil
284  //  //
285  bool EscriptDataset::saveVTK(string fileName)  bool EscriptDataset::saveVTK(string fileName)
286  {  {
287      if (numParts == 0)      if (meshBlocks.size() == 0)
288          return false;          return false;
289    
290      string meshName;      string meshName;
# Line 433  bool EscriptDataset::saveVTK(string file Line 374  bool EscriptDataset::saveVTK(string file
374          return false;          return false;
375      }      }
376    
     int myNumCells;  
     if (mpiSize > 1)  
         meshBlocks[0]->removeGhostZones(mpiRank);  
     ElementData_ptr elements = meshBlocks[0]->getElementsByName(meshName);  
     myNumCells = elements->getNumElements();  
     MPI_Reduce(&myNumCells, &gNumCells, 1, MPI_INT, MPI_SUM, 0, mpiComm);  
   
     int myCellSizeAndType[2];  
     myCellSizeAndType[0] = elements->getNodesPerElement();  
     myCellSizeAndType[1] = elements->getType();  
   
     // rank 0 needs to know element type and size but it's possible that this  
     // information is only available on other ranks (value=0) so retrieve it  
     MPI_Reduce(&myCellSizeAndType, &gCellSizeAndType, 2, MPI_INT, MPI_MAX, 0,  
             mpiComm);  
   
377      // these definitions make the code more readable and life easier -      // these definitions make the code more readable and life easier -
378      // WRITE_ORDERED causes all ranks to write the contents of the stream while      // WRITE_ORDERED causes all ranks to write the contents of the stream while
379      // WRITE_SHARED should only be called for rank 0.      // WRITE_SHARED should only be called for rank 0.
# Line 476  bool EscriptDataset::saveVTK(string file Line 401  bool EscriptDataset::saveVTK(string file
401    
402      ofstream ofs(fileName.c_str());      ofstream ofs(fileName.c_str());
403    
     int idx = 0;  
     for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {  
         if (numParts > 1)  
             (*meshIt)->removeGhostZones(idx);  
         ElementData_ptr elements = (*meshIt)->getElementsByName(meshName);  
         gNumCells += elements->getNumElements();  
         if (gCellSizeAndType[0] == 0)  
             gCellSizeAndType[0] = elements->getNodesPerElement();  
         if (gCellSizeAndType[1] == 0)  
             gCellSizeAndType[1] = elements->getType();  
     }  
   
404      // the non-MPI versions of WRITE_ORDERED and WRITE_SHARED are      // the non-MPI versions of WRITE_ORDERED and WRITE_SHARED are
405      // identical.      // identical.
406  #define WRITE_ORDERED(_stream) \  #define WRITE_ORDERED(_stream) \
# Line 502  bool EscriptDataset::saveVTK(string file Line 415  bool EscriptDataset::saveVTK(string file
415          _stream.str("");\          _stream.str("");\
416      } while(0)      } while(0)
417    
418    #endif // HAVE_MPI
419    
420        if (mpiSize > 1) {
421    #if HAVE_MPI
422            meshBlocks[0]->removeGhostZones(mpiRank);
423            ElementData_ptr elements = meshBlocks[0]->getElementsByName(meshName);
424            int myNumCells = elements->getNumElements();
425            MPI_Reduce(&myNumCells, &gNumCells, 1, MPI_INT, MPI_SUM, 0, mpiComm);
426    
427            int myCellSizeAndType[2];
428            myCellSizeAndType[0] = elements->getNodesPerElement();
429            myCellSizeAndType[1] = elements->getType();
430    
431            // rank 0 needs to know element type and size but it's possible that
432            // this information is only available on other ranks (value=0) so
433            // retrieve it
434            MPI_Reduce(&myCellSizeAndType, &gCellSizeAndType, 2, MPI_INT, MPI_MAX,
435                    0, mpiComm);
436  #endif  #endif
437        } else {
438            int idx = 0;
439            for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, idx++) {
440                if (meshBlocks.size() > 1)
441                    (*meshIt)->removeGhostZones(idx);
442                ElementData_ptr elements = (*meshIt)->getElementsByName(meshName);
443                gNumCells += elements->getNumElements();
444                if (gCellSizeAndType[0] == 0)
445                    gCellSizeAndType[0] = elements->getNodesPerElement();
446                if (gCellSizeAndType[1] == 0)
447                    gCellSizeAndType[1] = elements->getType();
448            }
449        }
450    
451      gNumPoints = meshBlocks[0]->getNodes()->getGlobalNumNodes();      gNumPoints = meshBlocks[0]->getNodes()->getGlobalNumNodes();
452    
# Line 541  bool EscriptDataset::saveVTK(string file Line 485  bool EscriptDataset::saveVTK(string file
485    
486      //      //
487      // coordinates (note that we are using the original nodes)      // coordinates (note that we are using the original nodes)
488      //      //
489      for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++) {      int blockNum = (mpiSize>1 ? mpiRank : 0);
490          (*meshIt)->getNodes()->writeCoordinatesVTK(oss, mpiRank);      for (meshIt = meshBlocks.begin(); meshIt != meshBlocks.end(); meshIt++, blockNum++) {
491            (*meshIt)->getNodes()->writeCoordinatesVTK(oss, blockNum);
492      }      }
493    
494      WRITE_ORDERED(oss);      WRITE_ORDERED(oss);
# Line 646  void EscriptDataset::writeVarToVTK(const Line 591  void EscriptDataset::writeVarToVTK(const
591              << "\" format=\"ascii\">" << endl;              << "\" format=\"ascii\">" << endl;
592      }      }
593    
594        // this is required in case we read a dataset with more than one chunk on
595        // one rank
596        int blockNum = (mpiSize>1 ? mpiRank : 0);
597      DataBlocks::const_iterator blockIt;      DataBlocks::const_iterator blockIt;
598      for (blockIt = varBlocks.begin(); blockIt != varBlocks.end(); blockIt++) {      for (blockIt = varBlocks.begin(); blockIt != varBlocks.end(); blockIt++, blockNum++) {
599          (*blockIt)->writeToVTK(os, mpiRank);          (*blockIt)->writeToVTK(os, blockNum);
600      }      }
601  }  }
602    
603  //  //
604    // Sets the domain using an escript::AbstractDomain instance.
605  //  //
606  //  bool EscriptDataset::setDomain(const escript::AbstractDomain* domain)
 bool EscriptDataset::loadMeshFromNetCDF()  
607  {  {
608      bool ok = true;  #ifndef VISIT_PLUGIN
609      char* str = new char[meshFmt.length()+10];      int myError = 0, gError;
610      if (mpiSize > 1) {  
611          FinleyMesh_ptr meshPart(new FinleyMesh());      // fail if the domain has already been set
612          sprintf(str, meshFmt.c_str(), mpiRank);      if (meshBlocks.size() > 0) {
613          string meshfile = str;          cerr << "Domain has already been set!" << endl;
614          if (meshPart->initFromNetCDF(meshfile)) {          myError = 1;
615              meshPart->reorderGhostZones(mpiRank);  
616              meshBlocks.push_back(meshPart);      } else {
617            FinleyMesh_ptr mesh(new FinleyMesh());
618            if (mesh->initFromEscript(domain)) {
619                if (mpiSize > 1)
620                    mesh->reorderGhostZones(mpiRank);
621                meshBlocks.push_back(mesh);
622          } else {          } else {
623              meshPart.reset();              cerr << "Error initializing domain!" << endl;
624              ok = false;              myError = 1;
625          }          }
626        }
627    
628        if (mpiSize > 1) {
629    #if HAVE_MPI
630            MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
631    #endif
632      } else {      } else {
633          for (int idx=0; idx < numParts; idx++) {          gError = myError;
634        }
635    
636        if (gError) {
637            meshBlocks.clear();
638        } else {
639            // Convert mesh data to variables
640            convertMeshVariables();
641        }
642        return !gError;
643    
644    #else // VISIT_PLUGIN
645        return false;
646    #endif
647    }
648    
649    //
650    // Sets the domain from dump files.
651    //
652    bool EscriptDataset::setDomain(const string filePattern, int nBlocks)
653    {
654        int myError = 0, gError;
655    
656        if (mpiSize > 1 && nBlocks != mpiSize) {
657            cerr << "Cannot load " << nBlocks << " chunks on " << mpiSize
658                << " MPI ranks!" << endl;
659            myError = 1;
660    
661        } else if (meshBlocks.size() > 0) {
662            cerr << "Domain has already been set!" << endl;
663            myError = 1;
664    
665        } else {
666            char* str = new char[filePattern.length()+10];
667            if (mpiSize > 1) {
668              FinleyMesh_ptr meshPart(new FinleyMesh());              FinleyMesh_ptr meshPart(new FinleyMesh());
669              sprintf(str, meshFmt.c_str(), idx);              sprintf(str, filePattern.c_str(), mpiRank);
670              string meshfile = str;              string meshfile = str;
671              if (meshPart->initFromNetCDF(meshfile)) {              if (meshPart->initFromNetCDF(meshfile)) {
672                  if (numParts > 1)                  meshPart->reorderGhostZones(mpiRank);
                     meshPart->reorderGhostZones(idx);  
673                  meshBlocks.push_back(meshPart);                  meshBlocks.push_back(meshPart);
674              } else {              } else {
675                  meshPart.reset();                  cerr << "Error initializing domain!" << endl;
676                  ok = false;                  myError = 1;
677                  break;              }
678            } else {
679                for (int idx=0; idx < nBlocks; idx++) {
680                    FinleyMesh_ptr meshPart(new FinleyMesh());
681                    sprintf(str, filePattern.c_str(), idx);
682                    string meshfile = str;
683                    if (meshPart->initFromNetCDF(meshfile)) {
684                        if (nBlocks > 1)
685                            meshPart->reorderGhostZones(idx);
686                        meshBlocks.push_back(meshPart);
687                    } else {
688                        cerr << "Error initializing domain block " << idx << endl;
689                        myError = 1;
690                        break;
691                    }
692              }              }
693          }          }
694            delete[] str;
695        }
696    
697        if (mpiSize > 1) {
698    #if HAVE_MPI
699            MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
700    #endif
701        } else {
702            gError = myError;
703        }
704    
705        if (gError) {
706            meshBlocks.clear();
707        } else {
708            // Convert mesh data to variables
709            convertMeshVariables();
710        }
711    
712        return !gError;
713    }
714    
715    //
716    // Sets an already converted domain.
717    //
718    bool EscriptDataset::setDomain(const MeshBlocks& domain)
719    {
720        int myError = 0, gError;
721    
722        if (mpiSize > 1 && domain.size() > 1) {
723            cerr << "Can only add one domain block per rank when using MPI!"
724                << endl;
725            myError = 1;
726    
727        } else if (meshBlocks.size() > 0) {
728            cerr << "Domain has already been set!" << endl;
729            myError = 1;
730    
731        }
732    
733        if (mpiSize > 1) {
734    #if HAVE_MPI
735            MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
736    #endif
737        } else {
738            gError = myError;
739        }
740    
741        if (!gError) {
742            externalMesh = true;
743            meshBlocks = domain;
744      }      }
745      delete[] str;  
746      return ok;      return !gError;
747  }  }
748    
749  //  //
750  //  //
751  //  //
752  bool EscriptDataset::loadVarFromNetCDF(const string& fileName,  bool EscriptDataset::addData(escript::Data& data, const string name)
                                        const string& varName)  
753  {  {
754      int myError = false;  #ifndef VISIT_PLUGIN
755      char* str = new char[fileName.length()+10];      bool success = true;
     VarInfo vi;  
756    
757      vi.varName = varName;      // fail if no domain has been set
758      vi.valid = true;      if (meshBlocks.size() == 0) {
759            success = false;
760    
761      // read all parts of the variable      } else {
762      MeshBlocks::iterator mIt;          // initialize variable
763      int idx = (mpiSize > 1) ? mpiRank : 0;          VarInfo vi;
764      for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++, idx++) {          vi.varName = name;
765          sprintf(str, fileName.c_str(), idx);  
766          string dfile = str;          DataVar_ptr var(new DataVar(vi.varName));
767          DataVar_ptr var(new DataVar(varName));          if (var->initFromEscript(data, meshBlocks[0])) {
         if (var->initFromNetCDF(dfile, *mIt))  
768              vi.dataBlocks.push_back(var);              vi.dataBlocks.push_back(var);
769          else {              updateSampleDistribution(vi);
770              cerr << "Error reading " << dfile << endl;              vi.valid = true;
771            } else {
772              var.reset();              var.reset();
773              myError = true;              vi.valid = false;
             break;  
774          }          }
775            variables.push_back(vi);
776      }      }
777      delete[] str;      return success;
778    
779      int gError = myError;  #else // VISIT_PLUGIN
780        return false;
781    #endif
782    }
783    
784      if (mpiSize > 1) {  //
785    //
786    //
787    bool EscriptDataset::addData(const string filePattern, const string name)
788    {
789        int myError = 0, gError;
790    
791        // fail if no domain has been set
792        if (meshBlocks.size() == 0) {
793            gError = 1;
794    
795        } else {
796            // initialize variable
797            VarInfo vi;
798            vi.varName = name;
799            vi.valid = true;
800            char* str = new char[filePattern.length()+10];
801    
802            // read all parts of the variable
803            MeshBlocks::iterator mIt;
804            int idx = (mpiSize > 1) ? mpiRank : 0;
805            for (mIt = meshBlocks.begin(); mIt != meshBlocks.end(); mIt++, idx++) {
806                sprintf(str, filePattern.c_str(), idx);
807                string dfile = str;
808                DataVar_ptr var(new DataVar(name));
809                if (var->initFromNetCDF(dfile, *mIt))
810                    vi.dataBlocks.push_back(var);
811                else {
812                    cerr << "Error reading " << dfile << endl;
813                    myError = 1;
814                    break;
815                }
816            }
817            delete[] str;
818    
819            if (mpiSize > 1) {
820  #if HAVE_MPI  #if HAVE_MPI
821          MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_LOR, mpiComm);              MPI_Allreduce(&myError, &gError, 1, MPI_INT, MPI_MAX, mpiComm);
822  #endif  #endif
823      }          } else {
824                gError = myError;
825            }
826    
827      if (gError) {          if (!gError) {
828          // at least one chunk was not read correctly              // only add variable if all chunks have been read without error
829          vi.dataBlocks.clear();              updateSampleDistribution(vi);
830          vi.valid = false;              variables.push_back(vi);
831            }
832      }      }
833    
834      updateSampleDistribution(vi);      return !gError;
     variables.push_back(vi);  
     return vi.valid;  
835  }  }
836    
837  //  //

Legend:
Removed from v.3090  
changed lines
  Added in v.3091

  ViewVC Help
Powered by ViewVC 1.1.26