ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/libmdtools/InitializeFromFile.cpp
Revision: 1334
Committed: Fri Jul 16 18:58:03 2004 UTC (19 years, 11 months ago) by gezelter
File size: 16156 byte(s)
Log Message:
Initial import of OOPSE-1.0 source tree

File Contents

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