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 1921 by gezelter, Thu Aug 1 18:23:07 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 302 | 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;
# Line 310 | Line 310 | namespace OpenMD {
310      SimInfo::MoleculeIterator mi;
311      Molecule::IntegrableObjectIterator ii;
312      RigidBody::AtomIterator ai;
313 –    Atom* atom;
313  
314   #ifndef IS_MPI
315      os << "  <Snapshot>\n";
# Line 318 | Line 317 | namespace OpenMD {
317      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
318  
319      os << "    <StuntDoubles>\n";
320 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
321 <
320 >    for (mol = info_->beginMolecule(mi); mol != NULL;
321 >         mol = info_->nextMolecule(mi)) {
322        
323        for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
324             sd = mol->nextIntegrableObject(ii)) {        
# Line 331 | Line 330 | namespace OpenMD {
330  
331      if (doSiteData_) {
332        os << "    <SiteData>\n";
333 <      for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
333 >      for (mol = info_->beginMolecule(mi); mol != NULL;
334 >           mol = info_->nextMolecule(mi)) {
335                
336          for (sd = mol->beginIntegrableObject(ii); sd != NULL;  
337 <             sd = mol->nextIntegrableObject(ii)) {      
338 <          
337 >           sd = mol->nextIntegrableObject(ii)) {        
338 >
339            int ioIndex = sd->getGlobalIntegrableObjectIndex();
340            // do one for the IO itself
341            os << prepareSiteLine(sd, ioIndex, 0);
# Line 344 | Line 344 | namespace OpenMD {
344              
345              RigidBody* rb = static_cast<RigidBody*>(sd);
346              int siteIndex = 0;
347 <            for (atom = rb->beginAtom(ai); atom != NULL;  
347 >            for (Atom* atom = rb->beginAtom(ai); atom != NULL;  
348                   atom = rb->nextAtom(ai)) {                                            
349                os << prepareSiteLine(atom, ioIndex, siteIndex);
350                siteIndex++;
# Line 358 | Line 358 | namespace OpenMD {
358  
359      os.flush();
360   #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)) {
361 >
362 >    const int masterNode = 0;
363 >    int worldRank = MPI::COMM_WORLD.Get_rank();
364 >    int nProc = MPI::COMM_WORLD.Get_size();
365  
366 +    if (worldRank == masterNode) {      
367 +      os << "  <Snapshot>\n";  
368 +      writeFrameProperties(os,
369 +                           info_->getSnapshotManager()->getCurrentSnapshot());
370 +      os << "    <StuntDoubles>\n";
371 +    }
372  
373 +    //every node prepares the dump lines for integrable objects belong to itself
374 +    std::string buffer;
375 +    for (mol = info_->beginMolecule(mi); mol != NULL;
376 +         mol = info_->nextMolecule(mi)) {
377        for (sd = mol->beginIntegrableObject(ii); sd != NULL;
378             sd = mol->nextIntegrableObject(ii)) {        
379 <          buffer += prepareDumpLine(sd);
379 >        buffer += prepareDumpLine(sd);
380        }
381      }
382      
372 –    const int masterNode = 0;
373 –    int nProc;
374 –    MPI_Comm_size(MPI_COMM_WORLD, &nProc);
383      if (worldRank == masterNode) {      
376 –      os << "  <Snapshot>\n";  
377 –      writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
378 –      os << "    <StuntDoubles>\n";
379 –        
384        os << buffer;
385 <
385 >      
386        for (int i = 1; i < nProc; ++i) {
387 +        // tell processor i to start sending us data:
388 +        MPI::COMM_WORLD.Bcast(&i, 1, MPI::INT, masterNode);
389  
390          // receive the length of the string buffer that was
391 <        // prepared by processor i
386 <
387 <        MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
391 >        // prepared by processor i:        
392          int recvLength;
393 <        MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
393 >        MPI::COMM_WORLD.Recv(&recvLength, 1, MPI::INT, i, MPI::ANY_TAG,
394 >                             istatus);
395 >
396 >        // create a buffer to receive the data
397          char* recvBuffer = new char[recvLength];
398          if (recvBuffer == NULL) {
399          } else {
400 <          MPI_Recv(recvBuffer, recvLength, MPI_CHAR, i, 0, MPI_COMM_WORLD, &istatus);
400 >          // receive the data:
401 >          MPI::COMM_WORLD.Recv(recvBuffer, recvLength, MPI::CHAR, i,
402 >                               MPI::ANY_TAG, istatus);
403 >          // send it to the file:
404            os << recvBuffer;
405 +          // get rid of the receive buffer:
406            delete [] recvBuffer;
407          }
408        }
398 –      os << "    </StuntDoubles>\n";
399 –      
400 –      os << "  </Snapshot>\n";
401 –      os.flush();
409      } else {
410        int sendBufferLength = buffer.size() + 1;
411        int myturn = 0;
412        for (int i = 1; i < nProc; ++i){
413 <        MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
413 >        // wait for the master node to call our number:
414 >        MPI::COMM_WORLD.Bcast(&myturn, 1, MPI::INT, masterNode);
415          if (myturn == worldRank){
416 <          MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
417 <          MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
416 >          // send the length of our buffer:
417 >          MPI::COMM_WORLD.Send(&sendBufferLength, 1, MPI::INT, masterNode, 0);
418 >
419 >          // send our buffer:
420 >          MPI::COMM_WORLD.Send((void *)buffer.c_str(), sendBufferLength,
421 >                               MPI::CHAR, masterNode, 0);
422 >
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 +          if (sd->isRigidBody()) {
447 +            
448 +            RigidBody* rb = static_cast<RigidBody*>(sd);
449 +            int siteIndex = 0;
450 +            for (Atom* 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 +      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) {
# Line 527 | Line 625 | namespace OpenMD {
625    }
626  
627    std::string DumpWriter::prepareSiteLine(StuntDouble* sd, int ioIndex, int siteIndex) {
628 <        
628 >    int storageLayout = info_->getSnapshotManager()->getStorageLayout();
629  
630      std::string id;
631      std::string type;
# Line 543 | Line 641 | namespace OpenMD {
641      }
642                
643      if (needFlucQ_) {
644 <      type += "cw";
645 <      RealType fqPos = sd->getFlucQPos();
646 <      if (isinf(fqPos) || isnan(fqPos) ) {      
647 <        sprintf( painCave.errMsg,
648 <                 "DumpWriter detected a numerical error writing the"
649 <                 " fluctuating charge for object %s", id.c_str());      
650 <        painCave.isFatal = 1;
651 <        simError();
652 <      }
653 <      sprintf(tempBuffer, " %13e ", fqPos);
654 <      line += tempBuffer;
655 <    
656 <      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;
644 >      if (storageLayout & DataStorage::dslFlucQPosition) {
645 >        type += "c";
646 >        RealType fqPos = sd->getFlucQPos();
647 >        if (isinf(fqPos) || isnan(fqPos) ) {      
648 >          sprintf( painCave.errMsg,
649 >                   "DumpWriter detected a numerical error writing the"
650 >                   " fluctuating charge for object %s", id.c_str());      
651 >          painCave.isFatal = 1;
652 >          simError();
653 >        }
654 >        sprintf(tempBuffer, " %13e ", fqPos);
655 >        line += tempBuffer;
656 >      }
657  
658 <      if (needForceVector_) {
659 <        type += "g";
660 <        RealType fqFrc = sd->getFlucQFrc();        
661 <        if (isinf(fqFrc) || isnan(fqFrc) ) {      
658 >      if (storageLayout & DataStorage::dslFlucQVelocity) {
659 >        type += "w";    
660 >        RealType fqVel = sd->getFlucQVel();
661 >        if (isinf(fqVel) || isnan(fqVel) ) {      
662            sprintf( painCave.errMsg,
663                     "DumpWriter detected a numerical error writing the"
664 <                   " fluctuating charge force for object %s", id.c_str());      
664 >                   " fluctuating charge velocity for object %s", id.c_str());      
665            painCave.isFatal = 1;
666            simError();
667          }
668 <        sprintf(tempBuffer, " %13e ", fqFrc);        
668 >        sprintf(tempBuffer, " %13e ", fqVel);
669          line += tempBuffer;
670        }
582 –    }
671  
672 +      if (needForceVector_) {
673 +        if (storageLayout & DataStorage::dslFlucQForce) {          
674 +          type += "g";
675 +          RealType fqFrc = sd->getFlucQFrc();        
676 +          if (isinf(fqFrc) || isnan(fqFrc) ) {      
677 +            sprintf( painCave.errMsg,
678 +                     "DumpWriter detected a numerical error writing the"
679 +                     " fluctuating charge force for object %s", id.c_str());      
680 +            painCave.isFatal = 1;
681 +            simError();
682 +          }
683 +          sprintf(tempBuffer, " %13e ", fqFrc);        
684 +          line += tempBuffer;
685 +        }
686 +      }
687 +    }
688 +    
689      if (needElectricField_) {
690 <      type += "e";
691 <      Vector3d eField= sd->getElectricField();
692 <      if (isinf(eField[0]) || isnan(eField[0]) ||
693 <          isinf(eField[1]) || isnan(eField[1]) ||
694 <          isinf(eField[2]) || isnan(eField[2]) ) {      
695 <        sprintf( painCave.errMsg,
696 <                 "DumpWriter detected a numerical error writing the electric"
697 <                 " field for object %s", id.c_str());      
698 <        painCave.isFatal = 1;
699 <        simError();
690 >      if (storageLayout & DataStorage::dslElectricField) {
691 >        type += "e";
692 >        Vector3d eField= sd->getElectricField();
693 >        if (isinf(eField[0]) || isnan(eField[0]) ||
694 >            isinf(eField[1]) || isnan(eField[1]) ||
695 >            isinf(eField[2]) || isnan(eField[2]) ) {      
696 >          sprintf( painCave.errMsg,
697 >                   "DumpWriter detected a numerical error writing the electric"
698 >                   " field for object %s", id.c_str());      
699 >          painCave.isFatal = 1;
700 >          simError();
701 >        }
702 >        sprintf(tempBuffer, " %13e %13e %13e",
703 >                eField[0], eField[1], eField[2]);
704 >        line += tempBuffer;
705        }
596 –      sprintf(tempBuffer, " %13e %13e %13e",
597 –              eField[0], eField[1], eField[2]);
598 –      line += tempBuffer;
706      }
707  
708  
709      if (needParticlePot_) {
710 <      type += "u";
711 <      RealType particlePot = sd->getParticlePot();
712 <      if (isinf(particlePot) || isnan(particlePot)) {      
713 <        sprintf( painCave.errMsg,
714 <                 "DumpWriter detected a numerical error writing the particle "
715 <                 " potential for object %s", id.c_str());      
716 <        painCave.isFatal = 1;
717 <        simError();
710 >      if (storageLayout & DataStorage::dslParticlePot) {
711 >        type += "u";
712 >        RealType particlePot = sd->getParticlePot();
713 >        if (isinf(particlePot) || isnan(particlePot)) {      
714 >          sprintf( painCave.errMsg,
715 >                   "DumpWriter detected a numerical error writing the particle "
716 >                   " potential for object %s", id.c_str());      
717 >          painCave.isFatal = 1;
718 >          simError();
719 >        }
720 >        sprintf(tempBuffer, " %13e", particlePot);
721 >        line += tempBuffer;
722        }
612 –      sprintf(tempBuffer, " %13e", particlePot);
613 –      line += tempBuffer;
723      }
724 <    
616 <
724 >  
725      sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
726      return std::string(tempBuffer);
727    }
# Line 628 | Line 736 | namespace OpenMD {
736   #ifdef IS_MPI
737      if (worldRank == 0) {
738   #endif // is_mpi
739 <
739 >      
740        eorStream = createOStream(eorFilename_);
741  
742   #ifdef IS_MPI
743      }
744 < #endif // is_mpi    
745 <
744 > #endif
745 >    
746      writeFrame(*eorStream);
747 <
747 >      
748   #ifdef IS_MPI
749      if (worldRank == 0) {
750 < #endif // is_mpi
750 > #endif
751 >      
752        writeClosing(*eorStream);
753        delete eorStream;
754 +      
755   #ifdef IS_MPI
756      }
757   #endif // is_mpi  

Diff Legend

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