ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/libmdtools/DumpWriter.cpp
Revision: 1421
Committed: Tue Jul 27 18:44:49 2004 UTC (19 years, 11 months ago) by tim
File size: 17466 byte(s)
Log Message: 
change parameter of calling sprintf

File Contents

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