NAMD: Patch.C Source File

00001 
00007 #include "InfoStream.h"
00008 #include "Patch.h"
00009 #include "PatchMap.h"
00010 #include "Compute.h"
00011 
00012 #include "AtomMap.h"
00013 #include "ComputeMap.h"
00014 #include "Node.h"
00015 #include "Molecule.h"
00016 #include "SimParameters.h"
00017 #include "ResizeArrayPrimIter.h"
00018 #include "ComputeNonbondedMICKernel.h"
00019 #include "Priorities.h"
00020 #include "ReductionMgr.h"
00021 
00022 #include "Sync.h"
00023 
00024 #include <algorithm>
00025 
00026 typedef ResizeArrayPrimIter<Compute*> ComputePtrListIter;
00027 
00028 struct cptr_sortop_priority {
00029   bool operator() (Compute *i, Compute *j) const {
00030     return ( i->priority() < j->priority() );
00031   }
00032 };
00033 
00034 //#define  DEBUGM
00035 #define MIN_DEBUG_LEVEL 4
00036 #include "Debug.h"
00037 
00038 Patch::~Patch() {
00039   delete atomMapper;
00040 #ifndef NAMD_CUDA
00041   delete reduction;
00042 #endif
00043 }
00044 
00045 Patch::Patch(PatchID pd) :
00046    lattice(flags.lattice),
00047    patchID(pd), numAtoms(0), numFixedAtoms(0),
00048    avgPositionPtrBegin(0), avgPositionPtrEnd(0),
00049    velocityPtrBegin(0), velocityPtrEnd(0),      // BEGIN LA, END LA
00050    positionBox(this,&Patch::positionBoxClosed,pd,0),
00051    avgPositionBox(this,&Patch::avgPositionBoxClosed,pd,3),
00052    velocityBox(this,&Patch::velocityBoxClosed,pd,4), // BEGIN LA, END LA
00053    psiSumBox(this,&Patch::psiSumBoxClosed,pd,5), // begin gbis
00054    intRadBox(this,&Patch::intRadBoxClosed,pd,6),
00055    bornRadBox(this,&Patch::bornRadBoxClosed,pd,7),
00056    dEdaSumBox(this,&Patch::dEdaSumBoxClosed,pd,8),
00057    dHdrPrefixBox(this,&Patch::dHdrPrefixBoxClosed,pd,9), //end gbis
00058    lcpoTypeBox(this,&Patch::lcpoTypeBoxClosed,pd,10),
00059    forceBox(this,&Patch::forceBoxClosed,pd,1),
00060    boxesOpen(0), _hasNewAtoms(0),
00061    computesSortedByPriority(0), firstHoldableCompute(0)
00062 
00063    // DMK - Atom Separation (water vs. non-water)
00064    #if NAMD_SeparateWaters != 0
00065      ,numWaterAtoms(-1)
00066    #endif
00067 {
00068   //CkPrintf("GBIS: PatchCreated\n");
00069 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00070     positionPtrBegin = 0;
00071     positionPtrEnd = 0;
00072 #endif
00073 
00074         nChild = 0;
00075         child = NULL;
00076 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00077   #ifdef USE_NODEPATCHMGR
00078   nodeChildren = NULL;
00079   numNodeChild = 0;
00080   #endif
00081 #endif
00082 
00083   lattice = Node::Object()->simParameters->lattice;
00084   atomMapper = new AtomMapper(pd);
00085 #ifndef NAMD_CUDA
00086   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00087 #endif
00088 
00089   // DMK
00090   #if defined(NAMD_MIC)
00091     cudaAtomPtr = NULL;
00092     #if MIC_SUBMIT_ATOMS_ON_ARRIVAL != 0
00093       pthread_mutex_init(&mic_atomData_mutex, NULL);
00094       mic_atomData = NULL;
00095       mic_atomData_seq = -1;
00096       mic_atomData_allocSize_host = 0;
00097       for (int i = 0; i < MIC_MAX_DEVICES_PER_NODE; i++) {
00098         mic_atomData_prev[i] = NULL;
00099         mic_atomData_deviceSeq[i] = -1;
00100         mic_atomData_devicePtr[i] = 0;
00101         mic_atomData_allocSize_device[i] = 0;
00102       }
00103     #endif
00104   #endif
00105 }
00106 
00107 Box<Patch,CompAtom>* Patch::registerPositionPickup(Compute *cid)
00108 {
00109    if ( computesSortedByPriority ) {
00110      positionComputeList.sort();
00111      computesSortedByPriority = 0;
00112    }
00113    //DebugM(4, "registerPositionPickupa("<<patchID<<") from " << cid->cid << "\n");
00114    if (positionComputeList.add(cid) < 0)
00115    {
00116      DebugM(7, "registerPositionPickup() failed for cid " << cid->cid << std::endl);
00117      return NULL;
00118    }
00119    return positionBox.checkOut(cid->cid);
00120 }
00121 
00122 void Patch::unregisterPositionPickup(Compute *cid, Box<Patch,CompAtom> **const box)
00123 {
00124    if ( computesSortedByPriority ) {
00125      positionComputeList.sort();
00126      computesSortedByPriority = 0;
00127    }
00128    DebugM(4, "UnregisterPositionPickup from " << cid->cid << "\n");
00129    positionComputeList.del(cid);
00130    positionBox.checkIn(*box);
00131    *box = 0;
00132 }
00133 
00134 Box<Patch,CompAtom>* Patch::registerAvgPositionPickup(Compute *cid)
00135 {
00136    //DebugM(4, "registerAvgPositionPickup("<<patchID<<") from " << cid->cid << "\n");
00137    return avgPositionBox.checkOut(cid->cid);
00138 }
00139 
00140 void Patch::unregisterAvgPositionPickup(Compute *cid, Box<Patch,CompAtom> **const box)
00141 {
00142    DebugM(4, "UnregisterAvgPositionPickup from " << cid->cid << "\n");
00143    avgPositionBox.checkIn(*box);
00144    *box = 0;
00145 }
00146 
00147 // BEGIN LA
00148 Box<Patch,CompAtom>* Patch::registerVelocityPickup(Compute *cid)
00149 {
00150    //DebugM(4, "registerVelocityPickup("<<patchID<<") from " << cid->cid << "\n");
00151    return velocityBox.checkOut(cid->cid);
00152 }
00153 
00154 void Patch::unregisterVelocityPickup(Compute *cid, Box<Patch,CompAtom> **const box)
00155 {
00156    DebugM(4, "UnregisterVelocityPickup from " << cid->cid << "\n");
00157    velocityBox.checkIn(*box);
00158    *box = 0;
00159 }
00160 // END LA
00161 
00162 //begin gbis
00163 //deposit, not pickup
00164 Box<Patch,GBReal>* Patch::registerPsiSumDeposit(Compute *cid) {
00165 
00166   if (psiSumComputeList.add(cid) < 0) {
00167     DebugM(7, "registerPsiSumDeposit() failed for cid " << cid->cid << std::endl);
00168     DebugM(7, "  size of psiSumCompueList " << psiSumComputeList.size() << std::endl);
00169      return NULL;
00170   }
00171   return psiSumBox.checkOut(cid->cid);
00172 }
00173 
00174 void Patch::unregisterPsiSumDeposit(Compute *cid,Box<Patch,GBReal> **const box) {
00175   psiSumComputeList.del(cid);
00176   psiSumBox.checkIn(*box);
00177   *box = 0;
00178 }
00179 Box<Patch,Real>* Patch::registerIntRadPickup(Compute *cid) {
00180   return intRadBox.checkOut(cid->cid);
00181 }
00182 void Patch::unregisterIntRadPickup(Compute *cid,Box<Patch,Real> **const box) {
00183   intRadBox.checkIn(*box);
00184   *box = 0;
00185 }
00186 
00187 //LCPO
00188 Box<Patch,int>* Patch::registerLcpoTypePickup(Compute *cid) {
00189   return lcpoTypeBox.checkOut(cid->cid);
00190 }
00191 void Patch::unregisterLcpoTypePickup(Compute *cid,Box<Patch,int> **const box) {
00192   lcpoTypeBox.checkIn(*box);
00193   *box = 0;
00194 }
00195 
00196 Box<Patch,Real>* Patch::registerBornRadPickup(Compute *cid) {
00197   return bornRadBox.checkOut(cid->cid);
00198 }
00199 void Patch::unregisterBornRadPickup(Compute *cid,Box<Patch,Real> **const box) {
00200   bornRadBox.checkIn(*box);
00201   *box = 0;
00202 }
00203 
00204 Box<Patch,GBReal>* Patch::registerDEdaSumDeposit(Compute *cid) {
00205   if (dEdaSumComputeList.add(cid) < 0) {
00206     DebugM(7, "registerDEdaSumDeposit() failed for cid " << cid->cid << std::endl);
00207     DebugM(7, "  size of dEdaSumCompueList " << dEdaSumComputeList.size() << std::endl);
00208      return NULL;
00209   }
00210   return dEdaSumBox.checkOut(cid->cid);
00211 }
00212 void Patch::unregisterDEdaSumDeposit(Compute *cid,Box<Patch,GBReal> **const box){
00213   dEdaSumComputeList.del(cid);
00214   dEdaSumBox.checkIn(*box);
00215   *box = 0;
00216 }
00217 
00218 Box<Patch,Real>* Patch::registerDHdrPrefixPickup(Compute *cid)
00219 {
00220   return dHdrPrefixBox.checkOut(cid->cid);
00221 }
00222 void Patch::unregisterDHdrPrefixPickup(Compute *cid,Box<Patch,Real> **const box) {
00223   dHdrPrefixBox.checkIn(*box);
00224   *box = 0;
00225 }
00226 //end gbis
00227 
00228 Box<Patch,Results>* Patch::registerForceDeposit(Compute *cid)
00229 {
00230    if (forceComputeList.add(cid) < 0)
00231    {
00232      DebugM(7, "registerForceDeposit() failed for cid " << cid->cid << std::endl);
00233      DebugM(7, "  size of forceCompueList " << forceComputeList.size() << std::endl);
00234      return NULL;
00235    }
00236    return forceBox.checkOut(cid->cid);
00237 }
00238 
00239 void Patch::unregisterForceDeposit(Compute *cid, Box<Patch,Results> **const box)
00240 {
00241    DebugM(4, "unregisterForceDeposit() computeID("<<cid<<")"<<std::endl);
00242    forceComputeList.del(cid);
00243    forceBox.checkIn(*box);
00244    *box = 0;
00245 }
00246 
00247 void Patch::positionBoxClosed(void)
00248 {
00249    //positionPtrBegin = 0;
00250    this->boxClosed(0);
00251 }
00252 
00253 void Patch::forceBoxClosed(void)
00254 {
00255    DebugM(4, "patchID("<<patchID<<") forceBoxClosed! call\n");
00256 #ifndef NAMD_CUDA
00257    // calculate direct nonbonded virial from segregated forces and aggregate forces
00258    const Vector center = lattice.unscale( PatchMap::Object()->center(patchID) );
00259 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00260    const CompAtom * const __restrict pd = positionPtrBegin;
00261 #else
00262    const CompAtom * const __restrict pd = p.begin();
00263 #endif
00264    const int n = numAtoms;
00265    int roff = 0;
00266    for ( int j1 = Results::nbond; j1 <= Results::slow; ++j1, roff += REDUCTION_VIRIAL_SLOW_XX - REDUCTION_VIRIAL_NBOND_XX ) {
00267       int j2 = j1 + ( Results::nbond_virial - Results::nbond );
00268       Force * __restrict f1 = results.f[j1];
00269       Force * __restrict f2 = results.f[j2];
00270       BigReal virial_xx = 0.;
00271       BigReal virial_xy = 0.;
00272       BigReal virial_xz = 0.;
00273       BigReal virial_yy = 0.;
00274       BigReal virial_yz = 0.;
00275       BigReal virial_zz = 0.;
00276 #pragma simd reduction(+:virial_xx,virial_xy,virial_xz,virial_yy,virial_yz,virial_zz)
00277 #pragma ivdep
00278       for ( int i=0; i<n; ++i ) {
00279         BigReal p_x = pd[i].position.x - center.x;
00280         BigReal p_y = pd[i].position.y - center.y;
00281         BigReal p_z = pd[i].position.z - center.z;
00282         BigReal f_x = f2[i].x;
00283         BigReal f_y = f2[i].y;
00284         BigReal f_z = f2[i].z;
00285         virial_xx += f_x * p_x;
00286         virial_xy += f_x * p_y;
00287         virial_xz += f_x * p_z;
00288         virial_yy += f_y * p_y;
00289         virial_yz += f_y * p_z;
00290         virial_zz += f_z * p_z;
00291         f1[i].x += f_x;
00292         f1[i].y += f_y;
00293         f1[i].z += f_z;
00294       }
00295       f[j2].resize(0);
00296       reduction->item(roff + REDUCTION_VIRIAL_NBOND_XX) += virial_xx;
00297       reduction->item(roff + REDUCTION_VIRIAL_NBOND_XY) += virial_xy;
00298       reduction->item(roff + REDUCTION_VIRIAL_NBOND_XZ) += virial_xz;
00299       reduction->item(roff + REDUCTION_VIRIAL_NBOND_YX) += virial_xy;
00300       reduction->item(roff + REDUCTION_VIRIAL_NBOND_YY) += virial_yy;
00301       reduction->item(roff + REDUCTION_VIRIAL_NBOND_YZ) += virial_yz;
00302       reduction->item(roff + REDUCTION_VIRIAL_NBOND_ZX) += virial_xz;
00303       reduction->item(roff + REDUCTION_VIRIAL_NBOND_ZY) += virial_yz;
00304       reduction->item(roff + REDUCTION_VIRIAL_NBOND_ZZ) += virial_zz;
00305    }
00306    reduction->submit();
00307 #endif
00308 
00309    for (int j = 0; j < Results::maxNumForces; ++j )
00310    {
00311      results.f[j] = 0;
00312    }
00313    this->boxClosed(1);
00314 }
00315 
00316 void Patch::avgPositionBoxClosed(void)
00317 {
00318    avgPositionPtrBegin = 0;
00319    this->boxClosed(3);
00320 }
00321 
00322 // BEGIN LA
00323 void Patch::velocityBoxClosed(void)
00324 {
00325    DebugM(4, "patchID("<<patchID<<") velocityBoxClosed! call\n");
00326    velocityPtrBegin = 0;
00327    this->boxClosed(4);  // ?? Don't know about number
00328 }
00329 // END LA
00330 
00331 // void Patch::boxClosed(int box) is virtual
00332 
00333 // begin gbis
00334 void Patch::psiSumBoxClosed(void) {
00335   this->boxClosed(5);
00336 }
00337 void Patch::intRadBoxClosed(void) {
00338    //dHdrPrefixPtr = 0;
00339    this->boxClosed(6);
00340 }
00341 void Patch::bornRadBoxClosed(void) {
00342    //bornRadPtr = 0;
00343    this->boxClosed(7);
00344 }
00345 void Patch::dEdaSumBoxClosed(void) {
00346    //dEdaSumPtr = 0;
00347    this->boxClosed(8);
00348 }
00349 void Patch::dHdrPrefixBoxClosed(void) {
00350    //dHdrPrefixPtr = 0;
00351    this->boxClosed(9);
00352 }
00353 // end gbis
00354 
00355 //LCPO
00356 void Patch::lcpoTypeBoxClosed(void) {
00357    this->boxClosed(10);
00358 }
00359 
00360 void Patch::positionsReady(int doneMigration)
00361 {
00362    DebugM(4,"Patch::positionsReady() - patchID(" << patchID <<")"<<std::endl );
00363 
00364    if ( doneMigration ){
00365 // #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00366 //       AtomMap::Object()->registerIDs(patchID,positionPtrBegin,positionPtrEnd);       
00367 // #else
00368        atomMapper->registerIDsCompAtomExt(pExt.begin(),pExt.end());
00369 // #endif
00370 
00371    }
00372 
00373 #ifdef NAMD_KNL
00374    {
00375      const Vector center = lattice.unscale( PatchMap::Object()->center(patchID) );
00376      const int n = numAtoms;
00377      pFlt.resize(n);
00378      CompAtomFlt * const pf = pFlt.begin();
00379 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00380      const CompAtom * const pd = positionPtrBegin;
00381 #else
00382      const CompAtom * const pd = p.begin();
00383 #endif
00384      for ( int i=0; i<n; ++i ) {
00385        // need to subtract center in double precision, then assign to float
00386        pf[i].position.x = pd[i].position.x - center.x;
00387        pf[i].position.y = pd[i].position.y - center.y;
00388        pf[i].position.z = pd[i].position.z - center.z;
00389        pf[i].vdwType = pd[i].vdwType;
00390      }
00391    }
00392 #endif
00393 
00394    boxesOpen = 2;
00395    if ( flags.doMolly ) boxesOpen++;
00396    // BEGIN LA
00397    if (flags.doLoweAndersen) {
00398        DebugM(4, "Patch::positionsReady, flags.doMolly = " << flags.doMolly << "\n");
00399        boxesOpen++;
00400    }
00401    // END LA
00402    _hasNewAtoms = (doneMigration != 0);
00403 
00404 #if CMK_BLUEGENEL
00405    CmiNetworkProgressAfter (0);
00406 #endif
00407 
00408    // Give all position pickup boxes access to positions
00409    //positionPtrBegin = p.begin();
00410 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00411    positionBox.open(positionPtrBegin);
00412 #else
00413    positionBox.open(p.begin());
00414 #endif
00415    if ( flags.doMolly ) {
00416      //avgPositionPtrBegin = p_avg.begin();
00417      avgPositionBox.open(avgPositionPtrBegin);
00418    }
00419    
00420    // BEGIN LA
00421    if (flags.doLoweAndersen) {
00422        velocityBox.open(velocityPtrBegin);
00423    }
00424    // END LA
00425    // begin gbis
00426     if (flags.doGBIS) {
00427       boxesOpen += 5;
00428       //intRad should already be taken care of
00429       intRadBox.open(intRad.begin());
00430       psiSum.resize(numAtoms);//resize array
00431       psiSum.setall(0);
00432       psiSumBox.open(psiSum.begin());
00433       psiFin.resize(numAtoms);//has no box
00434       psiFin.setall(0);
00435       bornRad.resize(numAtoms);
00436       bornRad.setall(0);
00437       bornRadBox.open(bornRad.begin());
00438       dEdaSum.resize(numAtoms);//resize array
00439       dEdaSum.setall(0);
00440       dEdaSumBox.open(dEdaSum.begin());
00441       dHdrPrefix.resize(numAtoms);
00442       dHdrPrefix.setall(0);
00443       dHdrPrefixBox.open(dHdrPrefix.begin());
00444     }
00445    // end gbis
00446 
00447   //LCPO
00448   if (flags.doLCPO) {
00449     boxesOpen++;
00450     lcpoTypeBox.open(lcpoType.begin());
00451   }
00452 
00453 #if CMK_BLUEGENEL
00454    CmiNetworkProgressAfter (0);
00455 #endif
00456    
00457    // Give all force deposit boxes access to forces
00458    Force *forcePtr;
00459    for ( int j = 0; j < Results::maxNumForces; ++j )
00460     {
00461       f[j].resize(numAtoms);
00462       forcePtr = f[j].begin();
00463       memset (forcePtr, 0, sizeof (Force) * numAtoms);
00464       results.f[j] = forcePtr;
00465     }
00466    forceBox.open(&results);
00467 
00468    if ( ! computesSortedByPriority ) {
00469      if (positionComputeList.size() == 0 && PatchMap::Object()->node(patchID) != CkMyPe()) {
00470        iout << iINFO << "PATCH_COUNT: Patch " << patchID 
00471             << " on PE " << CkMyPe() <<" home patch " 
00472             << PatchMap::Object()->node(patchID)
00473             << " does not have any computes\n" 
00474             << endi;
00475      }
00476 
00477      cptr_sortop_priority so;
00478      std::sort(positionComputeList.begin(), positionComputeList.end(), so);
00479      computesSortedByPriority = 1;
00480      int i;
00481      for ( i=0; i<positionComputeList.size(); ++i ) {
00482        if ( positionComputeList[i]->priority() > PME_PRIORITY ) break;
00483      }
00484      firstHoldableCompute = i;
00485    }
00486 
00487    int seq = flags.sequence;
00488 
00489    // Iterate over compute objects that need to be informed we are ready
00490    Compute **cid = positionComputeList.begin();
00491    for ( int i=0; i < firstHoldableCompute; ++i, ++cid ) {
00492      (*cid)->patchReady(patchID,doneMigration,seq);
00493    }
00494    Compute **cend = positionComputeList.end();
00495 
00496    // gzheng
00497    if (Sync::Object()->holdComputes(patchID, cid, cend, doneMigration, seq)) {
00498      return;
00499    }
00500 
00501    for( ; cid != cend; cid++ ) {
00502      (*cid)->patchReady(patchID,doneMigration,seq);
00503    }
00504 }
00505 
00506 // begin gbis
00507 
00508 void Patch::gbisP2Ready() {
00509  ComputePtrListIter cid(positionComputeList);
00510 
00511   int seq = flags.sequence;
00512   for(cid = cid.begin(); cid != cid.end(); cid++) {
00513     if ( (*cid)->type() == computeNonbondedSelfType ||
00514          (*cid)->type() == computeNonbondedPairType ||
00515          (*cid)->type() == computeNonbondedCUDAType
00516 #ifdef NAMD_CUDA
00517          || (*cid)->type() == computeNonbondedCUDA2Type
00518 #endif
00519          ) {
00520       (*cid)->gbisP2PatchReady(patchID,seq);
00521     }
00522   }
00523 }
00524 
00525 void Patch::gbisP3Ready() {
00526 
00527   ComputePtrListIter cid(positionComputeList);
00528 
00529   int seq = flags.sequence;
00530   for(cid = cid.begin(); cid != cid.end(); cid++) {
00531     if ( (*cid)->type() == computeNonbondedSelfType ||
00532          (*cid)->type() == computeNonbondedPairType ||
00533          (*cid)->type() == computeNonbondedCUDAType
00534 #ifdef NAMD_CUDA
00535          || (*cid)->type() == computeNonbondedCUDA2Type
00536 #endif
00537          ) {
00538       (*cid)->gbisP3PatchReady(patchID,seq);
00539     }
00540   }
00541 }
00542 
00543 //end gbis