ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 1000
Committed: Fri Jan 30 21:47:22 2004 UTC (20 years, 5 months ago) by tim
File size: 14799 byte(s)
Log Message: 
using class  Minimizer1D derived from MinimizerBase instead of a functor to do line seach

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     }
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 gezelter 916
225 mmeineke 787 #else //is_mpi
226     int nAtoms = entry_plug->n_atoms;
227     #endif //is_mpi
228    
229 mmeineke 377 double q[4];
230     DirectionalAtom* dAtom;
231     Atom** atoms = entry_plug->atoms;
232 mmeineke 670 double pos[3], vel[3];
233 tim 837
234 mmeineke 377 #ifndef IS_MPI
235 tim 934
236     for(k = 0; k < outFile.size(); k++){
237     *outFile[k] << nAtoms << "\n";
238 tim 837
239 tim 934 *outFile[k] << currentTime << ";\t"
240     << entry_plug->Hmat[0][0] << "\t"
241     << entry_plug->Hmat[1][0] << "\t"
242     << entry_plug->Hmat[2][0] << ";\t"
243    
244     << entry_plug->Hmat[0][1] << "\t"
245     << entry_plug->Hmat[1][1] << "\t"
246     << entry_plug->Hmat[2][1] << ";\t"
247 tim 837
248 tim 934 << entry_plug->Hmat[0][2] << "\t"
249     << entry_plug->Hmat[1][2] << "\t"
250     << entry_plug->Hmat[2][2] << ";";
251 mmeineke 572
252 tim 934 //write out additional parameters, such as chi and eta
253     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
254     }
255    
256 mmeineke 377 for( i=0; i<nAtoms; i++ ){
257 tim 837
258 mmeineke 670 atoms[i]->getPos(pos);
259     atoms[i]->getVel(vel);
260 mmeineke 377
261     sprintf( tempBuffer,
262     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
263     atoms[i]->getType(),
264 mmeineke 670 pos[0],
265     pos[1],
266     pos[2],
267     vel[0],
268     vel[1],
269     vel[2]);
270 mmeineke 377 strcpy( writeLine, tempBuffer );
271    
272     if( atoms[i]->isDirectional() ){
273 tim 837
274 mmeineke 377 dAtom = (DirectionalAtom *)atoms[i];
275     dAtom->getQ( q );
276 tim 837
277 mmeineke 377 sprintf( tempBuffer,
278     "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
279     q[0],
280     q[1],
281     q[2],
282     q[3],
283     dAtom->getJx(),
284     dAtom->getJy(),
285     dAtom->getJz());
286     strcat( writeLine, tempBuffer );
287     }
288     else
289     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
290 tim 837
291 tim 934 for(k = 0; k < outFile.size(); k++)
292     *outFile[k] << writeLine;
293 mmeineke 377 }
294    
295     #else // is_mpi
296 gezelter 416
297 chuckv 913 /* code to find maximum tag value */
298 tim 919
299 tim 920 int *tagub, flag, MAXTAG;
300 chuckv 913 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
301     if (flag) {
302 tim 920 MAXTAG = *tagub;
303 chuckv 913 } else {
304     MAXTAG = 32767;
305 gezelter 916 }
306 mmeineke 440
307     int haveError;
308    
309 mmeineke 447 MPI_Status istatus;
310 gezelter 416 int *AtomToProcMap = mpiSim->getAtomToProcMap();
311 tim 837
312 mmeineke 377 // write out header and node 0's coordinates
313 tim 837
314 mmeineke 377 if( worldRank == 0 ){
315 gezelter 916
316     // Node 0 needs a list of the magic potatoes for each processor;
317    
318     nProc = mpiSim->getNumberProcessors();
319     potatoes = new int[nProc];
320    
321 tim 934 //write out the comment lines
322 gezelter 916 for (i = 0; i < nProc; i++)
323     potatoes[i] = 0;
324    
325 tim 934 for(k = 0; k < outFile.size(); k++){
326     *outFile[k] << mpiSim->getTotAtoms() << "\n";
327 tim 837
328 tim 934 *outFile[k] << currentTime << ";\t"
329     << entry_plug->Hmat[0][0] << "\t"
330     << entry_plug->Hmat[1][0] << "\t"
331     << entry_plug->Hmat[2][0] << ";\t"
332 tim 837
333 tim 934 << entry_plug->Hmat[0][1] << "\t"
334     << entry_plug->Hmat[1][1] << "\t"
335     << entry_plug->Hmat[2][1] << ";\t"
336 tim 837
337 tim 934 << entry_plug->Hmat[0][2] << "\t"
338     << entry_plug->Hmat[1][2] << "\t"
339     << entry_plug->Hmat[2][2] << ";";
340    
341     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
342     }
343 tim 837
344 tim 929 currentIndex = 0;
345 tim 934
346 gezelter 416 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
347 chuckv 913
348 gezelter 417 // Get the Node number which has this atom;
349 chuckv 913
350 tim 837 which_node = AtomToProcMap[i];
351 chuckv 913
352 gezelter 907 if (which_node != 0) {
353 gezelter 916
354     if (potatoes[which_node] + 3 >= MAXTAG) {
355     // The potato was going to exceed the maximum value,
356     // so wrap this processor potato back to 0:
357    
358     potatoes[which_node] = 0;
359     MPI_Send(0, 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
360    
361     }
362    
363     myPotato = potatoes[which_node];
364 gezelter 907
365 gezelter 910 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
366 gezelter 916 myPotato, MPI_COMM_WORLD, &istatus);
367 gezelter 907
368 tim 920 atomTypeString = MPIatomTypeString;
369    
370 gezelter 916 myPotato++;
371    
372 gezelter 907 MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
373 gezelter 916 myPotato, MPI_COMM_WORLD, &istatus);
374    
375     myPotato++;
376 gezelter 907
377 gezelter 916 if (isDirectional) {
378     MPI_Recv(atomData13, 13, MPI_DOUBLE, which_node,
379     myPotato, MPI_COMM_WORLD, &istatus);
380     } else {
381     MPI_Recv(atomData6, 6, MPI_DOUBLE, which_node,
382     myPotato, MPI_COMM_WORLD, &istatus);
383 gezelter 907 }
384 gezelter 916
385     myPotato++;
386     potatoes[which_node] = myPotato;
387 gezelter 907
388     } else {
389    
390 tim 934 haveError = 0;
391 chuckv 949 which_atom = i;
392 gezelter 916
393 chuckv 949 local_index = indexArray[currentIndex].first;
394 gezelter 916
395 chuckv 949 if (which_atom == indexArray[currentIndex].second) {
396    
397     atomTypeString = atoms[local_index]->getType();
398    
399     atoms[local_index]->getPos(pos);
400     atoms[local_index]->getVel(vel);
401    
402 gezelter 916 atomData6[0] = pos[0];
403     atomData6[1] = pos[1];
404     atomData6[2] = pos[2];
405 tim 837
406 gezelter 916 atomData6[3] = vel[0];
407     atomData6[4] = vel[1];
408     atomData6[5] = vel[2];
409 gezelter 907
410     isDirectional = 0;
411    
412 chuckv 436 if( atoms[local_index]->isDirectional() ){
413 tim 837
414 gezelter 907 isDirectional = 1;
415    
416 chuckv 436 dAtom = (DirectionalAtom *)atoms[local_index];
417     dAtom->getQ( q );
418 gezelter 916
419     for (int j = 0; j < 6 ; j++)
420     atomData13[j] = atomData6[j];
421 gezelter 907
422 gezelter 916 atomData13[6] = q[0];
423     atomData13[7] = q[1];
424     atomData13[8] = q[2];
425     atomData13[9] = q[3];
426    
427     atomData13[10] = dAtom->getJx();
428     atomData13[11] = dAtom->getJy();
429     atomData13[12] = dAtom->getJz();
430 gezelter 907 }
431 gezelter 916
432 gezelter 907 } else {
433 chuckv 949 sprintf(painCave.errMsg,
434     "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
435     which_atom, worldRank, currentIndex, local_index );
436     haveError= 1;
437     simError();
438     }
439 gezelter 916
440 tim 934 if(haveError) DieDieDie();
441 gezelter 916
442 chuckv 949 currentIndex++;
443 tim 926 }
444     // If we've survived to here, format the line:
445    
446     if (!isDirectional) {
447    
448 tim 934 sprintf( writeLine,
449 chuckv 949 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
450     atomTypeString,
451     atomData6[0],
452     atomData6[1],
453     atomData6[2],
454     atomData6[3],
455     atomData6[4],
456     atomData6[5]);
457 tim 926
458 chuckv 949 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
459    
460 tim 926 } else {
461    
462 chuckv 949 sprintf( writeLine,
463     "%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",
464     atomTypeString,
465     atomData13[0],
466     atomData13[1],
467     atomData13[2],
468     atomData13[3],
469     atomData13[4],
470     atomData13[5],
471     atomData13[6],
472     atomData13[7],
473     atomData13[8],
474     atomData13[9],
475     atomData13[10],
476     atomData13[11],
477     atomData13[12]);
478 gezelter 907
479     }
480 tim 926
481 tim 934 for(k = 0; k < outFile.size(); k++)
482     *outFile[k] << writeLine;
483 mmeineke 377 }
484 tim 926
485 tim 934 for(k = 0; k < outFile.size(); k++)
486     outFile[k]->flush();
487    
488 gezelter 907 sprintf( checkPointMsg,
489     "Sucessfully took a dump.\n");
490 chuckv 949
491 gezelter 907 MPIcheckPoint();
492 chuckv 949
493 tim 919 delete[] potatoes;
494 chuckv 949
495 gezelter 415 } else {
496 tim 837
497 gezelter 907 // worldRank != 0, so I'm a remote node.
498 gezelter 916
499     // Set my magic potato to 0:
500    
501     myPotato = 0;
502 tim 929 currentIndex = 0;
503 gezelter 907
504     for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
505    
506     // Am I the node which has this atom?
507    
508     if (AtomToProcMap[i] == worldRank) {
509 tim 837
510 gezelter 916 if (myPotato + 3 >= MAXTAG) {
511 chuckv 949
512 gezelter 916 // The potato was going to exceed the maximum value,
513     // so wrap this processor potato back to 0 (and block until
514     // node 0 says we can go:
515 chuckv 949
516 gezelter 916 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
517    
518     }
519 tim 919 which_atom = i;
520 chuckv 949
521     local_index = indexArray[currentIndex].first;
522 tim 920
523 chuckv 949 if (which_atom == indexArray[currentIndex].second) {
524 gezelter 907
525     atomTypeString = atoms[local_index]->getType();
526 chuckv 949
527     atoms[local_index]->getPos(pos);
528     atoms[local_index]->getVel(vel);
529    
530 gezelter 916 atomData6[0] = pos[0];
531     atomData6[1] = pos[1];
532     atomData6[2] = pos[2];
533 tim 837
534 gezelter 916 atomData6[3] = vel[0];
535     atomData6[4] = vel[1];
536     atomData6[5] = vel[2];
537 gezelter 907
538     isDirectional = 0;
539 tim 837
540 gezelter 907 if( atoms[local_index]->isDirectional() ){
541 mmeineke 377
542 gezelter 907 isDirectional = 1;
543    
544     dAtom = (DirectionalAtom *)atoms[local_index];
545     dAtom->getQ( q );
546    
547 gezelter 916 for (int j = 0; j < 6 ; j++)
548     atomData13[j] = atomData6[j];
549    
550     atomData13[6] = q[0];
551     atomData13[7] = q[1];
552     atomData13[8] = q[2];
553     atomData13[9] = q[3];
554 tim 934
555 gezelter 916 atomData13[10] = dAtom->getJx();
556     atomData13[11] = dAtom->getJy();
557     atomData13[12] = dAtom->getJz();
558 gezelter 907 }
559 tim 837
560 gezelter 907 } else {
561 chuckv 949 sprintf(painCave.errMsg,
562     "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
563     which_atom, worldRank, currentIndex, local_index );
564     haveError= 1;
565     simError();
566     }
567    
568 gezelter 916 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
569 tim 837
570 gezelter 916 // null terminate the string before sending (just in case):
571     MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
572 tim 837
573 gezelter 910 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
574 tim 934 myPotato, MPI_COMM_WORLD);
575 gezelter 907
576 gezelter 916 myPotato++;
577    
578 gezelter 907 MPI_Send(&isDirectional, 1, MPI_INT, 0,
579 tim 934 myPotato, MPI_COMM_WORLD);
580 gezelter 907
581 gezelter 916 myPotato++;
582    
583 gezelter 907 if (isDirectional) {
584    
585 gezelter 916 MPI_Send(atomData13, 13, MPI_DOUBLE, 0,
586     myPotato, MPI_COMM_WORLD);
587 gezelter 907
588 gezelter 916 } else {
589    
590     MPI_Send(atomData6, 6, MPI_DOUBLE, 0,
591     myPotato, MPI_COMM_WORLD);
592 gezelter 907 }
593 gezelter 916
594 tim 929 myPotato++;
595     currentIndex++;
596 mmeineke 377 }
597 gezelter 907 }
598 mmeineke 440
599 gezelter 907 sprintf( checkPointMsg,
600 gezelter 916 "Sucessfully took a dump.\n");
601 gezelter 907 MPIcheckPoint();
602    
603 gezelter 916 }
604 gezelter 907
605 mmeineke 377 #endif // is_mpi
606     }
607 mmeineke 440
608     #ifdef IS_MPI
609    
610     // a couple of functions to let us escape the write loop
611    
612 gezelter 907 void dWrite::DieDieDie( void ){
613 tim 837
614 mmeineke 440 MPI_Finalize();
615     exit (0);
616     }
617    
618     #endif //is_mpi