ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 929
Committed: Tue Jan 13 15:46:49 2004 UTC (20 years, 5 months ago) by tim
File size: 12186 byte(s)
Log Message:
 Merge the code of writeFinal and writeDump;
 Adding sortingIndex into DumpWriter;
 Fix a bug of writing last frame twice in integrator

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