src/io/DumpReader.cpp

/* 
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. 
 * 
 * The University of Notre Dame grants you ("Licensee") a 
 * non-exclusive, royalty free, license to use, modify and 
 * redistribute this software in source and binary code form, provided 
 * that the following conditions are met: 
 * 
 * 1. Acknowledgement of the program authors must be made in any 
 *    publication of scientific results based in part on use of the 
 *    program.  An acceptable form of acknowledgement is citation of 
 *    the article in which the program was described (Matthew 
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher 
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented 
 *    Parallel Simulation Engine for Molecular Dynamics," 
 *    J. Comput. Chem. 26, pp. 252-271 (2005)) 
 * 
 * 2. Redistributions of source code must retain the above copyright 
 *    notice, this list of conditions and the following disclaimer. 
 * 
 * 3. Redistributions in binary form must reproduce the above copyright 
 *    notice, this list of conditions and the following disclaimer in the 
 *    documentation and/or other materials provided with the 
 *    distribution. 
 * 
 * This software is provided "AS IS," without a warranty of any 
 * kind. All express or implied conditions, representations and 
 * warranties, including any implied warranty of merchantability, 
 * fitness for a particular purpose or non-infringement, are hereby 
 * excluded.  The University of Notre Dame and its licensors shall not 
 * be liable for any damages suffered by licensee as a result of 
 * using, modifying or distributing the software or its 
 * derivatives. In no event will the University of Notre Dame or its 
 * licensors be liable for any lost revenue, profit or data, or for 
 * direct, indirect, special, consequential, incidental or punitive 
 * damages, however caused and regardless of the theory of liability, 
 * arising out of the use of or inability to use software, even if the 
 * University of Notre Dame has been advised of the possibility of 
 * such damages. 
 */ 
  
#define _LARGEFILE_SOURCE64 
#define _FILE_OFFSET_BITS 64 
 
#include <sys/types.h> 
#include <sys/stat.h> 
 
#include <iostream> 
#include <math.h> 
 
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 
 
#include "io/DumpReader.hpp" 
#include "primitives/Molecule.hpp" 
#include "utils/simError.h" 
#include "utils/MemoryUtils.hpp" 
#include "utils/StringTokenizer.hpp" 
 
#ifdef IS_MPI 
 
#include <mpi.h> 
#define TAKE_THIS_TAG_CHAR 0 
#define TAKE_THIS_TAG_INT 1 
 
#endif // is_mpi 
 
 
namespace oopse { 
   
  DumpReader::DumpReader(SimInfo* info, const std::string& filename) 
    : info_(info), filename_(filename), isScanned_(false), nframes_(0) { 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif 
      
      inFile_ = new std::ifstream(filename_.c_str()); 
      
      if (inFile_->fail()) { 
        sprintf(painCave.errMsg, 
                "DumpReader: Cannot open file: %s\n", 
                filename_.c_str()); 
        painCave.isFatal = 1; 
        simError(); 
      } 
      
#ifdef IS_MPI 
      
    } 
    
    strcpy(checkPointMsg, "Dump file opened for reading successfully."); 
    MPIcheckPoint(); 
    
#endif 
    
