ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/io/DumpReader.cpp
Revision: 1784
Committed: Wed Nov 24 22:12:12 2004 UTC (19 years, 9 months ago) by tim
File size: 16424 byte(s)
Log Message:
more get built

File Contents

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