ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 1000
Committed: Fri Jan 30 21:47:22 2004 UTC (20 years, 5 months ago) by tim
File size: 14799 byte(s)
Log Message:
using class  Minimizer1D derived from MinimizerBase instead of a functor to do line seach

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