    return; 
  } 
  
  DumpReader::~DumpReader() { 
    
#ifdef IS_MPI 
    
    if (worldRank == 0) { 
#endif 
      
      delete inFile_; 
      
#ifdef IS_MPI 
      
    } 
    
    strcpy(checkPointMsg, "Dump file closed successfully."); 
    MPIcheckPoint(); 
    
#endif 
    
    return; 
  } 
  
  int DumpReader::getNFrames(void) { 
     
    if (!isScanned_) 
      scanFile(); 
     
    return nframes_; 
  } 
   
  void DumpReader::scanFile(void) { 
    int i, j; 
    int lineNum = 0; 
    char readBuffer[maxBufferSize]; 
    std::streampos  currPos; 
     
#ifdef IS_MPI 
     
    if (worldRank == 0) { 
#endif // is_mpi 
       
      inFile_->seekg (0, std::ios::beg); 
       
 
      currPos = inFile_->tellg(); 
      inFile_->getline(readBuffer, sizeof(readBuffer)); 
      lineNum++; 
       
      if (inFile_->eof()) { 
        sprintf(painCave.errMsg, 
                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n", 
                filename_.c_str(), 
                lineNum); 
        painCave.isFatal = 1; 
        simError(); 
      } 
       
      while (!inFile_->eof()) { 
        framePos_.push_back(currPos); 
         
        i = atoi(readBuffer); 
         
        inFile_->getline(readBuffer, sizeof(readBuffer)); 
        lineNum++; 
         
        if (inFile_->eof()) { 
          sprintf(painCave.errMsg, 
                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n", 
                  filename_.c_str(), 
                  lineNum); 
          painCave.isFatal = 1; 
          simError(); 
        } 
         
        for(j = 0; j < i; j++) { 
          inFile_->getline(readBuffer, sizeof(readBuffer)); 
          lineNum++; 
           
          if (inFile_->eof()) { 
            sprintf(painCave.errMsg, 
                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d," 
                    " with atom %d\n", filename_.c_str(), 
                    lineNum, 
                    j); 
             
            painCave.isFatal = 1; 
            simError(); 
          } 
        } 
         
        currPos = inFile_->tellg(); 
        inFile_->getline(readBuffer, sizeof(readBuffer)); 
        lineNum++; 
      } 
 
      inFile_->seekg (0, std::ios::beg); 
       
      nframes_ = framePos_.size(); 
#ifdef IS_MPI 
    } 
     
    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD); 
     
    strcpy(checkPointMsg, "Successfully scanned DumpFile\n"); 
    MPIcheckPoint(); 
     
#endif // is_mpi 
     
    isScanned_ = true; 
  } 
   
  void DumpReader::readFrame(int whichFrame) { 
    if (!isScanned_) 
      scanFile(); 
         
    int storageLayout = info_->getSnapshotManager()->getStorageLayout(); 
     
    if (storageLayout & DataStorage::dslPosition) { 
      needPos_ = true; 
    } else { 
      needPos_ = false; 
    } 
     
    if (storageLayout & DataStorage::dslVelocity) { 
      needVel_ = true; 
    } else { 
      needVel_ = false; 
    } 
     
    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) { 
      needQuaternion_ = true; 
    } else { 
      needQuaternion_ = false; 
    } 
     
    if (storageLayout & DataStorage::dslAngularMomentum) { 
      needAngMom_ = true; 
    } else { 
      needAngMom_ = false;     
    } 
     
    readSet(whichFrame); 
  } 
   
  void DumpReader::readSet(int whichFrame) { 
    int i; 
    int nTotObjs;                  // the number of atoms 
    char read_buffer[maxBufferSize];  //the line buffer for reading 
    char * eof_test;               // ptr to see when we reach the end of the file 
     
    Molecule* mol; 
    StuntDouble* integrableObject; 
    SimInfo::MoleculeIterator mi; 
    Molecule::IntegrableObjectIterator ii; 
     
#ifndef IS_MPI 
    inFile_->clear();  
    inFile_->seekg(framePos_[whichFrame]); 
         
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
      sprintf(painCave.errMsg, 
              "DumpReader error: error reading 1st line of \"%s\"\n", 
              filename_.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 
     
    nTotObjs = atoi(read_buffer); 
     
    if (nTotObjs != info_->getNGlobalIntegrableObjects()) { 
      sprintf(painCave.errMsg, 
              "DumpReader error. %s nIntegrable, %d, " 
              "does not match the meta-data file's nIntegrable, %d.\n", 
              filename_.c_str(), 
              nTotObjs, 
              info_->getNGlobalIntegrableObjects()); 
       
      painCave.isFatal = 1; 
      simError(); 
    } 
     
    //read the box mat from the comment line 
     
     
    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n", 
              filename_.c_str()); 
      painCave.isFatal = 1; 
      simError(); 
    } 
     
    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot()); 
     
    //parse dump lines 
     
    i = 0; 
    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { 
       
      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
           integrableObject = mol->nextIntegrableObject(ii)) {            
         
         
         
        if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
          sprintf(painCave.errMsg, 
                  "DumpReader Error: error in reading file %s\n" 
                  "natoms  = %d; index = %d\n" 
                  "error reading the line from the file.\n", 
                  filename_.c_str(), 
                  nTotObjs, 
                  i); 
           
          painCave.isFatal = 1; 
          simError(); 
        } 
         
        parseDumpLine(read_buffer, integrableObject); 
        i++; 
      } 
    } 
     
    // MPI Section of code.......... 
     
