ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 947
Committed: Thu Jan 15 14:22:16 2004 UTC (20 years, 5 months ago) by gezelter
File size: 15023 byte(s)
Log Message:
Documented the Spud Toss

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