ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 934
Committed: Tue Jan 13 20:04:28 2004 UTC (20 years, 5 months ago) by tim
File size: 13174 byte(s)
Log Message:
change the interface of writeFrame

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