ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/io/DumpWriter.cpp
Revision: 1695
Committed: Mon Nov 1 22:52:57 2004 UTC (19 years, 10 months ago) by tim
File size: 18217 byte(s)
Log Message:
Molecule, Atom, DirectionalAtom, RigidBody and StuntDouble classes get compiled

File Contents

# User Rev Content
1 tim 1695 #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 "brains/mpiSimulation.hpp"
13    
14     namespace dWrite{
15     void DieDieDie( void );
16     }
17    
18     using namespace dWrite;
19     #endif //is_mpi
20    
21     #include "io/ReadWrite.hpp"
22     #include "utils/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     dumpFile.open(entry_plug->sampleName.c_str(), ios::out | ios::trunc );
33    
34     if( !dumpFile ){
35    
36     sprintf( painCave.errMsg,
37     "Could not open \"%s\" for dump output.\n",
38     entry_plug->sampleName.c_str());
39     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     finalOut.open( entry_plug->finalName.c_str(), 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.c_str() );
115     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     finalOut.open( entry_plug->finalName.c_str(), 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.c_str() );
148     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     Quat4d q;
229     Vector3d ji;
230     DirectionalAtom* dAtom;
231     Vector3d pos;
232     Vector3d vel;
233    
234     int nTotObjects;
235     StuntDouble* sd;
236     char* molName;
237     vector<StuntDouble*> integrableObjects;
238     vector<StuntDouble*>::iterator iter;
239     nTotObjects = entry_plug->getTotIntegrableObjects();
240     #ifndef IS_MPI
241    
242     for(k = 0; k < outFile.size(); k++){
243     *outFile[k] << nTotObjects << "\n";
244    
245     *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    
254     << entry_plug->Hmat[0][2] << "\t"
255     << entry_plug->Hmat[1][2] << "\t"
256     << entry_plug->Hmat[2][2] << ";";
257    
258     //write out additional parameters, such as chi and eta
259     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
260     }
261    
262     for( i=0; i< entry_plug->n_mol; i++ ){
263    
264     integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
265     molName = (entry_plug->compStamps[entry_plug->molecules[i].getStampID()])->getID();
266    
267     for( iter = integrableObjects.begin();iter != integrableObjects.end(); ++iter){
268     sd = *iter;
269     pos = sd->getPos();
270     vel = sd->getVel();
271    
272     sprintf( tempBuffer,
273     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
274     sd->getType(),
275     pos[0],
276     pos[1],
277     pos[2],
278     vel[0],
279     vel[1],
280     vel[2]);
281     strcpy( writeLine, tempBuffer );
282    
283     if( sd->isDirectional() ){
284    
285     q = sd->getQ();
286     ji = sd->getJ();
287    
288     sprintf( tempBuffer,
289     "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
290     q[0],
291     q[1],
292     q[2],
293     q[3],
294     ji[0],
295     ji[1],
296     ji[2]);
297     strcat( writeLine, tempBuffer );
298     }
299     else
300     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
301    
302     for(k = 0; k < outFile.size(); k++)
303     *outFile[k] << writeLine;
304     }
305    
306     }
307    
308     #else // is_mpi
309    
310     /* code to find maximum tag value */
311    
312     int *tagub, flag, MAXTAG;
313     MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
314     if (flag) {
315     MAXTAG = *tagub;
316     } else {
317     MAXTAG = 32767;
318     }
319    
320     int haveError;
321    
322     MPI_Status istatus;
323     int nCurObj;
324     int *MolToProcMap = mpiSim->getMolToProcMap();
325    
326     // write out header and node 0's coordinates
327    
328     if( worldRank == 0 ){
329    
330     // Node 0 needs a list of the magic potatoes for each processor;
331    
332     nProc = mpiSim->getNProcessors();
333     potatoes = new int[nProc];
334    
335     //write out the comment lines
336     for (i = 0; i < nProc; i++)
337     potatoes[i] = 0;
338    
339     for(k = 0; k < outFile.size(); k++){
340     *outFile[k] << nTotObjects << "\n";
341    
342     *outFile[k] << currentTime << ";\t"
343     << entry_plug->Hmat[0][0] << "\t"
344     << entry_plug->Hmat[1][0] << "\t"
345     << entry_plug->Hmat[2][0] << ";\t"
346    
347     << entry_plug->Hmat[0][1] << "\t"
348     << entry_plug->Hmat[1][1] << "\t"
349     << entry_plug->Hmat[2][1] << ";\t"
350    
351     << entry_plug->Hmat[0][2] << "\t"
352     << entry_plug->Hmat[1][2] << "\t"
353     << entry_plug->Hmat[2][2] << ";";
354    
355     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
356     }
357    
358     currentIndex = 0;
359    
360     for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
361    
362     // Get the Node number which has this atom;
363    
364     which_node = MolToProcMap[i];
365    
366     if (which_node != 0) {
367    
368     if (potatoes[which_node] + 1 >= MAXTAG) {
369     // The potato was going to exceed the maximum value,
370     // so wrap this processor potato back to 0:
371    
372     potatoes[which_node] = 0;
373     MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
374    
375     }
376    
377     myPotato = potatoes[which_node];
378    
379     //recieve the number of integrableObject in current molecule
380     MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
381     myPotato, MPI_COMM_WORLD, &istatus);
382     myPotato++;
383    
384     for(int l = 0; l < nCurObj; l++){
385    
386     if (potatoes[which_node] + 2 >= MAXTAG) {
387     // The potato was going to exceed the maximum value,
388     // so wrap this processor potato back to 0:
389    
390     potatoes[which_node] = 0;
391     MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
392    
393     }
394    
395     MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
396     myPotato, MPI_COMM_WORLD, &istatus);
397    
398     atomTypeString = MPIatomTypeString;
399    
400     myPotato++;
401    
402     MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
403     myPotato++;
404    
405     MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
406    
407     if(msgLen == 13)
408     isDirectional = 1;
409     else
410     isDirectional = 0;
411    
412     // If we've survived to here, format the line:
413    
414     if (!isDirectional) {
415    
416     sprintf( writeLine,
417     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
418     atomTypeString,
419     atomData[0],
420     atomData[1],
421     atomData[2],
422     atomData[3],
423     atomData[4],
424     atomData[5]);
425    
426     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
427    
428     }
429     else {
430    
431     sprintf( writeLine,
432     "%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",
433     atomTypeString,
434     atomData[0],
435     atomData[1],
436     atomData[2],
437     atomData[3],
438     atomData[4],
439     atomData[5],
440     atomData[6],
441     atomData[7],
442     atomData[8],
443     atomData[9],
444     atomData[10],
445     atomData[11],
446     atomData[12]);
447    
448     }
449    
450     for(k = 0; k < outFile.size(); k++)
451     *outFile[k] << writeLine;
452    
453     }// end for(int l =0)
454     potatoes[which_node] = myPotato;
455    
456     }
457     else {
458    
459     haveError = 0;
460    
461     local_index = indexArray[currentIndex].first;
462    
463     integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects();
464    
465     for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){
466     sd = *iter;
467     atomTypeString = sd->getType();
468    
469     pos = sd->getPos();
470     vel = sd->getVel();
471    
472     atomData[0] = pos[0];
473     atomData[1] = pos[1];
474     atomData[2] = pos[2];
475    
476     atomData[3] = vel[0];
477     atomData[4] = vel[1];
478     atomData[5] = vel[2];
479    
480     isDirectional = 0;
481    
482     if( sd->isDirectional() ){
483    
484     isDirectional = 1;
485    
486     q = sd->getQ();
487     ji = sd->getJ();
488    
489     for (int j = 0; j < 6 ; j++)
490     atomData[j] = atomData[j];
491    
492     atomData[6] = q[0];
493     atomData[7] = q[1];
494     atomData[8] = q[2];
495     atomData[9] = q[3];
496    
497     atomData[10] = ji[0];
498     atomData[11] = ji[1];
499     atomData[12] = ji[2];
500     }
501    
502     // If we've survived to here, format the line:
503    
504     if (!isDirectional) {
505    
506     sprintf( writeLine,
507     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
508     atomTypeString,
509     atomData[0],
510     atomData[1],
511     atomData[2],
512     atomData[3],
513     atomData[4],
514     atomData[5]);
515    
516     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
517    
518     }
519     else {
520    
521     sprintf( writeLine,
522     "%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",
523     atomTypeString,
524     atomData[0],
525     atomData[1],
526     atomData[2],
527     atomData[3],
528     atomData[4],
529     atomData[5],
530     atomData[6],
531     atomData[7],
532     atomData[8],
533     atomData[9],
534     atomData[10],
535     atomData[11],
536     atomData[12]);
537    
538     }
539    
540     for(k = 0; k < outFile.size(); k++)
541     *outFile[k] << writeLine;
542    
543    
544     }//end for(iter = integrableObject.begin())
545    
546     currentIndex++;
547     }
548    
549     }//end for(i = 0; i < mpiSim->getNmol())
550    
551     for(k = 0; k < outFile.size(); k++)
552     outFile[k]->flush();
553    
554     sprintf( checkPointMsg,
555     "Sucessfully took a dump.\n");
556    
557     MPIcheckPoint();
558    
559     delete[] potatoes;
560    
561     } else {
562    
563     // worldRank != 0, so I'm a remote node.
564    
565     // Set my magic potato to 0:
566    
567     myPotato = 0;
568     currentIndex = 0;
569    
570     for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
571    
572     // Am I the node which has this integrableObject?
573    
574     if (MolToProcMap[i] == worldRank) {
575    
576    
577     if (myPotato + 1 >= MAXTAG) {
578    
579     // The potato was going to exceed the maximum value,
580     // so wrap this processor potato back to 0 (and block until
581     // node 0 says we can go:
582    
583     MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
584    
585     }
586    
587     local_index = indexArray[currentIndex].first;
588     integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
589    
590     nCurObj = integrableObjects.size();
591    
592     MPI_Send(&nCurObj, 1, MPI_INT, 0,
593     myPotato, MPI_COMM_WORLD);
594     myPotato++;
595    
596     for( iter = integrableObjects.begin(); iter != integrableObjects.end(); iter++){
597    
598     if (myPotato + 2 >= MAXTAG) {
599    
600     // The potato was going to exceed the maximum value,
601     // so wrap this processor potato back to 0 (and block until
602     // node 0 says we can go:
603    
604     MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
605    
606     }
607    
608     sd = *iter;
609    
610     atomTypeString = sd->getType();
611    
612     pos = sd->getPos();
613     vel = sd->getVel();
614    
615     atomData[0] = pos[0];
616     atomData[1] = pos[1];
617     atomData[2] = pos[2];
618    
619     atomData[3] = vel[0];
620     atomData[4] = vel[1];
621     atomData[5] = vel[2];
622    
623     isDirectional = 0;
624    
625     if( sd->isDirectional() ){
626    
627     isDirectional = 1;
628    
629     q = sd->getQ();
630     ji = sd->getJ();
631    
632    
633     atomData[6] = q[0];
634     atomData[7] = q[1];
635     atomData[8] = q[2];
636     atomData[9] = q[3];
637    
638     atomData[10] = ji[0];
639     atomData[11] = ji[1];
640     atomData[12] = ji[2];
641     }
642    
643    
644     strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
645    
646     // null terminate the string before sending (just in case):
647     MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
648    
649     MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
650     myPotato, MPI_COMM_WORLD);
651    
652     myPotato++;
653    
654     if (isDirectional) {
655    
656     MPI_Send(atomData, 13, MPI_DOUBLE, 0,
657     myPotato, MPI_COMM_WORLD);
658    
659     } else {
660    
661     MPI_Send(atomData, 6, MPI_DOUBLE, 0,
662     myPotato, MPI_COMM_WORLD);
663     }
664    
665     myPotato++;
666    
667     }
668    
669     currentIndex++;
670    
671     }
672    
673     }
674    
675     sprintf( checkPointMsg,
676     "Sucessfully took a dump.\n");
677     MPIcheckPoint();
678    
679     }
680    
681    
682    
683     #endif // is_mpi
684     }
685    
686     #ifdef IS_MPI
687    
688     // a couple of functions to let us escape the write loop
689    
690     void dWrite::DieDieDie( void ){
691    
692     MPI_Finalize();
693     exit (0);
694     }
695    
696     #endif //is_mpi