ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/io/DumpReader.cpp
Revision: 1808
Committed: Tue Nov 30 23:14:29 2004 UTC (19 years, 9 months ago) by tim
File size: 16387 byte(s)
Log Message:
io get built, next is integrator

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