ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/io/DumpWriter.cpp
Revision: 2314
Committed: Tue Sep 20 21:22:38 2005 UTC (18 years, 9 months ago) by tim
File size: 17523 byte(s)
Log Message:
adding zlib support for DumpWriter

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