#else //IS_MPI 
     
    // first thing first, suspend fatalities. 
    int masterNode = 0; 
    int nCurObj; 
    painCave.isEventLoop = 1; 
     
    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone 
    int haveError; 
     
    MPI_Status istatus; 
    int nitems; 
     
    nTotObjs = info_->getNGlobalIntegrableObjects(); 
    haveError = 0; 
     
    if (worldRank == masterNode) { 
      inFile_->clear();             
      inFile_->seekg(framePos_[whichFrame]); 
       
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ", 
                filename_.c_str()); 
        painCave.isFatal = 1; 
        simError(); 
      } 
       
      nitems = atoi(read_buffer); 
       
      // Check to see that the number of integrable objects in the 
      // intial configuration file is the same as derived from the 
      // meta-data file. 
       
      if (nTotObjs != nitems) { 
        sprintf(painCave.errMsg, 
                "DumpReader Error. %s nIntegrable, %d, " 
                "does not match the meta-data file's nIntegrable, %d.\n", 
                filename_.c_str(), 
                nTotObjs, 
                info_->getNGlobalIntegrableObjects()); 
         
        painCave.isFatal = 1; 
        simError(); 
      } 
       
      //read the boxMat from the comment line 
       
       
       
      if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n", 
                filename_.c_str()); 
        painCave.isFatal = 1; 
        simError(); 
      } 
       
      //Every single processor will parse the comment line by itself 
      //By using this way, we might lose some efficiency, but if we want to add 
      //more parameters into comment line, we only need to modify function 
      //parseCommentLine 
       
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); 
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot()); 
       
      for(i = 0; i < info_->getNGlobalMolecules(); i++) { 
        int which_node = info_->getMolToProc(i); 
         
        if (which_node == masterNode) { 
          //molecules belong to master node 
           
          mol = info_->getMoleculeByGlobalIndex(i); 
           
          if (mol == NULL) { 
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank); 
            painCave.isFatal = 1; 
            simError(); 
          } 
           
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
               integrableObject = mol->nextIntegrableObject(ii)){ 
             
             
             
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
              sprintf(painCave.errMsg, 
                      "DumpReader Error: error in reading file %s\n" 
                      "natoms  = %d; index = %d\n" 
                      "error reading the line from the file.\n", 
                      filename_.c_str(), 
                      nTotObjs, 
                      i); 
               
              painCave.isFatal = 1; 
              simError(); 
            } 
             
            parseDumpLine(read_buffer, integrableObject); 
          } 
        } else { 
          //molecule belongs to slave nodes 
           
          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT, 
                   MPI_COMM_WORLD, &istatus); 
           
          for(int j = 0; j < nCurObj; j++) { 
             
             
            if (!inFile_->getline(read_buffer, sizeof(read_buffer))) { 
              sprintf(painCave.errMsg, 
                      "DumpReader Error: error in reading file %s\n" 
                      "natoms  = %d; index = %d\n" 
                      "error reading the line from the file.\n", 
                      filename_.c_str(), 
                      nTotObjs, 
                      i); 
               
              painCave.isFatal = 1; 
              simError(); 
            } 
             
            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node, 
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD); 
          } 
        } 
      } 
    } else { 
      //actions taken at slave nodes 
      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); 
       
      /**@todo*/ 
      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot()); 
       
      for(i = 0; i < info_->getNGlobalMolecules(); i++) { 
        int which_node = info_->getMolToProc(i); 
         
        if (which_node == worldRank) { 
          //molecule with global index i belongs to this processor 
           
          mol = info_->getMoleculeByGlobalIndex(i); 
          if (mol == NULL) { 
            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank); 
            painCave.isFatal = 1; 
            simError(); 
          } 
           
          nCurObj = mol->getNIntegrableObjects(); 
           
          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT, 
                   MPI_COMM_WORLD); 
           
          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
               integrableObject = mol->nextIntegrableObject(ii)){ 
             
            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode, 
                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus); 
             
            parseDumpLine(read_buffer, integrableObject); 
          } 
           
        } 
         
      } 
       
    } 
     
