ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 1231
Committed: Thu Jun 3 21:06:51 2004 UTC (20 years ago) by tim
File size: 17403 byte(s)
Log Message: 
fixed a bug which only writes out the first atom of a molecule

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 tim 1152 double atomData[13];
220 gezelter 907 int isDirectional;
221     char* atomTypeString;
222 gezelter 910 char MPIatomTypeString[MINIBUFFERSIZE];
223 tim 1108 int nObjects;
224 tim 1152 int msgLen; // the length of message actually recieved at master nodes
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 tim 1231
298     for(k = 0; k < outFile.size(); k++)
299     *outFile[k] << writeLine;
300 mmeineke 377 }
301 tim 837
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 gezelter 1203 nProc = mpiSim->getNProcessors();
329 gezelter 916 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 gezelter 1203 for (i = 0 ; i < mpiSim->getNMolGlobal(); 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 1152 if (potatoes[which_node] + 2 >= MAXTAG) {
383 tim 1108 // 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 tim 1152 MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
399 tim 1108 myPotato++;
400 gezelter 907
401 tim 1152 MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
402 tim 1108
403 tim 1152 if(msgLen == 13)
404     isDirectional = 1;
405     else
406     isDirectional = 0;
407 tim 1231
408     // If we've survived to here, format the line:
409 tim 1152
410 tim 1231 if (!isDirectional) {
411    
412     sprintf( writeLine,
413     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
414     atomTypeString,
415     atomData[0],
416     atomData[1],
417     atomData[2],
418     atomData[3],
419     atomData[4],
420     atomData[5]);
421    
422     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
423    
424     }
425     else {
426    
427     sprintf( writeLine,
428     "%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",
429     atomTypeString,
430     atomData[0],
431     atomData[1],
432     atomData[2],
433     atomData[3],
434     atomData[4],
435     atomData[5],
436     atomData[6],
437     atomData[7],
438     atomData[8],
439     atomData[9],
440     atomData[10],
441     atomData[11],
442     atomData[12]);
443    
444     }
445    
446     for(k = 0; k < outFile.size(); k++)
447     *outFile[k] << writeLine;
448    
449     }// end for(int l =0)
450 gezelter 916 potatoes[which_node] = myPotato;
451 gezelter 907
452 tim 1231 }
453     else {
454 gezelter 907
455 tim 934 haveError = 0;
456 gezelter 916
457 tim 1108 local_index = indexArray[currentIndex].first;
458 tim 837
459 tim 1108 integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects();
460 gezelter 907
461 tim 1108 for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){
462     sd = *iter;
463     atomTypeString = sd->getType();
464    
465     sd->getPos(pos);
466     sd->getVel(vel);
467    
468 tim 1152 atomData[0] = pos[0];
469     atomData[1] = pos[1];
470     atomData[2] = pos[2];
471 tim 837
472 tim 1152 atomData[3] = vel[0];
473     atomData[4] = vel[1];
474     atomData[5] = vel[2];
475 tim 1108
476     isDirectional = 0;
477 gezelter 916
478 tim 1108 if( sd->isDirectional() ){
479    
480     isDirectional = 1;
481    
482     sd->getQ( q );
483     sd->getJ( ji );
484    
485     for (int j = 0; j < 6 ; j++)
486 tim 1152 atomData[j] = atomData[j];
487 tim 1108
488 tim 1152 atomData[6] = q[0];
489     atomData[7] = q[1];
490     atomData[8] = q[2];
491     atomData[9] = q[3];
492 tim 1108
493 tim 1152 atomData[10] = ji[0];
494     atomData[11] = ji[1];
495     atomData[12] = ji[2];
496 tim 1108 }
497 gezelter 907
498 tim 1231 // If we've survived to here, format the line:
499    
500     if (!isDirectional) {
501    
502     sprintf( writeLine,
503     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
504     atomTypeString,
505     atomData[0],
506     atomData[1],
507     atomData[2],
508     atomData[3],
509     atomData[4],
510     atomData[5]);
511    
512     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
513    
514     }
515     else {
516    
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    
536     for(k = 0; k < outFile.size(); k++)
537     *outFile[k] << writeLine;
538    
539    
540     }//end for(iter = integrableObject.begin())
541 gezelter 916
542 tim 1108 currentIndex++;
543 tim 926 }
544 tim 1231
545     }//end for(i = 0; i < mpiSim->getNmol())
546 tim 926
547 tim 934 for(k = 0; k < outFile.size(); k++)
548     outFile[k]->flush();
549    
550 gezelter 907 sprintf( checkPointMsg,
551     "Sucessfully took a dump.\n");
552 chuckv 949
553 gezelter 907 MPIcheckPoint();
554 chuckv 949
555 tim 919 delete[] potatoes;
556 chuckv 949
557 gezelter 415 } else {
558 tim 837
559 gezelter 907 // worldRank != 0, so I'm a remote node.
560 gezelter 916
561     // Set my magic potato to 0:
562    
563     myPotato = 0;
564 tim 929 currentIndex = 0;
565 gezelter 907
566 gezelter 1203 for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
567 gezelter 907
568 tim 1108 // Am I the node which has this integrableObject?
569 gezelter 907
570 tim 1108 if (MolToProcMap[i] == worldRank) {
571 tim 837
572 tim 1108
573     if (myPotato + 1 >= MAXTAG) {
574 chuckv 949
575 gezelter 916 // The potato was going to exceed the maximum value,
576     // so wrap this processor potato back to 0 (and block until
577     // node 0 says we can go:
578 chuckv 949
579 gezelter 916 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
580    
581     }
582 chuckv 949
583 tim 1108 local_index = indexArray[currentIndex].first;
584     integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
585    
586     nCurObj = integrableObjects.size();
587    
588     MPI_Send(&nCurObj, 1, MPI_INT, 0,
589     myPotato, MPI_COMM_WORLD);
590 tim 1129 myPotato++;
591 tim 1108
592     for( iter = integrableObjects.begin(); iter != integrableObjects.end(); iter++){
593    
594 tim 1152 if (myPotato + 2 >= MAXTAG) {
595 tim 1108
596     // The potato was going to exceed the maximum value,
597     // so wrap this processor potato back to 0 (and block until
598     // node 0 says we can go:
599    
600     MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
601    
602     }
603    
604     sd = *iter;
605    
606     atomTypeString = sd->getType();
607    
608     sd->getPos(pos);
609     sd->getVel(vel);
610    
611 tim 1152 atomData[0] = pos[0];
612     atomData[1] = pos[1];
613     atomData[2] = pos[2];
614 tim 1108
615 tim 1152 atomData[3] = vel[0];
616     atomData[4] = vel[1];
617     atomData[5] = vel[2];
618 tim 1108
619     isDirectional = 0;
620    
621     if( sd->isDirectional() ){
622    
623     isDirectional = 1;
624 tim 920
625 tim 1108 sd->getQ( q );
626     sd->getJ( ji );
627    
628    
629 tim 1152 atomData[6] = q[0];
630     atomData[7] = q[1];
631     atomData[8] = q[2];
632     atomData[9] = q[3];
633 tim 1108
634 tim 1152 atomData[10] = ji[0];
635     atomData[11] = ji[1];
636     atomData[12] = ji[2];
637 tim 1108 }
638 tim 837
639 tim 1108
640     strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
641 tim 837
642 tim 1108 // null terminate the string before sending (just in case):
643     MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
644 mmeineke 377
645 tim 1108 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
646     myPotato, MPI_COMM_WORLD);
647 gezelter 907
648 tim 1108 myPotato++;
649 gezelter 907
650 tim 1108 if (isDirectional) {
651 tim 837
652 tim 1152 MPI_Send(atomData, 13, MPI_DOUBLE, 0,
653 tim 1108 myPotato, MPI_COMM_WORLD);
654    
655     } else {
656 tim 837
657 tim 1152 MPI_Send(atomData, 6, MPI_DOUBLE, 0,
658 tim 1108 myPotato, MPI_COMM_WORLD);
659     }
660 tim 837
661 tim 1108 myPotato++;
662 gezelter 916
663 tim 1108 }
664 gezelter 907
665 tim 1108 currentIndex++;
666 gezelter 907
667     }
668 tim 1108
669 mmeineke 377 }
670 mmeineke 440
671 gezelter 907 sprintf( checkPointMsg,
672 gezelter 916 "Sucessfully took a dump.\n");
673 tim 1129 MPIcheckPoint();
674 gezelter 907
675 tim 1129 }
676    
677    
678 gezelter 907
679 mmeineke 377 #endif // is_mpi
680     }
681 mmeineke 440
682     #ifdef IS_MPI
683    
684     // a couple of functions to let us escape the write loop
685    
686 gezelter 907 void dWrite::DieDieDie( void ){
687 tim 837
688 mmeineke 440 MPI_Finalize();
689     exit (0);
690     }
691    
692     #endif //is_mpi