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_ = fopen(filename_.c_str(), "r");
      
        if (inFile_ == NULL) {
          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
      
      int error;
      error = fclose(inFile_);
      
      if (error) {
        sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
        painCave.isFatal = 1;            
        simError();
      }
      
      MemoryUtils::deletePointers(framePos_);
      
#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];
    fpos_t * currPos;
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      
      rewind(inFile_);
      
      currPos = new fpos_t;
      fgetpos(inFile_, currPos);
      fgets(readBuffer, sizeof(readBuffer), inFile_);
      lineNum++;
      
      if (feof(inFile_)) {
        sprintf(painCave.errMsg,
                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
                filename_.c_str(),
                lineNum);
        painCave.isFatal = 1;
        simError();
      }
      
      while (!feof(inFile_)) {
        framePos_.push_back(currPos);
        
        i = atoi(readBuffer);
        
        fgets(readBuffer, sizeof(readBuffer), inFile_);
        lineNum++;
        
        if (feof(inFile_)) {
          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++) {
          fgets(readBuffer, sizeof(readBuffer), inFile_);
          lineNum++;
          
          if (feof(inFile_)) {
            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 = new fpos_t;
        fgetpos(inFile_, currPos);
        fgets(readBuffer, sizeof(readBuffer), inFile_);
        lineNum++;
      }
      
      delete currPos;
      rewind(inFile_);
      
      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
    
    fsetpos(inFile_, framePos_[whichFrame]);
    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
    
    if (eof_test == NULL) {
      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
    
    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
    
    if (eof_test == NULL) {
      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)) {           
        
        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
        
        if (eof_test == NULL) {
          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) {
      fsetpos(inFile_, framePos_[whichFrame]);
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
      
      if (eof_test == NULL) {
        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
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
      
      if (eof_test == NULL) {
        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)){
            
            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
            
            if (eof_test == NULL) {
              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++) {
            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
            
            if (eof_test == NULL) {
              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();
      
      double 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) {
    double currTime;
    Mat3x3d hmat;
    double chi;
    double 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:	601
Committed:	Thu Sep 15 19:17:04 2005 UTC (19 years, 7 months ago) by tim
File size:	19198 byte(s)
Log Message:	Fix a bug in DumpReader in case readFrame is called without calling getNFrames
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* 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	gezelter	2	#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	gezelter	246	#include "io/DumpReader.hpp"
56			#include "primitives/Molecule.hpp"
57	tim	3	#include "utils/simError.h"
58	gezelter	246	#include "utils/MemoryUtils.hpp"
59			#include "utils/StringTokenizer.hpp"
60	gezelter	2
61			#ifdef IS_MPI
62	gezelter	246
63	gezelter	2	#include <mpi.h>
64			#define TAKE_THIS_TAG_CHAR 0
65			#define TAKE_THIS_TAG_INT 1
66	gezelter	246
67	gezelter	2	#endif // is_mpi
68
69
70	gezelter	246	namespace oopse {
71	chrisfen	417
72			DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73	gezelter	507	: info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74	chrisfen	417
75	gezelter	2	#ifdef IS_MPI
76	chrisfen	417
77	gezelter	507	if (worldRank == 0) {
78	gezelter	2	#endif
79	chrisfen	417
80	gezelter	507	inFile_ = fopen(filename_.c_str(), "r");
81	chrisfen	417
82	gezelter	507	if (inFile_ == NULL) {
83			sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84			painCave.isFatal = 1;
85			simError();
86			}
87	chrisfen	417
88	gezelter	2	#ifdef IS_MPI
89	chrisfen	417
90	gezelter	507	}
91	chrisfen	417
92	gezelter	507	strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93			MPIcheckPoint();
94	chrisfen	417
95	gezelter	2	#endif
96	chrisfen	417
97	gezelter	507	return;
98			}
99	chrisfen	417
100			DumpReader::~DumpReader() {
101
102	gezelter	246	#ifdef IS_MPI
103	chrisfen	417
104	gezelter	246	if (worldRank == 0) {
105			#endif
106	chrisfen	417
107			int error;
108			error = fclose(inFile_);
109
110			if (error) {
111			sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
112			painCave.isFatal = 1;
113			simError();
114			}
115
116			MemoryUtils::deletePointers(framePos_);
117
118	gezelter	2	#ifdef IS_MPI
119	chrisfen	417
120	gezelter	2	}
121	chrisfen	417
122	gezelter	246	strcpy(checkPointMsg, "Dump file closed successfully.");
123			MPIcheckPoint();
124	chrisfen	417
125	gezelter	246	#endif
126	chrisfen	417
127	gezelter	246	return;
128	chrisfen	417	}
129
130			int DumpReader::getNFrames(void) {
131
132	gezelter	246	if (!isScanned_)
133	chrisfen	417	scanFile();
134
135	gezelter	246	return nframes_;
136	chrisfen	417	}
137
138			void DumpReader::scanFile(void) {
139			int i, j;
140			int lineNum = 0;
141			char readBuffer[maxBufferSize];
142			fpos_t * currPos;
143
144	gezelter	2	#ifdef IS_MPI
145	chrisfen	417
146	gezelter	246	if (worldRank == 0) {
147			#endif // is_mpi
148	chrisfen	417
149			rewind(inFile_);
150
151			currPos = new fpos_t;
152			fgetpos(inFile_, currPos);
153			fgets(readBuffer, sizeof(readBuffer), inFile_);
154			lineNum++;
155
156			if (feof(inFile_)) {
157			sprintf(painCave.errMsg,
158			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
159			filename_.c_str(),
160			lineNum);
161			painCave.isFatal = 1;
162			simError();
163			}
164
165			while (!feof(inFile_)) {
166			framePos_.push_back(currPos);
167
168			i = atoi(readBuffer);
169
170	gezelter	246	fgets(readBuffer, sizeof(readBuffer), inFile_);
171			lineNum++;
172	chrisfen	417
173	gezelter	246	if (feof(inFile_)) {
174	chrisfen	417	sprintf(painCave.errMsg,
175			"DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
176			filename_.c_str(),
177			lineNum);
178			painCave.isFatal = 1;
179			simError();
180			}
181
182			for(j = 0; j < i; j++) {
183			fgets(readBuffer, sizeof(readBuffer), inFile_);
184			lineNum++;
185
186			if (feof(inFile_)) {
187	gezelter	246	sprintf(painCave.errMsg,
188	chrisfen	417	"DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
189			" with atom %d\n", filename_.c_str(),
190			lineNum,
191			j);
192
193	gezelter	246	painCave.isFatal = 1;
194			simError();
195	chrisfen	417	}
196	gezelter	246	}
197
198	chrisfen	417	currPos = new fpos_t;
199			fgetpos(inFile_, currPos);
200			fgets(readBuffer, sizeof(readBuffer), inFile_);
201			lineNum++;
202			}
203
204			delete currPos;
205			rewind(inFile_);
206
207			nframes_ = framePos_.size();
208	gezelter	246	#ifdef IS_MPI
209	gezelter	2	}
210	chrisfen	417
211	gezelter	246	MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
212	chrisfen	417
213	gezelter	246	strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
214			MPIcheckPoint();
215	chrisfen	417
216	gezelter	246	#endif // is_mpi
217	chrisfen	417
218	gezelter	246	isScanned_ = true;
219	chrisfen	417	}
220
221			void DumpReader::readFrame(int whichFrame) {
222	tim	601	if (!isScanned_)
223			scanFile();
224
225	tim	318	int storageLayout = info_->getSnapshotManager()->getStorageLayout();
226	chrisfen	417
227	tim	318	if (storageLayout & DataStorage::dslPosition) {
228	chrisfen	417	needPos_ = true;
229	tim	318	} else {
230	chrisfen	417	needPos_ = false;
231	tim	318	}
232	chrisfen	417
233	tim	318	if (storageLayout & DataStorage::dslVelocity) {
234	chrisfen	417	needVel_ = true;
235	tim	318	} else {
236	chrisfen	417	needVel_ = false;
237	tim	318	}
238	chrisfen	417
239	tim	318	if (storageLayout & DataStorage::dslAmat \|\| storageLayout & DataStorage::dslElectroFrame) {
240	chrisfen	417	needQuaternion_ = true;
241	tim	318	} else {
242	chrisfen	417	needQuaternion_ = false;
243	tim	318	}
244	chrisfen	417
245	tim	318	if (storageLayout & DataStorage::dslAngularMomentum) {
246	chrisfen	417	needAngMom_ = true;
247	tim	318	} else {
248	chrisfen	417	needAngMom_ = false;
249	tim	318	}
250	chrisfen	417
251	gezelter	246	readSet(whichFrame);
252	chrisfen	417	}
253
254			void DumpReader::readSet(int whichFrame) {
255			int i;
256			int nTotObjs; // the number of atoms
257			char read_buffer[maxBufferSize]; //the line buffer for reading
258			char * eof_test; // ptr to see when we reach the end of the file
259
260			Molecule* mol;
261			StuntDouble* integrableObject;
262			SimInfo::MoleculeIterator mi;
263			Molecule::IntegrableObjectIterator ii;
264
265	gezelter	246	#ifndef IS_MPI
266	chrisfen	417
267	gezelter	246	fsetpos(inFile_, framePos_[whichFrame]);
268			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
269	chrisfen	417
270	gezelter	246	if (eof_test == NULL) {
271	chrisfen	417	sprintf(painCave.errMsg,
272			"DumpReader error: error reading 1st line of \"%s\"\n",
273			filename_.c_str());
274			painCave.isFatal = 1;
275			simError();
276	gezelter	2	}
277	chrisfen	417
278	gezelter	246	nTotObjs = atoi(read_buffer);
279	chrisfen	417
280	gezelter	246	if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
281	chrisfen	417	sprintf(painCave.errMsg,
282			"DumpReader error. %s nIntegrable, %d, "
283			"does not match the meta-data file's nIntegrable, %d.\n",
284			filename_.c_str(),
285			nTotObjs,
286			info_->getNGlobalIntegrableObjects());
287
288			painCave.isFatal = 1;
289			simError();
290	gezelter	2	}
291	chrisfen	417
292	gezelter	246	//read the box mat from the comment line
293	chrisfen	417
294	gezelter	246	eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
295	chrisfen	417
296	gezelter	246	if (eof_test == NULL) {
297	chrisfen	417	sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
298			filename_.c_str());
299			painCave.isFatal = 1;
300			simError();
301	gezelter	246	}
302	chrisfen	417
303	gezelter	246	parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
304	chrisfen	417
305	gezelter	246	//parse dump lines
306	chrisfen	417
307			i = 0;
308	gezelter	246	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
309	chrisfen	417
310			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
311			integrableObject = mol->nextIntegrableObject(ii)) {
312
313			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
314
315			if (eof_test == NULL) {
316			sprintf(painCave.errMsg,
317			"DumpReader Error: error in reading file %s\n"
318			"natoms = %d; index = %d\n"
319			"error reading the line from the file.\n",
320			filename_.c_str(),
321			nTotObjs,
322			i);
323
324			painCave.isFatal = 1;
325			simError();
326			}
327
328			parseDumpLine(read_buffer, integrableObject);
329			i++;
330			}
331	gezelter	2	}
332	chrisfen	417
333	gezelter	246	// MPI Section of code..........
334	chrisfen	417
335	gezelter	246	#else //IS_MPI
336	chrisfen	417
337	gezelter	246	// first thing first, suspend fatalities.
338			int masterNode = 0;
339			int nCurObj;
340			painCave.isEventLoop = 1;
341	chrisfen	417
342	gezelter	246	int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
343			int haveError;
344	chrisfen	417
345	gezelter	246	MPI_Status istatus;
346			int nitems;
347	chrisfen	417
348	gezelter	246	nTotObjs = info_->getNGlobalIntegrableObjects();
349			haveError = 0;
350	chrisfen	417
351	gezelter	246	if (worldRank == masterNode) {
352	chrisfen	417	fsetpos(inFile_, framePos_[whichFrame]);
353
354			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
355
356			if (eof_test == NULL) {
357			sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
358			filename_.c_str());
359			painCave.isFatal = 1;
360			simError();
361			}
362
363			nitems = atoi(read_buffer);
364
365			// Check to see that the number of integrable objects in the
366			// intial configuration file is the same as derived from the
367			// meta-data file.
368
369			if (nTotObjs != nitems) {
370			sprintf(painCave.errMsg,
371			"DumpReader Error. %s nIntegrable, %d, "
372			"does not match the meta-data file's nIntegrable, %d.\n",
373			filename_.c_str(),
374			nTotObjs,
375			info_->getNGlobalIntegrableObjects());
376
377			painCave.isFatal = 1;
378			simError();
379			}
380
381			//read the boxMat from the comment line
382
383			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
384
385			if (eof_test == NULL) {
386			sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
387			filename_.c_str());
388			painCave.isFatal = 1;
389			simError();
390			}
391
392			//Every single processor will parse the comment line by itself
393			//By using this way, we might lose some efficiency, but if we want to add
394			//more parameters into comment line, we only need to modify function
395			//parseCommentLine
396
397			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
398			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
399
400			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
401			int which_node = info_->getMolToProc(i);
402
403			if (which_node == masterNode) {
404			//molecules belong to master node
405
406			mol = info_->getMoleculeByGlobalIndex(i);
407
408			if (mol == NULL) {
409			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
410	gezelter	246	painCave.isFatal = 1;
411			simError();
412	chrisfen	417	}
413
414			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
415			integrableObject = mol->nextIntegrableObject(ii)){
416
417			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
418
419			if (eof_test == NULL) {
420			sprintf(painCave.errMsg,
421			"DumpReader Error: error in reading file %s\n"
422			"natoms = %d; index = %d\n"
423			"error reading the line from the file.\n",
424			filename_.c_str(),
425			nTotObjs,
426			i);
427
428			painCave.isFatal = 1;
429			simError();
430			}
431
432			parseDumpLine(read_buffer, integrableObject);
433			}
434			} else {
435			//molecule belongs to slave nodes
436
437			MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
438			MPI_COMM_WORLD, &istatus);
439
440			for(int j = 0; j < nCurObj; j++) {
441			eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
442
443			if (eof_test == NULL) {
444			sprintf(painCave.errMsg,
445			"DumpReader Error: error in reading file %s\n"
446			"natoms = %d; index = %d\n"
447			"error reading the line from the file.\n",
448			filename_.c_str(),
449			nTotObjs,
450			i);
451
452			painCave.isFatal = 1;
453			simError();
454			}
455
456			MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
457			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
458			}
459	gezelter	2	}
460	chrisfen	417	}
461			} else {
462			//actions taken at slave nodes
463			MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
464
465			/*@todo/
466			parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
467
468			for(i = 0; i < info_->getNGlobalMolecules(); i++) {
469			int which_node = info_->getMolToProc(i);
470
471			if (which_node == worldRank) {
472			//molecule with global index i belongs to this processor
473
474			mol = info_->getMoleculeByGlobalIndex(i);
475			if (mol == NULL) {
476			sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
477	gezelter	246	painCave.isFatal = 1;
478			simError();
479	chrisfen	417	}
480
481			nCurObj = mol->getNIntegrableObjects();
482
483			MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
484			MPI_COMM_WORLD);
485
486			for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
487			integrableObject = mol->nextIntegrableObject(ii)){
488	gezelter	246
489	chrisfen	417	MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
490			TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
491
492			parseDumpLine(read_buffer, integrableObject);
493			}
494
495	gezelter	246	}
496
497	chrisfen	417	}
498
499	gezelter	2	}
500	chrisfen	417
501	gezelter	246	#endif
502	chrisfen	417
503			}
504
505			void DumpReader::parseDumpLine(char line, StuntDouble integrableObject) {
506
507	gezelter	246	Vector3d pos; // position place holders
508			Vector3d vel; // velocity placeholders
509			Quat4d q; // the quaternions
510			Vector3d ji; // angular velocity placeholders;
511			StringTokenizer tokenizer(line);
512			int nTokens;
513
514			nTokens = tokenizer.countTokens();
515	chrisfen	417
516	gezelter	246	if (nTokens < 14) {
517	chrisfen	417	sprintf(painCave.errMsg,
518			"DumpReader Error: Not enough Tokens.\n%s\n", line);
519			painCave.isFatal = 1;
520			simError();
521	gezelter	2	}
522	chrisfen	417
523	gezelter	246	std::string name = tokenizer.nextToken();
524	chrisfen	417
525	gezelter	246	if (name != integrableObject->getType()) {
526	chrisfen	417
527			sprintf(painCave.errMsg,
528			"DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
529			name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
530			painCave.isFatal = 1;
531			simError();
532	gezelter	2	}
533	chrisfen	417
534	gezelter	246	pos[0] = tokenizer.nextTokenAsDouble();
535			pos[1] = tokenizer.nextTokenAsDouble();
536			pos[2] = tokenizer.nextTokenAsDouble();
537	tim	318	if (needPos_) {
538	chrisfen	417	integrableObject->setPos(pos);
539	tim	318	}
540	gezelter	246
541			vel[0] = tokenizer.nextTokenAsDouble();
542			vel[1] = tokenizer.nextTokenAsDouble();
543			vel[2] = tokenizer.nextTokenAsDouble();
544	tim	318	if (needVel_) {
545	chrisfen	417	integrableObject->setVel(vel);
546	tim	318	}
547
548	gezelter	246	if (integrableObject->isDirectional()) {
549	chrisfen	417
550			q[0] = tokenizer.nextTokenAsDouble();
551			q[1] = tokenizer.nextTokenAsDouble();
552			q[2] = tokenizer.nextTokenAsDouble();
553			q[3] = tokenizer.nextTokenAsDouble();
554
555			double qlen = q.length();
556			if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
557	gezelter	246
558	chrisfen	417	sprintf(painCave.errMsg,
559			"DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
560			painCave.isFatal = 1;
561			simError();
562
563			}
564
565			q.normalize();
566			if (needQuaternion_) {
567			integrableObject->setQ(q);
568			}
569
570			ji[0] = tokenizer.nextTokenAsDouble();
571			ji[1] = tokenizer.nextTokenAsDouble();
572			ji[2] = tokenizer.nextTokenAsDouble();
573			if (needAngMom_) {
574			integrableObject->setJ(ji);
575			}
576	gezelter	2	}
577	chrisfen	417
578			}
579
580
581			void DumpReader::parseCommentLine(char* line, Snapshot* s) {
582	gezelter	246	double currTime;
583			Mat3x3d hmat;
584			double chi;
585			double integralOfChiDt;
586			Mat3x3d eta;
587	chrisfen	417
588	gezelter	246	StringTokenizer tokenizer(line);
589			int nTokens;
590	chrisfen	417
591	gezelter	246	nTokens = tokenizer.countTokens();
592	chrisfen	417
593	gezelter	246	//comment line should at least contain 10 tokens: current time(1 token) and h-matrix(9 tokens)
594			if (nTokens < 10) {
595	chrisfen	417	sprintf(painCave.errMsg,
596			"DumpReader Error: Not enough tokens in comment line: %s", line);
597			painCave.isFatal = 1;
598			simError();
599	gezelter	2	}
600	chrisfen	417
601	gezelter	246	//read current time
602			currTime = tokenizer.nextTokenAsDouble();
603			s->setTime(currTime);
604
605			//read h-matrix
606			hmat(0, 0) = tokenizer.nextTokenAsDouble();
607			hmat(0, 1) = tokenizer.nextTokenAsDouble();
608			hmat(0, 2) = tokenizer.nextTokenAsDouble();
609			hmat(1, 0) = tokenizer.nextTokenAsDouble();
610			hmat(1, 1) = tokenizer.nextTokenAsDouble();
611			hmat(1, 2) = tokenizer.nextTokenAsDouble();
612			hmat(2, 0) = tokenizer.nextTokenAsDouble();
613			hmat(2, 1) = tokenizer.nextTokenAsDouble();
614			hmat(2, 2) = tokenizer.nextTokenAsDouble();
615			s->setHmat(hmat);
616
617	tim	318	//read chi and integralOfChidt, they should apprear in pair
618	gezelter	246	if (tokenizer.countTokens() >= 2) {
619	chrisfen	417	chi = tokenizer.nextTokenAsDouble();
620			integralOfChiDt = tokenizer.nextTokenAsDouble();
621
622			s->setChi(chi);
623			s->setIntegralOfChiDt(integralOfChiDt);
624	gezelter	2	}
625
626	gezelter	246	//read eta (eta is 3x3 matrix)
627			if (tokenizer.countTokens() >= 9) {
628	chrisfen	417	eta(0, 0) = tokenizer.nextTokenAsDouble();
629			eta(0, 1) = tokenizer.nextTokenAsDouble();
630			eta(0, 2) = tokenizer.nextTokenAsDouble();
631			eta(1, 0) = tokenizer.nextTokenAsDouble();
632			eta(1, 1) = tokenizer.nextTokenAsDouble();
633			eta(1, 2) = tokenizer.nextTokenAsDouble();
634			eta(2, 0) = tokenizer.nextTokenAsDouble();
635			eta(2, 1) = tokenizer.nextTokenAsDouble();
636			eta(2, 2) = tokenizer.nextTokenAsDouble();
637
638			s->setEta(eta);
639	gezelter	2	}
640
641	chrisfen	417
642			}
643
644	gezelter	246	}//end namespace oopse