ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/InitializeFromFile.cpp
Revision: 1113
Committed: Thu Apr 15 16:18:26 2004 UTC (20 years, 2 months ago) by tim
File size: 15756 byte(s)
Log Message:
fix whole bunch of bugs :-)

File Contents

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