ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/io/DumpReader.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/io/DumpReader.cpp (file contents):
Revision 2101 by chrisfen, Thu Mar 10 15:10:24 2005 UTC vs.
Revision 2318 by tim, Wed Sep 21 20:59:31 2005 UTC

# Line 1 | Line 1
1 < /*
2 < * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 < *
4 < * The University of Notre Dame grants you ("Licensee") a
5 < * non-exclusive, royalty free, license to use, modify and
6 < * redistribute this software in source and binary code form, provided
7 < * that the following conditions are met:
8 < *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
19 < *    notice, this list of conditions and the following disclaimer.
20 < *
21 < * 3. Redistributions in binary form must reproduce the above copyright
22 < *    notice, this list of conditions and the following disclaimer in the
23 < *    documentation and/or other materials provided with the
24 < *    distribution.
25 < *
26 < * This software is provided "AS IS," without a warranty of any
27 < * kind. All express or implied conditions, representations and
28 < * warranties, including any implied warranty of merchantability,
29 < * fitness for a particular purpose or non-infringement, are hereby
30 < * excluded.  The University of Notre Dame and its licensors shall not
31 < * be liable for any damages suffered by licensee as a result of
32 < * using, modifying or distributing the software or its
33 < * derivatives. In no event will the University of Notre Dame or its
34 < * licensors be liable for any lost revenue, profit or data, or for
35 < * direct, indirect, special, consequential, incidental or punitive
36 < * damages, however caused and regardless of the theory of liability,
37 < * arising out of the use of or inability to use software, even if the
38 < * University of Notre Dame has been advised of the possibility of
39 < * such damages.
40 < */
41 <
42 < #define _LARGEFILE_SOURCE64
43 < #define _FILE_OFFSET_BITS 64
44 <
45 < #include <sys/types.h>
46 < #include <sys/stat.h>
47 <
48 < #include <iostream>
49 < #include <math.h>
50 <
51 < #include <stdio.h>
52 < #include <stdlib.h>
53 < #include <string.h>
54 <
55 < #include "io/DumpReader.hpp"
56 < #include "primitives/Molecule.hpp"
57 < #include "utils/simError.h"
58 < #include "utils/MemoryUtils.hpp"
59 < #include "utils/StringTokenizer.hpp"
60 <
61 < #ifdef IS_MPI
62 <
63 < #include <mpi.h>
64 < #define TAKE_THIS_TAG_CHAR 0
65 < #define TAKE_THIS_TAG_INT 1
66 <
67 < #endif // is_mpi
68 <
69 <
70 < namespace oopse {
71 <  
72 <  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73 <  : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74 <    
75 < #ifdef IS_MPI
76 <    
77 <    if (worldRank == 0) {
78 < #endif
79 <      
80 <      inFile_ = fopen(filename_.c_str(), "r");
81 <      
82 <      if (inFile_ == NULL) {
83 <        sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84 <        painCave.isFatal = 1;
85 <        simError();
86 <      }
87 <      
88 < #ifdef IS_MPI
89 <      
90 <    }
91 <    
92 <    strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93 <    MPIcheckPoint();
94 <    
95 < #endif
96 <    
97 <    return;
98 <  }
99 <  
100 <  DumpReader::~DumpReader() {
101 <    
102 < #ifdef IS_MPI
103 <    
104 <    if (worldRank == 0) {
105 < #endif
106 <      
107 <      int error;
108 <      error = fclose(inFile_);
109 <      
110 <      if (error) {
111 <        sprintf(painCave.errMsg, "DumpReader Error: Error closing %s\n", filename_.c_str());
112 <        painCave.isFatal = 1;            
113 <        simError();
114 <      }
115 <      
116 <      MemoryUtils::deletePointers(framePos_);
117 <      
118 < #ifdef IS_MPI
119 <      
120 <    }
121 <    
122 <    strcpy(checkPointMsg, "Dump file closed successfully.");
123 <    MPIcheckPoint();
124 <    
125 < #endif
126 <    
127 <    return;
128 <  }
129 <  
130 <  int DumpReader::getNFrames(void) {
131 <    
132 <    if (!isScanned_)
133 <      scanFile();
134 <    
135 <    return nframes_;
136 <  }
137 <  
138 <  void DumpReader::scanFile(void) {
139 <    int i, j;
140 <    int lineNum = 0;
141 <    char readBuffer[maxBufferSize];
142 <    fpos_t * currPos;
143 <    
144 < #ifdef IS_MPI
145 <    
146 <    if (worldRank == 0) {
147 < #endif // is_mpi
148 <      
149 <      rewind(inFile_);
150 <      
151 <      currPos = new fpos_t;
152 <      fgetpos(inFile_, currPos);
153 <      fgets(readBuffer, sizeof(readBuffer), inFile_);
154 <      lineNum++;
155 <      
156 <      if (feof(inFile_)) {
157 <        sprintf(painCave.errMsg,
158 <                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
159 <                filename_.c_str(),
160 <                lineNum);
161 <        painCave.isFatal = 1;
162 <        simError();
163 <      }
164 <      
165 <      while (!feof(inFile_)) {
166 <        framePos_.push_back(currPos);
167 <        
168 <        i = atoi(readBuffer);
169 <        
170 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
171 <        lineNum++;
172 <        
173 <        if (feof(inFile_)) {
174 <          sprintf(painCave.errMsg,
175 <                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
176 <                  filename_.c_str(),
177 <                  lineNum);
178 <          painCave.isFatal = 1;
179 <          simError();
180 <        }
181 <        
182 <        for(j = 0; j < i; j++) {
183 <          fgets(readBuffer, sizeof(readBuffer), inFile_);
184 <          lineNum++;
185 <          
186 <          if (feof(inFile_)) {
187 <            sprintf(painCave.errMsg,
188 <                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
189 <                    " with atom %d\n", filename_.c_str(),
190 <                    lineNum,
191 <                    j);
192 <            
193 <            painCave.isFatal = 1;
194 <            simError();
195 <          }
196 <        }
197 <        
198 <        currPos = new fpos_t;
199 <        fgetpos(inFile_, currPos);
200 <        fgets(readBuffer, sizeof(readBuffer), inFile_);
201 <        lineNum++;
202 <      }
203 <      
204 <      delete currPos;
205 <      rewind(inFile_);
206 <      
207 <      nframes_ = framePos_.size();
208 < #ifdef IS_MPI
209 <    }
210 <    
211 <    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
212 <    
213 <    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
214 <    MPIcheckPoint();
215 <    
216 < #endif // is_mpi
217 <    
218 <    isScanned_ = true;
219 <  }
220 <  
221 <  void DumpReader::readFrame(int whichFrame) {
222 <    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
223 <    
224 <    if (storageLayout & DataStorage::dslPosition) {
225 <      needPos_ = true;
226 <    } else {
227 <      needPos_ = false;
228 <    }
229 <    
230 <    if (storageLayout & DataStorage::dslVelocity) {
231 <      needVel_ = true;
232 <    } else {
233 <      needVel_ = false;
234 <    }
235 <    
236 <    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
237 <      needQuaternion_ = true;
238 <    } else {
239 <      needQuaternion_ = false;
240 <    }
241 <    
242 <    if (storageLayout & DataStorage::dslAngularMomentum) {
243 <      needAngMom_ = true;
244 <    } else {
245 <      needAngMom_ = false;    
246 <    }
247 <    
248 <    readSet(whichFrame);
249 <  }
250 <  
251 <  void DumpReader::readSet(int whichFrame) {
252 <    int i;
253 <    int nTotObjs;                  // the number of atoms
254 <    char read_buffer[maxBufferSize];  //the line buffer for reading
255 <    char * eof_test;               // ptr to see when we reach the end of the file
256 <    
257 <    Molecule* mol;
258 <    StuntDouble* integrableObject;
259 <    SimInfo::MoleculeIterator mi;
260 <    Molecule::IntegrableObjectIterator ii;
261 <    
262 < #ifndef IS_MPI
263 <    
264 <    fsetpos(inFile_, framePos_[whichFrame]);
265 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
266 <    
267 <    if (eof_test == NULL) {
268 <      sprintf(painCave.errMsg,
269 <              "DumpReader error: error reading 1st line of \"%s\"\n",
270 <              filename_.c_str());
271 <      painCave.isFatal = 1;
272 <      simError();
273 <    }
274 <    
275 <    nTotObjs = atoi(read_buffer);
276 <    
277 <    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
278 <      sprintf(painCave.errMsg,
279 <              "DumpReader error. %s nIntegrable, %d, "
280 <              "does not match the meta-data file's nIntegrable, %d.\n",
281 <              filename_.c_str(),
282 <              nTotObjs,
283 <              info_->getNGlobalIntegrableObjects());
284 <      
285 <      painCave.isFatal = 1;
286 <      simError();
287 <    }
288 <    
289 <    //read the box mat from the comment line
290 <    
291 <    eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
292 <    
293 <    if (eof_test == NULL) {
294 <      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
295 <              filename_.c_str());
296 <      painCave.isFatal = 1;
297 <      simError();
298 <    }
299 <    
300 <    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
301 <    
302 <    //parse dump lines
303 <    
304 <    i = 0;
305 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
306 <      
307 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
308 <           integrableObject = mol->nextIntegrableObject(ii)) {          
309 <        
310 <        eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
311 <        
312 <        if (eof_test == NULL) {
313 <          sprintf(painCave.errMsg,
314 <                  "DumpReader Error: error in reading file %s\n"
315 <                  "natoms  = %d; index = %d\n"
316 <                  "error reading the line from the file.\n",
317 <                  filename_.c_str(),
318 <                  nTotObjs,
319 <                  i);
320 <          
321 <          painCave.isFatal = 1;
322 <          simError();
323 <        }
324 <        
325 <        parseDumpLine(read_buffer, integrableObject);
326 <        i++;
327 <      }
328 <    }
329 <    
330 <    // MPI Section of code..........
331 <    
332 < #else //IS_MPI
333 <    
334 <    // first thing first, suspend fatalities.
335 <    int masterNode = 0;
336 <    int nCurObj;
337 <    painCave.isEventLoop = 1;
338 <    
339 <    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
340 <    int haveError;
341 <    
342 <    MPI_Status istatus;
343 <    int nitems;
344 <    
345 <    nTotObjs = info_->getNGlobalIntegrableObjects();
346 <    haveError = 0;
347 <    
348 <    if (worldRank == masterNode) {
349 <      fsetpos(inFile_, framePos_[whichFrame]);
350 <      
351 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
352 <      
353 <      if (eof_test == NULL) {
354 <        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
355 <                filename_.c_str());
356 <        painCave.isFatal = 1;
357 <        simError();
358 <      }
359 <      
360 <      nitems = atoi(read_buffer);
361 <      
362 <      // Check to see that the number of integrable objects in the
363 <      // intial configuration file is the same as derived from the
364 <      // meta-data file.
365 <      
366 <      if (nTotObjs != nitems) {
367 <        sprintf(painCave.errMsg,
368 <                "DumpReader Error. %s nIntegrable, %d, "
369 <                "does not match the meta-data file's nIntegrable, %d.\n",
370 <                filename_.c_str(),
371 <                nTotObjs,
372 <                info_->getNGlobalIntegrableObjects());
373 <        
374 <        painCave.isFatal = 1;
375 <        simError();
376 <      }
377 <      
378 <      //read the boxMat from the comment line
379 <      
380 <      eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
381 <      
382 <      if (eof_test == NULL) {
383 <        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
384 <                filename_.c_str());
385 <        painCave.isFatal = 1;
386 <        simError();
387 <      }
388 <      
389 <      //Every single processor will parse the comment line by itself
390 <      //By using this way, we might lose some efficiency, but if we want to add
391 <      //more parameters into comment line, we only need to modify function
392 <      //parseCommentLine
393 <      
394 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
395 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
396 <      
397 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
398 <        int which_node = info_->getMolToProc(i);
399 <        
400 <        if (which_node == masterNode) {
401 <          //molecules belong to master node
402 <          
403 <          mol = info_->getMoleculeByGlobalIndex(i);
404 <          
405 <          if (mol == NULL) {
406 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
407 <            painCave.isFatal = 1;
408 <            simError();
409 <          }
410 <          
411 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
412 <               integrableObject = mol->nextIntegrableObject(ii)){
413 <            
414 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
415 <            
416 <            if (eof_test == NULL) {
417 <              sprintf(painCave.errMsg,
418 <                      "DumpReader Error: error in reading file %s\n"
419 <                      "natoms  = %d; index = %d\n"
420 <                      "error reading the line from the file.\n",
421 <                      filename_.c_str(),
422 <                      nTotObjs,
423 <                      i);
424 <              
425 <              painCave.isFatal = 1;
426 <              simError();
427 <            }
428 <            
429 <            parseDumpLine(read_buffer, integrableObject);
430 <          }
431 <        } else {
432 <          //molecule belongs to slave nodes
433 <          
434 <          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
435 <                   MPI_COMM_WORLD, &istatus);
436 <          
437 <          for(int j = 0; j < nCurObj; j++) {
438 <            eof_test = fgets(read_buffer, sizeof(read_buffer), inFile_);
439 <            
440 <            if (eof_test == NULL) {
441 <              sprintf(painCave.errMsg,
442 <                      "DumpReader Error: error in reading file %s\n"
443 <                      "natoms  = %d; index = %d\n"
444 <                      "error reading the line from the file.\n",
445 <                      filename_.c_str(),
446 <                      nTotObjs,
447 <                      i);
448 <              
449 <              painCave.isFatal = 1;
450 <              simError();
451 <            }
452 <            
453 <            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
454 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
455 <          }
456 <        }
457 <      }
458 <    } else {
459 <      //actions taken at slave nodes
460 <      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
461 <      
462 <      /**@todo*/
463 <      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
464 <      
465 <      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
466 <        int which_node = info_->getMolToProc(i);
467 <        
468 <        if (which_node == worldRank) {
469 <          //molecule with global index i belongs to this processor
470 <          
471 <          mol = info_->getMoleculeByGlobalIndex(i);
472 <          if (mol == NULL) {
473 <            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
474 <            painCave.isFatal = 1;
475 <            simError();
476 <          }
477 <          
478 <          nCurObj = mol->getNIntegrableObjects();
479 <          
480 <          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
481 <                   MPI_COMM_WORLD);
482 <          
483 <          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
484 <               integrableObject = mol->nextIntegrableObject(ii)){
485 <            
486 <            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
487 <                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
488 <            
489 <            parseDumpLine(read_buffer, integrableObject);
490 <          }
491 <          
492 <        }
493 <        
494 <      }
495 <      
496 <    }
497 <    
498 < #endif
499 <    
500 <  }
501 <  
502 <  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
503 <    
504 <    Vector3d pos;  // position place holders
505 <    Vector3d vel;  // velocity placeholders
506 <    Quat4d q;    // the quaternions
507 <    Vector3d ji;   // angular velocity placeholders;
508 <    StringTokenizer tokenizer(line);
509 <    int nTokens;
510 <    
511 <    nTokens = tokenizer.countTokens();
512 <    
513 <    if (nTokens < 14) {
514 <      sprintf(painCave.errMsg,
515 <              "DumpReader Error: Not enough Tokens.\n%s\n", line);
516 <      painCave.isFatal = 1;
517 <      simError();
518 <    }
519 <    
520 <    std::string name = tokenizer.nextToken();
521 <    
522 <    if (name != integrableObject->getType()) {
523 <      
524 <      sprintf(painCave.errMsg,
525 <              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
526 <              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
527 <      painCave.isFatal = 1;
528 <      simError();        
529 <    }
530 <    
531 <    pos[0] = tokenizer.nextTokenAsDouble();
532 <    pos[1] = tokenizer.nextTokenAsDouble();
533 <    pos[2] = tokenizer.nextTokenAsDouble();
534 <    if (needPos_) {
535 <      integrableObject->setPos(pos);
536 <    }
537 <    
538 <    vel[0] = tokenizer.nextTokenAsDouble();
539 <    vel[1] = tokenizer.nextTokenAsDouble();
540 <    vel[2] = tokenizer.nextTokenAsDouble();
541 <    if (needVel_) {
542 <      integrableObject->setVel(vel);
543 <    }
544 <    
545 <    if (integrableObject->isDirectional()) {
546 <      
547 <      q[0] = tokenizer.nextTokenAsDouble();
548 <      q[1] = tokenizer.nextTokenAsDouble();
549 <      q[2] = tokenizer.nextTokenAsDouble();
550 <      q[3] = tokenizer.nextTokenAsDouble();
551 <      
552 <      double qlen = q.length();
553 <      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
554 <        
555 <        sprintf(painCave.errMsg,
556 <                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
557 <        painCave.isFatal = 1;
558 <        simError();
559 <        
560 <      }
561 <      
562 <      q.normalize();
563 <      if (needQuaternion_) {          
564 <        integrableObject->setQ(q);
565 <      }
566 <      
567 <      ji[0] = tokenizer.nextTokenAsDouble();
568 <      ji[1] = tokenizer.nextTokenAsDouble();
569 <      ji[2] = tokenizer.nextTokenAsDouble();
570 <      if (needAngMom_) {
571 <        integrableObject->setJ(ji);
572 <      }
573 <    }
574 <    
575 <  }
576 <  
577 <  
578 <  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
579 <    double currTime;
580 <    Mat3x3d hmat;
581 <    double chi;
582 <    double integralOfChiDt;
583 <    Mat3x3d eta;
584 <    
585 <    StringTokenizer tokenizer(line);
586 <    int nTokens;
587 <    
588 <    nTokens = tokenizer.countTokens();
589 <    
590 <    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
591 <    if (nTokens < 10) {
592 <      sprintf(painCave.errMsg,
593 <              "DumpReader Error: Not enough tokens in comment line: %s", line);
594 <      painCave.isFatal = 1;
595 <      simError();  
596 <    }
597 <    
598 <    //read current time
599 <    currTime = tokenizer.nextTokenAsDouble();
600 <    s->setTime(currTime);
601 <    
602 <    //read h-matrix
603 <    hmat(0, 0) = tokenizer.nextTokenAsDouble();
604 <    hmat(0, 1) = tokenizer.nextTokenAsDouble();
605 <    hmat(0, 2) = tokenizer.nextTokenAsDouble();
606 <    hmat(1, 0) = tokenizer.nextTokenAsDouble();
607 <    hmat(1, 1) = tokenizer.nextTokenAsDouble();
608 <    hmat(1, 2) = tokenizer.nextTokenAsDouble();
609 <    hmat(2, 0) = tokenizer.nextTokenAsDouble();
610 <    hmat(2, 1) = tokenizer.nextTokenAsDouble();
611 <    hmat(2, 2) = tokenizer.nextTokenAsDouble();
612 <    s->setHmat(hmat);
613 <    
614 <    //read chi and integralOfChidt, they should apprear in pair
615 <    if (tokenizer.countTokens() >= 2) {
616 <      chi = tokenizer.nextTokenAsDouble();
617 <      integralOfChiDt = tokenizer.nextTokenAsDouble();            
618 <      
619 <      s->setChi(chi);
620 <      s->setIntegralOfChiDt(integralOfChiDt);
621 <    }
622 <    
623 <    //read eta (eta is 3x3 matrix)
624 <    if (tokenizer.countTokens() >= 9) {
625 <      eta(0, 0) = tokenizer.nextTokenAsDouble();
626 <      eta(0, 1) = tokenizer.nextTokenAsDouble();
627 <      eta(0, 2) = tokenizer.nextTokenAsDouble();
628 <      eta(1, 0) = tokenizer.nextTokenAsDouble();
629 <      eta(1, 1) = tokenizer.nextTokenAsDouble();
630 <      eta(1, 2) = tokenizer.nextTokenAsDouble();
631 <      eta(2, 0) = tokenizer.nextTokenAsDouble();
632 <      eta(2, 1) = tokenizer.nextTokenAsDouble();
633 <      eta(2, 2) = tokenizer.nextTokenAsDouble();      
634 <      
635 <      s->setEta(eta);
636 <    }
637 <    
638 <    
639 <  }
640 <  
641 < }//end namespace oopse
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Acknowledgement of the program authors must be made in any
10 > *    publication of scientific results based in part on use of the
11 > *    program.  An acceptable form of acknowledgement is citation of
12 > *    the article in which the program was described (Matthew
13 > *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 > *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 > *    Parallel Simulation Engine for Molecular Dynamics,"
16 > *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 > *
18 > * 2. Redistributions of source code must retain the above copyright
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 > #define _LARGEFILE_SOURCE64
43 > #define _FILE_OFFSET_BITS 64
44 >
45 > #include <sys/types.h>
46 > #include <sys/stat.h>
47 >
48 > #include <iostream>
49 > #include <math.h>
50 >
51 > #include <stdio.h>
52 > #include <stdlib.h>
53 > #include <string.h>
54 >
55 > #include "io/DumpReader.hpp"
56 > #include "primitives/Molecule.hpp"
57 > #include "utils/simError.h"
58 > #include "utils/MemoryUtils.hpp"
59 > #include "utils/StringTokenizer.hpp"
60 >
61 > #ifdef IS_MPI
62 >
63 > #include <mpi.h>
64 > #define TAKE_THIS_TAG_CHAR 0
65 > #define TAKE_THIS_TAG_INT 1
66 >
67 > #endif // is_mpi
68 >
69 >
70 > namespace oopse {
71 >  
72 >  DumpReader::DumpReader(SimInfo* info, const std::string& filename)
73 >    : info_(info), filename_(filename), isScanned_(false), nframes_(0) {
74 >    
75 > #ifdef IS_MPI
76 >    
77 >      if (worldRank == 0) {
78 > #endif
79 >      
80 >      inFile_ = new std::ifstream(filename_.c_str());
81 >      
82 >        if (inFile_->fail()) {
83 >          sprintf(painCave.errMsg, "DumpReader: Cannot open file: %s\n", filename_.c_str());
84 >          painCave.isFatal = 1;
85 >          simError();
86 >        }
87 >      
88 > #ifdef IS_MPI
89 >      
90 >      }
91 >    
92 >      strcpy(checkPointMsg, "Dump file opened for reading successfully.");
93 >      MPIcheckPoint();
94 >    
95 > #endif
96 >    
97 >      return;
98 >    }
99 >  
100 >  DumpReader::~DumpReader() {
101 >    
102 > #ifdef IS_MPI
103 >    
104 >    if (worldRank == 0) {
105 > #endif
106 >      
107 >      delete inFile_;
108 >      
109 > #ifdef IS_MPI
110 >      
111 >    }
112 >    
113 >    strcpy(checkPointMsg, "Dump file closed successfully.");
114 >    MPIcheckPoint();
115 >    
116 > #endif
117 >    
118 >    return;
119 >  }
120 >  
121 >  int DumpReader::getNFrames(void) {
122 >    
123 >    if (!isScanned_)
124 >      scanFile();
125 >    
126 >    return nframes_;
127 >  }
128 >  
129 >  void DumpReader::scanFile(void) {
130 >    int i, j;
131 >    int lineNum = 0;
132 >    char readBuffer[maxBufferSize];
133 >    std::streampos  currPos;
134 >    
135 > #ifdef IS_MPI
136 >    
137 >    if (worldRank == 0) {
138 > #endif // is_mpi
139 >      
140 >      inFile_->seekg (0, std::ios::beg);
141 >      
142 >
143 >      currPos = inFile_->tellg();
144 >      inFile_->getline(readBuffer, sizeof(readBuffer));
145 >      lineNum++;
146 >      
147 >      if (inFile_->eof()) {
148 >        sprintf(painCave.errMsg,
149 >                "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
150 >                filename_.c_str(),
151 >                lineNum);
152 >        painCave.isFatal = 1;
153 >        simError();
154 >      }
155 >      
156 >      while (!inFile_->eof()) {
157 >        framePos_.push_back(currPos);
158 >        
159 >        i = atoi(readBuffer);
160 >        
161 >        inFile_->getline(readBuffer, sizeof(readBuffer));
162 >        lineNum++;
163 >        
164 >        if (inFile_->eof()) {
165 >          sprintf(painCave.errMsg,
166 >                  "DumpReader Error: File \"%s\" ended unexpectedly at line %d\n",
167 >                  filename_.c_str(),
168 >                  lineNum);
169 >          painCave.isFatal = 1;
170 >          simError();
171 >        }
172 >        
173 >        for(j = 0; j < i; j++) {
174 >          inFile_->getline(readBuffer, sizeof(readBuffer));
175 >          lineNum++;
176 >          
177 >          if (inFile_->eof()) {
178 >            sprintf(painCave.errMsg,
179 >                    "DumpReader Error: File \"%s\" ended unexpectedly at line %d,"
180 >                    " with atom %d\n", filename_.c_str(),
181 >                    lineNum,
182 >                    j);
183 >            
184 >            painCave.isFatal = 1;
185 >            simError();
186 >          }
187 >        }
188 >        
189 >        currPos = inFile_->tellg();
190 >        inFile_->getline(readBuffer, sizeof(readBuffer));
191 >        lineNum++;
192 >      }
193 >
194 >      inFile_->seekg (0, std::ios::beg);
195 >      
196 >      nframes_ = framePos_.size();
197 > #ifdef IS_MPI
198 >    }
199 >    
200 >    MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
201 >    
202 >    strcpy(checkPointMsg, "Successfully scanned DumpFile\n");
203 >    MPIcheckPoint();
204 >    
205 > #endif // is_mpi
206 >    
207 >    isScanned_ = true;
208 >  }
209 >  
210 >  void DumpReader::readFrame(int whichFrame) {
211 >    if (!isScanned_)
212 >      scanFile();
213 >        
214 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
215 >    
216 >    if (storageLayout & DataStorage::dslPosition) {
217 >      needPos_ = true;
218 >    } else {
219 >      needPos_ = false;
220 >    }
221 >    
222 >    if (storageLayout & DataStorage::dslVelocity) {
223 >      needVel_ = true;
224 >    } else {
225 >      needVel_ = false;
226 >    }
227 >    
228 >    if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
229 >      needQuaternion_ = true;
230 >    } else {
231 >      needQuaternion_ = false;
232 >    }
233 >    
234 >    if (storageLayout & DataStorage::dslAngularMomentum) {
235 >      needAngMom_ = true;
236 >    } else {
237 >      needAngMom_ = false;    
238 >    }
239 >    
240 >    readSet(whichFrame);
241 >  }
242 >  
243 >  void DumpReader::readSet(int whichFrame) {
244 >    int i;
245 >    int nTotObjs;                  // the number of atoms
246 >    char read_buffer[maxBufferSize];  //the line buffer for reading
247 >    char * eof_test;               // ptr to see when we reach the end of the file
248 >    
249 >    Molecule* mol;
250 >    StuntDouble* integrableObject;
251 >    SimInfo::MoleculeIterator mi;
252 >    Molecule::IntegrableObjectIterator ii;
253 >    
254 > #ifndef IS_MPI
255 >    inFile_->clear();
256 >    inFile_->seekg(framePos_[whichFrame]);
257 >        
258 >    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
259 >      sprintf(painCave.errMsg,
260 >              "DumpReader error: error reading 1st line of \"%s\"\n",
261 >              filename_.c_str());
262 >      painCave.isFatal = 1;
263 >      simError();
264 >    }
265 >    
266 >    nTotObjs = atoi(read_buffer);
267 >    
268 >    if (nTotObjs != info_->getNGlobalIntegrableObjects()) {
269 >      sprintf(painCave.errMsg,
270 >              "DumpReader error. %s nIntegrable, %d, "
271 >              "does not match the meta-data file's nIntegrable, %d.\n",
272 >              filename_.c_str(),
273 >              nTotObjs,
274 >              info_->getNGlobalIntegrableObjects());
275 >      
276 >      painCave.isFatal = 1;
277 >      simError();
278 >    }
279 >    
280 >    //read the box mat from the comment line
281 >    
282 >    
283 >    if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
284 >      sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
285 >              filename_.c_str());
286 >      painCave.isFatal = 1;
287 >      simError();
288 >    }
289 >    
290 >    parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
291 >    
292 >    //parse dump lines
293 >    
294 >    i = 0;
295 >    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
296 >      
297 >      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
298 >           integrableObject = mol->nextIntegrableObject(ii)) {          
299 >        
300 >        
301 >        
302 >        if (!inFile_->getline(read_buffer, sizeof(read_buffer))) {
303 >          sprintf(painCave.errMsg,
304 >                  "DumpReader Error: error in reading file %s\n"
305 >                  "natoms  = %d; index = %d\n"
306 >                  "error reading the line from the file.\n",
307 >                  filename_.c_str(),
308 >                  nTotObjs,
309 >                  i);
310 >          
311 >          painCave.isFatal = 1;
312 >          simError();
313 >        }
314 >        
315 >        parseDumpLine(read_buffer, integrableObject);
316 >        i++;
317 >      }
318 >    }
319 >    
320 >    // MPI Section of code..........
321 >    
322 > #else //IS_MPI
323 >    
324 >    // first thing first, suspend fatalities.
325 >    int masterNode = 0;
326 >    int nCurObj;
327 >    painCave.isEventLoop = 1;
328 >    
329 >    int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
330 >    int haveError;
331 >    
332 >    MPI_Status istatus;
333 >    int nitems;
334 >    
335 >    nTotObjs = info_->getNGlobalIntegrableObjects();
336 >    haveError = 0;
337 >    
338 >    if (worldRank == masterNode) {
339 >      inFile_->clear();            
340 >      inFile_->seekg(framePos_[whichFrame]);
341 >      
342 >      if (!inFile_->getline(readBuffer, sizeof(readBuffer))) {
343 >        sprintf(painCave.errMsg, "DumpReader Error: Error reading 1st line of %s \n ",
344 >                filename_.c_str());
345 >        painCave.isFatal = 1;
346 >        simError();
347 >      }
348 >      
349 >      nitems = atoi(read_buffer);
350 >      
351 >      // Check to see that the number of integrable objects in the
352 >      // intial configuration file is the same as derived from the
353 >      // meta-data file.
354 >      
355 >      if (nTotObjs != nitems) {
356 >        sprintf(painCave.errMsg,
357 >                "DumpReader Error. %s nIntegrable, %d, "
358 >                "does not match the meta-data file's nIntegrable, %d.\n",
359 >                filename_.c_str(),
360 >                nTotObjs,
361 >                info_->getNGlobalIntegrableObjects());
362 >        
363 >        painCave.isFatal = 1;
364 >        simError();
365 >      }
366 >      
367 >      //read the boxMat from the comment line
368 >      
369 >      
370 >      
371 >      if (!inFile_->getline(readBuffer, sizeof(readBuffer))) {
372 >        sprintf(painCave.errMsg, "DumpReader Error: error in reading commment in %s\n",
373 >                filename_.c_str());
374 >        painCave.isFatal = 1;
375 >        simError();
376 >      }
377 >      
378 >      //Every single processor will parse the comment line by itself
379 >      //By using this way, we might lose some efficiency, but if we want to add
380 >      //more parameters into comment line, we only need to modify function
381 >      //parseCommentLine
382 >      
383 >      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
384 >      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
385 >      
386 >      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
387 >        int which_node = info_->getMolToProc(i);
388 >        
389 >        if (which_node == masterNode) {
390 >          //molecules belong to master node
391 >          
392 >          mol = info_->getMoleculeByGlobalIndex(i);
393 >          
394 >          if (mol == NULL) {
395 >            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
396 >            painCave.isFatal = 1;
397 >            simError();
398 >          }
399 >          
400 >          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
401 >               integrableObject = mol->nextIntegrableObject(ii)){
402 >            
403 >            
404 >            
405 >            if (!inFile_->getline(readBuffer, sizeof(readBuffer))) {
406 >              sprintf(painCave.errMsg,
407 >                      "DumpReader Error: error in reading file %s\n"
408 >                      "natoms  = %d; index = %d\n"
409 >                      "error reading the line from the file.\n",
410 >                      filename_.c_str(),
411 >                      nTotObjs,
412 >                      i);
413 >              
414 >              painCave.isFatal = 1;
415 >              simError();
416 >            }
417 >            
418 >            parseDumpLine(read_buffer, integrableObject);
419 >          }
420 >        } else {
421 >          //molecule belongs to slave nodes
422 >          
423 >          MPI_Recv(&nCurObj, 1, MPI_INT, which_node, TAKE_THIS_TAG_INT,
424 >                   MPI_COMM_WORLD, &istatus);
425 >          
426 >          for(int j = 0; j < nCurObj; j++) {
427 >            
428 >            
429 >            if (!inFile_->getline(readBuffer, sizeof(readBuffer))) {
430 >              sprintf(painCave.errMsg,
431 >                      "DumpReader Error: error in reading file %s\n"
432 >                      "natoms  = %d; index = %d\n"
433 >                      "error reading the line from the file.\n",
434 >                      filename_.c_str(),
435 >                      nTotObjs,
436 >                      i);
437 >              
438 >              painCave.isFatal = 1;
439 >              simError();
440 >            }
441 >            
442 >            MPI_Send(read_buffer, maxBufferSize, MPI_CHAR, which_node,
443 >                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD);
444 >          }
445 >        }
446 >      }
447 >    } else {
448 >      //actions taken at slave nodes
449 >      MPI_Bcast(read_buffer, maxBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
450 >      
451 >      /**@todo*/
452 >      parseCommentLine(read_buffer, info_->getSnapshotManager()->getCurrentSnapshot());
453 >      
454 >      for(i = 0; i < info_->getNGlobalMolecules(); i++) {
455 >        int which_node = info_->getMolToProc(i);
456 >        
457 >        if (which_node == worldRank) {
458 >          //molecule with global index i belongs to this processor
459 >          
460 >          mol = info_->getMoleculeByGlobalIndex(i);
461 >          if (mol == NULL) {
462 >            sprintf(painCave.errMsg, "DumpReader Error: Molecule not found on node %d!", worldRank);
463 >            painCave.isFatal = 1;
464 >            simError();
465 >          }
466 >          
467 >          nCurObj = mol->getNIntegrableObjects();
468 >          
469 >          MPI_Send(&nCurObj, 1, MPI_INT, masterNode, TAKE_THIS_TAG_INT,
470 >                   MPI_COMM_WORLD);
471 >          
472 >          for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
473 >               integrableObject = mol->nextIntegrableObject(ii)){
474 >            
475 >            MPI_Recv(read_buffer, maxBufferSize, MPI_CHAR, masterNode,
476 >                     TAKE_THIS_TAG_CHAR, MPI_COMM_WORLD, &istatus);
477 >            
478 >            parseDumpLine(read_buffer, integrableObject);
479 >          }
480 >          
481 >        }
482 >        
483 >      }
484 >      
485 >    }
486 >    
487 > #endif
488 >    
489 >  }
490 >  
491 >  void DumpReader::parseDumpLine(char *line, StuntDouble *integrableObject) {
492 >    
493 >    Vector3d pos;  // position place holders
494 >    Vector3d vel;  // velocity placeholders
495 >    Quat4d q;    // the quaternions
496 >    Vector3d ji;   // angular velocity placeholders;
497 >    StringTokenizer tokenizer(line);
498 >    int nTokens;
499 >    
500 >    nTokens = tokenizer.countTokens();
501 >    
502 >    if (nTokens < 14) {
503 >      sprintf(painCave.errMsg,
504 >              "DumpReader Error: Not enough Tokens.\n%s\n", line);
505 >      painCave.isFatal = 1;
506 >      simError();
507 >    }
508 >    
509 >    std::string name = tokenizer.nextToken();
510 >    
511 >    if (name != integrableObject->getType()) {
512 >      
513 >      sprintf(painCave.errMsg,
514 >              "DumpReader Error: Atom type [%s] in %s does not match Atom Type [%s] in .md file.\n",
515 >              name.c_str(), filename_.c_str(), integrableObject->getType().c_str());
516 >      painCave.isFatal = 1;
517 >      simError();        
518 >    }
519 >    
520 >    pos[0] = tokenizer.nextTokenAsDouble();
521 >    pos[1] = tokenizer.nextTokenAsDouble();
522 >    pos[2] = tokenizer.nextTokenAsDouble();
523 >    if (needPos_) {
524 >      integrableObject->setPos(pos);
525 >    }
526 >    
527 >    vel[0] = tokenizer.nextTokenAsDouble();
528 >    vel[1] = tokenizer.nextTokenAsDouble();
529 >    vel[2] = tokenizer.nextTokenAsDouble();
530 >    if (needVel_) {
531 >      integrableObject->setVel(vel);
532 >    }
533 >    
534 >    if (integrableObject->isDirectional()) {
535 >      
536 >      q[0] = tokenizer.nextTokenAsDouble();
537 >      q[1] = tokenizer.nextTokenAsDouble();
538 >      q[2] = tokenizer.nextTokenAsDouble();
539 >      q[3] = tokenizer.nextTokenAsDouble();
540 >      
541 >      double qlen = q.length();
542 >      if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
543 >        
544 >        sprintf(painCave.errMsg,
545 >                "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2 ~ 0).\n");
546 >        painCave.isFatal = 1;
547 >        simError();
548 >        
549 >      }
550 >      
551 >      q.normalize();
552 >      if (needQuaternion_) {          
553 >        integrableObject->setQ(q);
554 >      }
555 >      
556 >      ji[0] = tokenizer.nextTokenAsDouble();
557 >      ji[1] = tokenizer.nextTokenAsDouble();
558 >      ji[2] = tokenizer.nextTokenAsDouble();
559 >      if (needAngMom_) {
560 >        integrableObject->setJ(ji);
561 >      }
562 >    }
563 >    
564 >  }
565 >  
566 >  
567 >  void DumpReader::parseCommentLine(char* line, Snapshot* s) {
568 >    double currTime;
569 >    Mat3x3d hmat;
570 >    double chi;
571 >    double integralOfChiDt;
572 >    Mat3x3d eta;
573 >    
574 >    StringTokenizer tokenizer(line);
575 >    int nTokens;
576 >    
577 >    nTokens = tokenizer.countTokens();
578 >    
579 >    //comment line should at least contain 10 tokens: current time(1 token) and  h-matrix(9 tokens)
580 >    if (nTokens < 10) {
581 >      sprintf(painCave.errMsg,
582 >              "DumpReader Error: Not enough tokens in comment line: %s", line);
583 >      painCave.isFatal = 1;
584 >      simError();  
585 >    }
586 >    
587 >    //read current time
588 >    currTime = tokenizer.nextTokenAsDouble();
589 >    s->setTime(currTime);
590 >    
591 >    //read h-matrix
592 >    hmat(0, 0) = tokenizer.nextTokenAsDouble();
593 >    hmat(0, 1) = tokenizer.nextTokenAsDouble();
594 >    hmat(0, 2) = tokenizer.nextTokenAsDouble();
595 >    hmat(1, 0) = tokenizer.nextTokenAsDouble();
596 >    hmat(1, 1) = tokenizer.nextTokenAsDouble();
597 >    hmat(1, 2) = tokenizer.nextTokenAsDouble();
598 >    hmat(2, 0) = tokenizer.nextTokenAsDouble();
599 >    hmat(2, 1) = tokenizer.nextTokenAsDouble();
600 >    hmat(2, 2) = tokenizer.nextTokenAsDouble();
601 >    s->setHmat(hmat);
602 >    
603 >    //read chi and integralOfChidt, they should apprear in pair
604 >    if (tokenizer.countTokens() >= 2) {
605 >      chi = tokenizer.nextTokenAsDouble();
606 >      integralOfChiDt = tokenizer.nextTokenAsDouble();            
607 >      
608 >      s->setChi(chi);
609 >      s->setIntegralOfChiDt(integralOfChiDt);
610 >    }
611 >    
612 >    //read eta (eta is 3x3 matrix)
613 >    if (tokenizer.countTokens() >= 9) {
614 >      eta(0, 0) = tokenizer.nextTokenAsDouble();
615 >      eta(0, 1) = tokenizer.nextTokenAsDouble();
616 >      eta(0, 2) = tokenizer.nextTokenAsDouble();
617 >      eta(1, 0) = tokenizer.nextTokenAsDouble();
618 >      eta(1, 1) = tokenizer.nextTokenAsDouble();
619 >      eta(1, 2) = tokenizer.nextTokenAsDouble();
620 >      eta(2, 0) = tokenizer.nextTokenAsDouble();
621 >      eta(2, 1) = tokenizer.nextTokenAsDouble();
622 >      eta(2, 2) = tokenizer.nextTokenAsDouble();      
623 >      
624 >      s->setEta(eta);
625 >    }
626 >    
627 >    
628 >  }
629 >  
630 > }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines