ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 949
Committed: Thu Jan 15 21:57:10 2004 UTC (20 years, 5 months ago) by chuckv
File size: 15449 byte(s)
Log Message:
Fixes for Dumps

File Contents

# Content
1 #define _FILE_OFFSET_BITS 64
2
3 #include <string.h>
4 #include <iostream>
5 #include <fstream>
6 #include <algorithm>
7 #include <utility>
8
9 #ifdef IS_MPI
10 #include <mpi.h>
11 #include "mpiSimulation.hpp"
12
13 namespace dWrite{
14 void DieDieDie( void );
15 }
16
17 using namespace dWrite;
18 #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
31 dumpFile.open(entry_plug->sampleName, ios::out | ios::trunc );
32
33 if( !dumpFile ){
34
35 sprintf( painCave.errMsg,
36 "Could not open \"%s\" for dump output.\n",
37 entry_plug->sampleName);
38 painCave.isFatal = 1;
39 simError();
40 }
41
42 #ifdef IS_MPI
43 }
44
45 //sort the local atoms by global index
46 sortByGlobalIndex();
47
48 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 dumpFile.close();
61
62 #ifdef IS_MPI
63 }
64 #endif // is_mpi
65 }
66
67 #ifdef IS_MPI
68
69 /**
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 for(int i = 0; i < mpiSim->getMyNlocal();i++)
95 indexArray.push_back(make_pair(i, atoms[i]->getGlobalIndex()));
96
97 sort(indexArray.begin(), indexArray.end(), indexSortingCriterion);
98
99 //for (int i = 0; i < mpiSim->getMyNlocal(); i++) {
100 // printf("node %d has global %d at local %d\n", worldRank, indexArray[i].second, indexArray[i].first);
101 //}
102
103 }
104
105 #endif
106
107 void DumpWriter::writeDump(double currentTime){
108
109 ofstream finalOut;
110 vector<ofstream*> fileStreams;
111
112 #ifdef IS_MPI
113 printf("Hello from node %d\n", worldRank);
114 sortByGlobalIndex();
115 if(worldRank == 0 ){
116
117 finalOut.open( entry_plug->finalName, ios::out | ios::trunc );
118 if( !finalOut ){
119 sprintf( painCave.errMsg,
120 "Could not open \"%s\" for final dump output.\n",
121 entry_plug->finalName );
122 painCave.isFatal = 1;
123 simError();
124 }
125 }
126 #endif // is_mpi
127
128 fileStreams.push_back(&finalOut);
129 fileStreams.push_back(&dumpFile);
130
131 writeFrame(fileStreams, currentTime);
132
133 #ifdef IS_MPI
134 finalOut.close();
135 #endif
136
137 }
138
139 void DumpWriter::writeFinal(double currentTime){
140
141 ofstream finalOut;
142 vector<ofstream*> fileStreams;
143
144 #ifdef IS_MPI
145 if(worldRank == 0 ){
146
147 finalOut.open( entry_plug->finalName, ios::out | ios::trunc );
148
149 if( !finalOut ){
150 sprintf( painCave.errMsg,
151 "Could not open \"%s\" for final dump output.\n",
152 entry_plug->finalName );
153 painCave.isFatal = 1;
154 simError();
155 }
156
157 }
158 #endif // is_mpi
159
160 fileStreams.push_back(&finalOut);
161 writeFrame(fileStreams, currentTime);
162
163 #ifdef IS_MPI
164 finalOut.close();
165 #endif
166
167 }
168
169 void DumpWriter::writeFrame( vector<ofstream*>& outFile, double currentTime ){
170
171 const int BUFFERSIZE = 2000;
172 const int MINIBUFFERSIZE = 100;
173
174 char tempBuffer[BUFFERSIZE];
175 char writeLine[BUFFERSIZE];
176
177 int i, k;
178
179 #ifdef IS_MPI
180
181 /*********************************************************************
182 * Documentation? You want DOCUMENTATION?
183 *
184 * Why all the potatoes below?
185 *
186 * To make a long story short, the original version of DumpWriter
187 * worked in the most inefficient way possible. Node 0 would
188 * poke each of the node for an individual atom's formatted data
189 * as node 0 worked its way down the global index. This was particularly
190 * inefficient since the method blocked all processors at every atom
191 * (and did it twice!).
192 *
193 * An intermediate version of DumpWriter could be described from Node
194 * zero's perspective as follows:
195 *
196 * 1) Have 100 of your friends stand in a circle.
197 * 2) When you say go, have all of them start tossing potatoes at
198 * you (one at a time).
199 * 3) Catch the potatoes.
200 *
201 * It was an improvement, but MPI has buffers and caches that could
202 * best be described in this analogy as "potato nets", so there's no
203 * need to block the processors atom-by-atom.
204 *
205 * This new and improved DumpWriter works in an even more efficient
206 * way:
207 *
208 * 1) Have 100 of your friend stand in a circle.
209 * 2) When you say go, have them start tossing 5-pound bags of
210 * potatoes at you.
211 * 3) Once you've caught a friend's bag of potatoes,
212 * toss them a spud to let them know they can toss another bag.
213 *
214 * How's THAT for documentation?
215 *
216 *********************************************************************/
217
218 int *potatoes;
219 int myPotato;
220
221 int nProc;
222 int j, which_node, done, which_atom, local_index, currentIndex;
223 double atomData6[6];
224 double atomData13[13];
225 int isDirectional;
226 char* atomTypeString;
227 char MPIatomTypeString[MINIBUFFERSIZE];
228
229 #else //is_mpi
230 int nAtoms = entry_plug->n_atoms;
231 #endif //is_mpi
232
233 double q[4];
234 DirectionalAtom* dAtom;
235 Atom** atoms = entry_plug->atoms;
236 double pos[3], vel[3];
237
238 #ifndef IS_MPI
239
240 for(k = 0; k < outFile.size(); k++){
241 *outFile[k] << nAtoms << "\n";
242
243 *outFile[k] << currentTime << ";\t"
244 << entry_plug->Hmat[0][0] << "\t"
245 << entry_plug->Hmat[1][0] << "\t"
246 << entry_plug->Hmat[2][0] << ";\t"
247
248 << entry_plug->Hmat[0][1] << "\t"
249 << entry_plug->Hmat[1][1] << "\t"
250 << entry_plug->Hmat[2][1] << ";\t"
251
252 << entry_plug->Hmat[0][2] << "\t"
253 << entry_plug->Hmat[1][2] << "\t"
254 << entry_plug->Hmat[2][2] << ";";
255
256 //write out additional parameters, such as chi and eta
257 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
258 }
259
260 for( i=0; i<nAtoms; i++ ){
261
262 atoms[i]->getPos(pos);
263 atoms[i]->getVel(vel);
264
265 sprintf( tempBuffer,
266 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
267 atoms[i]->getType(),
268 pos[0],
269 pos[1],
270 pos[2],
271 vel[0],
272 vel[1],
273 vel[2]);
274 strcpy( writeLine, tempBuffer );
275
276 if( atoms[i]->isDirectional() ){
277
278 dAtom = (DirectionalAtom *)atoms[i];
279 dAtom->getQ( q );
280
281 sprintf( tempBuffer,
282 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
283 q[0],
284 q[1],
285 q[2],
286 q[3],
287 dAtom->getJx(),
288 dAtom->getJy(),
289 dAtom->getJz());
290 strcat( writeLine, tempBuffer );
291 }
292 else
293 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
294
295 for(k = 0; k < outFile.size(); k++)
296 *outFile[k] << writeLine;
297 }
298
299 #else // is_mpi
300
301 /* code to find maximum tag value */
302
303 int *tagub, flag, MAXTAG;
304 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
305 if (flag) {
306 MAXTAG = *tagub;
307 } else {
308 MAXTAG = 32767;
309 }
310
311 int haveError;
312
313 MPI_Status istatus;
314 int *AtomToProcMap = mpiSim->getAtomToProcMap();
315
316 // write out header and node 0's coordinates
317
318 if( worldRank == 0 ){
319
320 // Node 0 needs a list of the magic potatoes for each processor;
321
322 nProc = mpiSim->getNumberProcessors();
323 potatoes = new int[nProc];
324
325 //write out the comment lines
326 for (i = 0; i < nProc; i++)
327 potatoes[i] = 0;
328
329 for(k = 0; k < outFile.size(); k++){
330 *outFile[k] << mpiSim->getTotAtoms() << "\n";
331
332 *outFile[k] << currentTime << ";\t"
333 << entry_plug->Hmat[0][0] << "\t"
334 << entry_plug->Hmat[1][0] << "\t"
335 << entry_plug->Hmat[2][0] << ";\t"
336
337 << entry_plug->Hmat[0][1] << "\t"
338 << entry_plug->Hmat[1][1] << "\t"
339 << entry_plug->Hmat[2][1] << ";\t"
340
341 << entry_plug->Hmat[0][2] << "\t"
342 << entry_plug->Hmat[1][2] << "\t"
343 << entry_plug->Hmat[2][2] << ";";
344
345 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
346 }
347
348 currentIndex = 0;
349
350 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
351
352 // Get the Node number which has this atom;
353
354 which_node = AtomToProcMap[i];
355
356 if (which_node != 0) {
357
358 if (potatoes[which_node] + 3 >= MAXTAG) {
359 // The potato was going to exceed the maximum value,
360 // so wrap this processor potato back to 0:
361
362 potatoes[which_node] = 0;
363 MPI_Send(0, 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
364
365 }
366
367 myPotato = potatoes[which_node];
368
369 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
370 myPotato, MPI_COMM_WORLD, &istatus);
371
372 atomTypeString = MPIatomTypeString;
373
374 myPotato++;
375
376 MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
377 myPotato, MPI_COMM_WORLD, &istatus);
378
379 myPotato++;
380
381 if (isDirectional) {
382 MPI_Recv(atomData13, 13, MPI_DOUBLE, which_node,
383 myPotato, MPI_COMM_WORLD, &istatus);
384 } else {
385 MPI_Recv(atomData6, 6, MPI_DOUBLE, which_node,
386 myPotato, MPI_COMM_WORLD, &istatus);
387 }
388
389 myPotato++;
390 potatoes[which_node] = myPotato;
391
392 } else {
393
394 haveError = 0;
395 which_atom = i;
396
397 //local_index = -1;
398
399 //for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
400 // if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
401 //}
402
403 //if (local_index != -1) {
404
405 local_index = indexArray[currentIndex].first;
406
407 if (which_atom == indexArray[currentIndex].second) {
408
409 atomTypeString = atoms[local_index]->getType();
410
411 atoms[local_index]->getPos(pos);
412 atoms[local_index]->getVel(vel);
413
414 atomData6[0] = pos[0];
415 atomData6[1] = pos[1];
416 atomData6[2] = pos[2];
417
418 atomData6[3] = vel[0];
419 atomData6[4] = vel[1];
420 atomData6[5] = vel[2];
421
422 isDirectional = 0;
423
424 if( atoms[local_index]->isDirectional() ){
425
426 isDirectional = 1;
427
428 dAtom = (DirectionalAtom *)atoms[local_index];
429 dAtom->getQ( q );
430
431 for (int j = 0; j < 6 ; j++)
432 atomData13[j] = atomData6[j];
433
434 atomData13[6] = q[0];
435 atomData13[7] = q[1];
436 atomData13[8] = q[2];
437 atomData13[9] = q[3];
438
439 atomData13[10] = dAtom->getJx();
440 atomData13[11] = dAtom->getJy();
441 atomData13[12] = dAtom->getJz();
442 }
443
444 } else {
445 sprintf(painCave.errMsg,
446 "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
447 which_atom, worldRank, currentIndex, local_index );
448 haveError= 1;
449 simError();
450 }
451
452 if(haveError) DieDieDie();
453
454 currentIndex++;
455 }
456 // If we've survived to here, format the line:
457
458 if (!isDirectional) {
459
460 sprintf( writeLine,
461 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
462 atomTypeString,
463 atomData6[0],
464 atomData6[1],
465 atomData6[2],
466 atomData6[3],
467 atomData6[4],
468 atomData6[5]);
469
470 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
471
472 } else {
473
474 sprintf( writeLine,
475 "%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",
476 atomTypeString,
477 atomData13[0],
478 atomData13[1],
479 atomData13[2],
480 atomData13[3],
481 atomData13[4],
482 atomData13[5],
483 atomData13[6],
484 atomData13[7],
485 atomData13[8],
486 atomData13[9],
487 atomData13[10],
488 atomData13[11],
489 atomData13[12]);
490
491 }
492
493 for(k = 0; k < outFile.size(); k++)
494 *outFile[k] << writeLine;
495 }
496
497 for(k = 0; k < outFile.size(); k++)
498 outFile[k]->flush();
499
500 sprintf( checkPointMsg,
501 "Sucessfully took a dump.\n");
502
503 MPIcheckPoint();
504
505 delete[] potatoes;
506
507 } else {
508
509 // worldRank != 0, so I'm a remote node.
510
511 // Set my magic potato to 0:
512
513 myPotato = 0;
514 currentIndex = 0;
515
516 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
517
518 // Am I the node which has this atom?
519
520 if (AtomToProcMap[i] == worldRank) {
521
522 if (myPotato + 3 >= MAXTAG) {
523
524 // The potato was going to exceed the maximum value,
525 // so wrap this processor potato back to 0 (and block until
526 // node 0 says we can go:
527
528 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
529
530 }
531 which_atom = i;
532
533 //local_index = -1;
534
535 //for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
536 // if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
537 //}
538
539 //if (local_index != -1) {
540
541 local_index = indexArray[currentIndex].first;
542
543 if (which_atom == indexArray[currentIndex].second) {
544
545 atomTypeString = atoms[local_index]->getType();
546
547 atoms[local_index]->getPos(pos);
548 atoms[local_index]->getVel(vel);
549
550 atomData6[0] = pos[0];
551 atomData6[1] = pos[1];
552 atomData6[2] = pos[2];
553
554 atomData6[3] = vel[0];
555 atomData6[4] = vel[1];
556 atomData6[5] = vel[2];
557
558 isDirectional = 0;
559
560 if( atoms[local_index]->isDirectional() ){
561
562 isDirectional = 1;
563
564 dAtom = (DirectionalAtom *)atoms[local_index];
565 dAtom->getQ( q );
566
567 for (int j = 0; j < 6 ; j++)
568 atomData13[j] = atomData6[j];
569
570 atomData13[6] = q[0];
571 atomData13[7] = q[1];
572 atomData13[8] = q[2];
573 atomData13[9] = q[3];
574
575 atomData13[10] = dAtom->getJx();
576 atomData13[11] = dAtom->getJy();
577 atomData13[12] = dAtom->getJz();
578 }
579
580 } else {
581 sprintf(painCave.errMsg,
582 "Atom %d not found on processor %d, currentIndex = %d, local_index = %d\n",
583 which_atom, worldRank, currentIndex, local_index );
584 haveError= 1;
585 simError();
586 }
587
588 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
589
590 // null terminate the string before sending (just in case):
591 MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
592
593 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
594 myPotato, MPI_COMM_WORLD);
595
596 myPotato++;
597
598 MPI_Send(&isDirectional, 1, MPI_INT, 0,
599 myPotato, MPI_COMM_WORLD);
600
601 myPotato++;
602
603 if (isDirectional) {
604
605 MPI_Send(atomData13, 13, MPI_DOUBLE, 0,
606 myPotato, MPI_COMM_WORLD);
607
608 } else {
609
610 MPI_Send(atomData6, 6, MPI_DOUBLE, 0,
611 myPotato, MPI_COMM_WORLD);
612 }
613
614 myPotato++;
615 currentIndex++;
616 }
617 }
618
619 sprintf( checkPointMsg,
620 "Sucessfully took a dump.\n");
621 MPIcheckPoint();
622
623 }
624
625 #endif // is_mpi
626 }
627
628 #ifdef IS_MPI
629
630 // a couple of functions to let us escape the write loop
631
632 void dWrite::DieDieDie( void ){
633
634 MPI_Finalize();
635 exit (0);
636 }
637
638 #endif //is_mpi