ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/io/DumpWriter.cpp
Revision: 1437
Committed: Wed Apr 21 14:59:18 2010 UTC (15 years ago) by gezelter
File size: 15426 byte(s)
Log Message:
Adding nan/inf detection to DumpWriter, fixing one message

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 createDumpFile_ = true;
64 #ifdef HAVE_LIBZ
65 if (needCompression_) {
66 filename_ += ".gz";
67 eorFilename_ += ".gz";
68 }
69 #endif
70
71 #ifdef IS_MPI
72
73 if (worldRank == 0) {
74 #endif // is_mpi
75
76 dumpFile_ = createOStream(filename_);
77
78 if (!dumpFile_) {
79 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
80 filename_.c_str());
81 painCave.isFatal = 1;
82 simError();
83 }
84
85 #ifdef IS_MPI
86
87 }
88
89 #endif // is_mpi
90
91 }
92
93
94 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename)
95 : info_(info), filename_(filename){
96
97 Globals* simParams = info->getSimParams();
98 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
99
100 needCompression_ = simParams->getCompressDumpFile();
101 needForceVector_ = simParams->getOutputForceVector();
102 createDumpFile_ = true;
103 #ifdef HAVE_LIBZ
104 if (needCompression_) {
105 filename_ += ".gz";
106 eorFilename_ += ".gz";
107 }
108 #endif
109
110 #ifdef IS_MPI
111
112 if (worldRank == 0) {
113 #endif // is_mpi
114
115
116 dumpFile_ = createOStream(filename_);
117
118 if (!dumpFile_) {
119 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
120 filename_.c_str());
121 painCave.isFatal = 1;
122 simError();
123 }
124
125 #ifdef IS_MPI
126
127 }
128
129 #endif // is_mpi
130
131 }
132
133 DumpWriter::DumpWriter(SimInfo* info, const std::string& filename, bool writeDumpFile)
134 : info_(info), filename_(filename){
135
136 Globals* simParams = info->getSimParams();
137 eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
138
139 needCompression_ = simParams->getCompressDumpFile();
140 needForceVector_ = simParams->getOutputForceVector();
141
142 #ifdef HAVE_LIBZ
143 if (needCompression_) {
144 filename_ += ".gz";
145 eorFilename_ += ".gz";
146 }
147 #endif
148
149 #ifdef IS_MPI
150
151 if (worldRank == 0) {
152 #endif // is_mpi
153
154 createDumpFile_ = writeDumpFile;
155 if (createDumpFile_) {
156 dumpFile_ = createOStream(filename_);
157
158 if (!dumpFile_) {
159 sprintf(painCave.errMsg, "Could not open \"%s\" for dump output.\n",
160 filename_.c_str());
161 painCave.isFatal = 1;
162 simError();
163 }
164 }
165 #ifdef IS_MPI
166
167 }
168
169
170 #endif // is_mpi
171
172 }
173
174 DumpWriter::~DumpWriter() {
175
176 #ifdef IS_MPI
177
178 if (worldRank == 0) {
179 #endif // is_mpi
180 if (createDumpFile_){
181 writeClosing(*dumpFile_);
182 delete dumpFile_;
183 }
184 #ifdef IS_MPI
185
186 }
187
188 #endif // is_mpi
189
190 }
191
192 void DumpWriter::writeFrameProperties(std::ostream& os, Snapshot* s) {
193
194 char buffer[1024];
195
196 os << " <FrameData>\n";
197
198 RealType currentTime = s->getTime();
199
200 if (isinf(currentTime) || isnan(currentTime)) {
201 sprintf( painCave.errMsg,
202 "DumpWriter detected a numerical error writing the time");
203 painCave.isFatal = 1;
204 simError();
205 }
206
207 sprintf(buffer, " Time: %.10g\n", currentTime);
208 os << buffer;
209
210 Mat3x3d hmat;
211 hmat = s->getHmat();
212
213 for (unsigned int i = 0; i < 3; i++) {
214 for (unsigned int j = 0; j < 3; j++) {
215 if (isinf(hmat(i,j)) || isnan(hmat(i,j))) {
216 sprintf( painCave.errMsg,
217 "DumpWriter detected a numerical error writing the box");
218 painCave.isFatal = 1;
219 simError();
220 }
221 }
222 }
223
224 sprintf(buffer, " Hmat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
225 hmat(0, 0), hmat(1, 0), hmat(2, 0),
226 hmat(0, 1), hmat(1, 1), hmat(2, 1),
227 hmat(0, 2), hmat(1, 2), hmat(2, 2));
228 os << buffer;
229
230 RealType chi = s->getChi();
231 RealType integralOfChiDt = s->getIntegralOfChiDt();
232 if (isinf(chi) || isnan(chi) ||
233 isinf(integralOfChiDt) || isnan(integralOfChiDt)) {
234 sprintf( painCave.errMsg,
235 "DumpWriter detected a numerical error writing the thermostat");
236 painCave.isFatal = 1;
237 simError();
238 }
239 sprintf(buffer, " Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
240 os << buffer;
241
242 Mat3x3d eta;
243 eta = s->getEta();
244
245 for (unsigned int i = 0; i < 3; i++) {
246 for (unsigned int j = 0; j < 3; j++) {
247 if (isinf(eta(i,j)) || isnan(eta(i,j))) {
248 sprintf( painCave.errMsg,
249 "DumpWriter detected a numerical error writing the barostat");
250 painCave.isFatal = 1;
251 simError();
252 }
253 }
254 }
255
256 sprintf(buffer, " Barostat: {{ %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }, { %.10g, %.10g, %.10g }}\n",
257 eta(0, 0), eta(1, 0), eta(2, 0),
258 eta(0, 1), eta(1, 1), eta(2, 1),
259 eta(0, 2), eta(1, 2), eta(2, 2));
260 os << buffer;
261
262 os << " </FrameData>\n";
263 }
264
265 void DumpWriter::writeFrame(std::ostream& os) {
266
267 #ifdef IS_MPI
268 MPI_Status istatus;
269 #endif
270
271 Molecule* mol;
272 StuntDouble* integrableObject;
273 SimInfo::MoleculeIterator mi;
274 Molecule::IntegrableObjectIterator ii;
275
276 #ifndef IS_MPI
277 os << " <Snapshot>\n";
278
279 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
280
281 os << " <StuntDoubles>\n";
282 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
283
284
285 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
286 integrableObject = mol->nextIntegrableObject(ii)) {
287 os << prepareDumpLine(integrableObject);
288
289 }
290 }
291 os << " </StuntDoubles>\n";
292
293 os << " </Snapshot>\n";
294
295 os.flush();
296 #else
297 //every node prepares the dump lines for integrable objects belong to itself
298 std::string buffer;
299 for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
300
301
302 for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
303 integrableObject = mol->nextIntegrableObject(ii)) {
304 buffer += prepareDumpLine(integrableObject);
305 }
306 }
307
308 const int masterNode = 0;
309
310 if (worldRank == masterNode) {
311 os << " <Snapshot>\n";
312 writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
313 os << " <StuntDoubles>\n";
314
315 os << buffer;
316
317 int nProc;
318 MPI_Comm_size(MPI_COMM_WORLD, &nProc);
319 for (int i = 1; i < nProc; ++i) {
320
321 // receive the length of the string buffer that was
322 // prepared by processor i
323
324 int recvLength;
325 MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
326 char* recvBuffer = new char[recvLength];
327 if (recvBuffer == NULL) {
328 } else {
329 MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
330 os << recvBuffer;
331 delete [] recvBuffer;
332 }
333 }
334 os << " </StuntDoubles>\n";
335
336 os << " </Snapshot>\n";
337 os.flush();
338 } else {
339 int sendBufferLength = buffer.size() + 1;
340 MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
341 MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
342 }
343
344 #endif // is_mpi
345
346 }
347
348 std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
349
350 int index = integrableObject->getGlobalIntegrableObjectIndex();
351 std::string type("pv");
352 std::string line;
353 char tempBuffer[4096];
354
355 Vector3d pos;
356 Vector3d vel;
357 pos = integrableObject->getPos();
358
359 if (isinf(pos[0]) || isnan(pos[0]) ||
360 isinf(pos[1]) || isnan(pos[1]) ||
361 isinf(pos[2]) || isnan(pos[2]) ) {
362 sprintf( painCave.errMsg,
363 "DumpWriter detected a numerical error writing the position"
364 " for object %d", index);
365 painCave.isFatal = 1;
366 simError();
367 }
368
369 vel = integrableObject->getVel();
370
371 if (isinf(vel[0]) || isnan(vel[0]) ||
372 isinf(vel[1]) || isnan(vel[1]) ||
373 isinf(vel[2]) || isnan(vel[2]) ) {
374 sprintf( painCave.errMsg,
375 "DumpWriter detected a numerical error writing the velocity"
376 " for object %d", index);
377 painCave.isFatal = 1;
378 simError();
379 }
380
381 sprintf(tempBuffer, "%18.10g %18.10g %18.10g %13e %13e %13e",
382 pos[0], pos[1], pos[2],
383 vel[0], vel[1], vel[2]);
384 line += tempBuffer;
385
386 if (integrableObject->isDirectional()) {
387 type += "qj";
388 Quat4d q;
389 Vector3d ji;
390 q = integrableObject->getQ();
391
392 if (isinf(q[0]) || isnan(q[0]) ||
393 isinf(q[1]) || isnan(q[1]) ||
394 isinf(q[2]) || isnan(q[2]) ||
395 isinf(q[3]) || isnan(q[3]) ) {
396 sprintf( painCave.errMsg,
397 "DumpWriter detected a numerical error writing the quaternion"
398 " for object %d", index);
399 painCave.isFatal = 1;
400 simError();
401 }
402
403 ji = integrableObject->getJ();
404
405 if (isinf(ji[0]) || isnan(ji[0]) ||
406 isinf(ji[1]) || isnan(ji[1]) ||
407 isinf(ji[2]) || isnan(ji[2]) ) {
408 sprintf( painCave.errMsg,
409 "DumpWriter detected a numerical error writing the angular"
410 " momentum for object %d", index);
411 painCave.isFatal = 1;
412 simError();
413 }
414
415 sprintf(tempBuffer, " %13e %13e %13e %13e %13e %13e %13e",
416 q[0], q[1], q[2], q[3],
417 ji[0], ji[1], ji[2]);
418 line += tempBuffer;
419 }
420
421 if (needForceVector_) {
422 type += "f";
423 Vector3d frc;
424
425 frc = integrableObject->getFrc();
426
427 if (isinf(frc[0]) || isnan(frc[0]) ||
428 isinf(frc[1]) || isnan(frc[1]) ||
429 isinf(frc[2]) || isnan(frc[2]) ) {
430 sprintf( painCave.errMsg,
431 "DumpWriter detected a numerical error writing the force"
432 " for object %d", index);
433 painCave.isFatal = 1;
434 simError();
435 }
436 sprintf(tempBuffer, " %13e %13e %13e",
437 frc[0], frc[1], frc[2]);
438 line += tempBuffer;
439
440 if (integrableObject->isDirectional()) {
441 type += "t";
442 Vector3d trq;
443
444 trq = integrableObject->getTrq();
445
446 if (isinf(trq[0]) || isnan(trq[0]) ||
447 isinf(trq[1]) || isnan(trq[1]) ||
448 isinf(trq[2]) || isnan(trq[2]) ) {
449 sprintf( painCave.errMsg,
450 "DumpWriter detected a numerical error writing the torque"
451 " for object %d", index);
452 painCave.isFatal = 1;
453 simError();
454 }
455
456 sprintf(tempBuffer, " %13e %13e %13e",
457 trq[0], trq[1], trq[2]);
458 line += tempBuffer;
459 }
460 }
461
462 sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
463 return std::string(tempBuffer);
464 }
465
466 void DumpWriter::writeDump() {
467 writeFrame(*dumpFile_);
468 }
469
470 void DumpWriter::writeEor() {
471 std::ostream* eorStream;
472
473 #ifdef IS_MPI
474 if (worldRank == 0) {
475 #endif // is_mpi
476
477 eorStream = createOStream(eorFilename_);
478
479 #ifdef IS_MPI
480 }
481 #endif // is_mpi
482
483 writeFrame(*eorStream);
484
485 #ifdef IS_MPI
486 if (worldRank == 0) {
487 #endif // is_mpi
488 writeClosing(*eorStream);
489 delete eorStream;
490 #ifdef IS_MPI
491 }
492 #endif // is_mpi
493
494 }
495
496
497 void DumpWriter::writeDumpAndEor() {
498 std::vector<std::streambuf*> buffers;
499 std::ostream* eorStream;
500 #ifdef IS_MPI
501 if (worldRank == 0) {
502 #endif // is_mpi
503
504 buffers.push_back(dumpFile_->rdbuf());
505
506 eorStream = createOStream(eorFilename_);
507
508 buffers.push_back(eorStream->rdbuf());
509
510 #ifdef IS_MPI
511 }
512 #endif // is_mpi
513
514 TeeBuf tbuf(buffers.begin(), buffers.end());
515 std::ostream os(&tbuf);
516
517 writeFrame(os);
518
519 #ifdef IS_MPI
520 if (worldRank == 0) {
521 #endif // is_mpi
522 writeClosing(*eorStream);
523 delete eorStream;
524 #ifdef IS_MPI
525 }
526 #endif // is_mpi
527
528 }
529
530 std::ostream* DumpWriter::createOStream(const std::string& filename) {
531
532 std::ostream* newOStream;
533 #ifdef HAVE_LIBZ
534 if (needCompression_) {
535 newOStream = new ogzstream(filename.c_str());
536 } else {
537 newOStream = new std::ofstream(filename.c_str());
538 }
539 #else
540 newOStream = new std::ofstream(filename.c_str());
541 #endif
542 //write out MetaData first
543 (*newOStream) << "<OpenMD version=1>" << std::endl;
544 (*newOStream) << " <MetaData>" << std::endl;
545 (*newOStream) << info_->getRawMetaData();
546 (*newOStream) << " </MetaData>" << std::endl;
547 return newOStream;
548 }
549
550 void DumpWriter::writeClosing(std::ostream& os) {
551
552 os << "</OpenMD>\n";
553 os.flush();
554 }
555
556 }//end namespace OpenMD