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

Comparing branches/development/src/rnemd/RNEMD.cpp (file contents):
Revision 1866 by gezelter, Thu Apr 25 14:32:56 2013 UTC vs.
Revision 1876 by gezelter, Fri May 17 17:10:11 2013 UTC

# Line 71 | Line 71 | namespace OpenMD {
71    RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
72                                  evaluatorA_(info), seleManA_(info),
73                                  commonA_(info), evaluatorB_(info),
74 <                                seleManB_(info), commonB_(info),
74 >                                seleManB_(info), commonB_(info),
75 >                                hasData_(false), hasDividingArea_(false),
76                                  usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
77  
78      trialCount_ = 0;
# Line 603 | Line 604 | namespace OpenMD {
604   #ifdef IS_MPI
605      }
606   #endif
607 +
608 +    // delete all of the objects we created:
609 +    delete areaAccumulator_;    
610 +    data_.clear();
611    }
612    
613    void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) {
# Line 1146 | Line 1151 | namespace OpenMD {
1151            //if w is in the right range, so should be x, y, z.
1152            vector<StuntDouble*>::iterator sdi;
1153            Vector3d vel;
1154 <          for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1154 >          for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1155              if (rnemdFluxType_ == rnemdFullKE) {
1156                vel = (*sdi)->getVel() * c;
1157                (*sdi)->setVel(vel);
# Line 1157 | Line 1162 | namespace OpenMD {
1162              }
1163            }
1164            w = sqrt(w);
1165 <          for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1165 >          for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1166              if (rnemdFluxType_ == rnemdFullKE) {
1167                vel = (*sdi)->getVel();
1168                vel.x() *= x;
# Line 1276 | Line 1281 | namespace OpenMD {
1281        vector<RealType>::iterator ri;
1282        RealType r1, r2, alpha0;
1283        vector<pair<RealType,RealType> > rps;
1284 <      for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
1284 >      for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) {
1285          r2 = *ri;
1286          //check if FindRealRoots() give the right answer
1287          if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
# Line 1308 | Line 1313 | namespace OpenMD {
1313          RealType diff;
1314          pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
1315          vector<pair<RealType,RealType> >::iterator rpi;
1316 <        for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
1316 >        for (rpi = rps.begin(); rpi != rps.end(); ++rpi) {
1317            r1 = (*rpi).first;
1318            r2 = (*rpi).second;
1319            switch(rnemdFluxType_) {
# Line 1375 | Line 1380 | namespace OpenMD {
1380          }
1381          vector<StuntDouble*>::iterator sdi;
1382          Vector3d vel;
1383 <        for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1383 >        for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1384            vel = (*sdi)->getVel();
1385            vel.x() *= x;
1386            vel.y() *= y;
# Line 1386 | Line 1391 | namespace OpenMD {
1391          x = 1.0 + px * (1.0 - x);
1392          y = 1.0 + py * (1.0 - y);
1393          z = 1.0 + pz * (1.0 - z);
1394 <        for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1394 >        for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1395            vel = (*sdi)->getVel();
1396            vel.x() *= x;
1397            vel.y() *= y;
# Line 1592 | Line 1597 | namespace OpenMD {
1597  
1598      Vector3d ac, acrec, bc, bcrec;
1599      Vector3d ah, ahrec, bh, bhrec;
1595 –    RealType cNumerator, cDenominator;
1596 –    RealType hNumerator, hDenominator;
1600  
1598 –
1601      bool successfulExchange = false;
1602      if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1603        Vector3d vc = Pc / Mc;
# Line 1609 | Line 1611 | namespace OpenMD {
1611        bc  = -(Ici * angularMomentumTarget_) + omegac;
1612        bcrec = bc - omegac;
1613        
1614 <      cNumerator = Kc - kineticTarget_;
1614 >      RealType cNumerator = Kc - kineticTarget_;
1615        if (doLinearPart)
1616          cNumerator -= 0.5 * Mc * ac.lengthSquare();
1617        
# Line 1618 | Line 1620 | namespace OpenMD {
1620  
1621        if (cNumerator > 0.0) {
1622          
1623 <        cDenominator = Kc;
1623 >        RealType cDenominator = Kc;
1624  
1625          if (doLinearPart)
1626            cDenominator -= 0.5 * Mc * vc.lengthSquare();
# Line 1641 | Line 1643 | namespace OpenMD {
1643              bh  = (Ihi * angularMomentumTarget_) + omegah;
1644              bhrec = bh - omegah;
1645              
1646 <            hNumerator = Kh + kineticTarget_;
1646 >            RealType hNumerator = Kh + kineticTarget_;
1647              if (doLinearPart)
1648                hNumerator -= 0.5 * Mh * ah.lengthSquare();
1649              
# Line 1650 | Line 1652 | namespace OpenMD {
1652                
1653              if (hNumerator > 0.0) {
1654                
1655 <              hDenominator = Kh;
1655 >              RealType hDenominator = Kh;
1656                if (doLinearPart)
1657                  hDenominator -= 0.5 * Mh * vh.lengthSquare();
1658                if (doAngularPart)
# Line 1664 | Line 1666 | namespace OpenMD {
1666                    Vector3d vel;
1667                    Vector3d rPos;
1668                    
1669 <                  for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1669 >                  for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1670                      //vel = (*sdi)->getVel();
1671                      rPos = (*sdi)->getPos() - coordinateOrigin_;
1672                      if (doLinearPart)
# Line 1680 | Line 1682 | namespace OpenMD {
1682                        }
1683                      }
1684                    }
1685 <                  for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1685 >                  for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1686                      //vel = (*sdi)->getVel();
1687                      rPos = (*sdi)->getPos() - coordinateOrigin_;
1688                      if (doLinearPart)
# Line 1730 | Line 1732 | namespace OpenMD {
1732        int isd;
1733        StuntDouble* sd;
1734        vector<StuntDouble*> aSites;
1733 –      ConvexHull* surfaceMeshA = new ConvexHull();
1735        seleManA_.setSelectionSet(evaluatorA_.evaluate());
1736        for (sd = seleManA_.beginSelected(isd); sd != NULL;
1737             sd = seleManA_.nextSelected(isd)) {
1738          aSites.push_back(sd);
1739        }
1740 <
1740 > #if defined(HAVE_QHULL)
1741 >      ConvexHull* surfaceMeshA = new ConvexHull();
1742        surfaceMeshA->computeHull(aSites);
1743        areaA = surfaceMeshA->getArea();
1744 +      delete surfaceMeshA;
1745 + #else
1746 +      sprintf( painCave.errMsg,
1747 +               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1748 +               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1749 +      painCave.severity = OPENMD_ERROR;
1750 +      painCave.isFatal = 1;
1751 +      simError();
1752 + #endif
1753 +
1754      } else {
1755        if (usePeriodicBoundaryConditions_) {
1756          // in periodic boundaries, the surface area is twice the x-y
# Line 1752 | Line 1764 | namespace OpenMD {
1764        }
1765      }
1766  
1755 –
1756 –
1767      if (hasSelectionB_) {
1768        int isd;
1769        StuntDouble* sd;
1770        vector<StuntDouble*> bSites;
1761 –
1771        seleManB_.setSelectionSet(evaluatorB_.evaluate());
1763 –    
1772        for (sd = seleManB_.beginSelected(isd); sd != NULL;
1773             sd = seleManB_.nextSelected(isd)) {
1774          bSites.push_back(sd);
1775        }
1776  
1777 <      ConvexHull* surfaceMeshB = new ConvexHull();
1778 <
1777 > #if defined(HAVE_QHULL)
1778 >      ConvexHull* surfaceMeshB = new ConvexHull();    
1779        surfaceMeshB->computeHull(bSites);
1780        areaB = surfaceMeshB->getArea();
1781 +      delete surfaceMeshB;
1782 + #else
1783 +      sprintf( painCave.errMsg,
1784 +               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1785 +               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1786 +      painCave.severity = OPENMD_ERROR;
1787 +      painCave.isFatal = 1;
1788 +      simError();
1789 + #endif
1790 +
1791 +
1792      } else {
1793        if (usePeriodicBoundaryConditions_) {
1794          // in periodic boundaries, the surface area is twice the x-y
# Line 2009 | Line 2028 | namespace OpenMD {
2028          }
2029        }
2030      }
2031 +    hasData_ = true;
2032    }
2033  
2034    void RNEMD::getStarted() {
# Line 2041 | Line 2061 | namespace OpenMD {
2061    
2062    void RNEMD::writeOutputFile() {
2063      if (!doRNEMD_) return;
2064 +    if (!hasData_) return;
2065      
2066   #ifdef IS_MPI
2067      // If we're the root node, should we print out the results
# Line 2062 | Line 2083 | namespace OpenMD {
2083        RealType time = currentSnap_->getTime();
2084        RealType avgArea;
2085        areaAccumulator_->getAverage(avgArea);
2065 –      RealType Jz = kineticExchange_ / (time * avgArea)
2066 –        / PhysicalConstants::energyConvert;
2067 –      Vector3d JzP = momentumExchange_ / (time * avgArea);      
2068 –      Vector3d JzL = angularMomentumExchange_ / (time * avgArea);      
2086  
2087 +      RealType Jz(0.0);
2088 +      Vector3d JzP(V3Zero);
2089 +      Vector3d JzL(V3Zero);
2090 +      if (time >= info_->getSimParams()->getDt()) {
2091 +        Jz = kineticExchange_ / (time * avgArea)
2092 +          / PhysicalConstants::energyConvert;
2093 +        JzP = momentumExchange_ / (time * avgArea);
2094 +        JzL = angularMomentumExchange_ / (time * avgArea);
2095 +      }
2096 +
2097        rnemdFile_ << "#######################################################\n";
2098        rnemdFile_ << "# RNEMD {\n";
2099  
# Line 2219 | Line 2246 | namespace OpenMD {
2246        rnemdFile_ << "\t" << s;
2247      } else{
2248        sprintf( painCave.errMsg,
2249 <               "RNEMD detected a numerical error writing: %s for bin %d",
2249 >               "RNEMD detected a numerical error writing: %s for bin %u",
2250                 data_[index].title.c_str(), bin);
2251        painCave.isFatal = 1;
2252        simError();
# Line 2242 | Line 2269 | namespace OpenMD {
2269          isinf(s[1]) || isnan(s[1]) ||
2270          isinf(s[2]) || isnan(s[2]) ) {      
2271        sprintf( painCave.errMsg,
2272 <               "RNEMD detected a numerical error writing: %s for bin %d",
2272 >               "RNEMD detected a numerical error writing: %s for bin %u",
2273                 data_[index].title.c_str(), bin);
2274        painCave.isFatal = 1;
2275        simError();
# Line 2267 | Line 2294 | namespace OpenMD {
2294        rnemdFile_ << "\t" << s;
2295      } else{
2296        sprintf( painCave.errMsg,
2297 <               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
2297 >               "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2298                 data_[index].title.c_str(), bin);
2299        painCave.isFatal = 1;
2300        simError();
# Line 2289 | Line 2316 | namespace OpenMD {
2316          isinf(s[1]) || isnan(s[1]) ||
2317          isinf(s[2]) || isnan(s[2]) ) {      
2318        sprintf( painCave.errMsg,
2319 <               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
2319 >               "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2320                 data_[index].title.c_str(), bin);
2321        painCave.isFatal = 1;
2322        simError();

Diff Legend

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