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

Comparing trunk/src/io/DumpWriter.cpp (file contents):
Revision 1793 by gezelter, Fri Aug 31 21:16:10 2012 UTC vs.
Revision 1983 by gezelter, Tue Apr 15 20:36:19 2014 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   */
42 +
43 + #include "config.h"
44 +
45 + #ifdef IS_MPI
46 + #include <mpi.h>
47 + #endif
48  
49   #include "io/DumpWriter.hpp"
50   #include "primitives/Molecule.hpp"
# Line 52 | Line 58
58   #ifdef _MSC_VER
59   #define isnan(x) _isnan((x))
60   #define isinf(x) (!_finite(x) && !_isnan(x))
55 – #endif
56 –
57 – #ifdef IS_MPI
58 – #include <mpi.h>
61   #endif
62  
63   using namespace std;
# Line 310 | Line 312 | namespace OpenMD {
312      SimInfo::MoleculeIterator mi;
313      Molecule::IntegrableObjectIterator ii;
314      RigidBody::AtomIterator ai;
313 –    Atom* atom;
315  
316   #ifndef IS_MPI
317      os << "  <Snapshot>\n";
# Line 318 | Line 319 | namespace OpenMD {
319      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
320  
321      os << "    <StuntDoubles>\n";
322 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
323 <
322 >    for (mol = info_->beginMolecule(mi); mol != NULL;
323 >         mol = info_->nextMolecule(mi)) {
324        
325        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
326             sd = mol->nextIntegrableObject(ii)) {        
# Line 331 | Line 332 | namespace OpenMD {
332  
333      if (doSiteData_) {
334        os << "    <SiteData>\n";
335 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
335 >      for (mol = info_->beginMolecule(mi); mol != NULL;
336 >           mol = info_->nextMolecule(mi)) {
337                
338          for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
339 <             sd = mol->nextIntegrableObject(ii)) {      
340 <          
339 >           sd = mol->nextIntegrableObject(ii)) {        
340 >
341            int ioIndex = sd->getGlobalIntegrableObjectIndex();
342            // do one for the IO itself
343            os << prepareSiteLine(sd, ioIndex, 0);
# Line 344 | Line 346 | namespace OpenMD {
346              
347              RigidBody* rb = static_cast<RigidBody*>(sd);
348              int siteIndex = 0;
349 <            for (atom = rb->beginAtom(ai); atom != NULL;  
349 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
350                   atom = rb->nextAtom(ai)) {                                            
351                os << prepareSiteLine(atom, ioIndex, siteIndex);
352                siteIndex++;
# Line 358 | Line 360 | namespace OpenMD {
360  
361      os.flush();
362   #else
361 –    //every node prepares the dump lines for integrable objects belong to itself
362 –    std::string buffer;
363 –    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
363  
365 –
366 –      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
367 –           sd = mol->nextIntegrableObject(ii)) {        
368 –          buffer += prepareDumpLine(sd);
369 –      }
370 –    }
371 –    
364      const int masterNode = 0;
365 +    int worldRank;
366      int nProc;
367 <    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
367 >
368 >    MPI_Comm_size( MPI_COMM_WORLD, &nProc);
369 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
370 >
371 >
372      if (worldRank == masterNode) {      
373        os << "  <Snapshot>\n";  
374 <      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
374 >      writeFrameProperties(os,
375 >                           info_->getSnapshotManager()->getCurrentSnapshot());
376        os << "    <StuntDoubles>\n";
377 <        
380 <      os << buffer;
377 >    }
378  
379 +    //every node prepares the dump lines for integrable objects belong to itself
380 +    std::string buffer;
381 +    for (mol = info_->beginMolecule(mi); mol != NULL;
382 +         mol = info_->nextMolecule(mi)) {
383 +      for (sd = mol->beginIntegrableObject(ii); sd != NULL;
384 +           sd = mol->nextIntegrableObject(ii)) {        
385 +        buffer += prepareDumpLine(sd);
386 +      }
387 +    }
388 +    
389 +    if (worldRank == masterNode) {      
390 +      os << buffer;
391 +      
392        for (int i = 1; i < nProc; ++i) {
393 +        // tell processor i to start sending us data:
394 +        MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
395  
396          // receive the length of the string buffer that was
397 <        // prepared by processor i
386 <
387 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
397 >        // prepared by processor i:        
398          int recvLength;
399 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
399 >        MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
400 >                 &istatus);
401 >
402 >        // create a buffer to receive the data
403          char* recvBuffer = new char[recvLength];
404          if (recvBuffer == NULL) {
405          } else {
406 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
406 >          // receive the data:
407 >          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
408 >                               MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
409 >          // send it to the file:
410            os << recvBuffer;
411 +          // get rid of the receive buffer:
412            delete [] recvBuffer;
413          }
414        }
398 –      os << "    </StuntDoubles>\n";
399 –      
400 –      os << "  </Snapshot>\n";
401 –      os.flush();
415      } else {
416        int sendBufferLength = buffer.size() + 1;
417        int myturn = 0;
418        for (int i = 1; i < nProc; ++i){
419 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
419 >        // wait for the master node to call our number:
420 >        MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
421          if (myturn == worldRank){
422 +          // send the length of our buffer:
423            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
424 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
424 >
425 >          // send our buffer:
426 >          MPI_Send((void *)buffer.c_str(), sendBufferLength,
427 >                   MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
428 >
429          }
430        }
431      }
432 +    
433 +    if (worldRank == masterNode) {      
434 +      os << "    </StuntDoubles>\n";
435 +    }
436  
437 < #endif // is_mpi
437 >    if (doSiteData_) {
438 >      if (worldRank == masterNode) {
439 >        os << "    <SiteData>\n";
440 >      }
441 >      buffer.clear();
442 >      for (mol = info_->beginMolecule(mi); mol != NULL;
443 >           mol = info_->nextMolecule(mi)) {
444 >              
445 >        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
446 >             sd = mol->nextIntegrableObject(ii)) {      
447 >          
448 >          int ioIndex = sd->getGlobalIntegrableObjectIndex();
449 >          // do one for the IO itself
450 >          buffer += prepareSiteLine(sd, ioIndex, 0);
451  
452 +          if (sd->isRigidBody()) {
453 +            
454 +            RigidBody* rb = static_cast<RigidBody*>(sd);
455 +            int siteIndex = 0;
456 +            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
457 +                 atom = rb->nextAtom(ai)) {                                            
458 +              buffer += prepareSiteLine(atom, ioIndex, siteIndex);
459 +              siteIndex++;
460 +            }
461 +          }
462 +        }
463 +      }
464 +
465 +      if (worldRank == masterNode) {    
466 +        os << buffer;
467 +        
468 +        for (int i = 1; i < nProc; ++i) {
469 +          
470 +          // tell processor i to start sending us data:
471 +          MPI_Bcast(&i, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
472 +          
473 +          // receive the length of the string buffer that was
474 +          // prepared by processor i:        
475 +          int recvLength;
476 +          MPI_Recv(&recvLength, 1, MPI_INT, i, MPI_ANY_TAG, MPI_COMM_WORLD,
477 +                   &istatus);
478 +          
479 +          // create a buffer to receive the data
480 +          char* recvBuffer = new char[recvLength];
481 +          if (recvBuffer == NULL) {
482 +          } else {
483 +            // receive the data:
484 +            MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i,
485 +                     MPI_ANY_TAG, MPI_COMM_WORLD, &istatus);
486 +            // send it to the file:
487 +            os << recvBuffer;
488 +            // get rid of the receive buffer:
489 +            delete [] recvBuffer;
490 +          }
491 +        }      
492 +      } else {
493 +        int sendBufferLength = buffer.size() + 1;
494 +        int myturn = 0;
495 +        for (int i = 1; i < nProc; ++i){
496 +          // wait for the master node to call our number:
497 +          MPI_Bcast(&myturn, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
498 +          if (myturn == worldRank){
499 +            // send the length of our buffer:
500 +            MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
501 +            // send our buffer:
502 +            MPI_Send((void *)buffer.c_str(), sendBufferLength,
503 +                     MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
504 +          }
505 +        }
506 +      }
507 +      
508 +      if (worldRank == masterNode) {    
509 +        os << "    </SiteData>\n";
510 +      }
511 +    }
512 +    
513 +    if (worldRank == masterNode) {
514 +      os << "  </Snapshot>\n";
515 +      os.flush();
516 +    }
517 +    
518 + #endif // is_mpi
519 +    
520    }
521  
522    std::string DumpWriter::prepareDumpLine(StuntDouble* sd) {
# Line 527 | Line 631 | namespace OpenMD {
631    }
632  
633    std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
634 <        
634 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
635  
636      std::string id;
637      std::string type;
# Line 543 | Line 647 | namespace OpenMD {
647      }
648                
649      if (needFlucQ_) {
650 <      type += "cw";
651 <      RealType fqPos = sd->getFlucQPos();
652 <      if (isinf(fqPos) || isnan(fqPos) ) {      
653 <        sprintf( painCave.errMsg,
654 <                 "DumpWriter detected a numerical error writing the"
655 <                 " fluctuating charge for object %s", id.c_str());      
656 <        painCave.isFatal = 1;
657 <        simError();
658 <      }
659 <      sprintf(tempBuffer, " %13e ", fqPos);
660 <      line += tempBuffer;
661 <    
662 <      RealType fqVel = sd->getFlucQVel();
559 <      if (isinf(fqVel) || isnan(fqVel) ) {      
560 <        sprintf( painCave.errMsg,
561 <                 "DumpWriter detected a numerical error writing the"
562 <                 " fluctuating charge velocity for object %s", id.c_str());      
563 <        painCave.isFatal = 1;
564 <        simError();
565 <      }
566 <      sprintf(tempBuffer, " %13e ", fqVel);
567 <      line += tempBuffer;
650 >      if (storageLayout & DataStorage::dslFlucQPosition) {
651 >        type += "c";
652 >        RealType fqPos = sd->getFlucQPos();
653 >        if (isinf(fqPos) || isnan(fqPos) ) {      
654 >          sprintf( painCave.errMsg,
655 >                   "DumpWriter detected a numerical error writing the"
656 >                   " fluctuating charge for object %s", id.c_str());      
657 >          painCave.isFatal = 1;
658 >          simError();
659 >        }
660 >        sprintf(tempBuffer, " %13e ", fqPos);
661 >        line += tempBuffer;
662 >      }
663  
664 <      if (needForceVector_) {
665 <        type += "g";
666 <        RealType fqFrc = sd->getFlucQFrc();        
667 <        if (isinf(fqFrc) || isnan(fqFrc) ) {      
664 >      if (storageLayout & DataStorage::dslFlucQVelocity) {
665 >        type += "w";    
666 >        RealType fqVel = sd->getFlucQVel();
667 >        if (isinf(fqVel) || isnan(fqVel) ) {      
668            sprintf( painCave.errMsg,
669                     "DumpWriter detected a numerical error writing the"
670 <                   " fluctuating charge force for object %s", id.c_str());      
670 >                   " fluctuating charge velocity for object %s", id.c_str());      
671            painCave.isFatal = 1;
672            simError();
673          }
674 <        sprintf(tempBuffer, " %13e ", fqFrc);        
674 >        sprintf(tempBuffer, " %13e ", fqVel);
675          line += tempBuffer;
676        }
582 –    }
677  
678 +      if (needForceVector_) {
679 +        if (storageLayout & DataStorage::dslFlucQForce) {          
680 +          type += "g";
681 +          RealType fqFrc = sd->getFlucQFrc();        
682 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
683 +            sprintf( painCave.errMsg,
684 +                     "DumpWriter detected a numerical error writing the"
685 +                     " fluctuating charge force for object %s", id.c_str());      
686 +            painCave.isFatal = 1;
687 +            simError();
688 +          }
689 +          sprintf(tempBuffer, " %13e ", fqFrc);        
690 +          line += tempBuffer;
691 +        }
692 +      }
693 +    }
694 +    
695      if (needElectricField_) {
696 <      type += "e";
697 <      Vector3d eField= sd->getElectricField();
698 <      if (isinf(eField[0]) || isnan(eField[0]) ||
699 <          isinf(eField[1]) || isnan(eField[1]) ||
700 <          isinf(eField[2]) || isnan(eField[2]) ) {      
701 <        sprintf( painCave.errMsg,
702 <                 "DumpWriter detected a numerical error writing the electric"
703 <                 " field for object %s", id.c_str());      
704 <        painCave.isFatal = 1;
705 <        simError();
696 >      if (storageLayout & DataStorage::dslElectricField) {
697 >        type += "e";
698 >        Vector3d eField= sd->getElectricField();
699 >        if (isinf(eField[0]) || isnan(eField[0]) ||
700 >            isinf(eField[1]) || isnan(eField[1]) ||
701 >            isinf(eField[2]) || isnan(eField[2]) ) {      
702 >          sprintf( painCave.errMsg,
703 >                   "DumpWriter detected a numerical error writing the electric"
704 >                   " field for object %s", id.c_str());      
705 >          painCave.isFatal = 1;
706 >          simError();
707 >        }
708 >        sprintf(tempBuffer, " %13e %13e %13e",
709 >                eField[0], eField[1], eField[2]);
710 >        line += tempBuffer;
711        }
596 –      sprintf(tempBuffer, " %13e %13e %13e",
597 –              eField[0], eField[1], eField[2]);
598 –      line += tempBuffer;
712      }
713  
714  
715      if (needParticlePot_) {
716 <      type += "u";
717 <      RealType particlePot = sd->getParticlePot();
718 <      if (isinf(particlePot) || isnan(particlePot)) {      
719 <        sprintf( painCave.errMsg,
720 <                 "DumpWriter detected a numerical error writing the particle "
721 <                 " potential for object %s", id.c_str());      
722 <        painCave.isFatal = 1;
723 <        simError();
716 >      if (storageLayout & DataStorage::dslParticlePot) {
717 >        type += "u";
718 >        RealType particlePot = sd->getParticlePot();
719 >        if (isinf(particlePot) || isnan(particlePot)) {      
720 >          sprintf( painCave.errMsg,
721 >                   "DumpWriter detected a numerical error writing the particle "
722 >                   " potential for object %s", id.c_str());      
723 >          painCave.isFatal = 1;
724 >          simError();
725 >        }
726 >        sprintf(tempBuffer, " %13e", particlePot);
727 >        line += tempBuffer;
728        }
612 –      sprintf(tempBuffer, " %13e", particlePot);
613 –      line += tempBuffer;
729      }
730 <    
616 <
730 >  
731      sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
732      return std::string(tempBuffer);
733    }
# Line 623 | Line 737 | namespace OpenMD {
737    }
738  
739    void DumpWriter::writeEor() {
740 <    std::ostream* eorStream;
741 <    
740 >
741 >    std::ostream* eorStream = NULL;
742 >
743   #ifdef IS_MPI
744      if (worldRank == 0) {
745   #endif // is_mpi
746 <
746 >      
747        eorStream = createOStream(eorFilename_);
748  
749   #ifdef IS_MPI
750      }
751 < #endif // is_mpi    
752 <
751 > #endif
752 >    
753      writeFrame(*eorStream);
754 <
754 >      
755   #ifdef IS_MPI
756      if (worldRank == 0) {
757 < #endif // is_mpi
757 > #endif
758 >      
759        writeClosing(*eorStream);
760        delete eorStream;
761 +      
762   #ifdef IS_MPI
763      }
764   #endif // is_mpi  
# Line 655 | Line 772 | namespace OpenMD {
772   #ifdef IS_MPI
773      if (worldRank == 0) {
774   #endif // is_mpi
658 –
775        buffers.push_back(dumpFile_->rdbuf());
660 –
776        eorStream = createOStream(eorFilename_);
662 –
777        buffers.push_back(eorStream->rdbuf());
664 –        
778   #ifdef IS_MPI
779      }
780   #endif // is_mpi    
781  
782      TeeBuf tbuf(buffers.begin(), buffers.end());
783      std::ostream os(&tbuf);
671 –
784      writeFrame(os);
785  
786   #ifdef IS_MPI
# Line 678 | Line 790 | namespace OpenMD {
790        delete eorStream;
791   #ifdef IS_MPI
792      }
793 < #endif // is_mpi  
682 <    
793 > #endif // is_mpi      
794    }
795  
796    std::ostream* DumpWriter::createOStream(const std::string& filename) {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines