ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 1231
Committed: Thu Jun 3 21:06:51 2004 UTC (20 years, 1 month ago) by tim
File size: 17403 byte(s)
Log Message:
fixed a bug which only writes out the first atom of a molecule

File Contents

# Content
1 #define _LARGEFILE_SOURCE64
2 #define _FILE_OFFSET_BITS 64
3
4 #include <string.h>
5 #include <iostream>
6 #include <fstream>
7 #include <algorithm>
8 #include <utility>
9
10 #ifdef IS_MPI
11 #include <mpi.h>
12 #include "mpiSimulation.hpp"
13
14 namespace dWrite{
15 void DieDieDie( void );
16 }
17
18 using namespace dWrite;
19 #endif //is_mpi
20
21 #include "ReadWrite.hpp"
22 #include "simError.h"
23
24 DumpWriter::DumpWriter( SimInfo* the_entry_plug ){
25
26 entry_plug = the_entry_plug;
27
28 #ifdef IS_MPI
29 if(worldRank == 0 ){
30 #endif // is_mpi
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 Molecule* mols = entry_plug->molecules;
92 indexArray.clear();
93
94 for(int i = 0; i < entry_plug->n_mol;i++)
95 indexArray.push_back(make_pair(i, mols[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 atomData[13];
220 int isDirectional;
221 char* atomTypeString;
222 char MPIatomTypeString[MINIBUFFERSIZE];
223 int nObjects;
224 int msgLen; // the length of message actually recieved at master nodes
225 #endif //is_mpi
226
227 double q[4], ji[3];
228 DirectionalAtom* dAtom;
229 double pos[3], vel[3];
230 int nTotObjects;
231 StuntDouble* sd;
232 char* molName;
233 vector<StuntDouble*> integrableObjects;
234 vector<StuntDouble*>::iterator iter;
235 nTotObjects = entry_plug->getTotIntegrableObjects();
236 #ifndef IS_MPI
237
238 for(k = 0; k < outFile.size(); k++){
239 *outFile[k] << nTotObjects << "\n";
240
241 *outFile[k] << currentTime << ";\t"
242 << entry_plug->Hmat[0][0] << "\t"
243 << entry_plug->Hmat[1][0] << "\t"
244 << entry_plug->Hmat[2][0] << ";\t"
245
246 << entry_plug->Hmat[0][1] << "\t"
247 << entry_plug->Hmat[1][1] << "\t"
248 << entry_plug->Hmat[2][1] << ";\t"
249
250 << entry_plug->Hmat[0][2] << "\t"
251 << entry_plug->Hmat[1][2] << "\t"
252 << entry_plug->Hmat[2][2] << ";";
253
254 //write out additional parameters, such as chi and eta
255 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
256 }
257
258 for( i=0; i< entry_plug->n_mol; i++ ){
259
260 integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
261 molName = (entry_plug->compStamps[entry_plug->molecules[i].getStampID()])->getID();
262
263 for( iter = integrableObjects.begin();iter != integrableObjects.end(); ++iter){
264 sd = *iter;
265 sd->getPos(pos);
266 sd->getVel(vel);
267
268 sprintf( tempBuffer,
269 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
270 sd->getType(),
271 pos[0],
272 pos[1],
273 pos[2],
274 vel[0],
275 vel[1],
276 vel[2]);
277 strcpy( writeLine, tempBuffer );
278
279 if( sd->isDirectional() ){
280
281 sd->getQ( q );
282 sd->getJ( ji );
283
284 sprintf( tempBuffer,
285 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
286 q[0],
287 q[1],
288 q[2],
289 q[3],
290 ji[0],
291 ji[1],
292 ji[2]);
293 strcat( writeLine, tempBuffer );
294 }
295 else
296 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
297
298 for(k = 0; k < outFile.size(); k++)
299 *outFile[k] << writeLine;
300 }
301
302 }
303
304 #else // is_mpi
305
306 /* code to find maximum tag value */
307
308 int *tagub, flag, MAXTAG;
309 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
310 if (flag) {
311 MAXTAG = *tagub;
312 } else {
313 MAXTAG = 32767;
314 }
315
316 int haveError;
317
318 MPI_Status istatus;
319 int nCurObj;
320 int *MolToProcMap = mpiSim->getMolToProcMap();
321
322 // write out header and node 0's coordinates
323
324 if( worldRank == 0 ){
325
326 // Node 0 needs a list of the magic potatoes for each processor;
327
328 nProc = mpiSim->getNProcessors();
329 potatoes = new int[nProc];
330
331 //write out the comment lines
332 for (i = 0; i < nProc; i++)
333 potatoes[i] = 0;
334
335 for(k = 0; k < outFile.size(); k++){
336 *outFile[k] << nTotObjects << "\n";
337
338 *outFile[k] << currentTime << ";\t"
339 << entry_plug->Hmat[0][0] << "\t"
340 << entry_plug->Hmat[1][0] << "\t"
341 << entry_plug->Hmat[2][0] << ";\t"
342
343 << entry_plug->Hmat[0][1] << "\t"
344 << entry_plug->Hmat[1][1] << "\t"
345 << entry_plug->Hmat[2][1] << ";\t"
346
347 << entry_plug->Hmat[0][2] << "\t"
348 << entry_plug->Hmat[1][2] << "\t"
349 << entry_plug->Hmat[2][2] << ";";
350
351 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
352 }
353
354 currentIndex = 0;
355
356 for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
357
358 // Get the Node number which has this atom;
359
360 which_node = MolToProcMap[i];
361
362 if (which_node != 0) {
363
364 if (potatoes[which_node] + 1 >= MAXTAG) {
365 // The potato was going to exceed the maximum value,
366 // so wrap this processor potato back to 0:
367
368 potatoes[which_node] = 0;
369 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
370
371 }
372
373 myPotato = potatoes[which_node];
374
375 //recieve the number of integrableObject in current molecule
376 MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
377 myPotato, MPI_COMM_WORLD, &istatus);
378 myPotato++;
379
380 for(int l = 0; l < nCurObj; l++){
381
382 if (potatoes[which_node] + 2 >= MAXTAG) {
383 // The potato was going to exceed the maximum value,
384 // so wrap this processor potato back to 0:
385
386 potatoes[which_node] = 0;
387 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
388
389 }
390
391 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
392 myPotato, MPI_COMM_WORLD, &istatus);
393
394 atomTypeString = MPIatomTypeString;
395
396 myPotato++;
397
398 MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
399 myPotato++;
400
401 MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
402
403 if(msgLen == 13)
404 isDirectional = 1;
405 else
406 isDirectional = 0;
407
408 // If we've survived to here, format the line:
409
410 if (!isDirectional) {
411
412 sprintf( writeLine,
413 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
414 atomTypeString,
415 atomData[0],
416 atomData[1],
417 atomData[2],
418 atomData[3],
419 atomData[4],
420 atomData[5]);
421
422 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
423
424 }
425 else {
426
427 sprintf( writeLine,
428 "%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",
429 atomTypeString,
430 atomData[0],
431 atomData[1],
432 atomData[2],
433 atomData[3],
434 atomData[4],
435 atomData[5],
436 atomData[6],
437 atomData[7],
438 atomData[8],
439 atomData[9],
440 atomData[10],
441 atomData[11],
442 atomData[12]);
443
444 }
445
446 for(k = 0; k < outFile.size(); k++)
447 *outFile[k] << writeLine;
448
449 }// end for(int l =0)
450 potatoes[which_node] = myPotato;
451
452 }
453 else {
454
455 haveError = 0;
456
457 local_index = indexArray[currentIndex].first;
458
459 integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects();
460
461 for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){
462 sd = *iter;
463 atomTypeString = sd->getType();
464
465 sd->getPos(pos);
466 sd->getVel(vel);
467
468 atomData[0] = pos[0];
469 atomData[1] = pos[1];
470 atomData[2] = pos[2];
471
472 atomData[3] = vel[0];
473 atomData[4] = vel[1];
474 atomData[5] = vel[2];
475
476 isDirectional = 0;
477
478 if( sd->isDirectional() ){
479
480 isDirectional = 1;
481
482 sd->getQ( q );
483 sd->getJ( ji );
484
485 for (int j = 0; j < 6 ; j++)
486 atomData[j] = atomData[j];
487
488 atomData[6] = q[0];
489 atomData[7] = q[1];
490 atomData[8] = q[2];
491 atomData[9] = q[3];
492
493 atomData[10] = ji[0];
494 atomData[11] = ji[1];
495 atomData[12] = ji[2];
496 }
497
498 // If we've survived to here, format the line:
499
500 if (!isDirectional) {
501
502 sprintf( writeLine,
503 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
504 atomTypeString,
505 atomData[0],
506 atomData[1],
507 atomData[2],
508 atomData[3],
509 atomData[4],
510 atomData[5]);
511
512 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
513
514 }
515 else {
516
517 sprintf( writeLine,
518 "%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",
519 atomTypeString,
520 atomData[0],
521 atomData[1],
522 atomData[2],
523 atomData[3],
524 atomData[4],
525 atomData[5],
526 atomData[6],
527 atomData[7],
528 atomData[8],
529 atomData[9],
530 atomData[10],
531 atomData[11],
532 atomData[12]);
533
534 }
535
536 for(k = 0; k < outFile.size(); k++)
537 *outFile[k] << writeLine;
538
539
540 }//end for(iter = integrableObject.begin())
541
542 currentIndex++;
543 }
544
545 }//end for(i = 0; i < mpiSim->getNmol())
546
547 for(k = 0; k < outFile.size(); k++)
548 outFile[k]->flush();
549
550 sprintf( checkPointMsg,
551 "Sucessfully took a dump.\n");
552
553 MPIcheckPoint();
554
555 delete[] potatoes;
556
557 } else {
558
559 // worldRank != 0, so I'm a remote node.
560
561 // Set my magic potato to 0:
562
563 myPotato = 0;
564 currentIndex = 0;
565
566 for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
567
568 // Am I the node which has this integrableObject?
569
570 if (MolToProcMap[i] == worldRank) {
571
572
573 if (myPotato + 1 >= MAXTAG) {
574
575 // The potato was going to exceed the maximum value,
576 // so wrap this processor potato back to 0 (and block until
577 // node 0 says we can go:
578
579 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
580
581 }
582
583 local_index = indexArray[currentIndex].first;
584 integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
585
586 nCurObj = integrableObjects.size();
587
588 MPI_Send(&nCurObj, 1, MPI_INT, 0,
589 myPotato, MPI_COMM_WORLD);
590 myPotato++;
591
592 for( iter = integrableObjects.begin(); iter != integrableObjects.end(); iter++){
593
594 if (myPotato + 2 >= MAXTAG) {
595
596 // The potato was going to exceed the maximum value,
597 // so wrap this processor potato back to 0 (and block until
598 // node 0 says we can go:
599
600 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
601
602 }
603
604 sd = *iter;
605
606 atomTypeString = sd->getType();
607
608 sd->getPos(pos);
609 sd->getVel(vel);
610
611 atomData[0] = pos[0];
612 atomData[1] = pos[1];
613 atomData[2] = pos[2];
614
615 atomData[3] = vel[0];
616 atomData[4] = vel[1];
617 atomData[5] = vel[2];
618
619 isDirectional = 0;
620
621 if( sd->isDirectional() ){
622
623 isDirectional = 1;
624
625 sd->getQ( q );
626 sd->getJ( ji );
627
628
629 atomData[6] = q[0];
630 atomData[7] = q[1];
631 atomData[8] = q[2];
632 atomData[9] = q[3];
633
634 atomData[10] = ji[0];
635 atomData[11] = ji[1];
636 atomData[12] = ji[2];
637 }
638
639
640 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
641
642 // null terminate the string before sending (just in case):
643 MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
644
645 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
646 myPotato, MPI_COMM_WORLD);
647
648 myPotato++;
649
650 if (isDirectional) {
651
652 MPI_Send(atomData, 13, MPI_DOUBLE, 0,
653 myPotato, MPI_COMM_WORLD);
654
655 } else {
656
657 MPI_Send(atomData, 6, MPI_DOUBLE, 0,
658 myPotato, MPI_COMM_WORLD);
659 }
660
661 myPotato++;
662
663 }
664
665 currentIndex++;
666
667 }
668
669 }
670
671 sprintf( checkPointMsg,
672 "Sucessfully took a dump.\n");
673 MPIcheckPoint();
674
675 }
676
677
678
679 #endif // is_mpi
680 }
681
682 #ifdef IS_MPI
683
684 // a couple of functions to let us escape the write loop
685
686 void dWrite::DieDieDie( void ){
687
688 MPI_Finalize();
689 exit (0);
690 }
691
692 #endif //is_mpi