ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 936
Committed: Tue Jan 13 20:35:25 2004 UTC (20 years, 5 months ago) by tim
File size: 13504 byte(s)
Log Message:
open and close the eor file whenever it is used instead of rewinding it

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    
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     Atom** atoms = entry_plug->atoms;
92    
93     indexArray.clear();
94    
95     for(int i = 0; i < mpiSim->getMyNlocal();i++)
96     indexArray.push_back(make_pair(i, atoms[i]->getGlobalIndex()));
97    
98     sort(indexArray.begin(), indexArray.end(), indexSortingCriterion);
99     }
100     #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 936
110     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 934 }
119     #endif // is_mpi
120    
121     fileStreams.push_back(&finalOut);
122     fileStreams.push_back(&dumpFile);
123    
124     writeFrame(fileStreams, currentTime);
125 tim 936
126     #ifdef IS_MPI
127     finalOut.close();
128     #endif
129 tim 929
130     }
131    
132     void DumpWriter::writeFinal(double currentTime){
133    
134 tim 936 ofstream finalOut;
135 tim 934 vector<ofstream*> fileStreams;
136    
137 tim 929 #ifdef IS_MPI
138     if(worldRank == 0 ){
139 tim 936
140     finalOut.open( entry_plug->finalName, ios::out | ios::trunc );
141    
142     if( !finalOut ){
143     sprintf( painCave.errMsg,
144     "Could not open \"%s\" for final dump output.\n",
145     entry_plug->finalName );
146     painCave.isFatal = 1;
147     simError();
148     }
149    
150 tim 934 }
151 tim 929 #endif // is_mpi
152    
153 tim 934 fileStreams.push_back(&finalOut);
154     writeFrame(fileStreams, currentTime);
155 tim 936
156     #ifdef IS_MPI
157     finalOut.close();
158     #endif
159 tim 929
160     }
161    
162 tim 934 void DumpWriter::writeFrame( vector<ofstream*>& outFile, double currentTime ){
163 tim 929
164 mmeineke 377 const int BUFFERSIZE = 2000;
165 gezelter 912 const int MINIBUFFERSIZE = 100;
166 gezelter 907
167 tim 936 char tempBuffer[BUFFERSIZE];
168 mmeineke 377 char writeLine[BUFFERSIZE];
169    
170 tim 934 int i, k;
171 gezelter 916
172 mmeineke 787 #ifdef IS_MPI
173 gezelter 916
174     int *potatoes;
175     int myPotato;
176    
177     int nProc;
178 tim 929 int j, which_node, done, which_atom, local_index, currentIndex;
179 gezelter 916 double atomData6[6];
180     double atomData13[13];
181 gezelter 907 int isDirectional;
182     char* atomTypeString;
183 gezelter 910 char MPIatomTypeString[MINIBUFFERSIZE];
184 gezelter 916
185 mmeineke 787 #else //is_mpi
186     int nAtoms = entry_plug->n_atoms;
187     #endif //is_mpi
188    
189 mmeineke 377 double q[4];
190     DirectionalAtom* dAtom;
191     Atom** atoms = entry_plug->atoms;
192 mmeineke 670 double pos[3], vel[3];
193 tim 837
194 mmeineke 377 #ifndef IS_MPI
195 tim 934
196     for(k = 0; k < outFile.size(); k++){
197     *outFile[k] << nAtoms << "\n";
198 tim 837
199 tim 934 *outFile[k] << currentTime << ";\t"
200     << entry_plug->Hmat[0][0] << "\t"
201     << entry_plug->Hmat[1][0] << "\t"
202     << entry_plug->Hmat[2][0] << ";\t"
203    
204     << entry_plug->Hmat[0][1] << "\t"
205     << entry_plug->Hmat[1][1] << "\t"
206     << entry_plug->Hmat[2][1] << ";\t"
207 tim 837
208 tim 934 << entry_plug->Hmat[0][2] << "\t"
209     << entry_plug->Hmat[1][2] << "\t"
210     << entry_plug->Hmat[2][2] << ";";
211 mmeineke 572
212 tim 934 //write out additional parameters, such as chi and eta
213     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
214     }
215    
216 mmeineke 377 for( i=0; i<nAtoms; i++ ){
217 tim 837
218 mmeineke 670 atoms[i]->getPos(pos);
219     atoms[i]->getVel(vel);
220 mmeineke 377
221     sprintf( tempBuffer,
222     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
223     atoms[i]->getType(),
224 mmeineke 670 pos[0],
225     pos[1],
226     pos[2],
227     vel[0],
228     vel[1],
229     vel[2]);
230 mmeineke 377 strcpy( writeLine, tempBuffer );
231    
232     if( atoms[i]->isDirectional() ){
233 tim 837
234 mmeineke 377 dAtom = (DirectionalAtom *)atoms[i];
235     dAtom->getQ( q );
236 tim 837
237 mmeineke 377 sprintf( tempBuffer,
238     "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
239     q[0],
240     q[1],
241     q[2],
242     q[3],
243     dAtom->getJx(),
244     dAtom->getJy(),
245     dAtom->getJz());
246     strcat( writeLine, tempBuffer );
247     }
248     else
249     strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
250 tim 837
251 tim 934 for(k = 0; k < outFile.size(); k++)
252     *outFile[k] << writeLine;
253 mmeineke 377 }
254    
255     #else // is_mpi
256 gezelter 416
257 chuckv 913 /* code to find maximum tag value */
258 tim 919
259 tim 920 int *tagub, flag, MAXTAG;
260 chuckv 913 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
261     if (flag) {
262 tim 920 MAXTAG = *tagub;
263 chuckv 913 } else {
264     MAXTAG = 32767;
265 gezelter 916 }
266 mmeineke 440
267     int haveError;
268    
269 mmeineke 447 MPI_Status istatus;
270 gezelter 416 int *AtomToProcMap = mpiSim->getAtomToProcMap();
271 tim 837
272 mmeineke 377 // write out header and node 0's coordinates
273 tim 837
274 mmeineke 377 if( worldRank == 0 ){
275 gezelter 916
276     // Node 0 needs a list of the magic potatoes for each processor;
277    
278     nProc = mpiSim->getNumberProcessors();
279     potatoes = new int[nProc];
280    
281 tim 934 //write out the comment lines
282 gezelter 916 for (i = 0; i < nProc; i++)
283     potatoes[i] = 0;
284    
285 tim 934 for(k = 0; k < outFile.size(); k++){
286     *outFile[k] << mpiSim->getTotAtoms() << "\n";
287 tim 837
288 tim 934 *outFile[k] << currentTime << ";\t"
289     << entry_plug->Hmat[0][0] << "\t"
290     << entry_plug->Hmat[1][0] << "\t"
291     << entry_plug->Hmat[2][0] << ";\t"
292 tim 837
293 tim 934 << entry_plug->Hmat[0][1] << "\t"
294     << entry_plug->Hmat[1][1] << "\t"
295     << entry_plug->Hmat[2][1] << ";\t"
296 tim 837
297 tim 934 << entry_plug->Hmat[0][2] << "\t"
298     << entry_plug->Hmat[1][2] << "\t"
299     << entry_plug->Hmat[2][2] << ";";
300    
301     *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
302     }
303 tim 837
304 tim 929 currentIndex = 0;
305 tim 934
306 gezelter 416 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
307 chuckv 913
308 gezelter 417 // Get the Node number which has this atom;
309 chuckv 913
310 tim 837 which_node = AtomToProcMap[i];
311 chuckv 913
312 gezelter 907 if (which_node != 0) {
313 gezelter 916
314     if (potatoes[which_node] + 3 >= MAXTAG) {
315     // The potato was going to exceed the maximum value,
316     // so wrap this processor potato back to 0:
317    
318     potatoes[which_node] = 0;
319     MPI_Send(0, 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
320    
321     }
322    
323     myPotato = potatoes[which_node];
324 gezelter 907
325 gezelter 910 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
326 gezelter 916 myPotato, MPI_COMM_WORLD, &istatus);
327 gezelter 907
328 tim 920 atomTypeString = MPIatomTypeString;
329    
330 gezelter 916 myPotato++;
331    
332 gezelter 907 MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
333 gezelter 916 myPotato, MPI_COMM_WORLD, &istatus);
334    
335     myPotato++;
336 gezelter 907
337 gezelter 916 if (isDirectional) {
338     MPI_Recv(atomData13, 13, MPI_DOUBLE, which_node,
339     myPotato, MPI_COMM_WORLD, &istatus);
340     } else {
341     MPI_Recv(atomData6, 6, MPI_DOUBLE, which_node,
342     myPotato, MPI_COMM_WORLD, &istatus);
343 gezelter 907 }
344 gezelter 916
345     myPotato++;
346     potatoes[which_node] = myPotato;
347 gezelter 907
348     } else {
349    
350 tim 934 haveError = 0;
351     which_atom = i;
352 gezelter 916
353 tim 929 local_index = indexArray[currentIndex].first;
354    
355     if (which_atom == indexArray[currentIndex].second) {
356 gezelter 916
357 gezelter 907 atomTypeString = atoms[local_index]->getType();
358    
359 tim 934 atoms[local_index]->getPos(pos);
360     atoms[local_index]->getVel(vel);
361 mmeineke 670
362 gezelter 916 atomData6[0] = pos[0];
363     atomData6[1] = pos[1];
364     atomData6[2] = pos[2];
365 tim 837
366 gezelter 916 atomData6[3] = vel[0];
367     atomData6[4] = vel[1];
368     atomData6[5] = vel[2];
369 gezelter 907
370     isDirectional = 0;
371    
372 chuckv 436 if( atoms[local_index]->isDirectional() ){
373 tim 837
374 gezelter 907 isDirectional = 1;
375    
376 chuckv 436 dAtom = (DirectionalAtom *)atoms[local_index];
377     dAtom->getQ( q );
378 gezelter 916
379     for (int j = 0; j < 6 ; j++)
380     atomData13[j] = atomData6[j];
381 gezelter 907
382 gezelter 916 atomData13[6] = q[0];
383     atomData13[7] = q[1];
384     atomData13[8] = q[2];
385     atomData13[9] = q[3];
386    
387     atomData13[10] = dAtom->getJx();
388     atomData13[11] = dAtom->getJy();
389     atomData13[12] = dAtom->getJz();
390 gezelter 907 }
391 gezelter 916
392 gezelter 907 } else {
393 tim 934 sprintf(painCave.errMsg,
394     "Atom %d not found on processor %d\n",
395     i, worldRank );
396     haveError= 1;
397     simError();
398     }
399 gezelter 916
400 tim 934 if(haveError) DieDieDie();
401 gezelter 916
402 tim 929 currentIndex ++;
403 tim 926 }
404     // If we've survived to here, format the line:
405    
406     if (!isDirectional) {
407    
408 tim 934 sprintf( writeLine,
409     "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
410     atomTypeString,
411     atomData6[0],
412     atomData6[1],
413     atomData6[2],
414     atomData6[3],
415     atomData6[4],
416     atomData6[5]);
417 tim 929
418 tim 934 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
419 tim 926
420     } else {
421    
422 tim 934 sprintf( writeLine,
423     "%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",
424     atomTypeString,
425     atomData13[0],
426     atomData13[1],
427     atomData13[2],
428     atomData13[3],
429     atomData13[4],
430     atomData13[5],
431     atomData13[6],
432     atomData13[7],
433     atomData13[8],
434     atomData13[9],
435     atomData13[10],
436     atomData13[11],
437     atomData13[12]);
438 gezelter 907
439     }
440 tim 926
441 tim 934 for(k = 0; k < outFile.size(); k++)
442     *outFile[k] << writeLine;
443 mmeineke 377 }
444 tim 926
445 tim 934 for(k = 0; k < outFile.size(); k++)
446     outFile[k]->flush();
447    
448 gezelter 907 sprintf( checkPointMsg,
449     "Sucessfully took a dump.\n");
450 tim 934
451 gezelter 907 MPIcheckPoint();
452 tim 934
453 tim 919 delete[] potatoes;
454 tim 934
455 gezelter 415 } else {
456 tim 837
457 gezelter 907 // worldRank != 0, so I'm a remote node.
458 gezelter 916
459     // Set my magic potato to 0:
460    
461     myPotato = 0;
462 tim 929 currentIndex = 0;
463 gezelter 907
464     for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
465    
466     // Am I the node which has this atom?
467    
468     if (AtomToProcMap[i] == worldRank) {
469 tim 837
470 gezelter 916 if (myPotato + 3 >= MAXTAG) {
471    
472     // The potato was going to exceed the maximum value,
473     // so wrap this processor potato back to 0 (and block until
474     // node 0 says we can go:
475    
476     MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
477    
478     }
479 tim 919 which_atom = i;
480 tim 929 local_index = indexArray[currentIndex].first;
481 tim 920
482 tim 929 if (which_atom == indexArray[currentIndex].second) {
483 gezelter 907
484     atomTypeString = atoms[local_index]->getType();
485 tim 837
486 tim 934 atoms[local_index]->getPos(pos);
487     atoms[local_index]->getVel(vel);
488 tim 837
489 gezelter 916 atomData6[0] = pos[0];
490     atomData6[1] = pos[1];
491     atomData6[2] = pos[2];
492 tim 837
493 gezelter 916 atomData6[3] = vel[0];
494     atomData6[4] = vel[1];
495     atomData6[5] = vel[2];
496 gezelter 907
497     isDirectional = 0;
498 tim 837
499 gezelter 907 if( atoms[local_index]->isDirectional() ){
500 mmeineke 377
501 gezelter 907 isDirectional = 1;
502    
503     dAtom = (DirectionalAtom *)atoms[local_index];
504     dAtom->getQ( q );
505    
506 gezelter 916 for (int j = 0; j < 6 ; j++)
507     atomData13[j] = atomData6[j];
508    
509     atomData13[6] = q[0];
510     atomData13[7] = q[1];
511     atomData13[8] = q[2];
512     atomData13[9] = q[3];
513 tim 934
514 gezelter 916 atomData13[10] = dAtom->getJx();
515     atomData13[11] = dAtom->getJy();
516     atomData13[12] = dAtom->getJz();
517 gezelter 907 }
518 tim 837
519 gezelter 907 } else {
520 tim 934 sprintf(painCave.errMsg,
521     "Atom %d not found on processor %d\n",
522     i, worldRank );
523     haveError= 1;
524     simError();
525     }
526 tim 837
527 gezelter 916 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
528 tim 837
529 gezelter 916 // null terminate the string before sending (just in case):
530     MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
531 tim 837
532 gezelter 910 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
533 tim 934 myPotato, MPI_COMM_WORLD);
534 gezelter 907
535 gezelter 916 myPotato++;
536    
537 gezelter 907 MPI_Send(&isDirectional, 1, MPI_INT, 0,
538 tim 934 myPotato, MPI_COMM_WORLD);
539 gezelter 907
540 gezelter 916 myPotato++;
541    
542 gezelter 907 if (isDirectional) {
543    
544 gezelter 916 MPI_Send(atomData13, 13, MPI_DOUBLE, 0,
545     myPotato, MPI_COMM_WORLD);
546 gezelter 907
547 gezelter 916 } else {
548    
549     MPI_Send(atomData6, 6, MPI_DOUBLE, 0,
550     myPotato, MPI_COMM_WORLD);
551 gezelter 907 }
552 gezelter 916
553 tim 929 myPotato++;
554     currentIndex++;
555 mmeineke 377 }
556 gezelter 907 }
557 mmeineke 440
558 gezelter 907 sprintf( checkPointMsg,
559 gezelter 916 "Sucessfully took a dump.\n");
560 gezelter 907 MPIcheckPoint();
561    
562 gezelter 916 }
563 gezelter 907
564 mmeineke 377 #endif // is_mpi
565     }
566 mmeineke 440
567     #ifdef IS_MPI
568    
569     // a couple of functions to let us escape the write loop
570    
571 gezelter 907 void dWrite::DieDieDie( void ){
572 tim 837
573 mmeineke 440 MPI_Finalize();
574     exit (0);
575     }
576    
577     #endif //is_mpi