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

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1791 by gezelter, Fri Aug 31 16:48:50 2012 UTC vs.
Revision 1804 by gezelter, Wed Oct 17 19:04:44 2012 UTC

# Line 281 | Line 281 | namespace OpenMD {
281      // do some sanity checking
282  
283      int selectionCount = seleMan_.getSelectionCount();
284 +
285      int nIntegrable = info->getNGlobalIntegrableObjects();
286  
287      if (selectionCount > nIntegrable) {
# Line 451 | Line 452 | namespace OpenMD {
452  
453      int selei;
454      StuntDouble* sd;
454 –    int idx;
455  
456      RealType min_val;
457      bool min_found = false;  
# Line 464 | Line 464 | namespace OpenMD {
464      for (sd = seleMan_.beginSelected(selei); sd != NULL;
465           sd = seleMan_.nextSelected(selei)) {
466  
467 –      idx = sd->getLocalIndex();
468 –
467        Vector3d pos = sd->getPos();
468  
469        // wrap the stuntdouble's position back into the box:
# Line 543 | Line 541 | namespace OpenMD {
541        }
542      }
543      
544 < #ifdef IS_MPI
545 <    int nProc, worldRank;
544 > #ifdef IS_MPI    
545 >    int worldRank = MPI::COMM_WORLD.Get_rank();
546      
549 –    nProc = MPI::COMM_WORLD.Get_size();
550 –    worldRank = MPI::COMM_WORLD.Get_rank();
551 –
547      bool my_min_found = min_found;
548      bool my_max_found = max_found;
549  
# Line 781 | Line 776 | namespace OpenMD {
776  
777      int selei;
778      StuntDouble* sd;
784 –    int idx;
779  
780      vector<StuntDouble*> hotBin, coldBin;
781  
# Line 802 | Line 796 | namespace OpenMD {
796  
797      for (sd = seleMan_.beginSelected(selei); sd != NULL;
798           sd = seleMan_.nextSelected(selei)) {
805 –
806 –      idx = sd->getLocalIndex();
799  
800        Vector3d pos = sd->getPos();
801  
# Line 917 | Line 909 | namespace OpenMD {
909  
910        if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients
911          c = sqrt(c);
912 <        //std::cerr << "cold slab scaling coefficient: " << c << endl;
921 <        //now convert to hotBin coefficient
912 >
913          RealType w = 0.0;
914          if (rnemdFluxType_ ==  rnemdFullKE) {
915            x = 1.0 + px * (1.0 - c);
# Line 956 | Line 947 | namespace OpenMD {
947              }
948            }
949            w = sqrt(w);
959 –          // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
960 –          //           << "\twh= " << w << endl;
950            for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
951              if (rnemdFluxType_ == rnemdFullKE) {
952                vel = (*sdi)->getVel();
# Line 1231 | Line 1220 | namespace OpenMD {
1220  
1221      int selei;
1222      StuntDouble* sd;
1234 –    int idx;
1223  
1224      vector<StuntDouble*> hotBin, coldBin;
1225  
# Line 1246 | Line 1234 | namespace OpenMD {
1234      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1235           sd = seleMan_.nextSelected(selei)) {
1236  
1249 –      idx = sd->getLocalIndex();
1250 –
1237        Vector3d pos = sd->getPos();
1238  
1239        // wrap the stuntdouble's position back into the box:
# Line 1266 | Line 1252 | namespace OpenMD {
1252        
1253          if (inA) {
1254            hotBin.push_back(sd);
1269 –          //std::cerr << "before, velocity = " << vel << endl;
1255            Ph += mass * vel;
1271 –          //std::cerr << "after, velocity = " << vel << endl;
1256            Mh += mass;
1257            Kh += mass * vel.lengthSquare();
1258            if (rnemdFluxType_ == rnemdFullKE) {
# Line 1316 | Line 1300 | namespace OpenMD {
1300      
1301      Kh *= 0.5;
1302      Kc *= 0.5;
1319 –
1320 –    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1321 –    //        << "\tKc= " << Kc << endl;
1322 –    // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1303      
1304   #ifdef IS_MPI
1305      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
# Line 1351 | Line 1331 | namespace OpenMD {
1331                if (hDenominator > 0.0) {
1332                  RealType h = sqrt(hNumerator / hDenominator);
1333                  if ((h > 0.9) && (h < 1.1)) {
1334 <                  // std::cerr << "cold slab scaling coefficient: " << c << "\n";
1355 <                  // std::cerr << "hot slab scaling coefficient: " << h <<  "\n";
1334 >
1335                    vector<StuntDouble*>::iterator sdi;
1336                    Vector3d vel;
1337                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
# Line 1427 | Line 1406 | namespace OpenMD {
1406  
1407      seleMan_.setSelectionSet(evaluator_.evaluate());
1408  
1409 <    int selei;
1409 >    int selei(0);
1410      StuntDouble* sd;
1432 –    int idx;
1411  
1412      vector<RealType> binMass(nBins_, 0.0);
1413      vector<RealType> binPx(nBins_, 0.0);
# Line 1453 | Line 1431 | namespace OpenMD {
1431               sd != NULL;
1432               sd = mol->nextIntegrableObject(iiter))
1433      */
1434 +
1435      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1436 <         sd = seleMan_.nextSelected(selei)) {
1437 <      
1459 <      idx = sd->getLocalIndex();
1460 <      
1436 >         sd = seleMan_.nextSelected(selei)) {    
1437 >    
1438        Vector3d pos = sd->getPos();
1439  
1440        // wrap the stuntdouble's position back into the box:
# Line 1472 | Line 1449 | namespace OpenMD {
1449        // The modulo operator is used to wrap the case when we are
1450        // beyond the end of the bins back to the beginning.
1451        int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1452 <    
1452 >
1453        RealType mass = sd->getMass();
1454        Vector3d vel = sd->getVel();
1455  
# Line 1503 | Line 1480 | namespace OpenMD {
1480        }
1481      }
1482      
1506 –
1483   #ifdef IS_MPI
1484      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1485                                nBins_, MPI::INT, MPI::SUM);
# Line 1534 | Line 1510 | namespace OpenMD {
1510        den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1511          / currentSnap_->getVolume() ;
1512  
1513 <      temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1514 <                               PhysicalConstants::energyConvert);
1515 <
1516 <      for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1517 <        if(outputMask_[j]) {
1518 <          switch(j) {
1519 <          case Z:
1520 <            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1521 <            break;
1522 <          case TEMPERATURE:
1523 <            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1524 <            break;
1525 <          case VELOCITY:
1526 <            dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1527 <            break;
1528 <          case DENSITY:
1529 <            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1530 <            break;
1513 >      if (binCount[i] > 0) {
1514 >        // only add values if there are things to add
1515 >        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1516 >                                 PhysicalConstants::energyConvert);
1517 >        
1518 >        for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1519 >          if(outputMask_[j]) {
1520 >            switch(j) {
1521 >            case Z:
1522 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1523 >              break;
1524 >            case TEMPERATURE:
1525 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1526 >              break;
1527 >            case VELOCITY:
1528 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1529 >              break;
1530 >            case DENSITY:
1531 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1532 >              break;
1533 >            }
1534            }
1535          }
1536        }

Diff Legend

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