ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/oopse-1.0/libmdtools/DumpReader.cpp
Revision: 1447
Committed: Fri Jul 30 21:01:35 2004 UTC (19 years, 11 months ago) by gezelter
File size: 18157 byte(s)
Log Message:
Initial import of OOPSE sources into cvs tree

File Contents

# Content
1 #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
15 #include "ReadWrite.hpp"
16 #include "simError.h"
17
18 #ifdef IS_MPI
19 #include <mpi.h>
20 #include "mpiSimulation.hpp"
21 #define TAKE_THIS_TAG_CHAR 0
22 #define TAKE_THIS_TAG_INT 1
23 #endif // is_mpi
24
25
26 DumpReader :: DumpReader(const char *in_name ){
27
28 isScanned = false;
29
30 #ifdef IS_MPI
31 if (worldRank == 0) {
32 #endif
33
34 inFile = fopen(in_name, "r");
35 if(inFile == NULL){
36 sprintf(painCave.errMsg,
37 "Cannot open file: %s\n", in_name);
38 painCave.isFatal = 1;
39 simError();
40 }
41
42 inFileName = in_name;
43 #ifdef IS_MPI
44 }
45 strcpy( checkPointMsg, "Dump file opened for reading successfully." );
46 MPIcheckPoint();
47 #endif
48 return;
49 }
50
51 DumpReader :: ~DumpReader( ){
52 #ifdef IS_MPI
53 if (worldRank == 0) {
54 #endif
55 vector<fpos_t*>::iterator i;
56
57 int error;
58 error = fclose( inFile );
59 if( error ){
60 sprintf( painCave.errMsg,
61 "Error closing %s\n", inFileName.c_str());
62 simError();
63 }
64
65 for(i = framePos.begin(); i != framePos.end(); ++i)
66 delete *i;
67 framePos.clear();
68
69 #ifdef IS_MPI
70 }
71 strcpy( checkPointMsg, "Dump file closed successfully." );
72 MPIcheckPoint();
73 #endif
74
75 return;
76 }
77
78 int DumpReader::getNframes( void ){
79
80 if( !isScanned )
81 scanFile();
82 return framePos.size();
83 }
84
85 void DumpReader::scanFile( void ){
86
87 int i, j;
88 int lineNum = 0;
89 char readBuffer[2000];
90 fpos_t *currPos;
91
92 #ifdef IS_MPI
93 if( worldRank == 0 ){
94 #endif // is_mpi
95
96 rewind( inFile );
97
98 currPos = new fpos_t;
99 fgetpos( inFile, currPos );
100 fgets( readBuffer, sizeof( readBuffer ), inFile );
101 lineNum++;
102 if( feof( inFile ) ){
103 sprintf( painCave.errMsg,
104 "File \"%s\" ended unexpectedly at line %d\n",
105 inFileName.c_str(),
106 lineNum );
107 painCave.isFatal = 1;
108 simError();
109 }
110
111 while( !feof( inFile ) ){
112
113 framePos.push_back(currPos);
114
115 i = atoi(readBuffer);
116
117 fgets( readBuffer, sizeof( readBuffer ), inFile );
118 lineNum++;
119 if( feof( inFile ) ){
120 sprintf( painCave.errMsg,
121 "File \"%s\" ended unexpectedly at line %d\n",
122 inFileName.c_str(),
123 lineNum );
124 painCave.isFatal = 1;
125 simError();
126 }
127
128 for(j=0; j<i; j++){
129
130 fgets( readBuffer, sizeof( readBuffer ), inFile );
131 lineNum++;
132 if( feof( inFile ) ){
133 sprintf( painCave.errMsg,
134 "File \"%s\" ended unexpectedly at line %d,"
135 " with atom %d\n",
136 inFileName.c_str(),
137 lineNum,
138 j );
139 painCave.isFatal = 1;
140 simError();
141 }
142
143 }
144
145 currPos = new fpos_t;
146 fgetpos( inFile, currPos );
147 fgets( readBuffer, sizeof( readBuffer ), inFile );
148 lineNum++;
149 }
150
151 delete currPos;
152 rewind( inFile );
153
154 isScanned = true;
155
156 #ifdef IS_MPI
157 }
158 strcpy( checkPointMsg, "Successfully scanned DumpFile\n" );
159 MPIcheckPoint();
160 #endif // is_mpi
161 }
162
163 void DumpReader :: readFrame( SimInfo* the_simnfo, int whichFrame){
164
165 simnfo = the_simnfo;
166
167 this->readSet( whichFrame );
168 }
169
170
171
172 void DumpReader :: readSet( int whichFrame ){
173
174 int i;
175 unsigned int j;
176
177 #ifdef IS_MPI
178 int done, which_node, which_atom; // loop counter
179 #endif //is_mpi
180
181 const int BUFFERSIZE = 2000; // size of the read buffer
182 int nTotObjs; // the number of atoms
183 char read_buffer[BUFFERSIZE]; //the line buffer for reading
184
185 char *eof_test; // ptr to see when we reach the end of the file
186 char *parseErr;
187
188 vector<StuntDouble*> integrableObjects;
189
190
191 #ifndef IS_MPI
192
193 fsetpos(inFile, framePos[whichFrame]);
194 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
195 if( eof_test == NULL ){
196 sprintf( painCave.errMsg,
197 "DumpReader error: error reading 1st line of \"%s\"\n",
198 inFileName.c_str() );
199 painCave.isFatal = 1;
200 simError();
201 }
202
203 nTotObjs = atoi( read_buffer );
204
205 if( nTotObjs != simnfo->getTotIntegrableObjects() ){
206 sprintf( painCave.errMsg,
207 "DumpReader error. %s nIntegrable, %d, "
208 "does not match the meta-data file's nIntegrable, %d.\n",
209 inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
210 painCave.isFatal = 1;
211 simError();
212 }
213
214 //read the box mat from the comment line
215
216 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
217 if(eof_test == NULL){
218 sprintf( painCave.errMsg,
219 "error in reading commment in %s\n", inFileName.c_str());
220 painCave.isFatal = 1;
221 simError();
222 }
223
224 parseErr = parseCommentLine( read_buffer, simnfo);
225 if( parseErr != NULL ){
226 strcpy( painCave.errMsg, parseErr );
227 painCave.isFatal = 1;
228 simError();
229 }
230
231 //parse dump lines
232
233 for( i=0; i < simnfo->n_mol; i++){
234
235 integrableObjects = (simnfo->molecules[i]).getIntegrableObjects();
236
237 for(j = 0; j < integrableObjects.size(); j++){
238
239 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
240 if(eof_test == NULL){
241 sprintf(painCave.errMsg,
242 "error in reading file %s\n"
243 "natoms = %d; index = %d\n"
244 "error reading the line from the file.\n",
245 inFileName.c_str(), nTotObjs, i );
246 painCave.isFatal = 1;
247 simError();
248 }
249
250 parseErr = parseDumpLine( read_buffer, integrableObjects[j]);
251 if( parseErr != NULL ){
252 strcpy( painCave.errMsg, parseErr );
253 painCave.isFatal = 1;
254 simError();
255 }
256 }
257 }
258
259 // MPI Section of code..........
260 #else //IS_MPI
261
262 // first thing first, suspend fatalities.
263 painCave.isEventLoop = 1;
264
265 int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
266 int haveError;
267
268 MPI_Status istatus;
269 int *MolToProcMap = mpiSim->getMolToProcMap();
270 int localIndex;
271 int nCurObj;
272 int nitems;
273
274 nTotObjs = simnfo->getTotIntegrableObjects();
275 haveError = 0;
276 if (worldRank == 0) {
277 fsetpos(inFile, framePos[whichFrame]);
278
279 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
280 if( eof_test == NULL ){
281 sprintf( painCave.errMsg,
282 "Error reading 1st line of %s \n ",inFileName.c_str());
283 haveError = 1;
284 simError();
285 }
286
287 nitems = atoi( read_buffer );
288
289 // 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
293 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 inFileName.c_str(), nTotObjs, simnfo->getTotIntegrableObjects());
298 haveError= 1;
299 simError();
300 }
301
302 //read the boxMat from the comment line
303
304 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
305 if(eof_test == NULL){
306 sprintf( painCave.errMsg,
307 "error in reading commment in %s\n", inFileName.c_str());
308 haveError = 1;
309 simError();
310 }
311
312 //Every single processor will parse the comment line by itself
313 //By using this way, we might lose some efficiency, but if we want to add
314 //more parameters into comment line, we only need to modify function
315 //parseCommentLine
316
317 MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
318
319 parseErr = parseCommentLine( read_buffer, simnfo);
320
321 if( parseErr != NULL ){
322 strcpy( painCave.errMsg, parseErr );
323 haveError = 1;
324 simError();
325 }
326
327 for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
328 which_node = MolToProcMap[i];
329 if(which_node == 0){
330 //molecules belong to master node
331
332 localIndex = mpiSim->getGlobalToLocalMol(i);
333
334 if(localIndex == -1) {
335 strcpy(painCave.errMsg, "Molecule not found on node 0!");
336 haveError = 1;
337 simError();
338 }
339
340 integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
341 for(j=0; j < integrableObjects.size(); j++){
342
343 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
344 if(eof_test == NULL){
345 sprintf(painCave.errMsg,
346 "error in reading file %s\n"
347 "natoms = %d; index = %d\n"
348 "error reading the line from the file.\n",
349 inFileName.c_str(), nTotObjs, i );
350 haveError= 1;
351 simError();
352 }
353
354 if(haveError) nodeZeroError();
355
356 parseDumpLine(read_buffer, integrableObjects[j]);
357
358 }
359
360
361 }
362 else{
363 //molecule belongs to slave nodes
364
365 MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
366 TAKE_THIS_TAG_INT, MPI_COMM_WORLD, &istatus);
367
368 for(j=0; j < nCurObj; j++){
369
370 eof_test = fgets(read_buffer, sizeof(read_buffer), inFile);
371 if(eof_test == NULL){
372 sprintf(painCave.errMsg,
373 "error in reading file %s\n"
374 "natoms = %d; index = %d\n"
375 "error reading the line from the file.\n",
376 inFileName.c_str(), nTotObjs, i );
377 haveError= 1;
378 simError();
379 }
380
381 if(haveError) nodeZeroError();
382
383 MPI_Send(read_buffer, BUFFERSIZE, MPI_CHAR, which_node,
384 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
385
386 }
387
388 }
389
390 }
391
392 }
393 else{
394 //actions taken at slave nodes
395 MPI_Bcast(read_buffer, BUFFERSIZE, MPI_CHAR, 0, MPI_COMM_WORLD);
396
397 parseErr = parseCommentLine( read_buffer, simnfo);
398
399 if( parseErr != NULL ){
400 strcpy( painCave.errMsg, parseErr );
401 haveError = 1;
402 simError();
403 }
404
405 for (i=0 ; i < mpiSim->getNMolGlobal(); i++) {
406 which_node = MolToProcMap[i];
407
408 if(which_node == worldRank){
409 //molecule with global index i belongs to this processor
410
411 localIndex = mpiSim->getGlobalToLocalMol(i);
412
413 if(localIndex == -1) {
414 sprintf(painCave.errMsg, "Molecule not found on node %d\n", worldRank);
415 haveError = 1;
416 simError();
417 }
418
419 integrableObjects = (simnfo->molecules[localIndex]).getIntegrableObjects();
420
421 nCurObj = integrableObjects.size();
422
423 MPI_Send(&nCurObj, 1, MPI_INT, 0,
424 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
425
426 for(j = 0; j < integrableObjects.size(); j++){
427
428 MPI_Recv(read_buffer, BUFFERSIZE, MPI_CHAR, 0,
429 TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
430
431 parseErr = parseDumpLine(read_buffer, integrableObjects[j]);
432
433 if( parseErr != NULL ){
434 strcpy( painCave.errMsg, parseErr );
435 simError();
436 }
437
438 }
439
440 }
441
442 }
443
444 }
445
446 #endif
447 }
448
449 char* DumpReader::parseDumpLine(char* readLine, StuntDouble* sd){
450
451 char *foo; // the pointer to the current string token
452
453 double pos[3]; // position place holders
454 double vel[3]; // velocity placeholders
455 double q[4]; // the quaternions
456 double ji[3]; // angular velocity placeholders;
457 double qSqr, qLength; // needed to normalize the quaternion vector.
458
459
460 // set the string tokenizer
461
462 foo = strtok(readLine, " ,;\t");
463
464 // check the atom name to the current atom
465
466 if( strcmp( foo, sd->getType() ) ){
467 sprintf( painCave.errMsg,
468 "DumpReader error. Does not"
469 " match the meta-data atom %s.\n",
470 sd->getType() );
471 return strdup( painCave.errMsg );
472 }
473
474 // get the positions
475
476 foo = strtok(NULL, " ,;\t");
477 if(foo == NULL){
478 sprintf( painCave.errMsg,
479 "error in reading postition x from %s\n",
480 inFileName.c_str());
481 return strdup( painCave.errMsg );
482 }
483 pos[0] = atof( foo );
484
485 foo = strtok(NULL, " ,;\t");
486 if(foo == NULL){
487 sprintf( painCave.errMsg,
488 "error in reading postition y from %s\n",
489 inFileName.c_str());
490 return strdup( painCave.errMsg );
491 }
492 pos[1] = atof( foo );
493
494 foo = strtok(NULL, " ,;\t");
495 if(foo == NULL){
496 sprintf( painCave.errMsg,
497 "error in reading postition z from %s\n",
498 inFileName.c_str());
499 return strdup( painCave.errMsg );
500 }
501 pos[2] = atof( foo );
502
503
504 // get the velocities
505
506 foo = strtok(NULL, " ,;\t");
507 if(foo == NULL){
508 sprintf( painCave.errMsg,
509 "error in reading velocity x from %s\n",
510 inFileName.c_str() );
511 return strdup( painCave.errMsg );
512 }
513 vel[0] = atof( foo );
514
515 foo = strtok(NULL, " ,;\t");
516 if(foo == NULL){
517 sprintf( painCave.errMsg,
518 "error in reading velocity x from %s\n",
519 inFileName.c_str() );
520 return strdup( painCave.errMsg );
521 }
522 vel[1] = atof( foo );
523
524 foo = strtok(NULL, " ,;\t");
525 if(foo == NULL){
526 sprintf( painCave.errMsg,
527 "error in reading velocity x from %s\n",
528 inFileName.c_str() );
529 return strdup( painCave.errMsg );
530 }
531 vel[2] = atof( foo );
532
533
534 // add the positions and velocities to the atom
535
536 sd->setPos( pos );
537 sd->setVel( vel );
538
539 if (!sd->isDirectional())
540 return NULL;
541
542 // get the quaternions
543
544 if( sd->isDirectional() ){
545
546 foo = strtok(NULL, " ,;\t");
547 if(foo == NULL){
548 sprintf( painCave.errMsg,
549 "error in reading velocity x from %s\n",
550 inFileName.c_str() );
551 return strdup( painCave.errMsg );
552 }
553 q[0] = atof( foo );
554
555 foo = strtok(NULL, " ,;\t");
556 if(foo == NULL){
557 sprintf( painCave.errMsg,
558 "error in reading velocity x from %s\n",
559 inFileName.c_str() );
560 return strdup( painCave.errMsg );
561 }
562 q[1] = atof( foo );
563
564 foo = strtok(NULL, " ,;\t");
565 if(foo == NULL){
566 sprintf( painCave.errMsg,
567 "error in reading velocity x from %s\n",
568 inFileName.c_str() );
569 return strdup( painCave.errMsg );
570 }
571 q[2] = atof( foo );
572
573 foo = strtok(NULL, " ,;\t");
574 if(foo == NULL){
575 sprintf( painCave.errMsg,
576 "error in reading velocity x from %s\n",
577 inFileName.c_str() );
578 return strdup( painCave.errMsg );
579 }
580 q[3] = atof( foo );
581
582 // get the angular velocities
583
584 foo = strtok(NULL, " ,;\t");
585 if(foo == NULL){
586 sprintf( painCave.errMsg,
587 "error in reading velocity x from %s\n",
588 inFileName.c_str() );
589 return strdup( painCave.errMsg );
590 }
591 ji[0] = atof( foo );
592
593 foo = strtok(NULL, " ,;\t");
594 if(foo == NULL){
595 sprintf( painCave.errMsg,
596 "error in reading velocity x from %s\n",
597 inFileName.c_str() );
598 return strdup( painCave.errMsg );
599 }
600 ji[1] = atof(foo );
601
602 foo = strtok(NULL, " ,;\t");
603 if(foo == NULL){
604 sprintf( painCave.errMsg,
605 "error in reading velocity x from %s\n",
606 inFileName.c_str() );
607 return strdup( painCave.errMsg );
608 }
609 ji[2] = atof( foo );
610
611
612 // check that the quaternion vector is normalized
613
614 qSqr = (q[0] * q[0]) + (q[1] * q[1]) + (q[2] * q[2]) + (q[3] * q[3]);
615
616 if (fabs(qSqr) < 1e-6) {
617 sprintf(painCave.errMsg,
618 "initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
619 return strdup(painCave.errMsg);
620 }
621
622 qLength = sqrt( qSqr );
623 q[0] = q[0] / qLength;
624 q[1] = q[1] / qLength;
625 q[2] = q[2] / qLength;
626 q[3] = q[3] / qLength;
627
628 // add quaternion and angular velocities
629
630 sd->setQ( q );
631 sd->setJ( ji );
632 }
633
634
635
636 return NULL;
637 }
638
639
640 char* DumpReader::parseCommentLine(char* readLine, SimInfo* entry_plug){
641
642 double currTime;
643 double boxMat[9];
644 double theBoxMat3[3][3];
645 double chi;
646 double integralOfChidt;
647 double eta[9];
648
649 char *foo; // the pointer to the current string token
650
651 // set the string tokenizer
652
653 foo = strtok(readLine, " ,;\t");
654 // set the timeToken.
655
656 if(foo == NULL){
657 sprintf( painCave.errMsg,
658 "error in reading Time from %s\n",
659 inFileName.c_str() );
660 return strdup( painCave.errMsg );
661 }
662
663 currTime = atof( foo );
664 entry_plug->setTime( currTime );
665
666 //get H-Matrix
667
668 for(int i = 0 ; i < 9; i++){
669 foo = strtok(NULL, " ,;\t");
670 if(foo == NULL){
671 sprintf( painCave.errMsg,
672 "error in reading H[%d] from %s\n", i, inFileName.c_str() );
673 return strdup( painCave.errMsg );
674 }
675 boxMat[i] = atof( foo );
676 }
677
678 for(int i=0;i<3;i++)
679 for(int j=0;j<3;j++) theBoxMat3[i][j] = boxMat[3*j+i];
680
681 //set H-Matrix
682 entry_plug->setBoxM( theBoxMat3 );
683
684 //get chi and integralOfChidt, they should appear by pair
685
686 if( entry_plug->useInitXSstate ){
687 foo = strtok(NULL, " ,;\t\n");
688 if(foo != NULL){
689 chi = atof(foo);
690
691 foo = strtok(NULL, " ,;\t\n");
692 if(foo == NULL){
693 sprintf( painCave.errMsg,
694 "chi and integralOfChidt should appear by pair in %s\n", inFileName.c_str() );
695 return strdup( painCave.errMsg );
696 }
697 integralOfChidt = atof( foo );
698
699 //push chi and integralOfChidt into SimInfo::properties which can be
700 //retrieved by integrator later
701 DoubleData* chiValue = new DoubleData();
702 chiValue->setID(CHIVALUE_ID);
703 chiValue->setData(chi);
704 entry_plug->addProperty(chiValue);
705
706 DoubleData* integralOfChidtValue = new DoubleData();
707 integralOfChidtValue->setID(INTEGRALOFCHIDT_ID);
708 integralOfChidtValue->setData(integralOfChidt);
709 entry_plug->addProperty(integralOfChidtValue);
710
711 }
712 else
713 return NULL;
714
715 //get eta
716 foo = strtok(NULL, " ,;\t\n");
717 if(foo != NULL ){
718
719 for(int i = 0 ; i < 9; i++){
720
721 if(foo == NULL){
722 sprintf( painCave.errMsg,
723 "error in reading eta[%d] from %s\n", i, inFileName.c_str() );
724 return strdup( painCave.errMsg );
725 }
726 eta[i] = atof( foo );
727 foo = strtok(NULL, " ,;\t\n");
728 }
729 }
730 else
731 return NULL;
732
733 //push eta into SimInfo::properties which can be
734 //retrieved by integrator later
735 //entry_plug->setBoxM( theBoxMat3 );
736 DoubleArrayData* etaValue = new DoubleArrayData();
737 etaValue->setID(ETAVALUE_ID);
738 etaValue->setData(eta, 9);
739 entry_plug->addProperty(etaValue);
740 }
741
742 return NULL;
743 }
744
745 #ifdef IS_MPI
746 void DumpReader::nodeZeroError( void ){
747 int j, myStatus;
748
749 myStatus = 0;
750 for (j = 0; j < mpiSim->getNProcessors(); j++) {
751 MPI_Send( &myStatus, 1, MPI_INT, j,
752 TAKE_THIS_TAG_INT, MPI_COMM_WORLD);
753 }
754
755
756 MPI_Finalize();
757 exit (0);
758
759 }
760
761 void DumpReader::anonymousNodeDie( void ){
762
763 MPI_Finalize();
764 exit (0);
765 }
766 #endif