#endif 
     
  } 
   
  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) { 
     
    Vector3d pos;  // position place holders 
    Vector3d vel;  // velocity placeholders 
    Quat4d q;    // the quaternions 
    Vector3d ji;   // angular velocity placeholders; 
    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    if (nTokens < 14) { 
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough Tokens.\n%s\n", line); 
      painCave.isFatal = 1; 
      simError(); 
    } 
     
    std::string name = tokenizer.nextToken(); 
     
    if (name != integrableObject->getType()) { 
       
      sprintf(painCave.errMsg, 
              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n", 
              name.c_str(), filename_.c_str(), integrableObject->getType().c_str()); 
      painCave.isFatal = 1; 
      simError();         
    } 
     
    pos[0] = tokenizer.nextTokenAsDouble(); 
    pos[1] = tokenizer.nextTokenAsDouble(); 
    pos[2] = tokenizer.nextTokenAsDouble(); 
    if (needPos_) { 
      integrableObject->setPos(pos); 
    } 
     
    vel[0] = tokenizer.nextTokenAsDouble(); 
    vel[1] = tokenizer.nextTokenAsDouble(); 
    vel[2] = tokenizer.nextTokenAsDouble(); 
    if (needVel_) { 
      integrableObject->setVel(vel); 
    } 
     
    if (integrableObject->isDirectional()) { 
       
      q[0] = tokenizer.nextTokenAsDouble(); 
      q[1] = tokenizer.nextTokenAsDouble(); 
      q[2] = tokenizer.nextTokenAsDouble(); 
      q[3] = tokenizer.nextTokenAsDouble(); 
       
      RealType qlen = q.length(); 
      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0 
         
        sprintf(painCave.errMsg, 
                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n"); 
        painCave.isFatal = 1; 
        simError(); 
         
      }  
       
      q.normalize(); 
      if (needQuaternion_) {            
        integrableObject->setQ(q); 
      } 
       
      ji[0] = tokenizer.nextTokenAsDouble(); 
      ji[1] = tokenizer.nextTokenAsDouble(); 
      ji[2] = tokenizer.nextTokenAsDouble(); 
      if (needAngMom_) { 
        integrableObject->setJ(ji); 
      } 
    } 
     
  } 
   
   
  void DumpReader::parseCommentLine(char* line, Snapshot* s) { 
    RealType currTime; 
    Mat3x3d hmat; 
    RealType chi; 
    RealType integralOfChiDt; 
    Mat3x3d eta; 
     
    StringTokenizer tokenizer(line); 
    int nTokens; 
     
    nTokens = tokenizer.countTokens(); 
     
    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens) 
    if (nTokens < 10) { 
      sprintf(painCave.errMsg, 
              "DumpReader Error: Not enough tokens in comment line: %s", line); 
      painCave.isFatal = 1; 
      simError();    
    } 
     
    //read current time 
    currTime = tokenizer.nextTokenAsDouble(); 
    s->setTime(currTime); 
     
    //read h-matrix 
    hmat(0, 0) = tokenizer.nextTokenAsDouble(); 
    hmat(0, 1) = tokenizer.nextTokenAsDouble(); 
    hmat(0, 2) = tokenizer.nextTokenAsDouble(); 
    hmat(1, 0) = tokenizer.nextTokenAsDouble(); 
    hmat(1, 1) = tokenizer.nextTokenAsDouble(); 
    hmat(1, 2) = tokenizer.nextTokenAsDouble(); 
    hmat(2, 0) = tokenizer.nextTokenAsDouble(); 
    hmat(2, 1) = tokenizer.nextTokenAsDouble(); 
    hmat(2, 2) = tokenizer.nextTokenAsDouble(); 
    s->setHmat(hmat); 
     
    //read chi and integralOfChidt, they should apprear in pair 
    if (tokenizer.countTokens() >= 2) { 
      chi = tokenizer.nextTokenAsDouble(); 
      integralOfChiDt = tokenizer.nextTokenAsDouble();             
       
      s->setChi(chi); 
      s->setIntegralOfChiDt(integralOfChiDt); 
    } 
     
    //read eta (eta is 3x3 matrix) 
    if (tokenizer.countTokens() >= 9) { 
      eta(0, 0) = tokenizer.nextTokenAsDouble(); 
      eta(0, 1) = tokenizer.nextTokenAsDouble(); 
      eta(0, 2) = tokenizer.nextTokenAsDouble(); 
      eta(1, 0) = tokenizer.nextTokenAsDouble(); 
      eta(1, 1) = tokenizer.nextTokenAsDouble(); 
      eta(1, 2) = tokenizer.nextTokenAsDouble(); 
      eta(2, 0) = tokenizer.nextTokenAsDouble(); 
      eta(2, 1) = tokenizer.nextTokenAsDouble(); 
      eta(2, 2) = tokenizer.nextTokenAsDouble();       
       
      s->setEta(eta); 
    } 
     
     
  } 
   
}//end namespace oopse 
Revision:	996
Committed:	Wed Jun 28 14:35:14 2006 UTC (18 years, 10 months ago) by chrisfen
File size:	19385 byte(s)
Log Message:	more efficient file reading for thermodynamic integration, and fixed some formatting
#	User	Rev	Content
1	chrisfen	721	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			#define _LARGEFILE_SOURCE64
43			#define _FILE_OFFSET_BITS 64
44
45			#include <sys/types.h>
46			#include <sys/stat.h>
47
48			#include <iostream>
49			#include <math.h>
50
51			#include <stdio.h>
52			#include <stdlib.h>
53			#include <string.h>
54
55			#include "io/DumpReader.hpp"
56			#include "primitives/Molecule.hpp"
57			#include "utils/simError.h"
58			#include "utils/MemoryUtils.hpp"
59			#include "utils/StringTokenizer.hpp"
60
61			#ifdef IS_MPI
62
63			#include <mpi.h>
64			#define TAKE_THIS_TAG_CHAR 0
65			#define TAKE_THIS_TAG_INT 1
66
67			#endif // is_mpi
68
69
70			namespace oopse {
71
72			DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73			: info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74	chrisfen	996
75	chrisfen	721	#ifdef IS_MPI
76	chrisfen	996
77			if (worldRank == 0) {
78	chrisfen	721	#endif
79	chrisfen	996
80	chrisfen	721	inFile_ = new std::ifstream(filename_.c_str());
81	chrisfen	996
82			if (inFile_->fail()) {
83			sprintf(painCave.errMsg,
84			"DumpReader: Cannot open file: %s\n",
85			filename_.c_str());
86			painCave.isFatal = 1;
87			simError();
88			}
89
90	chrisfen	721	#ifdef IS_MPI
91	chrisfen	996
92			}
93
94			strcpy(checkPointMsg, "Dump file opened for reading successfully.");
95			MPIcheckPoint();
96
97	chrisfen	721	#endif
98	chrisfen	996
99			return;
100			}
101
102	chrisfen	721	DumpReader::~DumpReader() {
103	chrisfen	996
104	chrisfen	721	#ifdef IS_MPI
105	chrisfen	996
106	chrisfen	721	if (worldRank == 0) {
107			#endif
108	chrisfen	996
109	chrisfen	721	delete inFile_;
110	chrisfen	996
111	chrisfen	721	#ifdef IS_MPI
112	chrisfen	996
113	chrisfen	721	}
114	chrisfen	996
115	chrisfen	721	strcpy(checkPointMsg, "Dump file closed successfully.");
116			MPIcheckPoint();
117	chrisfen	996
118	chrisfen	721	#endif
119	chrisfen	996
120	chrisfen	721	return;
121			}
122	chrisfen	996
123	chrisfen	721	int DumpReader::getNFrames(void) {
124
125			if (!isScanned_)
126			scanFile();
127
128			return nframes_;
129			}
130
131			void DumpReader::scanFile(void) {
132			int i, j;
133			int lineNum = 0;
134			char readBuffer[maxBufferSize];
135			std::streampos currPos;
136
137			#ifdef IS_MPI
138
139			if (worldRank == 0) {
140			#endif // is_mpi
141
142			inFile_->seekg (0, std::ios::beg);
143
144
145			currPos = inFile_->tellg();
146			inFile_->getline(readBuffer, sizeof(readBuffer));
147			lineNum++;
148
149			if (inFile_->eof()) {
150			sprintf(painCave.errMsg,
151			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
152			filename_.c_str(),
153			lineNum);
154			painCave.isFatal = 1;
155			simError();
156			}
157
158			while (!inFile_->eof()) {
159			framePos_.push_back(currPos);
160
161			i = atoi(readBuffer);
162
163			inFile_->getline(readBuffer, sizeof(readBuffer));
164			lineNum++;
165
166			if (inFile_->eof()) {
167			sprintf(painCave.errMsg,
168			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
169			filename_.c_str(),
170			lineNum);
171			painCave.isFatal = 1;
172			simError();
173			}
174
175			for(j = 0; j < i; j++) {
176			inFile_->getline(readBuffer, sizeof(readBuffer));
177			lineNum++;
178
179			if (inFile_->eof()) {
180			sprintf(painCave.errMsg,
181			"DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
182			" with atom %d\n", filename_.c_str(),
183			lineNum,
184			j);
185
186			painCave.isFatal = 1;
187			simError();
188			}
189			}
190
191			currPos = inFile_->tellg();
192			inFile_->getline(readBuffer, sizeof(readBuffer));
193			lineNum++;
194			}
195
196			inFile_->seekg (0, std::ios::beg);
197
198			nframes_ = framePos_.size();
199			#ifdef IS_MPI
200			}
201
202			MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
203
204			strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
205			MPIcheckPoint();
206
207			#endif // is_mpi
208
209			isScanned_ = true;
210			}
211
212			void DumpReader::readFrame(int whichFrame) {
213			if (!isScanned_)
214			scanFile();
215
216			int storageLayout = info_->getSnapshotManager()->getStorageLayout();
217
218			if (storageLayout & DataStorage::dslPosition) {
219			needPos_ = true;
220			} else {
221			needPos_ = false;
222			}
223
224			if (storageLayout & DataStorage::dslVelocity) {
225			needVel_ = true;
226			} else {
227			needVel_ = false;
228			}
229
230			if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
231			needQuaternion_ = true;
232			} else {
233			needQuaternion_ = false;
234			}
235
236			if (storageLayout & DataStorage::dslAngularMomentum) {
237			needAngMom_ = true;
238			} else {
239			needAngMom_ = false;
240			}
241
242			readSet(whichFrame);
243			}
244
245			void DumpReader::readSet(int whichFrame) {
246			int i;
247			int nTotObjs; // the number of atoms
248			char read_buffer[maxBufferSize]; //the line buffer for reading
249			char * eof_test; // ptr to see when we reach the end of the file
250
251			Molecule* mol;
252			StuntDouble* integrableObject;
253			SimInfo::MoleculeIterator mi;
254			Molecule::IntegrableObjectIterator ii;
255
256			#ifndef IS_MPI
257			inFile_->clear();
258			inFile_->seekg(framePos_[whichFrame]);
259
260			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
261			sprintf(painCave.errMsg,
262			"DumpReader error: error reading 1st line of \"%s\"\n",
263			filename_.c_str());
264			painCave.isFatal = 1;
265			simError();
266			}
267
268			nTotObjs = atoi(read_buffer);
269
270			if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
271			sprintf(painCave.errMsg,
272			"DumpReader error. %s nIntegrable, %d, "
273			"does not match the meta-data file's nIntegrable, %d.\n",
274			filename_.c_str(),
275			nTotObjs,
276			info_->getNGlobalIntegrableObjects());
277
278			painCave.isFatal = 1;
279			simError();
280			}
281
282			//read the box mat from the comment line
283
284
285			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
286			sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
287			filename_.c_str());
288			painCave.isFatal = 1;
289			simError();
290			}
291
292			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
293
294			//parse dump lines
295
296			i = 0;
297			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
298
299			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
300			integrableObject = mol->nextIntegrableObject(ii)) {
301
302
303
304			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
305			sprintf(painCave.errMsg,
306			"DumpReader Error: error in reading file %s\n"
307			"natoms = %d; index = %d\n"
308			"error reading the line from the file.\n",
309			filename_.c_str(),
310			nTotObjs,
311			i);
312
313			painCave.isFatal = 1;
314			simError();
315			}
316
317			parseDumpLine(read_buffer, integrableObject);
318			i++;
319			}
320			}
321
322			// MPI Section of code..........
323
324			#else //IS_MPI
325
326			// first thing first, suspend fatalities.
327			int masterNode = 0;
328			int nCurObj;
329			painCave.isEventLoop = 1;
330
331			int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
332			int haveError;
333
334			MPI_Status istatus;
335			int nitems;
336
337			nTotObjs = info_->getNGlobalIntegrableObjects();
338			haveError = 0;
339
340			if (worldRank == masterNode) {
341			inFile_->clear();
342			inFile_->seekg(framePos_[whichFrame]);
343
344			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
345			sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
346			filename_.c_str());
347			painCave.isFatal = 1;
348			simError();
349			}
350
351			nitems = atoi(read_buffer);
352
353			// Check to see that the number of integrable objects in the
354			// intial configuration file is the same as derived from the
355			// meta-data file.
356
357			if (nTotObjs != nitems) {
358			sprintf(painCave.errMsg,
359			"DumpReader Error. %s nIntegrable, %d, "
360			"does not match the meta-data file's nIntegrable, %d.\n",
361			filename_.c_str(),
362			nTotObjs,
363			info_->getNGlobalIntegrableObjects());
364
365			painCave.isFatal = 1;
366			simError();
367			}
368
369			//read the boxMat from the comment line
370
371
372
373			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
374			sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
375			filename_.c_str());
376			painCave.isFatal = 1;
377			simError();
378			}
379
380			//Every single processor will parse the comment line by itself
381			//By using this way, we might lose some efficiency, but if we want to add
382			//more parameters into comment line, we only need to modify function
383			//parseCommentLine
384
385			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
386			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
387
388			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
389			int which_node = info_->getMolToProc(i);
390
391			if (which_node == masterNode) {
392			//molecules belong to master node
393
394			mol = info_->getMoleculeByGlobalIndex(i);
395
396			if (mol == NULL) {
397			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
398			painCave.isFatal = 1;
399			simError();
400			}
401
402			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
403			integrableObject = mol->nextIntegrableObject(ii)){
404
405
406
407			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
408			sprintf(painCave.errMsg,
409			"DumpReader Error: error in reading file %s\n"
410			"natoms = %d; index = %d\n"
411			"error reading the line from the file.\n",
412			filename_.c_str(),
413			nTotObjs,
414			i);
415
416			painCave.isFatal = 1;
417			simError();
418			}
419
420			parseDumpLine(read_buffer, integrableObject);
421			}
422			} else {
423			//molecule belongs to slave nodes
424
425			MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
426			MPI_COMM_WORLD, &istatus);
427
428			for(int j = 0; j < nCurObj; j++) {
429
430
431			if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
432			sprintf(painCave.errMsg,
433			"DumpReader Error: error in reading file %s\n"
434			"natoms = %d; index = %d\n"
435			"error reading the line from the file.\n",
436			filename_.c_str(),
437			nTotObjs,
438			i);
439
440			painCave.isFatal = 1;
441			simError();
442			}
443
444			MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
445			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
446			}
447			}
448			}
449			} else {
450			//actions taken at slave nodes
451			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
452
453			/*@todo/
454			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
455
456			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
457			int which_node = info_->getMolToProc(i);
458
459			if (which_node == worldRank) {
460			//molecule with global index i belongs to this processor
461
462			mol = info_->getMoleculeByGlobalIndex(i);
463			if (mol == NULL) {
464			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
465			painCave.isFatal = 1;
466			simError();
467			}
468
469			nCurObj = mol->getNIntegrableObjects();
470
471			MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
472			MPI_COMM_WORLD);
473
474			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
475			integrableObject = mol->nextIntegrableObject(ii)){
476
477			MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
478			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
479
480			parseDumpLine(read_buffer, integrableObject);
481			}
482
483			}
484
485			}
486
487			}
488
489			#endif
490
491			}
492
493			void DumpReader::parseDumpLine(char line, StuntDouble integrableObject) {
494
495			Vector3d pos; // position place holders
496			Vector3d vel; // velocity placeholders
497			Quat4d q; // the quaternions
498			Vector3d ji; // angular velocity placeholders;
499			StringTokenizer tokenizer(line);
500			int nTokens;
501
502			nTokens = tokenizer.countTokens();
503
504			if (nTokens < 14) {
505			sprintf(painCave.errMsg,
506			"DumpReader Error: Not enough Tokens.\n%s\n", line);
507			painCave.isFatal = 1;
508			simError();
509			}
510
511			std::string name = tokenizer.nextToken();
512
513			if (name != integrableObject->getType()) {
514
515			sprintf(painCave.errMsg,
516			"DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
517			name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
518			painCave.isFatal = 1;
519			simError();
520			}
521
522			pos[0] = tokenizer.nextTokenAsDouble();
523			pos[1] = tokenizer.nextTokenAsDouble();
524			pos[2] = tokenizer.nextTokenAsDouble();
525			if (needPos_) {
526			integrableObject->setPos(pos);
527			}
528
529			vel[0] = tokenizer.nextTokenAsDouble();
530			vel[1] = tokenizer.nextTokenAsDouble();
531			vel[2] = tokenizer.nextTokenAsDouble();
532			if (needVel_) {
533			integrableObject->setVel(vel);
534			}
535
536			if (integrableObject->isDirectional()) {
537
538			q[0] = tokenizer.nextTokenAsDouble();
539			q[1] = tokenizer.nextTokenAsDouble();
540			q[2] = tokenizer.nextTokenAsDouble();
541			q[3] = tokenizer.nextTokenAsDouble();
542
543	tim	963	RealType qlen = q.length();
544	chrisfen	721	if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
545
546			sprintf(painCave.errMsg,
547	chrisfen	984	"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
548	chrisfen	721	painCave.isFatal = 1;
549			simError();
550
551			}
552
553			q.normalize();
554			if (needQuaternion_) {
555			integrableObject->setQ(q);
556			}
557
558			ji[0] = tokenizer.nextTokenAsDouble();
559			ji[1] = tokenizer.nextTokenAsDouble();
560			ji[2] = tokenizer.nextTokenAsDouble();
561			if (needAngMom_) {
562			integrableObject->setJ(ji);
563			}
564			}
565
566			}
567
568
569			void DumpReader::parseCommentLine(char* line, Snapshot* s) {
570	tim	963	RealType currTime;
571	chrisfen	721	Mat3x3d hmat;
572	tim	963	RealType chi;
573			RealType integralOfChiDt;
574	chrisfen	721	Mat3x3d eta;
575
576			StringTokenizer tokenizer(line);
577			int nTokens;
578
579			nTokens = tokenizer.countTokens();
580
581			//comment line should at least contain 10 tokens: current time(1 token) and h-matrix(9 tokens)
582			if (nTokens < 10) {
583			sprintf(painCave.errMsg,
584			"DumpReader Error: Not enough tokens in comment line: %s", line);
585			painCave.isFatal = 1;
586			simError();
587			}
588
589			//read current time
590			currTime = tokenizer.nextTokenAsDouble();
591			s->setTime(currTime);
592
593			//read h-matrix
594			hmat(0, 0) = tokenizer.nextTokenAsDouble();
595			hmat(0, 1) = tokenizer.nextTokenAsDouble();
596			hmat(0, 2) = tokenizer.nextTokenAsDouble();
597			hmat(1, 0) = tokenizer.nextTokenAsDouble();
598			hmat(1, 1) = tokenizer.nextTokenAsDouble();
599			hmat(1, 2) = tokenizer.nextTokenAsDouble();
600			hmat(2, 0) = tokenizer.nextTokenAsDouble();
601			hmat(2, 1) = tokenizer.nextTokenAsDouble();
602			hmat(2, 2) = tokenizer.nextTokenAsDouble();
603			s->setHmat(hmat);
604
605			//read chi and integralOfChidt, they should apprear in pair
606			if (tokenizer.countTokens() >= 2) {
607			chi = tokenizer.nextTokenAsDouble();
608			integralOfChiDt = tokenizer.nextTokenAsDouble();
609
610			s->setChi(chi);
611			s->setIntegralOfChiDt(integralOfChiDt);
612			}
613
614			//read eta (eta is 3x3 matrix)
615			if (tokenizer.countTokens() >= 9) {
616			eta(0, 0) = tokenizer.nextTokenAsDouble();
617			eta(0, 1) = tokenizer.nextTokenAsDouble();
618			eta(0, 2) = tokenizer.nextTokenAsDouble();
619			eta(1, 0) = tokenizer.nextTokenAsDouble();
620			eta(1, 1) = tokenizer.nextTokenAsDouble();
621			eta(1, 2) = tokenizer.nextTokenAsDouble();
622			eta(2, 0) = tokenizer.nextTokenAsDouble();
623			eta(2, 1) = tokenizer.nextTokenAsDouble();
624			eta(2, 2) = tokenizer.nextTokenAsDouble();
625
626			s->setEta(eta);
627			}
628
629
630			}
631
632			}//end namespace oopse