ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/io/DumpWriter.cpp
Revision: 1695
Committed: Mon Nov 1 22:52:57 2004 UTC (19 years, 8 months ago) by tim
File size: 18217 byte(s)
Log Message:
Molecule, Atom, DirectionalAtom, RigidBody and StuntDouble classes get compiled

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 "brains/mpiSimulation.hpp"
13
14 namespace dWrite{
15 void DieDieDie( void );
16 }
17
18 using namespace dWrite;
19 #endif //is_mpi
20
21 #include "io/ReadWrite.hpp"
22 #include "utils/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.c_str(), 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.c_str());
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.c_str(), 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.c_str() );
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.c_str(), 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.c_str() );
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;
174 unsigned int k;
175
176 #ifdef IS_MPI
177
178 /*********************************************************************
179 * Documentation? You want DOCUMENTATION?
180 *
181 * Why all the potatoes below?
182 *
183 * To make a long story short, the original version of DumpWriter
184 * worked in the most inefficient way possible. Node 0 would
185 * poke each of the node for an individual atom's formatted data
186 * as node 0 worked its way down the global index. This was particularly
187 * inefficient since the method blocked all processors at every atom
188 * (and did it twice!).
189 *
190 * An intermediate version of DumpWriter could be described from Node
191 * zero's perspective as follows:
192 *
193 * 1) Have 100 of your friends stand in a circle.
194 * 2) When you say go, have all of them start tossing potatoes at
195 * you (one at a time).
196 * 3) Catch the potatoes.
197 *
198 * It was an improvement, but MPI has buffers and caches that could
199 * best be described in this analogy as "potato nets", so there's no
200 * need to block the processors atom-by-atom.
201 *
202 * This new and improved DumpWriter works in an even more efficient
203 * way:
204 *
205 * 1) Have 100 of your friend stand in a circle.
206 * 2) When you say go, have them start tossing 5-pound bags of
207 * potatoes at you.
208 * 3) Once you've caught a friend's bag of potatoes,
209 * toss them a spud to let them know they can toss another bag.
210 *
211 * How's THAT for documentation?
212 *
213 *********************************************************************/
214
215 int *potatoes;
216 int myPotato;
217
218 int nProc;
219 int j, which_node, done, which_atom, local_index, currentIndex;
220 double atomData[13];
221 int isDirectional;
222 char* atomTypeString;
223 char MPIatomTypeString[MINIBUFFERSIZE];
224 int nObjects;
225 int msgLen; // the length of message actually recieved at master nodes
226 #endif //is_mpi
227
228 Quat4d q;
229 Vector3d ji;
230 DirectionalAtom* dAtom;
231 Vector3d pos;
232 Vector3d vel;
233
234 int nTotObjects;
235 StuntDouble* sd;
236 char* molName;
237 vector<StuntDouble*> integrableObjects;
238 vector<StuntDouble*>::iterator iter;
239 nTotObjects = entry_plug->getTotIntegrableObjects();
240 #ifndef IS_MPI
241
242 for(k = 0; k < outFile.size(); k++){
243 *outFile[k] << nTotObjects << "\n";
244
245 *outFile[k] << currentTime << ";\t"
246 << entry_plug->Hmat[0][0] << "\t"
247 << entry_plug->Hmat[1][0] << "\t"
248 << entry_plug->Hmat[2][0] << ";\t"
249
250 << entry_plug->Hmat[0][1] << "\t"
251 << entry_plug->Hmat[1][1] << "\t"
252 << entry_plug->Hmat[2][1] << ";\t"
253
254 << entry_plug->Hmat[0][2] << "\t"
255 << entry_plug->Hmat[1][2] << "\t"
256 << entry_plug->Hmat[2][2] << ";";
257
258 //write out additional parameters, such as chi and eta
259 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
260 }
261
262 for( i=0; i< entry_plug->n_mol; i++ ){
263
264 integrableObjects = entry_plug->molecules[i].getIntegrableObjects();
265 molName = (entry_plug->compStamps[entry_plug->molecules[i].getStampID()])->getID();
266
267 for( iter = integrableObjects.begin();iter != integrableObjects.end(); ++iter){
268 sd = *iter;
269 pos = sd->getPos();
270 vel = sd->getVel();
271
272 sprintf( tempBuffer,
273 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
274 sd->getType(),
275 pos[0],
276 pos[1],
277 pos[2],
278 vel[0],
279 vel[1],
280 vel[2]);
281 strcpy( writeLine, tempBuffer );
282
283 if( sd->isDirectional() ){
284
285 q = sd->getQ();
286 ji = sd->getJ();
287
288 sprintf( tempBuffer,
289 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
290 q[0],
291 q[1],
292 q[2],
293 q[3],
294 ji[0],
295 ji[1],
296 ji[2]);
297 strcat( writeLine, tempBuffer );
298 }
299 else
300 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
301
302 for(k = 0; k < outFile.size(); k++)
303 *outFile[k] << writeLine;
304 }
305
306 }
307
308 #else // is_mpi
309
310 /* code to find maximum tag value */
311
312 int *tagub, flag, MAXTAG;
313 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
314 if (flag) {
315 MAXTAG = *tagub;
316 } else {
317 MAXTAG = 32767;
318 }
319
320 int haveError;
321
322 MPI_Status istatus;
323 int nCurObj;
324 int *MolToProcMap = mpiSim->getMolToProcMap();
325
326 // write out header and node 0's coordinates
327
328 if( worldRank == 0 ){
329
330 // Node 0 needs a list of the magic potatoes for each processor;
331
332 nProc = mpiSim->getNProcessors();
333 potatoes = new int[nProc];
334
335 //write out the comment lines
336 for (i = 0; i < nProc; i++)
337 potatoes[i] = 0;
338
339 for(k = 0; k < outFile.size(); k++){
340 *outFile[k] << nTotObjects << "\n";
341
342 *outFile[k] << currentTime << ";\t"
343 << entry_plug->Hmat[0][0] << "\t"
344 << entry_plug->Hmat[1][0] << "\t"
345 << entry_plug->Hmat[2][0] << ";\t"
346
347 << entry_plug->Hmat[0][1] << "\t"
348 << entry_plug->Hmat[1][1] << "\t"
349 << entry_plug->Hmat[2][1] << ";\t"
350
351 << entry_plug->Hmat[0][2] << "\t"
352 << entry_plug->Hmat[1][2] << "\t"
353 << entry_plug->Hmat[2][2] << ";";
354
355 *outFile[k] << entry_plug->the_integrator->getAdditionalParameters() << endl;
356 }
357
358 currentIndex = 0;
359
360 for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
361
362 // Get the Node number which has this atom;
363
364 which_node = MolToProcMap[i];
365
366 if (which_node != 0) {
367
368 if (potatoes[which_node] + 1 >= MAXTAG) {
369 // The potato was going to exceed the maximum value,
370 // so wrap this processor potato back to 0:
371
372 potatoes[which_node] = 0;
373 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
374
375 }
376
377 myPotato = potatoes[which_node];
378
379 //recieve the number of integrableObject in current molecule
380 MPI_Recv(&nCurObj, 1, MPI_INT, which_node,
381 myPotato, MPI_COMM_WORLD, &istatus);
382 myPotato++;
383
384 for(int l = 0; l < nCurObj; l++){
385
386 if (potatoes[which_node] + 2 >= MAXTAG) {
387 // The potato was going to exceed the maximum value,
388 // so wrap this processor potato back to 0:
389
390 potatoes[which_node] = 0;
391 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0, MPI_COMM_WORLD);
392
393 }
394
395 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
396 myPotato, MPI_COMM_WORLD, &istatus);
397
398 atomTypeString = MPIatomTypeString;
399
400 myPotato++;
401
402 MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato, MPI_COMM_WORLD, &istatus);
403 myPotato++;
404
405 MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
406
407 if(msgLen == 13)
408 isDirectional = 1;
409 else
410 isDirectional = 0;
411
412 // If we've survived to here, format the line:
413
414 if (!isDirectional) {
415
416 sprintf( writeLine,
417 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
418 atomTypeString,
419 atomData[0],
420 atomData[1],
421 atomData[2],
422 atomData[3],
423 atomData[4],
424 atomData[5]);
425
426 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
427
428 }
429 else {
430
431 sprintf( writeLine,
432 "%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",
433 atomTypeString,
434 atomData[0],
435 atomData[1],
436 atomData[2],
437 atomData[3],
438 atomData[4],
439 atomData[5],
440 atomData[6],
441 atomData[7],
442 atomData[8],
443 atomData[9],
444 atomData[10],
445 atomData[11],
446 atomData[12]);
447
448 }
449
450 for(k = 0; k < outFile.size(); k++)
451 *outFile[k] << writeLine;
452
453 }// end for(int l =0)
454 potatoes[which_node] = myPotato;
455
456 }
457 else {
458
459 haveError = 0;
460
461 local_index = indexArray[currentIndex].first;
462
463 integrableObjects = (entry_plug->molecules[local_index]).getIntegrableObjects();
464
465 for(iter= integrableObjects.begin(); iter != integrableObjects.end(); ++iter){
466 sd = *iter;
467 atomTypeString = sd->getType();
468
469 pos = sd->getPos();
470 vel = sd->getVel();
471
472 atomData[0] = pos[0];
473 atomData[1] = pos[1];
474 atomData[2] = pos[2];
475
476 atomData[3] = vel[0];
477 atomData[4] = vel[1];
478 atomData[5] = vel[2];
479
480 isDirectional = 0;
481
482 if( sd->isDirectional() ){
483
484 isDirectional = 1;
485
486 q = sd->getQ();
487 ji = sd->getJ();
488
489 for (int j = 0; j < 6 ; j++)
490 atomData[j] = atomData[j];
491
492 atomData[6] = q[0];
493 atomData[7] = q[1];
494 atomData[8] = q[2];
495 atomData[9] = q[3];
496
497 atomData[10] = ji[0];
498 atomData[11] = ji[1];
499 atomData[12] = ji[2];
500 }
501
502 // If we've survived to here, format the line:
503
504 if (!isDirectional) {
505
506 sprintf( writeLine,
507 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
508 atomTypeString,
509 atomData[0],
510 atomData[1],
511 atomData[2],
512 atomData[3],
513 atomData[4],
514 atomData[5]);
515
516 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
517
518 }
519 else {
520
521 sprintf( writeLine,
522 "%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",
523 atomTypeString,
524 atomData[0],
525 atomData[1],
526 atomData[2],
527 atomData[3],
528 atomData[4],
529 atomData[5],
530 atomData[6],
531 atomData[7],
532 atomData[8],
533 atomData[9],
534 atomData[10],
535 atomData[11],
536 atomData[12]);
537
538 }
539
540 for(k = 0; k < outFile.size(); k++)
541 *outFile[k] << writeLine;
542
543
544 }//end for(iter = integrableObject.begin())
545
546 currentIndex++;
547 }
548
549 }//end for(i = 0; i < mpiSim->getNmol())
550
551 for(k = 0; k < outFile.size(); k++)
552 outFile[k]->flush();
553
554 sprintf( checkPointMsg,
555 "Sucessfully took a dump.\n");
556
557 MPIcheckPoint();
558
559 delete[] potatoes;
560
561 } else {
562
563 // worldRank != 0, so I'm a remote node.
564
565 // Set my magic potato to 0:
566
567 myPotato = 0;
568 currentIndex = 0;
569
570 for (i = 0 ; i < mpiSim->getNMolGlobal(); i++ ) {
571
572 // Am I the node which has this integrableObject?
573
574 if (MolToProcMap[i] == worldRank) {
575
576
577 if (myPotato + 1 >= MAXTAG) {
578
579 // The potato was going to exceed the maximum value,
580 // so wrap this processor potato back to 0 (and block until
581 // node 0 says we can go:
582
583 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
584
585 }
586
587 local_index = indexArray[currentIndex].first;
588 integrableObjects = entry_plug->molecules[local_index].getIntegrableObjects();
589
590 nCurObj = integrableObjects.size();
591
592 MPI_Send(&nCurObj, 1, MPI_INT, 0,
593 myPotato, MPI_COMM_WORLD);
594 myPotato++;
595
596 for( iter = integrableObjects.begin(); iter != integrableObjects.end(); iter++){
597
598 if (myPotato + 2 >= MAXTAG) {
599
600 // The potato was going to exceed the maximum value,
601 // so wrap this processor potato back to 0 (and block until
602 // node 0 says we can go:
603
604 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &istatus);
605
606 }
607
608 sd = *iter;
609
610 atomTypeString = sd->getType();
611
612 pos = sd->getPos();
613 vel = sd->getVel();
614
615 atomData[0] = pos[0];
616 atomData[1] = pos[1];
617 atomData[2] = pos[2];
618
619 atomData[3] = vel[0];
620 atomData[4] = vel[1];
621 atomData[5] = vel[2];
622
623 isDirectional = 0;
624
625 if( sd->isDirectional() ){
626
627 isDirectional = 1;
628
629 q = sd->getQ();
630 ji = sd->getJ();
631
632
633 atomData[6] = q[0];
634 atomData[7] = q[1];
635 atomData[8] = q[2];
636 atomData[9] = q[3];
637
638 atomData[10] = ji[0];
639 atomData[11] = ji[1];
640 atomData[12] = ji[2];
641 }
642
643
644 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
645
646 // null terminate the string before sending (just in case):
647 MPIatomTypeString[MINIBUFFERSIZE-1] = '\0';
648
649 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
650 myPotato, MPI_COMM_WORLD);
651
652 myPotato++;
653
654 if (isDirectional) {
655
656 MPI_Send(atomData, 13, MPI_DOUBLE, 0,
657 myPotato, MPI_COMM_WORLD);
658
659 } else {
660
661 MPI_Send(atomData, 6, MPI_DOUBLE, 0,
662 myPotato, MPI_COMM_WORLD);
663 }
664
665 myPotato++;
666
667 }
668
669 currentIndex++;
670
671 }
672
673 }
674
675 sprintf( checkPointMsg,
676 "Sucessfully took a dump.\n");
677 MPIcheckPoint();
678
679 }
680
681
682
683 #endif // is_mpi
684 }
685
686 #ifdef IS_MPI
687
688 // a couple of functions to let us escape the write loop
689
690 void dWrite::DieDieDie( void ){
691
692 MPI_Finalize();
693 exit (0);
694 }
695
696 #endif //is_mpi