ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
Revision: 1629
Committed: Wed Sep 14 21:15:17 2011 UTC (13 years, 7 months ago) by gezelter
File size: 16300 byte(s)
Log Message:
Merging changes from old branch into development branch

File Contents

# Content
1 /*
2 * Copyright (c) 2009 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41
42 #include "io/DumpWriter.hpp"
43 #include "primitives/Molecule.hpp"
44 #include "utils/simError.h"
45 #include "io/basic_teebuf.hpp"
46 #include "io/gzstream.hpp"
47 #include "io/Globals.hpp"
48
49
50 #ifdef IS_MPI
51 #include <mpi.h>
52 #endif //is_mpi
53
54 using namespace std;
55 namespace OpenMD {
56
57 DumpWriter::DumpWriter(SimInfo* info)
58 : info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
59
60 Globals* simParams = info->getSimParams();
61 needCompression_ = simParams->getCompressDumpFile();
62 needForceVector_ = simParams->getOutputForceVector();
63 needParticlePot_ = simParams->getOutputParticlePotential();
64 createDumpFile_ = true;
65 #ifdef HAVE_LIBZ
66 if (needCompression_) {
67 filename_ += ".gz";
68 eorFilename_ += ".gz";
69 }
70 #endif
71
72 #ifdef IS_MPI
73
74 if (worldRank == 0) {
75 #endif // is_mpi
76
77 dumpFile_ = createOStream(filename_);
78
79 if (!dumpFile_) {
80 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
81 filename_.c_str());
82 painCave.isFatal = 1;
83 simError();
84 }
85
86 #ifdef IS_MPI
87
88 }
89
90 #endif // is_mpi
91
92 }
93
94
95 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
96 : info_(info), filename_(filename){
97
98 Globals* simParams = info->getSimParams();
99 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
100
101 needCompression_ = simParams->getCompressDumpFile();
102 needForceVector_ = simParams->getOutputForceVector();
103 createDumpFile_ = true;
104 #ifdef HAVE_LIBZ
105 if (needCompression_) {
106 filename_ += ".gz";
107 eorFilename_ += ".gz";
108 }
109 #endif
110
111 #ifdef IS_MPI
112
113 if (worldRank == 0) {
114 #endif // is_mpi
115
116
117 dumpFile_ = createOStream(filename_);
118
119 if (!dumpFile_) {
120 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
121 filename_.c_str());
122 painCave.isFatal = 1;
123 simError();
124 }
125
126 #ifdef IS_MPI
127
128 }
129
130 #endif // is_mpi
131
132 }
133
134 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
135 : info_(info), filename_(filename){
136
137 Globals* simParams = info->getSimParams();
138 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
139
140 needCompression_ = simParams->getCompressDumpFile();
141 needForceVector_ = simParams->getOutputForceVector();
142 needParticlePot_ = simParams->getOutputParticlePotential();
143
144 #ifdef HAVE_LIBZ
145 if (needCompression_) {
146 filename_ += ".gz";
147 eorFilename_ += ".gz";
148 }
149 #endif
150
151 #ifdef IS_MPI
152
153 if (worldRank == 0) {
154 #endif // is_mpi
155
156 createDumpFile_ = writeDumpFile;
157 if (createDumpFile_) {
158 dumpFile_ = createOStream(filename_);
159
160 if (!dumpFile_) {
161 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
162 filename_.c_str());
163 painCave.isFatal = 1;
164 simError();
165 }
166 }
167 #ifdef IS_MPI
168
169 }
170
171
172 #endif // is_mpi
173
174 }
175
176 DumpWriter::~DumpWriter() {
177
178 #ifdef IS_MPI
179
180 if (worldRank == 0) {
181 #endif // is_mpi
182 if (createDumpFile_){
183 writeClosing(*dumpFile_);
184 delete dumpFile_;
185 }
186 #ifdef IS_MPI
187
188 }
189
190 #endif // is_mpi
191
192 }
193
194 void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
195
196 char buffer[1024];
197
198 os << " <FrameData>\n";
199
200 RealType currentTime = s->getTime();
201
202 if (isinf(currentTime) || isnan(currentTime)) {
203 sprintf( painCave.errMsg,
204 "DumpWriter detected a numerical error writing the time");
205 painCave.isFatal = 1;
206 simError();
207 }
208
209 sprintf(buffer, " Time: %.10g\n", currentTime);
210 os << buffer;
211
212 Mat3x3d hmat;
213 hmat = s->getHmat();
214
215 for (unsigned int i = 0; i < 3; i++) {
216 for (unsigned int j = 0; j < 3; j++) {
217 if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {
218 sprintf( painCave.errMsg,
219 "DumpWriter detected a numerical error writing the box");
220 painCave.isFatal = 1;
221 simError();
222 }
223 }
224 }
225
226 sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
227 hmat(0, 0), hmat(1, 0), hmat(2, 0),
228 hmat(0, 1), hmat(1, 1), hmat(2, 1),
229 hmat(0, 2), hmat(1, 2), hmat(2, 2));
230 os << buffer;
231
232 RealType chi = s->getChi();
233 RealType integralOfChiDt = s->getIntegralOfChiDt();
234 if (isinf(chi) || isnan(chi) ||
235 isinf(integralOfChiDt) || isnan(integralOfChiDt)) {
236 sprintf( painCave.errMsg,
237 "DumpWriter detected a numerical error writing the thermostat");
238 painCave.isFatal = 1;
239 simError();
240 }
241 sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
242 os << buffer;
243
244 Mat3x3d eta;
245 eta = s->getEta();
246
247 for (unsigned int i = 0; i < 3; i++) {
248 for (unsigned int j = 0; j < 3; j++) {
249 if (isinf(eta(i,j)) || isnan(eta(i,j))) {
250 sprintf( painCave.errMsg,
251 "DumpWriter detected a numerical error writing the barostat");
252 painCave.isFatal = 1;
253 simError();
254 }
255 }
256 }
257
258 sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
259 eta(0, 0), eta(1, 0), eta(2, 0),
260 eta(0, 1), eta(1, 1), eta(2, 1),
261 eta(0, 2), eta(1, 2), eta(2, 2));
262 os << buffer;
263
264 os << " </FrameData>\n";
265 }
266
267 void DumpWriter::writeFrame(std::ostream& os) {
268
269 #ifdef IS_MPI
270 MPI_Status istatus;
271 #endif
272
273 Molecule* mol;
274 StuntDouble* integrableObject;
275 SimInfo::MoleculeIterator mi;
276 Molecule::IntegrableObjectIterator ii;
277
278 #ifndef IS_MPI
279 os << " <Snapshot>\n";
280
281 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
282
283 os << " <StuntDoubles>\n";
284 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
285
286
287 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
288 integrableObject = mol->nextIntegrableObject(ii)) {
289 os << prepareDumpLine(integrableObject);
290
291 }
292 }
293 os << " </StuntDoubles>\n";
294
295 os << " </Snapshot>\n";
296
297 os.flush();
298 #else
299 //every node prepares the dump lines for integrable objects belong to itself
300 std::string buffer;
301 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
302
303
304 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
305 integrableObject = mol->nextIntegrableObject(ii)) {
306 buffer += prepareDumpLine(integrableObject);
307 }
308 }
309
310 const int masterNode = 0;
311 int nProc;
312 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
313 if (worldRank == masterNode) {
314 os << " <Snapshot>\n";
315 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
316 os << " <StuntDoubles>\n";
317
318 os << buffer;
319
320 for (int i = 1; i < nProc; ++i) {
321
322 // receive the length of the string buffer that was
323 // prepared by processor i
324
325 MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
326 int recvLength;
327 MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
328 char* recvBuffer = new char[recvLength];
329 if (recvBuffer == NULL) {
330 } else {
331 MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
332 os << recvBuffer;
333 delete [] recvBuffer;
334 }
335 }
336 os << " </StuntDoubles>\n";
337
338 os << " </Snapshot>\n";
339 os.flush();
340 } else {
341 int sendBufferLength = buffer.size() + 1;
342 int myturn = 0;
343 for (int i = 1; i < nProc; ++i){
344 MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
345 if (myturn == worldRank){
346 MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
347 MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
348 }
349 }
350 }
351
352 #endif // is_mpi
353
354 }
355
356 std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
357
358 int index = integrableObject->getGlobalIntegrableObjectIndex();
359 std::string type("pv");
360 std::string line;
361 char tempBuffer[4096];
362
363 Vector3d pos;
364 Vector3d vel;
365 pos = integrableObject->getPos();
366
367 if (isinf(pos[0]) || isnan(pos[0]) ||
368 isinf(pos[1]) || isnan(pos[1]) ||
369 isinf(pos[2]) || isnan(pos[2]) ) {
370 sprintf( painCave.errMsg,
371 "DumpWriter detected a numerical error writing the position"
372 " for object %d", index);
373 painCave.isFatal = 1;
374 simError();
375 }
376
377 vel = integrableObject->getVel();
378
379 if (isinf(vel[0]) || isnan(vel[0]) ||
380 isinf(vel[1]) || isnan(vel[1]) ||
381 isinf(vel[2]) || isnan(vel[2]) ) {
382 sprintf( painCave.errMsg,
383 "DumpWriter detected a numerical error writing the velocity"
384 " for object %d", index);
385 painCave.isFatal = 1;
386 simError();
387 }
388
389 sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
390 pos[0], pos[1], pos[2],
391 vel[0], vel[1], vel[2]);
392 line += tempBuffer;
393
394 if (integrableObject->isDirectional()) {
395 type += "qj";
396 Quat4d q;
397 Vector3d ji;
398 q = integrableObject->getQ();
399
400 if (isinf(q[0]) || isnan(q[0]) ||
401 isinf(q[1]) || isnan(q[1]) ||
402 isinf(q[2]) || isnan(q[2]) ||
403 isinf(q[3]) || isnan(q[3]) ) {
404 sprintf( painCave.errMsg,
405 "DumpWriter detected a numerical error writing the quaternion"
406 " for object %d", index);
407 painCave.isFatal = 1;
408 simError();
409 }
410
411 ji = integrableObject->getJ();
412
413 if (isinf(ji[0]) || isnan(ji[0]) ||
414 isinf(ji[1]) || isnan(ji[1]) ||
415 isinf(ji[2]) || isnan(ji[2]) ) {
416 sprintf( painCave.errMsg,
417 "DumpWriter detected a numerical error writing the angular"
418 " momentum for object %d", index);
419 painCave.isFatal = 1;
420 simError();
421 }
422
423 sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
424 q[0], q[1], q[2], q[3],
425 ji[0], ji[1], ji[2]);
426 line += tempBuffer;
427 }
428
429 if (needForceVector_) {
430 type += "f";
431 Vector3d frc;
432
433 frc = integrableObject->getFrc();
434
435 if (isinf(frc[0]) || isnan(frc[0]) ||
436 isinf(frc[1]) || isnan(frc[1]) ||
437 isinf(frc[2]) || isnan(frc[2]) ) {
438 sprintf( painCave.errMsg,
439 "DumpWriter detected a numerical error writing the force"
440 " for object %d", index);
441 painCave.isFatal = 1;
442 simError();
443 }
444 sprintf(tempBuffer, " %13e %13e %13e",
445 frc[0], frc[1], frc[2]);
446 line += tempBuffer;
447
448 if (integrableObject->isDirectional()) {
449 type += "t";
450 Vector3d trq;
451
452 trq = integrableObject->getTrq();
453
454 if (isinf(trq[0]) || isnan(trq[0]) ||
455 isinf(trq[1]) || isnan(trq[1]) ||
456 isinf(trq[2]) || isnan(trq[2]) ) {
457 sprintf( painCave.errMsg,
458 "DumpWriter detected a numerical error writing the torque"
459 " for object %d", index);
460 painCave.isFatal = 1;
461 simError();
462 }
463
464 sprintf(tempBuffer, " %13e %13e %13e",
465 trq[0], trq[1], trq[2]);
466 line += tempBuffer;
467 }
468 }
469 if (needParticlePot_) {
470 type += "u";
471 RealType particlePot;
472
473 particlePot = integrableObject->getParticlePot();
474
475 if (isinf(particlePot) || isnan(particlePot)) {
476 sprintf( painCave.errMsg,
477 "DumpWriter detected a numerical error writing the particle "
478 " potential for object %d", index);
479 painCave.isFatal = 1;
480 simError();
481 }
482 sprintf(tempBuffer, " %13e", particlePot);
483 line += tempBuffer;
484 }
485
486 sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
487 return std::string(tempBuffer);
488 }
489
490 void DumpWriter::writeDump() {
491 writeFrame(*dumpFile_);
492 }
493
494 void DumpWriter::writeEor() {
495 std::ostream* eorStream;
496
497 #ifdef IS_MPI
498 if (worldRank == 0) {
499 #endif // is_mpi
500
501 eorStream = createOStream(eorFilename_);
502
503 #ifdef IS_MPI
504 }
505 #endif // is_mpi
506
507 writeFrame(*eorStream);
508
509 #ifdef IS_MPI
510 if (worldRank == 0) {
511 #endif // is_mpi
512 writeClosing(*eorStream);
513 delete eorStream;
514 #ifdef IS_MPI
515 }
516 #endif // is_mpi
517
518 }
519
520
521 void DumpWriter::writeDumpAndEor() {
522 std::vector<std::streambuf*> buffers;
523 std::ostream* eorStream;
524 #ifdef IS_MPI
525 if (worldRank == 0) {
526 #endif // is_mpi
527
528 buffers.push_back(dumpFile_->rdbuf());
529
530 eorStream = createOStream(eorFilename_);
531
532 buffers.push_back(eorStream->rdbuf());
533
534 #ifdef IS_MPI
535 }
536 #endif // is_mpi
537
538 TeeBuf tbuf(buffers.begin(), buffers.end());
539 std::ostream os(&tbuf);
540
541 writeFrame(os);
542
543 #ifdef IS_MPI
544 if (worldRank == 0) {
545 #endif // is_mpi
546 writeClosing(*eorStream);
547 delete eorStream;
548 #ifdef IS_MPI
549 }
550 #endif // is_mpi
551
552 }
553
554 std::ostream* DumpWriter::createOStream(const std::string& filename) {
555
556 std::ostream* newOStream;
557 #ifdef HAVE_LIBZ
558 if (needCompression_) {
559 newOStream = new ogzstream(filename.c_str());
560 } else {
561 newOStream = new std::ofstream(filename.c_str());
562 }
563 #else
564 newOStream = new std::ofstream(filename.c_str());
565 #endif
566 //write out MetaData first
567 (*newOStream) << "<OpenMD version=1>" << std::endl;
568 (*newOStream) << " <MetaData>" << std::endl;
569 (*newOStream) << info_->getRawMetaData();
570 (*newOStream) << " </MetaData>" << std::endl;
571 return newOStream;
572 }
573
574 void DumpWriter::writeClosing(std::ostream& os) {
575
576 os << "</OpenMD>\n";
577 os.flush();
578 }
579
580 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date