ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpReader.cpp
Revision: 1249
Committed: Mon Jun 2 02:42:50 2008 UTC (16 years, 11 months ago) by gezelter
File size: 16222 byte(s)
Log Message:
Fixed a bug in BlockSnapshotManager.

File Contents

# 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 gezelter 1104 : info_(info), filename_(filename), isScanned_(false), nframes_(0), needCOMprops_(false) {
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 gezelter 1241 errorCheckPoint();
96 chrisfen 996
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 gezelter 1241 errorCheckPoint();
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 tim 1024 int lineNo = 0;
133     std::streampos prevPos;
134 chrisfen 721 std::streampos currPos;
135 tim 1024
136 chrisfen 721 #ifdef IS_MPI
137 tim 1024
138 chrisfen 721 if (worldRank == 0) {
139     #endif // is_mpi
140 tim 1024
141     currPos = inFile_->tellg();
142     prevPos = currPos;
143     bool foundOpenSnapshotTag = false;
144     bool foundClosedSnapshotTag = false;
145     while(inFile_->getline(buffer, bufferSize)) {
146     ++lineNo;
147    
148     std::string line = buffer;
149     currPos = inFile_->tellg();
150     if (line.find("<Snapshot>")!= std::string::npos) {
151     if (foundOpenSnapshotTag) {
152 chrisfen 721 sprintf(painCave.errMsg,
153 tim 1024 "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
154     filename_.c_str());
155 chrisfen 721 painCave.isFatal = 1;
156 tim 1024 simError();
157     }
158     foundOpenSnapshotTag = true;
159     foundClosedSnapshotTag = false;
160     framePos_.push_back(prevPos);
161    
162     } else if (line.find("</Snapshot>") != std::string::npos){
163     if (!foundOpenSnapshotTag) {
164     sprintf(painCave.errMsg,
165     "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
166     filename_.c_str());
167     painCave.isFatal = 1;
168 chrisfen 721 simError();
169 tim 1024 }
170    
171     if (foundClosedSnapshotTag) {
172     sprintf(painCave.errMsg,
173     "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
174     filename_.c_str());
175     painCave.isFatal = 1;
176     simError();
177     }
178     foundClosedSnapshotTag = true;
179     foundOpenSnapshotTag = false;
180     }
181     prevPos = currPos;
182     }
183    
184     // only found <Snapshot> for the last frame means the file is corrupted, we should discard
185     // it and give a warning message
186     if (foundOpenSnapshotTag) {
187     sprintf(painCave.errMsg,
188     "DumpReader: last frame in %s is invalid\n", filename_.c_str());
189     painCave.isFatal = 0;
190     simError();
191     framePos_.pop_back();
192     }
193    
194 chrisfen 721 nframes_ = framePos_.size();
195 tim 1024
196     if (nframes_ == 0) {
197     sprintf(painCave.errMsg,
198     "DumpReader: %s does not contain a valid frame\n", filename_.c_str());
199     painCave.isFatal = 1;
200     simError();
201     }
202 chrisfen 721 #ifdef IS_MPI
203     }
204    
205     MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
206 tim 1024
207 chrisfen 721 #endif // is_mpi
208 tim 1024
209 chrisfen 721 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 gezelter 1104
244     if (needCOMprops_) {
245     Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
246     Vector3d com;
247     Vector3d comvel;
248     Vector3d comw;
249 chuckv 1112 if (needPos_ && needVel_){
250     info_->getComAll(com, comvel);
251     comw = info_->getAngularMomentum();
252     }else{
253     com = info_->getCom();
254     comvel = 0.0;
255     comw = 0.0;
256     }
257 gezelter 1104 s->setCOMprops(com, comvel, comw);
258     }
259    
260 chrisfen 721 }
261    
262 tim 1024 void DumpReader::readSet(int whichFrame) {
263     std::string line;
264    
265 chrisfen 721 #ifndef IS_MPI
266     inFile_->clear();
267     inFile_->seekg(framePos_[whichFrame]);
268 tim 1024
269     std::istream& inputStream = *inFile_;
270    
271     #else
272     int masterNode = 0;
273     std::stringstream sstream;
274     if (worldRank == masterNode) {
275     std::string sendBuffer;
276    
277     inFile_->clear();
278     inFile_->seekg(framePos_[whichFrame]);
279    
280     while (inFile_->getline(buffer, bufferSize)) {
281    
282     line = buffer;
283     sendBuffer += line;
284     sendBuffer += '\n';
285     if (line.find("</Snapshot>") != std::string::npos) {
286     break;
287     }
288     }
289    
290     int sendBufferSize = sendBuffer.size();
291     MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
292     MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
293    
294     sstream.str(sendBuffer);
295     } else {
296     int sendBufferSize;
297     MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
298     char * recvBuffer = new char[sendBufferSize+1];
299     MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
300     sstream.str(recvBuffer);
301     }
302    
303     std::istream& inputStream = sstream;
304     #endif
305    
306     inputStream.getline(buffer, bufferSize);
307    
308     line = buffer;
309     if (line.find("<Snapshot>") == std::string::npos) {
310 chrisfen 721 sprintf(painCave.errMsg,
311 tim 1024 "DumpReader Error: can not find <Snapshot>\n");
312 chrisfen 721 painCave.isFatal = 1;
313     simError();
314     }
315 tim 1024
316     //read frameData
317     readFrameProperties(inputStream);
318    
319     //read StuntDoubles
320     readStuntDoubles(inputStream);
321    
322     inputStream.getline(buffer, bufferSize);
323     line = buffer;
324     if (line.find("</Snapshot>") == std::string::npos) {
325 chrisfen 721 sprintf(painCave.errMsg,
326 tim 1024 "DumpReader Error: can not find </Snapshot>\n");
327 chrisfen 721 painCave.isFatal = 1;
328     simError();
329 tim 1024 }
330    
331 chrisfen 721 }
332    
333 tim 1024 void DumpReader::parseDumpLine(const std::string& line) {
334    
335    
336 chrisfen 721 StringTokenizer tokenizer(line);
337     int nTokens;
338    
339     nTokens = tokenizer.countTokens();
340    
341 tim 1024 if (nTokens < 2) {
342 chrisfen 721 sprintf(painCave.errMsg,
343 tim 1024 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
344 chrisfen 721 painCave.isFatal = 1;
345     simError();
346     }
347 tim 1024
348     int index = tokenizer.nextTokenAsInt();
349    
350     StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
351    
352     if (integrableObject == NULL) {
353     return;
354     }
355     std::string type = tokenizer.nextToken();
356     int size = type.size();
357 gezelter 1037
358 tim 1024 for(int i = 0; i < size; ++i) {
359     switch(type[i]) {
360    
361     case 'p': {
362     Vector3d pos;
363     pos[0] = tokenizer.nextTokenAsDouble();
364     pos[1] = tokenizer.nextTokenAsDouble();
365     pos[2] = tokenizer.nextTokenAsDouble();
366     if (needPos_) {
367     integrableObject->setPos(pos);
368     }
369     break;
370     }
371     case 'v' : {
372     Vector3d vel;
373     vel[0] = tokenizer.nextTokenAsDouble();
374     vel[1] = tokenizer.nextTokenAsDouble();
375     vel[2] = tokenizer.nextTokenAsDouble();
376     if (needVel_) {
377     integrableObject->setVel(vel);
378     }
379     break;
380     }
381    
382     case 'q' : {
383     Quat4d q;
384     if (integrableObject->isDirectional()) {
385    
386     q[0] = tokenizer.nextTokenAsDouble();
387     q[1] = tokenizer.nextTokenAsDouble();
388     q[2] = tokenizer.nextTokenAsDouble();
389     q[3] = tokenizer.nextTokenAsDouble();
390    
391     RealType qlen = q.length();
392     if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
393    
394     sprintf(painCave.errMsg,
395     "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
396     painCave.isFatal = 1;
397     simError();
398    
399     }
400    
401     q.normalize();
402     if (needQuaternion_) {
403     integrableObject->setQ(q);
404     }
405     }
406     break;
407     }
408     case 'j' : {
409     Vector3d ji;
410     if (integrableObject->isDirectional()) {
411     ji[0] = tokenizer.nextTokenAsDouble();
412     ji[1] = tokenizer.nextTokenAsDouble();
413     ji[2] = tokenizer.nextTokenAsDouble();
414     if (needAngMom_) {
415     integrableObject->setJ(ji);
416     }
417     }
418 gezelter 1037 break;
419 tim 1024 }
420     case 'f': {
421    
422     Vector3d force;
423     force[0] = tokenizer.nextTokenAsDouble();
424     force[1] = tokenizer.nextTokenAsDouble();
425     force[2] = tokenizer.nextTokenAsDouble();
426     integrableObject->setFrc(force);
427     break;
428     }
429     case 't' : {
430    
431     Vector3d torque;
432     torque[0] = tokenizer.nextTokenAsDouble();
433     torque[1] = tokenizer.nextTokenAsDouble();
434     torque[2] = tokenizer.nextTokenAsDouble();
435     integrableObject->setTrq(torque);
436     break;
437     }
438     default: {
439     sprintf(painCave.errMsg,
440     "DumpReader Error: %s is an unrecognized type\n", type.c_str());
441     painCave.isFatal = 1;
442     simError();
443     break;
444     }
445    
446     }
447     }
448 chrisfen 721
449 tim 1024 }
450    
451    
452     void DumpReader::readStuntDoubles(std::istream& inputStream) {
453    
454     inputStream.getline(buffer, bufferSize);
455     std::string line(buffer);
456    
457     if (line.find("<StuntDoubles>") == std::string::npos) {
458 chrisfen 721 sprintf(painCave.errMsg,
459 tim 1024 "DumpReader Error: Missing <StuntDoubles>\n");
460 chrisfen 721 painCave.isFatal = 1;
461 tim 1024 simError();
462     }
463    
464     while(inputStream.getline(buffer, bufferSize)) {
465     line = buffer;
466    
467     if(line.find("</StuntDoubles>") != std::string::npos) {
468     break;
469     }
470    
471     parseDumpLine(line);
472     }
473    
474     }
475    
476     void DumpReader::readFrameProperties(std::istream& inputStream) {
477    
478     Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
479     inputStream.getline(buffer, bufferSize);
480     std::string line(buffer);
481    
482     if (line.find("<FrameData>") == std::string::npos) {
483     sprintf(painCave.errMsg,
484     "DumpReader Error: Missing <FrameData>\n");
485     painCave.isFatal = 1;
486     simError();
487     }
488    
489     while(inputStream.getline(buffer, bufferSize)) {
490     line = buffer;
491    
492     if(line.find("</FrameData>") != std::string::npos) {
493     break;
494     }
495    
496     StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
497     if (!tokenizer.hasMoreTokens()) {
498 chrisfen 721 sprintf(painCave.errMsg,
499 tim 1024 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
500 chrisfen 721 painCave.isFatal = 1;
501 tim 1024 simError();
502     }
503    
504     std::string propertyName = tokenizer.nextToken();
505     if (propertyName == "Time") {
506     RealType currTime = tokenizer.nextTokenAsDouble();
507     s->setTime(currTime);
508     } else if (propertyName == "Hmat"){
509     Mat3x3d hmat;
510     hmat(0, 0) = tokenizer.nextTokenAsDouble();
511     hmat(0, 1) = tokenizer.nextTokenAsDouble();
512     hmat(0, 2) = tokenizer.nextTokenAsDouble();
513     hmat(1, 0) = tokenizer.nextTokenAsDouble();
514     hmat(1, 1) = tokenizer.nextTokenAsDouble();
515     hmat(1, 2) = tokenizer.nextTokenAsDouble();
516     hmat(2, 0) = tokenizer.nextTokenAsDouble();
517     hmat(2, 1) = tokenizer.nextTokenAsDouble();
518     hmat(2, 2) = tokenizer.nextTokenAsDouble();
519     s->setHmat(hmat);
520     } else if (propertyName == "Thermostat") {
521     RealType chi = tokenizer.nextTokenAsDouble();
522     RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
523     s->setChi(chi);
524     s->setIntegralOfChiDt(integralOfChiDt);
525     } else if (propertyName == "Barostat") {
526     Mat3x3d eta;
527     eta(0, 0) = tokenizer.nextTokenAsDouble();
528     eta(0, 1) = tokenizer.nextTokenAsDouble();
529     eta(0, 2) = tokenizer.nextTokenAsDouble();
530     eta(1, 0) = tokenizer.nextTokenAsDouble();
531     eta(1, 1) = tokenizer.nextTokenAsDouble();
532     eta(1, 2) = tokenizer.nextTokenAsDouble();
533     eta(2, 0) = tokenizer.nextTokenAsDouble();
534     eta(2, 1) = tokenizer.nextTokenAsDouble();
535     eta(2, 2) = tokenizer.nextTokenAsDouble();
536     s->setEta(eta);
537     } else {
538     sprintf(painCave.errMsg,
539     "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
540     painCave.isFatal = 0;
541     simError();
542     }
543    
544     }
545    
546     }
547    
548 chrisfen 721
549     }//end namespace oopse