ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/DumpWriter.cpp
Revision: 907
Committed: Thu Jan 8 17:40:56 2004 UTC (20 years, 6 months ago) by gezelter
File size: 19093 byte(s)
Log Message:
First Stab at fixing DumpWriter

File Contents

# Content
1 #define _FILE_OFFSET_BITS 64
2
3 #include <string.h>
4 #include <iostream>
5 #include <fstream>
6
7 #ifdef IS_MPI
8 #include <mpi.h>
9 #include "mpiSimulation.hpp"
10
11 namespace dWrite{
12 void DieDieDie( void );
13 }
14
15 using namespace dWrite;
16 #endif //is_mpi
17
18 #include "ReadWrite.hpp"
19 #include "simError.h"
20
21 DumpWriter::DumpWriter( SimInfo* the_entry_plug ){
22
23 entry_plug = the_entry_plug;
24
25 #ifdef IS_MPI
26 if(worldRank == 0 ){
27 #endif // is_mpi
28
29 strcpy( outName, entry_plug->sampleName );
30
31 outFile.open(outName, ios::out | ios::trunc );
32
33 if( !outFile ){
34
35 sprintf( painCave.errMsg,
36 "Could not open \"%s\" for dump output.\n",
37 outName);
38 painCave.isFatal = 1;
39 simError();
40 }
41
42 //outFile.setf( ios::scientific );
43
44 #ifdef IS_MPI
45 }
46
47 sprintf( checkPointMsg,
48 "Sucessfully opened output file for dumping.\n");
49 MPIcheckPoint();
50 #endif // is_mpi
51 }
52
53 DumpWriter::~DumpWriter( ){
54
55 #ifdef IS_MPI
56 if(worldRank == 0 ){
57 #endif // is_mpi
58
59 outFile.close();
60
61 #ifdef IS_MPI
62 }
63 #endif // is_mpi
64 }
65
66 void DumpWriter::writeDump( double currentTime ){
67
68 const int BUFFERSIZE = 2000;
69 const int MINIBUFFERSIZE = 10;
70
71 char tempBuffer[BUFFERSIZE];
72 char writeLine[BUFFERSIZE];
73
74 int i;
75 #ifdef IS_MPI
76 int j, which_node, done, which_atom, local_index;
77 double atomTransData[6];
78 double atomOrientData[7];
79 int isDirectional;
80 char* atomTypeString;
81 int me;
82 int atomTypeTag;
83 int atomIsDirectionalTag;
84 int atomTransDataTag;
85 int atomOrientDataTag;
86 #else //is_mpi
87 int nAtoms = entry_plug->n_atoms;
88 #endif //is_mpi
89
90 double q[4];
91 DirectionalAtom* dAtom;
92 Atom** atoms = entry_plug->atoms;
93 double pos[3], vel[3];
94
95 // write current frame to the eor file
96
97 this->writeFinal( currentTime );
98
99 #ifndef IS_MPI
100
101 outFile << nAtoms << "\n";
102
103 outFile << currentTime << ";\t"
104 << entry_plug->Hmat[0][0] << "\t"
105 << entry_plug->Hmat[1][0] << "\t"
106 << entry_plug->Hmat[2][0] << ";\t"
107
108 << entry_plug->Hmat[0][1] << "\t"
109 << entry_plug->Hmat[1][1] << "\t"
110 << entry_plug->Hmat[2][1] << ";\t"
111
112 << entry_plug->Hmat[0][2] << "\t"
113 << entry_plug->Hmat[1][2] << "\t"
114 << entry_plug->Hmat[2][2] << ";";
115 //write out additional parameters, such as chi and eta
116 outFile << entry_plug->the_integrator->getAdditionalParameters();
117 outFile << endl;
118
119 for( i=0; i<nAtoms; i++ ){
120
121 atoms[i]->getPos(pos);
122 atoms[i]->getVel(vel);
123
124 sprintf( tempBuffer,
125 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
126 atoms[i]->getType(),
127 pos[0],
128 pos[1],
129 pos[2],
130 vel[0],
131 vel[1],
132 vel[2]);
133 strcpy( writeLine, tempBuffer );
134
135 if( atoms[i]->isDirectional() ){
136
137 dAtom = (DirectionalAtom *)atoms[i];
138 dAtom->getQ( q );
139
140 sprintf( tempBuffer,
141 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
142 q[0],
143 q[1],
144 q[2],
145 q[3],
146 dAtom->getJx(),
147 dAtom->getJy(),
148 dAtom->getJz());
149 strcat( writeLine, tempBuffer );
150 }
151 else
152 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
153
154 outFile << writeLine;
155 }
156 outFile.flush();
157
158 #else // is_mpi
159
160 // first thing first, suspend fatalities.
161 painCave.isEventLoop = 1;
162
163 int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
164 int haveError;
165
166 MPI_Status istatus;
167 int *AtomToProcMap = mpiSim->getAtomToProcMap();
168
169 // write out header and node 0's coordinates
170
171 if( worldRank == 0 ){
172 outFile << mpiSim->getTotAtoms() << "\n";
173
174 outFile << currentTime << ";\t"
175 << entry_plug->Hmat[0][0] << "\t"
176 << entry_plug->Hmat[1][0] << "\t"
177 << entry_plug->Hmat[2][0] << ";\t"
178
179 << entry_plug->Hmat[0][1] << "\t"
180 << entry_plug->Hmat[1][1] << "\t"
181 << entry_plug->Hmat[2][1] << ";\t"
182
183 << entry_plug->Hmat[0][2] << "\t"
184 << entry_plug->Hmat[1][2] << "\t"
185 << entry_plug->Hmat[2][2] << ";";
186
187 outFile << entry_plug->the_integrator->getAdditionalParameters();
188 outFile << endl;
189 outFile.flush();
190 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
191 // Get the Node number which has this atom;
192
193 which_node = AtomToProcMap[i];
194
195 if (which_node != 0) {
196
197 atomTypeTag = 4*i;
198 atomIsDirectionalTag = 4*i + 1;
199 atomTransDataTag = 4*i + 2;
200 atomOrientDataTag = 4*i + 3;
201
202 MPI_Recv(atomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
203 atomTypeTag, MPI_COMM_WORLD, &istatus);
204
205 MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
206 atomIsDirectionalTag, MPI_COMM_WORLD, &istatus);
207
208 MPI_Recv(atomTransData, 6, MPI_DOUBLE, which_node,
209 atomTransDataTag, MPI_COMM_WORLD, &istatus);
210
211 if (isDirectional) {
212
213 MPI_Recv(atomOrientData, 7, MPI_DOUBLE, which_node,
214 atomOrientDataTag, MPI_COMM_WORLD, &istatus);
215
216 }
217
218 } else {
219
220 haveError = 0;
221 which_atom = i;
222 local_index=-1;
223
224 for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
225 if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
226 }
227
228 if (local_index != -1) {
229
230 atomTypeString = atoms[local_index]->getType();
231
232 atoms[local_index]->getPos(pos);
233 atoms[local_index]->getVel(vel);
234
235 atomTransData[0] = pos[0];
236 atomTransData[1] = pos[1];
237 atomTransData[2] = pos[2];
238
239 atomTransData[3] = vel[0];
240 atomTransData[4] = vel[1];
241 atomTransData[5] = vel[2];
242
243 isDirectional = 0;
244
245 if( atoms[local_index]->isDirectional() ){
246
247 isDirectional = 1;
248
249 dAtom = (DirectionalAtom *)atoms[local_index];
250 dAtom->getQ( q );
251
252 atomOrientData[0] = q[0];
253 atomOrientData[1] = q[1];
254 atomOrientData[2] = q[2];
255 atomOrientData[3] = q[3];
256
257 atomOrientData[4] = dAtom->getJx();
258 atomOrientData[5] = dAtom->getJy();
259 atomOrientData[6] = dAtom->getJz();
260 }
261
262 } else {
263 sprintf(painCave.errMsg,
264 "Atom %d not found on processor %d\n",
265 i, worldRank );
266 haveError= 1;
267 simError();
268 }
269
270 if(haveError) DieDieDie();
271
272 // If we've survived to here, format the line:
273
274 sprintf( tempBuffer,
275 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
276 atomTypeString,
277 atomTransData[0],
278 atomTransData[1],
279 atomTransData[2],
280 atomTransData[3],
281 atomTransData[4],
282 atomTransData[5]);
283
284 strcpy( writeLine, tempBuffer );
285
286 if (isDirectional) {
287
288 sprintf( tempBuffer,
289 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
290 atomOrientData[0],
291 atomOrientData[1],
292 atomOrientData[2],
293 atomOrientData[3],
294 atomOrientData[4],
295 atomOrientData[5],
296 atomOrientData[6]);
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
303 outFile << writeLine;
304 outFile.flush();
305 }
306 }
307
308 outFile.flush();
309 sprintf( checkPointMsg,
310 "Sucessfully took a dump.\n");
311 MPIcheckPoint();
312
313 } else {
314
315 // worldRank != 0, so I'm a remote node.
316
317 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
318
319 // Am I the node which has this atom?
320
321 if (AtomToProcMap[i] == worldRank) {
322
323 local_index=-1;
324 for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
325 if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
326 }
327 if (local_index != -1) {
328
329 atomTypeString = atoms[local_index]->getType();
330
331 atoms[local_index]->getPos(pos);
332 atoms[local_index]->getVel(vel);
333
334 atomTransData[0] = pos[0];
335 atomTransData[1] = pos[1];
336 atomTransData[2] = pos[2];
337
338 atomTransData[3] = vel[0];
339 atomTransData[4] = vel[1];
340 atomTransData[5] = vel[2];
341
342 isDirectional = 0;
343
344 if( atoms[local_index]->isDirectional() ){
345
346 isDirectional = 1;
347
348 dAtom = (DirectionalAtom *)atoms[local_index];
349 dAtom->getQ( q );
350
351 atomOrientData[0] = q[0];
352 atomOrientData[1] = q[1];
353 atomOrientData[2] = q[2];
354 atomOrientData[3] = q[3];
355
356 atomOrientData[4] = dAtom->getJx();
357 atomOrientData[5] = dAtom->getJy();
358 atomOrientData[6] = dAtom->getJz();
359 }
360
361 } else {
362 sprintf(painCave.errMsg,
363 "Atom %d not found on processor %d\n",
364 i, worldRank );
365 haveError= 1;
366 simError();
367 }
368
369 // I've survived this far, so send off the data!
370
371 atomTypeTag = 4*i;
372 atomIsDirectionalTag = 4*i + 1;
373 atomTransDataTag = 4*i + 2;
374 atomOrientDataTag = 4*i + 3;
375
376 MPI_Send(atomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
377 atomTypeTag, MPI_COMM_WORLD);
378
379 MPI_Send(&isDirectional, 1, MPI_INT, 0,
380 atomIsDirectionalTag, MPI_COMM_WORLD);
381
382 MPI_Send(atomTransData, 6, MPI_DOUBLE, 0,
383 atomTransDataTag, MPI_COMM_WORLD);
384
385 if (isDirectional) {
386
387 MPI_Send(atomOrientData, 7, MPI_DOUBLE, 0,
388 atomOrientDataTag, MPI_COMM_WORLD);
389
390 }
391
392 }
393 }
394
395 sprintf( checkPointMsg,
396 "Sucessfully took a dump.\n");
397 MPIcheckPoint();
398
399 }
400
401 painCave.isEventLoop = 0;
402
403 #endif // is_mpi
404 }
405
406 void DumpWriter::writeFinal(double finalTime){
407
408 char finalName[500];
409 ofstream finalOut;
410
411 const int BUFFERSIZE = 2000;
412 const int MINIBUFFERSIZE = 10;
413 char tempBuffer[BUFFERSIZE];
414 char writeLine[BUFFERSIZE];
415
416 double q[4];
417 DirectionalAtom* dAtom;
418 Atom** atoms = entry_plug->atoms;
419 int i;
420 #ifdef IS_MPI
421 int j, which_node, done, which_atom, local_index;
422 double atomTransData[6];
423 double atomOrientData[7];
424 int isDirectional;
425 char* atomTypeString;
426 int atomTypeTag;
427 int atomIsDirectionalTag;
428 int atomTransDataTag;
429 int atomOrientDataTag;
430 #else //is_mpi
431 int nAtoms = entry_plug->n_atoms;
432 #endif //is_mpi
433
434 double pos[3], vel[3];
435
436 #ifdef IS_MPI
437 if(worldRank == 0 ){
438 #endif // is_mpi
439
440 strcpy( finalName, entry_plug->finalName );
441
442 finalOut.open( finalName, ios::out | ios::trunc );
443 if( !finalOut ){
444 sprintf( painCave.errMsg,
445 "Could not open \"%s\" for final dump output.\n",
446 finalName );
447 painCave.isFatal = 1;
448 simError();
449 }
450
451 // finalOut.setf( ios::scientific );
452
453 #ifdef IS_MPI
454 }
455
456 sprintf(checkPointMsg,"Opened file for final configuration\n");
457 MPIcheckPoint();
458
459 #endif //is_mpi
460
461
462 #ifndef IS_MPI
463
464 finalOut << nAtoms << "\n";
465
466 finalOut << finalTime << ";\t"
467 << entry_plug->Hmat[0][0] << "\t"
468 << entry_plug->Hmat[1][0] << "\t"
469 << entry_plug->Hmat[2][0] << ";\t"
470
471 << entry_plug->Hmat[0][1] << "\t"
472 << entry_plug->Hmat[1][1] << "\t"
473 << entry_plug->Hmat[2][1] << ";\t"
474
475 << entry_plug->Hmat[0][2] << "\t"
476 << entry_plug->Hmat[1][2] << "\t"
477 << entry_plug->Hmat[2][2] << ";";
478
479 //write out additional parameters, such as chi and eta
480 finalOut << entry_plug->the_integrator->getAdditionalParameters();
481 finalOut << endl;
482
483 for( i=0; i<nAtoms; i++ ){
484
485 atoms[i]->getPos(pos);
486 atoms[i]->getVel(vel);
487
488 sprintf( tempBuffer,
489 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
490 atoms[i]->getType(),
491 pos[0],
492 pos[1],
493 pos[2],
494 vel[0],
495 vel[1],
496 vel[2]);
497 strcpy( writeLine, tempBuffer );
498
499 if( atoms[i]->isDirectional() ){
500
501 dAtom = (DirectionalAtom *)atoms[i];
502 dAtom->getQ( q );
503
504 sprintf( tempBuffer,
505 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
506 q[0],
507 q[1],
508 q[2],
509 q[3],
510 dAtom->getJx(),
511 dAtom->getJy(),
512 dAtom->getJz());
513 strcat( writeLine, tempBuffer );
514 }
515 else
516 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
517
518 finalOut << writeLine;
519 }
520 finalOut.flush();
521 finalOut.close();
522
523 #else // is_mpi
524
525 // first thing first, suspend fatalities.
526 painCave.isEventLoop = 1;
527
528 int myStatus; // 1 = wakeup & success; 0 = error; -1 = AllDone
529 int haveError;
530
531 MPI_Status istatus;
532 int *AtomToProcMap = mpiSim->getAtomToProcMap();
533
534 // write out header and node 0's coordinates
535
536 if( worldRank == 0 ){
537 finalOut << mpiSim->getTotAtoms() << "\n";
538
539 finalOut << finalTime << ";\t"
540 << entry_plug->Hmat[0][0] << "\t"
541 << entry_plug->Hmat[1][0] << "\t"
542 << entry_plug->Hmat[2][0] << ";\t"
543
544 << entry_plug->Hmat[0][1] << "\t"
545 << entry_plug->Hmat[1][1] << "\t"
546 << entry_plug->Hmat[2][1] << ";\t"
547
548 << entry_plug->Hmat[0][2] << "\t"
549 << entry_plug->Hmat[1][2] << "\t"
550 << entry_plug->Hmat[2][2] << ";";
551
552 finalOut << entry_plug->the_integrator->getAdditionalParameters();
553 finalOut << endl;
554 finalOut.flush();
555 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
556 // Get the Node number which has this atom;
557
558 which_node = AtomToProcMap[i];
559
560 if (which_node != 0) {
561
562 atomTypeTag = 4*i;
563 atomIsDirectionalTag = 4*i + 1;
564 atomTransDataTag = 4*i + 2;
565 atomOrientDataTag = 4*i + 3;
566
567 MPI_Recv(atomTypeString, MINIBUFFERSIZE, MPI_CHAR, which_node,
568 atomTypeTag, MPI_COMM_WORLD, &istatus);
569
570 MPI_Recv(&isDirectional, 1, MPI_INT, which_node,
571 atomIsDirectionalTag, MPI_COMM_WORLD, &istatus);
572
573 MPI_Recv(atomTransData, 6, MPI_DOUBLE, which_node,
574 atomTransDataTag, MPI_COMM_WORLD, &istatus);
575
576 if (isDirectional) {
577
578 MPI_Recv(atomOrientData, 7, MPI_DOUBLE, which_node,
579 atomOrientDataTag, MPI_COMM_WORLD, &istatus);
580
581 }
582
583 } else {
584
585 haveError = 0;
586 which_atom = i;
587 local_index=-1;
588
589 for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
590 if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
591 }
592
593 if (local_index != -1) {
594
595 atomTypeString = atoms[local_index]->getType();
596
597 atoms[local_index]->getPos(pos);
598 atoms[local_index]->getVel(vel);
599
600 atomTransData[0] = pos[0];
601 atomTransData[1] = pos[1];
602 atomTransData[2] = pos[2];
603
604 atomTransData[3] = vel[0];
605 atomTransData[4] = vel[1];
606 atomTransData[5] = vel[2];
607
608 isDirectional = 0;
609
610 if( atoms[local_index]->isDirectional() ){
611
612 isDirectional = 1;
613
614 dAtom = (DirectionalAtom *)atoms[local_index];
615 dAtom->getQ( q );
616
617 atomOrientData[0] = q[0];
618 atomOrientData[1] = q[1];
619 atomOrientData[2] = q[2];
620 atomOrientData[3] = q[3];
621
622 atomOrientData[4] = dAtom->getJx();
623 atomOrientData[5] = dAtom->getJy();
624 atomOrientData[6] = dAtom->getJz();
625 }
626
627 } else {
628 sprintf(painCave.errMsg,
629 "Atom %d not found on processor %d\n",
630 i, worldRank );
631 haveError= 1;
632 simError();
633 }
634
635 if(haveError) DieDieDie();
636
637 // If we've survived to here, format the line:
638
639 sprintf( tempBuffer,
640 "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
641 atomTypeString,
642 atomTransData[0],
643 atomTransData[1],
644 atomTransData[2],
645 atomTransData[3],
646 atomTransData[4],
647 atomTransData[5]);
648
649 strcpy( writeLine, tempBuffer );
650
651 if (isDirectional) {
652
653 sprintf( tempBuffer,
654 "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
655 atomOrientData[0],
656 atomOrientData[1],
657 atomOrientData[2],
658 atomOrientData[3],
659 atomOrientData[4],
660 atomOrientData[5],
661 atomOrientData[6]);
662 strcat( writeLine, tempBuffer );
663
664 } else {
665 strcat( writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n" );
666 }
667
668 finalOut << writeLine;
669 finalOut.flush();
670 }
671 }
672
673 finalOut.flush();
674 sprintf( checkPointMsg,
675 "Sucessfully took a dump.\n");
676 MPIcheckPoint();
677
678 } else {
679
680 // worldRank != 0, so I'm a remote node.
681
682 for (i = 0 ; i < mpiSim->getTotAtoms(); i++ ) {
683
684 // Am I the node which has this atom?
685
686 if (AtomToProcMap[i] == worldRank) {
687
688 local_index=-1;
689 for (j=0; (j<mpiSim->getMyNlocal()) && (local_index < 0); j++) {
690 if (atoms[j]->getGlobalIndex() == which_atom) local_index = j;
691 }
692 if (local_index != -1) {
693
694 atomTypeString = atoms[local_index]->getType();
695
696 atoms[local_index]->getPos(pos);
697 atoms[local_index]->getVel(vel);
698
699 atomTransData[0] = pos[0];
700 atomTransData[1] = pos[1];
701 atomTransData[2] = pos[2];
702
703 atomTransData[3] = vel[0];
704 atomTransData[4] = vel[1];
705 atomTransData[5] = vel[2];
706
707 isDirectional = 0;
708
709 if( atoms[local_index]->isDirectional() ){
710
711 isDirectional = 1;
712
713 dAtom = (DirectionalAtom *)atoms[local_index];
714 dAtom->getQ( q );
715
716 atomOrientData[0] = q[0];
717 atomOrientData[1] = q[1];
718 atomOrientData[2] = q[2];
719 atomOrientData[3] = q[3];
720
721 atomOrientData[4] = dAtom->getJx();
722 atomOrientData[5] = dAtom->getJy();
723 atomOrientData[6] = dAtom->getJz();
724 }
725
726 } else {
727 sprintf(painCave.errMsg,
728 "Atom %d not found on processor %d\n",
729 i, worldRank );
730 haveError= 1;
731 simError();
732 }
733
734 // I've survived this far, so send off the data!
735
736 atomTypeTag = 4*i;
737 atomIsDirectionalTag = 4*i + 1;
738 atomTransDataTag = 4*i + 2;
739 atomOrientDataTag = 4*i + 3;
740
741 MPI_Send(atomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
742 atomTypeTag, MPI_COMM_WORLD);
743
744 MPI_Send(&isDirectional, 1, MPI_INT, 0,
745 atomIsDirectionalTag, MPI_COMM_WORLD);
746
747 MPI_Send(atomTransData, 6, MPI_DOUBLE, 0,
748 atomTransDataTag, MPI_COMM_WORLD);
749
750 if (isDirectional) {
751
752 MPI_Send(atomOrientData, 7, MPI_DOUBLE, 0,
753 atomOrientDataTag, MPI_COMM_WORLD);
754
755 }
756
757 }
758 }
759
760 sprintf( checkPointMsg,
761 "Sucessfully wrote final file.\n");
762 MPIcheckPoint();
763
764 }
765
766 painCave.isEventLoop = 0;
767
768 if( worldRank == 0 ) finalOut.close();
769 #endif // is_mpi
770 }
771
772
773
774 #ifdef IS_MPI
775
776 // a couple of functions to let us escape the write loop
777
778 void dWrite::DieDieDie( void ){
779
780 MPI_Finalize();
781 exit (0);
782 }
783
784 #endif //is_mpi