ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/io/DumpWriter.cpp
Revision: 1727
Committed: Thu Nov 11 16:41:58 2004 UTC (19 years, 8 months ago) by tim
File size: 19696 byte(s)
Log Message:
add Snapshot.cpp, remove useless mpiSimulation

File Contents

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