ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/io/DumpWriter.cpp
Revision: 2318
Committed: Wed Sep 21 20:59:31 2005 UTC (18 years, 9 months ago) by tim
File size: 17658 byte(s)
Log Message:
DumpReader using pure c++ io

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 #include "io/DumpWriter.hpp"
43 #include "primitives/Molecule.hpp"
44 #include "utils/simError.h"
45 #include "io/basic_teebuf.hpp"
46 #include "io/gzstream.hpp"
47 #include "io/Globals.hpp"
48
49 #ifdef IS_MPI
50 #include <mpi.h>
51 #endif //is_mpi
52
53 namespace oopse {
54
55 DumpWriter::DumpWriter(SimInfo* info)
56 : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
57
58 Globals* simParams = info->getSimParams();
59 needCompression_ = simParams->getCompressDumpFile();
60
61 #ifdef HAVE_LIBZ
62 if (needCompression_) {
63 filename_ += ".gz";
64 eorFilename_ += ".gz";
65 }
66 #endif
67
68 #ifdef IS_MPI
69
70 if (worldRank == 0) {
71 #endif // is_mpi
72
73
74 dumpFile_ = createOStream(filename_);
75
76 if (!dumpFile_) {
77 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
78 filename_.c_str());
79 painCave.isFatal = 1;
80 simError();
81 }
82
83 #ifdef IS_MPI
84
85 }
86
87 sprintf(checkPointMsg, "Sucessfully opened output file for dumping.\n");
88 MPIcheckPoint();
89
90 #endif // is_mpi
91
92 }
93
94
95 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
96 : info_(info), filename_(filename){
97
98 Globals* simParams = info->getSimParams();
99 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
100
101 needCompression_ = simParams->getCompressDumpFile();
102
103 #ifdef HAVE_LIBZ
104 if (needCompression_) {
105 filename_ += ".gz";
106 eorFilename_ += ".gz";
107 }
108 #endif
109
110 #ifdef IS_MPI
111
112 if (worldRank == 0) {
113 #endif // is_mpi
114
115
116 dumpFile_ = createOStream(filename_);
117
118 if (!dumpFile_) {
119 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
120 filename_.c_str());
121 painCave.isFatal = 1;
122 simError();
123 }
124
125 #ifdef IS_MPI
126
127 }
128
129 sprintf(checkPointMsg, "Sucessfully opened output file for dumping.\n");
130 MPIcheckPoint();
131
132 #endif // is_mpi
133
134 }
135
136 DumpWriter::~DumpWriter() {
137
138 #ifdef IS_MPI
139
140 if (worldRank == 0) {
141 #endif // is_mpi
142
143 delete dumpFile_;
144
145 #ifdef IS_MPI
146
147 }
148
149 #endif // is_mpi
150
151 }
152
153 void DumpWriter::writeCommentLine(std::ostream& os, Snapshot* s) {
154
155 double currentTime;
156 Mat3x3d hmat;
157 double chi;
158 double integralOfChiDt;
159 Mat3x3d eta;
160
161 currentTime = s->getTime();
162 hmat = s->getHmat();
163 chi = s->getChi();
164 integralOfChiDt = s->getIntegralOfChiDt();
165 eta = s->getEta();
166
167 os << currentTime << ";\t"
168 << hmat(0, 0) << "\t" << hmat(1, 0) << "\t" << hmat(2, 0) << ";\t"
169 << hmat(0, 1) << "\t" << hmat(1, 1) << "\t" << hmat(2, 1) << ";\t"
170 << hmat(0, 2) << "\t" << hmat(1, 2) << "\t" << hmat(2, 2) << ";\t";
171
172 //write out additional parameters, such as chi and eta
173
174 os << chi << "\t" << integralOfChiDt << "\t;";
175
176 os << eta(0, 0) << "\t" << eta(1, 0) << "\t" << eta(2, 0) << ";\t"
177 << eta(0, 1) << "\t" << eta(1, 1) << "\t" << eta(2, 1) << ";\t"
178 << eta(0, 2) << "\t" << eta(1, 2) << "\t" << eta(2, 2) << ";";
179
180 os << "\n";
181 }
182
183 void DumpWriter::writeFrame(std::ostream& os) {
184 const int BUFFERSIZE = 2000;
185 const int MINIBUFFERSIZE = 100;
186
187 char tempBuffer[BUFFERSIZE];
188 char writeLine[BUFFERSIZE];
189
190 Quat4d q;
191 Vector3d ji;
192 Vector3d pos;
193 Vector3d vel;
194
195 Molecule* mol;
196 StuntDouble* integrableObject;
197 SimInfo::MoleculeIterator mi;
198 Molecule::IntegrableObjectIterator ii;
199
200 int nTotObjects;
201 nTotObjects = info_->getNGlobalIntegrableObjects();
202
203 #ifndef IS_MPI
204
205
206 os << nTotObjects << "\n";
207
208 writeCommentLine(os, info_->getSnapshotManager()->getCurrentSnapshot());
209
210 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
211
212 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
213 integrableObject = mol->nextIntegrableObject(ii)) {
214
215
216 pos = integrableObject->getPos();
217 vel = integrableObject->getVel();
218
219 sprintf(tempBuffer, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
220 integrableObject->getType().c_str(),
221 pos[0], pos[1], pos[2],
222 vel[0], vel[1], vel[2]);
223
224 strcpy(writeLine, tempBuffer);
225
226 if (integrableObject->isDirectional()) {
227 q = integrableObject->getQ();
228 ji = integrableObject->getJ();
229
230 sprintf(tempBuffer, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
231 q[0], q[1], q[2], q[3],
232 ji[0], ji[1], ji[2]);
233 strcat(writeLine, tempBuffer);
234 } else {
235 strcat(writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
236 }
237
238 os << writeLine;
239
240 }
241 }
242
243 os.flush();
244 #else // is_mpi
245 /*********************************************************************
246 * Documentation? You want DOCUMENTATION?
247 *
248 * Why all the potatoes below?
249 *
250 * To make a long story short, the original version of DumpWriter
251 * worked in the most inefficient way possible. Node 0 would
252 * poke each of the node for an individual atom's formatted data
253 * as node 0 worked its way down the global index. This was particularly
254 * inefficient since the method blocked all processors at every atom
255 * (and did it twice!).
256 *
257 * An intermediate version of DumpWriter could be described from Node
258 * zero's perspective as follows:
259 *
260 * 1) Have 100 of your friends stand in a circle.
261 * 2) When you say go, have all of them start tossing potatoes at
262 * you (one at a time).
263 * 3) Catch the potatoes.
264 *
265 * It was an improvement, but MPI has buffers and caches that could
266 * best be described in this analogy as "potato nets", so there's no
267 * need to block the processors atom-by-atom.
268 *
269 * This new and improved DumpWriter works in an even more efficient
270 * way:
271 *
272 * 1) Have 100 of your friend stand in a circle.
273 * 2) When you say go, have them start tossing 5-pound bags of
274 * potatoes at you.
275 * 3) Once you've caught a friend's bag of potatoes,
276 * toss them a spud to let them know they can toss another bag.
277 *
278 * How's THAT for documentation?
279 *
280 *********************************************************************/
281 const int masterNode = 0;
282
283 int * potatoes;
284 int myPotato;
285 int nProc;
286 int which_node;
287 double atomData[13];
288 int isDirectional;
289 char MPIatomTypeString[MINIBUFFERSIZE];
290 int msgLen; // the length of message actually recieved at master nodes
291 int haveError;
292 MPI_Status istatus;
293 int nCurObj;
294
295 // code to find maximum tag value
296 int * tagub;
297 int flag;
298 int MAXTAG;
299 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
300
301 if (flag) {
302 MAXTAG = *tagub;
303 } else {
304 MAXTAG = 32767;
305 }
306
307 if (worldRank == masterNode) { //master node (node 0) is responsible for writing the dump file
308
309 // Node 0 needs a list of the magic potatoes for each processor;
310
311 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
312 potatoes = new int[nProc];
313
314 //write out the comment lines
315 for(int i = 0; i < nProc; i++) {
316 potatoes[i] = 0;
317 }
318
319
320 os << nTotObjects << "\n";
321 writeCommentLine(os, info_->getSnapshotManager()->getCurrentSnapshot());
322
323 for(int i = 0; i < info_->getNGlobalMolecules(); i++) {
324
325 // Get the Node number which has this atom;
326
327 which_node = info_->getMolToProc(i);
328
329 if (which_node != masterNode) { //current molecule is in slave node
330 if (potatoes[which_node] + 1 >= MAXTAG) {
331 // The potato was going to exceed the maximum value,
332 // so wrap this processor potato back to 0:
333
334 potatoes[which_node] = 0;
335 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0,
336 MPI_COMM_WORLD);
337 }
338
339 myPotato = potatoes[which_node];
340
341 //recieve the number of integrableObject in current molecule
342 MPI_Recv(&nCurObj, 1, MPI_INT, which_node, myPotato,
343 MPI_COMM_WORLD, &istatus);
344 myPotato++;
345
346 for(int l = 0; l < nCurObj; l++) {
347 if (potatoes[which_node] + 2 >= MAXTAG) {
348 // The potato was going to exceed the maximum value,
349 // so wrap this processor potato back to 0:
350
351 potatoes[which_node] = 0;
352 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node,
353 0, MPI_COMM_WORLD);
354 }
355
356 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR,
357 which_node, myPotato, MPI_COMM_WORLD,
358 &istatus);
359
360 myPotato++;
361
362 MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato,
363 MPI_COMM_WORLD, &istatus);
364 myPotato++;
365
366 MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
367
368 if (msgLen == 13)
369 isDirectional = 1;
370 else
371 isDirectional = 0;
372
373 // If we've survived to here, format the line:
374
375 if (!isDirectional) {
376 sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
377 MPIatomTypeString, atomData[0],
378 atomData[1], atomData[2],
379 atomData[3], atomData[4],
380 atomData[5]);
381
382 strcat(writeLine,
383 "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
384 } else {
385 sprintf(writeLine,
386 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
387 MPIatomTypeString,
388 atomData[0],
389 atomData[1],
390 atomData[2],
391 atomData[3],
392 atomData[4],
393 atomData[5],
394 atomData[6],
395 atomData[7],
396 atomData[8],
397 atomData[9],
398 atomData[10],
399 atomData[11],
400 atomData[12]);
401 }
402
403 os << writeLine;
404
405 } // end for(int l =0)
406
407 potatoes[which_node] = myPotato;
408 } else { //master node has current molecule
409
410 mol = info_->getMoleculeByGlobalIndex(i);
411
412 if (mol == NULL) {
413 sprintf(painCave.errMsg, "Molecule not found on node %d!", worldRank);
414 painCave.isFatal = 1;
415 simError();
416 }
417
418 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419 integrableObject = mol->nextIntegrableObject(ii)) {
420
421 pos = integrableObject->getPos();
422 vel = integrableObject->getVel();
423
424 atomData[0] = pos[0];
425 atomData[1] = pos[1];
426 atomData[2] = pos[2];
427
428 atomData[3] = vel[0];
429 atomData[4] = vel[1];
430 atomData[5] = vel[2];
431
432 isDirectional = 0;
433
434 if (integrableObject->isDirectional()) {
435 isDirectional = 1;
436
437 q = integrableObject->getQ();
438 ji = integrableObject->getJ();
439
440 for(int j = 0; j < 6; j++) {
441 atomData[j] = atomData[j];
442 }
443
444 atomData[6] = q[0];
445 atomData[7] = q[1];
446 atomData[8] = q[2];
447 atomData[9] = q[3];
448
449 atomData[10] = ji[0];
450 atomData[11] = ji[1];
451 atomData[12] = ji[2];
452 }
453
454 // If we've survived to here, format the line:
455
456 if (!isDirectional) {
457 sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
458 integrableObject->getType().c_str(), atomData[0],
459 atomData[1], atomData[2],
460 atomData[3], atomData[4],
461 atomData[5]);
462
463 strcat(writeLine,
464 "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
465 } else {
466 sprintf(writeLine,
467 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
468 integrableObject->getType().c_str(),
469 atomData[0],
470 atomData[1],
471 atomData[2],
472 atomData[3],
473 atomData[4],
474 atomData[5],
475 atomData[6],
476 atomData[7],
477 atomData[8],
478 atomData[9],
479 atomData[10],
480 atomData[11],
481 atomData[12]);
482 }
483
484
485 os << writeLine;
486
487 } //end for(iter = integrableObject.begin())
488 }
489 } //end for(i = 0; i < mpiSim->getNmol())
490
491 os.flush();
492
493 sprintf(checkPointMsg, "Sucessfully took a dump.\n");
494 MPIcheckPoint();
495
496 delete [] potatoes;
497 } else {
498
499 // worldRank != 0, so I'm a remote node.
500
501 // Set my magic potato to 0:
502
503 myPotato = 0;
504
505 for(int i = 0; i < info_->getNGlobalMolecules(); i++) {
506
507 // Am I the node which has this integrableObject?
508 int whichNode = info_->getMolToProc(i);
509 if (whichNode == worldRank) {
510 if (myPotato + 1 >= MAXTAG) {
511
512 // The potato was going to exceed the maximum value,
513 // so wrap this processor potato back to 0 (and block until
514 // node 0 says we can go:
515
516 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD,
517 &istatus);
518 }
519
520 mol = info_->getMoleculeByGlobalIndex(i);
521
522
523 nCurObj = mol->getNIntegrableObjects();
524
525 MPI_Send(&nCurObj, 1, MPI_INT, 0, myPotato, MPI_COMM_WORLD);
526 myPotato++;
527
528 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
529 integrableObject = mol->nextIntegrableObject(ii)) {
530
531 if (myPotato + 2 >= MAXTAG) {
532
533 // The potato was going to exceed the maximum value,
534 // so wrap this processor potato back to 0 (and block until
535 // node 0 says we can go:
536
537 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD,
538 &istatus);
539 }
540
541 pos = integrableObject->getPos();
542 vel = integrableObject->getVel();
543
544 atomData[0] = pos[0];
545 atomData[1] = pos[1];
546 atomData[2] = pos[2];
547
548 atomData[3] = vel[0];
549 atomData[4] = vel[1];
550 atomData[5] = vel[2];
551
552 isDirectional = 0;
553
554 if (integrableObject->isDirectional()) {
555 isDirectional = 1;
556
557 q = integrableObject->getQ();
558 ji = integrableObject->getJ();
559
560 atomData[6] = q[0];
561 atomData[7] = q[1];
562 atomData[8] = q[2];
563 atomData[9] = q[3];
564
565 atomData[10] = ji[0];
566 atomData[11] = ji[1];
567 atomData[12] = ji[2];
568 }
569
570 strncpy(MPIatomTypeString, integrableObject->getType().c_str(), MINIBUFFERSIZE);
571
572 // null terminate the std::string before sending (just in case):
573 MPIatomTypeString[MINIBUFFERSIZE - 1] = '\0';
574
575 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
576 myPotato, MPI_COMM_WORLD);
577
578 myPotato++;
579
580 if (isDirectional) {
581 MPI_Send(atomData, 13, MPI_DOUBLE, 0, myPotato,
582 MPI_COMM_WORLD);
583 } else {
584 MPI_Send(atomData, 6, MPI_DOUBLE, 0, myPotato,
585 MPI_COMM_WORLD);
586 }
587
588 myPotato++;
589 }
590
591 }
592
593 }
594 sprintf(checkPointMsg, "Sucessfully took a dump.\n");
595 MPIcheckPoint();
596 }
597
598 #endif // is_mpi
599
600 }
601
602 void DumpWriter::writeDump() {
603 writeFrame(*dumpFile_);
604 }
605
606 void DumpWriter::writeEor() {
607 std::ostream* eorStream;
608
609 #ifdef IS_MPI
610 if (worldRank == 0) {
611 #endif // is_mpi
612
613 eorStream = createOStream(eorFilename_);
614
615 #ifdef IS_MPI
616 }
617 #endif // is_mpi
618
619 writeFrame(*eorStream);
620
621 #ifdef IS_MPI
622 if (worldRank == 0) {
623 #endif // is_mpi
624 delete eorStream;
625
626 #ifdef IS_MPI
627 }
628 #endif // is_mpi
629
630 }
631
632
633 void DumpWriter::writeDumpAndEor() {
634 std::vector<std::streambuf*> buffers;
635 std::ostream* eorStream;
636 #ifdef IS_MPI
637 if (worldRank == 0) {
638 #endif // is_mpi
639
640 buffers.push_back(dumpFile_->rdbuf());
641
642 eorStream = createOStream(eorFilename_);
643
644 buffers.push_back(eorStream->rdbuf());
645
646 #ifdef IS_MPI
647 }
648 #endif // is_mpi
649
650 TeeBuf tbuf(buffers.begin(), buffers.end());
651 std::ostream os(&tbuf);
652
653 writeFrame(os);
654
655 #ifdef IS_MPI
656 if (worldRank == 0) {
657 #endif // is_mpi
658 delete eorStream;
659
660 #ifdef IS_MPI
661 }
662 #endif // is_mpi
663
664 }
665
666 std::ostream* DumpWriter::createOStream(const std::string& filename) {
667
668 std::ostream* newOStream;
669 #ifdef HAVE_LIBZ
670 if (needCompression_) {
671 newOStream = new ogzstream(filename.c_str());
672 } else {
673 newOStream = new std::ofstream(filename.c_str());
674 }
675 #else
676 newOStream = new std::ofstream(filename.c_str());
677 #endif
678 return newOStream;
679 }
680
681 }//end namespace oopse