ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 951
Committed: Fri Jan 16 21:51:51 2004 UTC (20 years, 5 months ago) by mmeineke
File size: 15480 byte(s)
Log Message: 
fixed a bug where only MPI jobs could write eor files

File Contents

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