ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/io/DumpWriter.cpp
Revision: 1867
Committed: Tue Dec 7 23:08:14 2004 UTC (19 years, 6 months ago) by tim
File size: 16805 byte(s)
Log Message:
NPT in progress

File Contents

# Content
1 #include "io/DumpWriter.hpp"
2 #include "primitives/Molecule.hpp"
3 #include "utils/simError.h"
4
5 #ifdef IS_MPI
6 #include <mpi.h>
7 #endif //is_mpi
8
9 namespace oopse {
10
11 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
12 : info_(info), filename_(filename){
13 #ifdef IS_MPI
14
15 if (worldRank == 0) {
16 #endif // is_mpi
17
18 dumpFile_.open(filename_.c_str(), std::ios::out | std::ios::trunc);
19
20 if (!dumpFile_) {
21 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
22 filename_.c_str());
23 painCave.isFatal = 1;
24 simError();
25 }
26
27 #ifdef IS_MPI
28
29 }
30
31 sprintf(checkPointMsg, "Sucessfully opened output file for dumping.\n");
32 MPIcheckPoint();
33
34 #endif // is_mpi
35
36 }
37
38 DumpWriter::~DumpWriter() {
39
40 #ifdef IS_MPI
41
42 if (worldRank == 0) {
43 #endif // is_mpi
44
45 dumpFile_.close();
46
47 #ifdef IS_MPI
48
49 }
50
51 #endif // is_mpi
52
53 }
54
55 void DumpWriter::writeCommentLine(std::ostream& os, Snapshot* s) {
56
57 double currentTime;
58 Mat3x3d hmat;
59 double chi;
60 double integralOfChiDt;
61 Mat3x3d eta;
62
63 currentTime = s->getTime();
64 hmat = s->getHmat();
65 chi = s->getChi();
66 integralOfChiDt = s->getIntegralOfChiDt();
67 eta = s->getEta();
68
69 os << currentTime << ";\t"
70 << hmat(0, 0) << "\t" << hmat(1, 0) << "\t" << hmat(2, 0) << ";\t"
71 << hmat(0, 1) << "\t" << hmat(1, 1) << "\t" << hmat(2, 1) << ";\t"
72 << hmat(0, 2) << "\t" << hmat(1, 2) << "\t" << hmat(2, 2) << ";\t";
73
74 //write out additional parameters, such as chi and eta
75
76 os << chi << "\t" << integralOfChiDt << "\t;";
77
78 os << eta(0, 0) << "\t" << eta(1, 0) << "\t" << eta(2, 0) << ";\t"
79 << eta(0, 1) << "\t" << eta(1, 1) << "\t" << eta(2, 1) << ";\t"
80 << eta(0, 2) << "\t" << eta(1, 2) << "\t" << eta(2, 2) << ";";
81
82 os << std::endl;
83 }
84
85 void DumpWriter::writeFrame(std::ostream& os) {
86 const int BUFFERSIZE = 2000;
87 const int MINIBUFFERSIZE = 100;
88
89 char tempBuffer[BUFFERSIZE];
90 char writeLine[BUFFERSIZE];
91
92 Quat4d q;
93 Vector3d ji;
94 Vector3d pos;
95 Vector3d vel;
96
97 Molecule* mol;
98 StuntDouble* integrableObject;
99 SimInfo::MoleculeIterator mi;
100 Molecule::IntegrableObjectIterator ii;
101
102 int nTotObjects;
103 nTotObjects = info_->getNGlobalIntegrableObjects();
104
105 #ifndef IS_MPI
106
107
108 os << nTotObjects << "\n";
109
110 writeCommentLine(os, info_->getSnapshotManager()->getCurrentSnapshot());
111
112 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
113
114 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
115 integrableObject = mol->nextIntegrableObject(ii)) {
116
117
118 pos = integrableObject->getPos();
119 vel = integrableObject->getVel();
120
121 sprintf(tempBuffer, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
122 integrableObject->getType().c_str(),
123 pos[0], pos[1], pos[2],
124 vel[0], vel[1], vel[2]);
125
126 strcpy(writeLine, tempBuffer);
127
128 if (integrableObject->isDirectional()) {
129 q = integrableObject->getQ();
130 ji = integrableObject->getJ();
131
132 sprintf(tempBuffer, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
133 q[0], q[1], q[2], q[3],
134 ji[0], ji[1], ji[2]);
135 strcat(writeLine, tempBuffer);
136 } else {
137 strcat(writeLine, "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
138 }
139
140 os << writeLine;
141
142 }
143 }
144
145 #else // is_mpi
146 /*********************************************************************
147 * Documentation? You want DOCUMENTATION?
148 *
149 * Why all the potatoes below?
150 *
151 * To make a long story short, the original version of DumpWriter
152 * worked in the most inefficient way possible. Node 0 would
153 * poke each of the node for an individual atom's formatted data
154 * as node 0 worked its way down the global index. This was particularly
155 * inefficient since the method blocked all processors at every atom
156 * (and did it twice!).
157 *
158 * An intermediate version of DumpWriter could be described from Node
159 * zero's perspective as follows:
160 *
161 * 1) Have 100 of your friends stand in a circle.
162 * 2) When you say go, have all of them start tossing potatoes at
163 * you (one at a time).
164 * 3) Catch the potatoes.
165 *
166 * It was an improvement, but MPI has buffers and caches that could
167 * best be described in this analogy as "potato nets", so there's no
168 * need to block the processors atom-by-atom.
169 *
170 * This new and improved DumpWriter works in an even more efficient
171 * way:
172 *
173 * 1) Have 100 of your friend stand in a circle.
174 * 2) When you say go, have them start tossing 5-pound bags of
175 * potatoes at you.
176 * 3) Once you've caught a friend's bag of potatoes,
177 * toss them a spud to let them know they can toss another bag.
178 *
179 * How's THAT for documentation?
180 *
181 *********************************************************************/
182 const int masterNode = 0;
183
184 int * potatoes;
185 int myPotato;
186 int nProc;
187 int which_node;
188 double atomData[13];
189 int isDirectional;
190 const char * atomTypeString;
191 char MPIatomTypeString[MINIBUFFERSIZE];
192 int msgLen; // the length of message actually recieved at master nodes
193 int haveError;
194 MPI_Status istatus;
195 int nCurObj;
196
197 // code to find maximum tag value
198 int * tagub;
199 int flag;
200 int MAXTAG;
201 MPI_Attr_get(MPI_COMM_WORLD, MPI_TAG_UB, &tagub, &flag);
202
203 if (flag) {
204 MAXTAG = *tagub;
205 } else {
206 MAXTAG = 32767;
207 }
208
209 if (worldRank == masterNode) { //master node (node 0) is responsible for writing the dump file
210
211 // Node 0 needs a list of the magic potatoes for each processor;
212
213 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
214 potatoes = new int[nProc];
215
216 //write out the comment lines
217 for(i = 0; i < nProc; i++) {
218 potatoes[i] = 0;
219 }
220
221
222 os << nTotObjects << "\n";
223 writeCommentLine(info_->getSnapshotManager()->getCurrentSnapshot());
224
225 for(i = 0; i < info_->getNGlobalMolecules(); i++) {
226
227 // Get the Node number which has this atom;
228
229 which_node = info_->getMolToProc(i);
230
231 if (which_node != masterNode) { //current molecule is in slave node
232 if (potatoes[which_node] + 1 >= MAXTAG) {
233 // The potato was going to exceed the maximum value,
234 // so wrap this processor potato back to 0:
235
236 potatoes[which_node] = 0;
237 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node, 0,
238 MPI_COMM_WORLD);
239 }
240
241 myPotato = potatoes[which_node];
242
243 //recieve the number of integrableObject in current molecule
244 MPI_Recv(&nCurObj, 1, MPI_INT, which_node, myPotato,
245 MPI_COMM_WORLD, &istatus);
246 myPotato++;
247
248 for(int l = 0; l < nCurObj; l++) {
249 if (potatoes[which_node] + 2 >= MAXTAG) {
250 // The potato was going to exceed the maximum value,
251 // so wrap this processor potato back to 0:
252
253 potatoes[which_node] = 0;
254 MPI_Send(&potatoes[which_node], 1, MPI_INT, which_node,
255 0, MPI_COMM_WORLD);
256 }
257
258 MPI_Recv(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR,
259 which_node, myPotato, MPI_COMM_WORLD,
260 &istatus);
261
262 atomTypeString = MPIatomTypeString;
263
264 myPotato++;
265
266 MPI_Recv(atomData, 13, MPI_DOUBLE, which_node, myPotato,
267 MPI_COMM_WORLD, &istatus);
268 myPotato++;
269
270 MPI_Get_count(&istatus, MPI_DOUBLE, &msgLen);
271
272 if (msgLen == 13)
273 isDirectional = 1;
274 else
275 isDirectional = 0;
276
277 // If we've survived to here, format the line:
278
279 if (!isDirectional) {
280 sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
281 atomTypeString, atomData[0],
282 atomData[1], atomData[2],
283 atomData[3], atomData[4],
284 atomData[5]);
285
286 strcat(writeLine,
287 "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
288 } else {
289 sprintf(writeLine,
290 "%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",
291 atomTypeString,
292 atomData[0],
293 atomData[1],
294 atomData[2],
295 atomData[3],
296 atomData[4],
297 atomData[5],
298 atomData[6],
299 atomData[7],
300 atomData[8],
301 atomData[9],
302 atomData[10],
303 atomData[11],
304 atomData[12]);
305 }
306
307 os << writeLine;
308
309 } // end for(int l =0)
310
311 potatoes[which_node] = myPotato;
312 } else { //master node has current molecule
313
314 mol = info_->getMoleculeByGlobalIndex(i);
315
316 if (mol == NULL) {
317 strcpy(painCave.errMsg, "Molecule not found on node %d!", worldRank);
318 painCave.isFatal = 1;
319 simError();
320 }
321
322 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
323 integrableObject = mol->nextIntegrableObject(ii)) {
324
325 atomTypeString = integrableObject->getType().c_str();
326
327 pos = integrableObject->getPos();
328 vel = integrableObject->getVel();
329
330 atomData[0] = pos[0];
331 atomData[1] = pos[1];
332 atomData[2] = pos[2];
333
334 atomData[3] = vel[0];
335 atomData[4] = vel[1];
336 atomData[5] = vel[2];
337
338 isDirectional = 0;
339
340 if (integrableObject->isDirectional()) {
341 isDirectional = 1;
342
343 q = integrableObject->getQ();
344 ji = integrableObject->getJ();
345
346 for(int j = 0; j < 6; j++) {
347 atomData[j] = atomData[j];
348 }
349
350 atomData[6] = q[0];
351 atomData[7] = q[1];
352 atomData[8] = q[2];
353 atomData[9] = q[3];
354
355 atomData[10] = ji[0];
356 atomData[11] = ji[1];
357 atomData[12] = ji[2];
358 }
359
360 // If we've survived to here, format the line:
361
362 if (!isDirectional) {
363 sprintf(writeLine, "%s\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t",
364 atomTypeString, atomData[0],
365 atomData[1], atomData[2],
366 atomData[3], atomData[4],
367 atomData[5]);
368
369 strcat(writeLine,
370 "0.0\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0\n");
371 } else {
372 sprintf(writeLine,
373 "%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",
374 atomTypeString,
375 atomData[0],
376 atomData[1],
377 atomData[2],
378 atomData[3],
379 atomData[4],
380 atomData[5],
381 atomData[6],
382 atomData[7],
383 atomData[8],
384 atomData[9],
385 atomData[10],
386 atomData[11],
387 atomData[12]);
388 }
389
390
391 os << writeLine;
392
393 } //end for(iter = integrableObject.begin())
394 }
395 } //end for(i = 0; i < mpiSim->getNmol())
396
397 os.flush();
398
399 sprintf(checkPointMsg, "Sucessfully took a dump.\n");
400
401 MPIcheckPoint();
402
403 delete [] potatoes;
404 } else {
405
406 // worldRank != 0, so I'm a remote node.
407
408 // Set my magic potato to 0:
409
410 myPotato = 0;
411
412 for(i = 0; i < info_->getNGlobalMolecules(); i++) {
413
414 // Am I the node which has this integrableObject?
415 whichNode = info_->getMolToProc(i);
416 if (whichNode == worldRank) {
417 if (myPotato + 1 >= MAXTAG) {
418
419 // The potato was going to exceed the maximum value,
420 // so wrap this processor potato back to 0 (and block until
421 // node 0 says we can go:
422
423 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD,
424 &istatus);
425 }
426
427 mol = info_->getMoleculeByGlobalIndex(i);
428
429
430 nCurObj = mol->getNIntegrableObjects();
431
432 MPI_Send(&nCurObj, 1, MPI_INT, 0, myPotato, MPI_COMM_WORLD);
433 myPotato++;
434
435 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
436 integrableObject = mol->nextIntegrableObject(ii)) {
437
438 if (myPotato + 2 >= MAXTAG) {
439
440 // The potato was going to exceed the maximum value,
441 // so wrap this processor potato back to 0 (and block until
442 // node 0 says we can go:
443
444 MPI_Recv(&myPotato, 1, MPI_INT, 0, 0, MPI_COMM_WORLD,
445 &istatus);
446 }
447
448 atomTypeString = integrableObject->getType().c_str();
449
450 pos = integrableObject->getPos();
451 vel = integrableObject->getVel();
452
453 atomData[0] = pos[0];
454 atomData[1] = pos[1];
455 atomData[2] = pos[2];
456
457 atomData[3] = vel[0];
458 atomData[4] = vel[1];
459 atomData[5] = vel[2];
460
461 isDirectional = 0;
462
463 if (integrableObject->isDirectional()) {
464 isDirectional = 1;
465
466 q = integrableObject->getQ();
467 ji = integrableObject->getJ();
468
469 atomData[6] = q[0];
470 atomData[7] = q[1];
471 atomData[8] = q[2];
472 atomData[9] = q[3];
473
474 atomData[10] = ji[0];
475 atomData[11] = ji[1];
476 atomData[12] = ji[2];
477 }
478
479 strncpy(MPIatomTypeString, atomTypeString, MINIBUFFERSIZE);
480
481 // null terminate the std::string before sending (just in case):
482 MPIatomTypeString[MINIBUFFERSIZE - 1] = '\0';
483
484 MPI_Send(MPIatomTypeString, MINIBUFFERSIZE, MPI_CHAR, 0,
485 myPotato, MPI_COMM_WORLD);
486
487 myPotato++;
488
489 if (isDirectional) {
490 MPI_Send(atomData, 13, MPI_DOUBLE, 0, myPotato,
491 MPI_COMM_WORLD);
492 } else {
493 MPI_Send(atomData, 6, MPI_DOUBLE, 0, myPotato,
494 MPI_COMM_WORLD);
495 }
496
497 myPotato++;
498 }
499
500 }
501
502 }
503
504 sprintf(checkPointMsg, "Sucessfully took a dump.\n");
505 MPIcheckPoint();
506 }
507
508 #endif // is_mpi
509
510 }
511
512 }//end namespace oopse