ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/io/DumpWriter.cpp
(Generate patch)

Comparing branches/development/src/io/DumpWriter.cpp (file contents):
Revision 1752 by gezelter, Sun Jun 10 14:05:02 2012 UTC vs.
Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 UTC

# Line 35 | Line 35
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 44 | Line 44
44   #include "primitives/Molecule.hpp"
45   #include "utils/simError.h"
46   #include "io/basic_teebuf.hpp"
47 + #ifdef HAVE_ZLIB
48   #include "io/gzstream.hpp"
49 + #endif
50   #include "io/Globals.hpp"
51  
52 + #ifdef _MSC_VER
53 + #define isnan(x) _isnan((x))
54 + #define isinf(x) (!_finite(x) && !_isnan(x))
55 + #endif
56  
57   #ifdef IS_MPI
58   #include <mpi.h>
59 < #endif //is_mpi
59 > #endif
60  
61   using namespace std;
62   namespace OpenMD {
# Line 257 | Line 263 | namespace OpenMD {
263              hmat(0, 2), hmat(1, 2), hmat(2, 2));
264      os << buffer;
265  
266 <    RealType chi = s->getChi();
267 <    RealType integralOfChiDt = s->getIntegralOfChiDt();
268 <    if (isinf(chi) || isnan(chi) ||
269 <        isinf(integralOfChiDt) || isnan(integralOfChiDt)) {      
266 >    pair<RealType, RealType> thermostat = s->getThermostat();
267 >
268 >    if (isinf(thermostat.first)  || isnan(thermostat.first) ||
269 >        isinf(thermostat.second) || isnan(thermostat.second)) {      
270        sprintf( painCave.errMsg,
271                 "DumpWriter detected a numerical error writing the thermostat");
272        painCave.isFatal = 1;
273        simError();
274      }
275 <    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", chi, integralOfChiDt);
275 >    sprintf(buffer, "  Thermostat: %.10g , %.10g\n", thermostat.first,
276 >            thermostat.second);
277      os << buffer;
278  
279      Mat3x3d eta;
280 <    eta = s->getEta();
280 >    eta = s->getBarostat();
281  
282      for (unsigned int i = 0; i < 3; i++) {
283        for (unsigned int j = 0; j < 3; j++) {
# Line 295 | Line 302 | namespace OpenMD {
302    void DumpWriter::writeFrame(std::ostream& os) {
303  
304   #ifdef IS_MPI
305 <    MPI_Status istatus;
305 >    MPI::Status istatus;
306   #endif
307  
308      Molecule* mol;
309 <    StuntDouble* integrableObject;
309 >    StuntDouble* sd;
310      SimInfo::MoleculeIterator mi;
311      Molecule::IntegrableObjectIterator ii;
312      RigidBody::AtomIterator ai;
# Line 314 | Line 321 | namespace OpenMD {
321      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
322  
323        
324 <      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
325 <           integrableObject = mol->nextIntegrableObject(ii)) {  
326 <          os << prepareDumpLine(integrableObject);
324 >      for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
325 >           sd = mol->nextIntegrableObject(ii)) {        
326 >          os << prepareDumpLine(sd);
327            
328        }
329      }    
# Line 324 | Line 331 | namespace OpenMD {
331  
332      if (doSiteData_) {
333        os << "    <SiteData>\n";
334 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
335 <              
336 <        for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  
337 <           integrableObject = mol->nextIntegrableObject(ii)) {  
334 >      for (mol = info_->beginMolecule(mi); mol != NULL;
335 >           mol = info_->nextMolecule(mi)) {
336 >              
337 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
338 >           sd = mol->nextIntegrableObject(ii)) {        
339  
340 <          int ioIndex = integrableObject->getGlobalIntegrableObjectIndex();
340 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
341            // do one for the IO itself
342 <          os << prepareSiteLine(integrableObject, ioIndex, 0);
342 >          os << prepareSiteLine(sd, ioIndex, 0);
343  
344 <          if (integrableObject->isRigidBody()) {
344 >          if (sd->isRigidBody()) {
345              
346 <            RigidBody* rb = static_cast<RigidBody*>(integrableObject);
346 >            RigidBody* rb = static_cast<RigidBody*>(sd);
347              int siteIndex = 0;
348              for (atom = rb->beginAtom(ai); atom != NULL;  
349                   atom = rb->nextAtom(ai)) {                                            
# Line 351 | Line 359 | namespace OpenMD {
359  
360      os.flush();
361   #else
354    //every node prepares the dump lines for integrable objects belong to itself
355    std::string buffer;
356    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
362  
358
359      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
360           integrableObject = mol->nextIntegrableObject(ii)) {  
361          buffer += prepareDumpLine(integrableObject);
362      }
363    }
364    
363      const int masterNode = 0;
364 <    int nProc;
365 <    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
364 >    int worldRank = MPI::COMM_WORLD.Get_rank();
365 >    int nProc = MPI::COMM_WORLD.Get_size();
366 >
367      if (worldRank == masterNode) {      
368        os << "  <Snapshot>\n";  
369 <      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
369 >      writeFrameProperties(os,
370 >                           info_->getSnapshotManager()->getCurrentSnapshot());
371        os << "    <StuntDoubles>\n";
372 <        
373 <      os << buffer;
372 >    }
373  
374 +    //every node prepares the dump lines for integrable objects belong to itself
375 +    std::string buffer;
376 +    for (mol = info_->beginMolecule(mi); mol != NULL;
377 +         mol = info_->nextMolecule(mi)) {
378 +      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
379 +           sd = mol->nextIntegrableObject(ii)) {        
380 +        buffer += prepareDumpLine(sd);
381 +      }
382 +    }
383 +    
384 +    if (worldRank == masterNode) {      
385 +      os << buffer;
386 +      
387        for (int i = 1; i < nProc; ++i) {
388 +        // tell processor i to start sending us data:
389 +        MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
390  
391          // receive the length of the string buffer that was
392 <        // prepared by processor i
379 <
380 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
392 >        // prepared by processor i:        
393          int recvLength;
394 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
394 >        MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
395 >                             istatus);
396 >
397 >        // create a buffer to receive the data
398          char* recvBuffer = new char[recvLength];
399          if (recvBuffer == NULL) {
400          } else {
401 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
401 >          // receive the data:
402 >          MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
403 >                               MPI::ANY_TAG, istatus);
404 >          // send it to the file:
405            os << recvBuffer;
406 +          // get rid of the receive buffer:
407            delete [] recvBuffer;
408          }
409        }
391      os << "    </StuntDoubles>\n";
392      
393      os << "  </Snapshot>\n";
394      os.flush();
410      } else {
411        int sendBufferLength = buffer.size() + 1;
412        int myturn = 0;
413        for (int i = 1; i < nProc; ++i){
414 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
414 >        // wait for the master node to call our number:
415 >        MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
416          if (myturn == worldRank){
417 <          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
418 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
417 >          // send the length of our buffer:
418 >          MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
419 >
420 >          // send our buffer:
421 >          MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
422 >                               MPI::CHAR, masterNode, 0);
423          }
424        }
425      }
426 +    
427 +    if (worldRank == masterNode) {      
428 +      os << "    </StuntDoubles>\n";
429 +    }
430  
431 < #endif // is_mpi
431 >    if (doSiteData_) {
432 >      if (worldRank == masterNode) {
433 >        os << "    <SiteData>\n";
434 >      }
435 >      buffer.clear();
436 >      for (mol = info_->beginMolecule(mi); mol != NULL;
437 >           mol = info_->nextMolecule(mi)) {
438 >              
439 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
440 >             sd = mol->nextIntegrableObject(ii)) {      
441 >          
442 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
443 >          // do one for the IO itself
444 >          buffer += prepareSiteLine(sd, ioIndex, 0);
445  
446 <  }
446 >          if (sd->isRigidBody()) {
447 >            
448 >            RigidBody* rb = static_cast<RigidBody*>(sd);
449 >            int siteIndex = 0;
450 >            for (atom = rb->beginAtom(ai); atom != NULL;  
451 >                 atom = rb->nextAtom(ai)) {                                            
452 >              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
453 >              siteIndex++;
454 >            }
455 >          }
456 >        }
457 >      }
458  
459 <  std::string DumpWriter::prepareDumpLine(StuntDouble* integrableObject) {
459 >      if (worldRank == masterNode) {    
460 >        os << buffer;
461 >        
462 >        for (int i = 1; i < nProc; ++i) {
463 >          
464 >          // tell processor i to start sending us data:
465 >          MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
466 >          
467 >          // receive the length of the string buffer that was
468 >          // prepared by processor i:        
469 >          int recvLength;
470 >          MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
471 >                               istatus);
472 >          
473 >          // create a buffer to receive the data
474 >          char* recvBuffer = new char[recvLength];
475 >          if (recvBuffer == NULL) {
476 >          } else {
477 >            // receive the data:
478 >            MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
479 >                                 MPI::ANY_TAG, istatus);
480 >            // send it to the file:
481 >            os << recvBuffer;
482 >            // get rid of the receive buffer:
483 >            delete [] recvBuffer;
484 >          }
485 >        }      
486 >      } else {
487 >        int sendBufferLength = buffer.size() + 1;
488 >        int myturn = 0;
489 >        for (int i = 1; i < nProc; ++i){
490 >          // wait for the master node to call our number:
491 >          MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
492 >          if (myturn == worldRank){
493 >            // send the length of our buffer:
494 >            MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
495 >            // send our buffer:
496 >            MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
497 >                                 MPI::CHAR, masterNode, 0);
498 >          }
499 >        }
500 >      }
501 >      
502 >      if (worldRank == masterNode) {    
503 >        os << "    </SiteData>\n";
504 >      }
505 >    }
506 >    
507 >    if (worldRank == masterNode) {
508 >      os << "  </Snapshot>\n";
509 >      os.flush();
510 >    }
511 >    
512 > #endif // is_mpi
513 >    
514 >  }
515 >
516 >  std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
517          
518 <    int index = integrableObject->getGlobalIntegrableObjectIndex();
518 >    int index = sd->getGlobalIntegrableObjectIndex();
519      std::string type("pv");
520      std::string line;
521      char tempBuffer[4096];
522  
523      Vector3d pos;
524      Vector3d vel;
525 <    pos = integrableObject->getPos();
525 >    pos = sd->getPos();
526  
527      if (isinf(pos[0]) || isnan(pos[0]) ||
528          isinf(pos[1]) || isnan(pos[1]) ||
# Line 429 | Line 534 | namespace OpenMD {
534        simError();
535      }
536  
537 <    vel = integrableObject->getVel();          
537 >    vel = sd->getVel();        
538  
539      if (isinf(vel[0]) || isnan(vel[0]) ||
540          isinf(vel[1]) || isnan(vel[1]) ||
# Line 446 | Line 551 | namespace OpenMD {
551              vel[0], vel[1], vel[2]);                    
552      line += tempBuffer;
553  
554 <    if (integrableObject->isDirectional()) {
554 >    if (sd->isDirectional()) {
555        type += "qj";
556        Quat4d q;
557        Vector3d ji;
558 <      q = integrableObject->getQ();
558 >      q = sd->getQ();
559  
560        if (isinf(q[0]) || isnan(q[0]) ||
561            isinf(q[1]) || isnan(q[1]) ||
# Line 463 | Line 568 | namespace OpenMD {
568          simError();
569        }
570  
571 <      ji = integrableObject->getJ();
571 >      ji = sd->getJ();
572  
573        if (isinf(ji[0]) || isnan(ji[0]) ||
574            isinf(ji[1]) || isnan(ji[1]) ||
# Line 483 | Line 588 | namespace OpenMD {
588  
589      if (needForceVector_) {
590        type += "f";
591 <      Vector3d frc = integrableObject->getFrc();
591 >      Vector3d frc = sd->getFrc();
592        if (isinf(frc[0]) || isnan(frc[0]) ||
593            isinf(frc[1]) || isnan(frc[1]) ||
594            isinf(frc[2]) || isnan(frc[2]) ) {      
# Line 497 | Line 602 | namespace OpenMD {
602                frc[0], frc[1], frc[2]);
603        line += tempBuffer;
604        
605 <      if (integrableObject->isDirectional()) {
605 >      if (sd->isDirectional()) {
606          type += "t";
607 <        Vector3d trq = integrableObject->getTrq();        
607 >        Vector3d trq = sd->getTrq();        
608          if (isinf(trq[0]) || isnan(trq[0]) ||
609              isinf(trq[1]) || isnan(trq[1]) ||
610              isinf(trq[2]) || isnan(trq[2]) ) {      
# Line 519 | Line 624 | namespace OpenMD {
624      return std::string(tempBuffer);
625    }
626  
627 <  std::string DumpWriter::prepareSiteLine(StuntDouble* integrableObject, int ioIndex, int siteIndex) {
627 >  std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
628          
629  
630      std::string id;
# Line 527 | Line 632 | namespace OpenMD {
632      std::string line;
633      char tempBuffer[4096];
634  
635 <    if (integrableObject->isRigidBody()) {
635 >    if (sd->isRigidBody()) {
636        sprintf(tempBuffer, "%10d           ", ioIndex);
637        id = std::string(tempBuffer);
638      } else {
# Line 537 | Line 642 | namespace OpenMD {
642                
643      if (needFlucQ_) {
644        type += "cw";
645 <      RealType fqPos = integrableObject->getFlucQPos();
645 >      RealType fqPos = sd->getFlucQPos();
646        if (isinf(fqPos) || isnan(fqPos) ) {      
647          sprintf( painCave.errMsg,
648                   "DumpWriter detected a numerical error writing the"
# Line 548 | Line 653 | namespace OpenMD {
653        sprintf(tempBuffer, " %13e ", fqPos);
654        line += tempBuffer;
655      
656 <      RealType fqVel = integrableObject->getFlucQVel();
656 >      RealType fqVel = sd->getFlucQVel();
657        if (isinf(fqVel) || isnan(fqVel) ) {      
658          sprintf( painCave.errMsg,
659                   "DumpWriter detected a numerical error writing the"
# Line 561 | Line 666 | namespace OpenMD {
666  
667        if (needForceVector_) {
668          type += "g";
669 <        RealType fqFrc = integrableObject->getFlucQFrc();        
669 >        RealType fqFrc = sd->getFlucQFrc();        
670          if (isinf(fqFrc) || isnan(fqFrc) ) {      
671            sprintf( painCave.errMsg,
672                     "DumpWriter detected a numerical error writing the"
# Line 576 | Line 681 | namespace OpenMD {
681  
682      if (needElectricField_) {
683        type += "e";
684 <      Vector3d eField= integrableObject->getElectricField();
684 >      Vector3d eField= sd->getElectricField();
685        if (isinf(eField[0]) || isnan(eField[0]) ||
686            isinf(eField[1]) || isnan(eField[1]) ||
687            isinf(eField[2]) || isnan(eField[2]) ) {      
# Line 594 | Line 699 | namespace OpenMD {
699  
700      if (needParticlePot_) {
701        type += "u";
702 <      RealType particlePot = integrableObject->getParticlePot();
702 >      RealType particlePot = sd->getParticlePot();
703        if (isinf(particlePot) || isnan(particlePot)) {      
704          sprintf( painCave.errMsg,
705                   "DumpWriter detected a numerical error writing the particle "
# Line 678 | Line 783 | namespace OpenMD {
783    std::ostream* DumpWriter::createOStream(const std::string& filename) {
784  
785      std::ostream* newOStream;
786 < #ifdef HAVE_LIBZ
786 > #ifdef HAVE_ZLIB
787      if (needCompression_) {
788        newOStream = new ogzstream(filename.c_str());
789      } else {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines