ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/io/DumpReader.cpp
Revision: 3102
Committed: Fri Dec 29 21:43:18 2006 UTC (18 years, 4 months ago) by gezelter
File size: 16121 byte(s)
Log Message:
Adding the ability to compute Center of Mass properties on
reading in a dump file.  DumpReader must be told if it needs
to do this, however.

File Contents

# Content
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), needCOMprops_(false) {
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,
84 "DumpReader: Cannot open file: %s\n",
85 filename_.c_str());
86 painCave.isFatal = 1;
87 simError();
88 }
89
90 #ifdef IS_MPI
91
92 }
93
94 strcpy(checkPointMsg, "Dump file opened for reading successfully.");
95 MPIcheckPoint();
96
97 #endif
98
99 return;
100 }
101
102 DumpReader::~DumpReader() {
103
104 #ifdef IS_MPI
105
106 if (worldRank == 0) {
107 #endif
108
109 delete inFile_;
110
111 #ifdef IS_MPI
112
113 }
114
115 strcpy(checkPointMsg, "Dump file closed successfully.");
116 MPIcheckPoint();
117
118 #endif
119
120 return;
121 }
122
123 int DumpReader::getNFrames(void) {
124
125 if (!isScanned_)
126 scanFile();
127
128 return nframes_;
129 }
130
131 void DumpReader::scanFile(void) {
132 int lineNo = 0;
133 std::streampos prevPos;
134 std::streampos currPos;
135
136 #ifdef IS_MPI
137
138 if (worldRank == 0) {
139 #endif // is_mpi
140
141 currPos = inFile_->tellg();
142 prevPos = currPos;
143 bool foundOpenSnapshotTag = false;
144 bool foundClosedSnapshotTag = false;
145 while(inFile_->getline(buffer, bufferSize)) {
146 ++lineNo;
147
148 std::string line = buffer;
149 currPos = inFile_->tellg();
150 if (line.find("<Snapshot>")!= std::string::npos) {
151 if (foundOpenSnapshotTag) {
152 sprintf(painCave.errMsg,
153 "DumpReader:<Snapshot> is multiply nested at line %d in %s \n", lineNo,
154 filename_.c_str());
155 painCave.isFatal = 1;
156 simError();
157 }
158 foundOpenSnapshotTag = true;
159 foundClosedSnapshotTag = false;
160 framePos_.push_back(prevPos);
161
162 } else if (line.find("</Snapshot>") != std::string::npos){
163 if (!foundOpenSnapshotTag) {
164 sprintf(painCave.errMsg,
165 "DumpReader:</Snapshot> appears before <Snapshot> at line %d in %s \n", lineNo,
166 filename_.c_str());
167 painCave.isFatal = 1;
168 simError();
169 }
170
171 if (foundClosedSnapshotTag) {
172 sprintf(painCave.errMsg,
173 "DumpReader:</Snapshot> appears multiply nested at line %d in %s \n", lineNo,
174 filename_.c_str());
175 painCave.isFatal = 1;
176 simError();
177 }
178 foundClosedSnapshotTag = true;
179 foundOpenSnapshotTag = false;
180 }
181 prevPos = currPos;
182 }
183
184 // only found <Snapshot> for the last frame means the file is corrupted, we should discard
185 // it and give a warning message
186 if (foundOpenSnapshotTag) {
187 sprintf(painCave.errMsg,
188 "DumpReader: last frame in %s is invalid\n", filename_.c_str());
189 painCave.isFatal = 0;
190 simError();
191 framePos_.pop_back();
192 }
193
194 nframes_ = framePos_.size();
195
196 if (nframes_ == 0) {
197 sprintf(painCave.errMsg,
198 "DumpReader: %s does not contain a valid frame\n", filename_.c_str());
199 painCave.isFatal = 1;
200 simError();
201 }
202 #ifdef IS_MPI
203 }
204
205 MPI_Bcast(&nframes_, 1, MPI_INT, 0, MPI_COMM_WORLD);
206
207 #endif // is_mpi
208
209 isScanned_ = true;
210 }
211
212 void DumpReader::readFrame(int whichFrame) {
213 if (!isScanned_)
214 scanFile();
215
216 int storageLayout = info_->getSnapshotManager()->getStorageLayout();
217
218 if (storageLayout & DataStorage::dslPosition) {
219 needPos_ = true;
220 } else {
221 needPos_ = false;
222 }
223
224 if (storageLayout & DataStorage::dslVelocity) {
225 needVel_ = true;
226 } else {
227 needVel_ = false;
228 }
229
230 if (storageLayout & DataStorage::dslAmat || storageLayout & DataStorage::dslElectroFrame) {
231 needQuaternion_ = true;
232 } else {
233 needQuaternion_ = false;
234 }
235
236 if (storageLayout & DataStorage::dslAngularMomentum) {
237 needAngMom_ = true;
238 } else {
239 needAngMom_ = false;
240 }
241
242 readSet(whichFrame);
243
244 if (needCOMprops_) {
245 Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
246 Vector3d com;
247 Vector3d comvel;
248 Vector3d comw;
249 info_->getComAll(com, comvel);
250 comw = info_->getAngularMomentum();
251 s->setCOMprops(com, comvel, comw);
252 }
253
254 }
255
256 void DumpReader::readSet(int whichFrame) {
257 std::string line;
258
259 #ifndef IS_MPI
260
261 inFile_->clear();
262 inFile_->seekg(framePos_[whichFrame]);
263
264 std::istream& inputStream = *inFile_;
265
266 #else
267 int masterNode = 0;
268 std::stringstream sstream;
269 if (worldRank == masterNode) {
270 std::string sendBuffer;
271
272 inFile_->clear();
273 inFile_->seekg(framePos_[whichFrame]);
274
275 while (inFile_->getline(buffer, bufferSize)) {
276
277 line = buffer;
278 sendBuffer += line;
279 sendBuffer += '\n';
280 if (line.find("</Snapshot>") != std::string::npos) {
281 break;
282 }
283 }
284
285 int sendBufferSize = sendBuffer.size();
286 MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
287 MPI_Bcast((void *)sendBuffer.c_str(), sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
288
289 sstream.str(sendBuffer);
290 } else {
291 int sendBufferSize;
292 MPI_Bcast(&sendBufferSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
293 char * recvBuffer = new char[sendBufferSize+1];
294 MPI_Bcast(recvBuffer, sendBufferSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
295 sstream.str(recvBuffer);
296 }
297
298 std::istream& inputStream = sstream;
299 #endif
300
301 inputStream.getline(buffer, bufferSize);
302
303 line = buffer;
304 if (line.find("<Snapshot>") == std::string::npos) {
305 sprintf(painCave.errMsg,
306 "DumpReader Error: can not find <Snapshot>\n");
307 painCave.isFatal = 1;
308 simError();
309 }
310
311 //read frameData
312 readFrameProperties(inputStream);
313
314 //read StuntDoubles
315 readStuntDoubles(inputStream);
316
317 inputStream.getline(buffer, bufferSize);
318 line = buffer;
319 if (line.find("</Snapshot>") == std::string::npos) {
320 sprintf(painCave.errMsg,
321 "DumpReader Error: can not find </Snapshot>\n");
322 painCave.isFatal = 1;
323 simError();
324 }
325
326 }
327
328 void DumpReader::parseDumpLine(const std::string& line) {
329
330
331 StringTokenizer tokenizer(line);
332 int nTokens;
333
334 nTokens = tokenizer.countTokens();
335
336 if (nTokens < 2) {
337 sprintf(painCave.errMsg,
338 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
339 painCave.isFatal = 1;
340 simError();
341 }
342
343 int index = tokenizer.nextTokenAsInt();
344
345 StuntDouble* integrableObject = info_->getIOIndexToIntegrableObject(index);
346
347 if (integrableObject == NULL) {
348 return;
349 }
350 std::string type = tokenizer.nextToken();
351 int size = type.size();
352
353 for(int i = 0; i < size; ++i) {
354 switch(type[i]) {
355
356 case 'p': {
357 Vector3d pos;
358 pos[0] = tokenizer.nextTokenAsDouble();
359 pos[1] = tokenizer.nextTokenAsDouble();
360 pos[2] = tokenizer.nextTokenAsDouble();
361 if (needPos_) {
362 integrableObject->setPos(pos);
363 }
364 break;
365 }
366 case 'v' : {
367 Vector3d vel;
368 vel[0] = tokenizer.nextTokenAsDouble();
369 vel[1] = tokenizer.nextTokenAsDouble();
370 vel[2] = tokenizer.nextTokenAsDouble();
371 if (needVel_) {
372 integrableObject->setVel(vel);
373 }
374 break;
375 }
376
377 case 'q' : {
378 Quat4d q;
379 if (integrableObject->isDirectional()) {
380
381 q[0] = tokenizer.nextTokenAsDouble();
382 q[1] = tokenizer.nextTokenAsDouble();
383 q[2] = tokenizer.nextTokenAsDouble();
384 q[3] = tokenizer.nextTokenAsDouble();
385
386 RealType qlen = q.length();
387 if (qlen < oopse::epsilon) { //check quaternion is not equal to 0
388
389 sprintf(painCave.errMsg,
390 "DumpReader Error: initial quaternion error (q0^2 + q1^2 + q2^2 + q3^2) ~ 0\n");
391 painCave.isFatal = 1;
392 simError();
393
394 }
395
396 q.normalize();
397 if (needQuaternion_) {
398 integrableObject->setQ(q);
399 }
400 }
401 break;
402 }
403 case 'j' : {
404 Vector3d ji;
405 if (integrableObject->isDirectional()) {
406 ji[0] = tokenizer.nextTokenAsDouble();
407 ji[1] = tokenizer.nextTokenAsDouble();
408 ji[2] = tokenizer.nextTokenAsDouble();
409 if (needAngMom_) {
410 integrableObject->setJ(ji);
411 }
412 }
413 break;
414 }
415 case 'f': {
416
417 Vector3d force;
418 force[0] = tokenizer.nextTokenAsDouble();
419 force[1] = tokenizer.nextTokenAsDouble();
420 force[2] = tokenizer.nextTokenAsDouble();
421 integrableObject->setFrc(force);
422 break;
423 }
424 case 't' : {
425
426 Vector3d torque;
427 torque[0] = tokenizer.nextTokenAsDouble();
428 torque[1] = tokenizer.nextTokenAsDouble();
429 torque[2] = tokenizer.nextTokenAsDouble();
430 integrableObject->setTrq(torque);
431 break;
432 }
433 default: {
434 sprintf(painCave.errMsg,
435 "DumpReader Error: %s is an unrecognized type\n", type.c_str());
436 painCave.isFatal = 1;
437 simError();
438 break;
439 }
440
441 }
442 }
443
444 }
445
446
447 void DumpReader::readStuntDoubles(std::istream& inputStream) {
448
449 inputStream.getline(buffer, bufferSize);
450 std::string line(buffer);
451
452 if (line.find("<StuntDoubles>") == std::string::npos) {
453 sprintf(painCave.errMsg,
454 "DumpReader Error: Missing <StuntDoubles>\n");
455 painCave.isFatal = 1;
456 simError();
457 }
458
459 while(inputStream.getline(buffer, bufferSize)) {
460 line = buffer;
461
462 if(line.find("</StuntDoubles>") != std::string::npos) {
463 break;
464 }
465
466 parseDumpLine(line);
467 }
468
469 }
470
471 void DumpReader::readFrameProperties(std::istream& inputStream) {
472
473 Snapshot* s = info_->getSnapshotManager()->getCurrentSnapshot();
474 inputStream.getline(buffer, bufferSize);
475 std::string line(buffer);
476
477 if (line.find("<FrameData>") == std::string::npos) {
478 sprintf(painCave.errMsg,
479 "DumpReader Error: Missing <FrameData>\n");
480 painCave.isFatal = 1;
481 simError();
482 }
483
484 while(inputStream.getline(buffer, bufferSize)) {
485 line = buffer;
486
487 if(line.find("</FrameData>") != std::string::npos) {
488 break;
489 }
490
491 StringTokenizer tokenizer(line, " ;\t\n\r{}:,");
492 if (!tokenizer.hasMoreTokens()) {
493 sprintf(painCave.errMsg,
494 "DumpReader Error: Not enough Tokens.\n%s\n", line.c_str());
495 painCave.isFatal = 1;
496 simError();
497 }
498
499 std::string propertyName = tokenizer.nextToken();
500 if (propertyName == "Time") {
501 RealType currTime = tokenizer.nextTokenAsDouble();
502 s->setTime(currTime);
503 } else if (propertyName == "Hmat"){
504 Mat3x3d hmat;
505 hmat(0, 0) = tokenizer.nextTokenAsDouble();
506 hmat(0, 1) = tokenizer.nextTokenAsDouble();
507 hmat(0, 2) = tokenizer.nextTokenAsDouble();
508 hmat(1, 0) = tokenizer.nextTokenAsDouble();
509 hmat(1, 1) = tokenizer.nextTokenAsDouble();
510 hmat(1, 2) = tokenizer.nextTokenAsDouble();
511 hmat(2, 0) = tokenizer.nextTokenAsDouble();
512 hmat(2, 1) = tokenizer.nextTokenAsDouble();
513 hmat(2, 2) = tokenizer.nextTokenAsDouble();
514 s->setHmat(hmat);
515 } else if (propertyName == "Thermostat") {
516 RealType chi = tokenizer.nextTokenAsDouble();
517 RealType integralOfChiDt = tokenizer.nextTokenAsDouble();
518 s->setChi(chi);
519 s->setIntegralOfChiDt(integralOfChiDt);
520 } else if (propertyName == "Barostat") {
521 Mat3x3d eta;
522 eta(0, 0) = tokenizer.nextTokenAsDouble();
523 eta(0, 1) = tokenizer.nextTokenAsDouble();
524 eta(0, 2) = tokenizer.nextTokenAsDouble();
525 eta(1, 0) = tokenizer.nextTokenAsDouble();
526 eta(1, 1) = tokenizer.nextTokenAsDouble();
527 eta(1, 2) = tokenizer.nextTokenAsDouble();
528 eta(2, 0) = tokenizer.nextTokenAsDouble();
529 eta(2, 1) = tokenizer.nextTokenAsDouble();
530 eta(2, 2) = tokenizer.nextTokenAsDouble();
531 s->setEta(eta);
532 } else {
533 sprintf(painCave.errMsg,
534 "DumpReader Error: %s is an invalid property in <FrameData>\n", propertyName.c_str());
535 painCave.isFatal = 0;
536 simError();
537 }
538
539 }
540
541 }
542
543
544 }//end namespace oopse