ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/io/DumpWriter.cpp
Revision: 1739
Committed: Mon Nov 15 18:02:15 2004 UTC (19 years, 8 months ago) by tim
File size: 16723 byte(s)
Log Message:
finish DumpReader, DumpWriter.Next Step is LJFF and integrators

File Contents

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