ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 1129
Committed: Thu Apr 22 03:29:30 2004 UTC (20 years, 2 months ago) by tim
File size: 16286 byte(s)
Log Message: 
fixed two bugs in MPI version of InitfromFile and one unmatch MPI command in DumpWriter

File Contents

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