ComputeConsForceCUDA.C

Go to the documentation of this file.

#include "ComputeConsForceCUDA.h"
#include "Molecule.h"
#include "Node.h"
#include "SimParameters.h"
#include "ComputeConsForceCUDAKernel.h"
#define MIN_DEBUG_LEVEL 4
//#define DEBUGM
#include "Debug.h"

#ifdef NODEGROUP_FORCE_REGISTER




ComputeConsForceCUDA::ComputeConsForceCUDA(
                       std::vector<HomePatch*> &patchList,
                       std::vector<AtomMap*> &atomMapList,
                       bool _mGpuOn)
{
  mGpuOn=_mGpuOn;
  resetVirial=true;
  Molecule *molecule       = Node::Object()->molecule;
  SimParameters *simParams = Node::Object()->simParameters;
  Vector *cf = molecule->consForce;  // Force array
  nConsForceAtoms = 0;
  int32 *index = molecule->consForceIndexes;  // Indexes into the force array
  int numAtoms = molecule->numAtoms;
  int forceID;
  for(int gid = 0; gid < numAtoms; gid++){
    if ((forceID=index[gid])!=-1)
      {
        consAtomsIndexMap[gid]=h_consForce.size();
        h_consForce.push_back(cf[forceID]);
        h_consForceID.push_back(gid); // Pushes back global ID vector
        DebugM(4, "ComputeConsForceCUDA::ComputeConsForcedCUDA gid " << gid <<" forceid "<< forceID << " cf[forceID] "<< cf[forceID] <<" nConsForceAtoms " << h_consForce.size() <<"\n"<< endi);
      }
  }
  nConsForceAtoms= h_consForce.size();

  h_consForceSOA.resize(nConsForceAtoms);
  allocate_device<unsigned int>(&d_tbcatomic, 1);
  allocate_device<int>(&d_consForceSOA, nConsForceAtoms);
  allocate_device<int>(&d_consForceID, nConsForceAtoms);
  allocate_device<double3>(&d_consForce, nConsForceAtoms);
  copy_HtoD_sync<int>(h_consForceID.data(), d_consForceID, nConsForceAtoms);
  copy_HtoD_sync<double3>(h_consForce.data(), d_consForce, nConsForceAtoms);
  cudaCheck(cudaMemset(d_tbcatomic, 0, sizeof(unsigned int))); // sets the scalar to zero

}

void ComputeConsForceCUDA::updateConsForceAtoms(
                            std::vector<AtomMap*> &atomMapsList,
                            std::vector<CudaLocalRecord> &localRecords,
                            const int* h_globalToLocalID)
{
  DebugM(4, "ComputeConsForceCUDA::updateConsForcedAtoms "<< nConsForceAtoms <<"\n"<< endi);
  consForceLocalAtomsIndex.clear();
  int gid;

  for(int i = 0; i < nConsForceAtoms; i++){
    // translates the global ID to the SOA ID.
    gid = h_consForceID[i];
    LocalID lid;
    // Search for a valid localID in all atoms
    for(int j = 0 ; j < atomMapsList.size(); j++){
      lid = atomMapsList[j]->localID(gid);
      if( lid.pid != -1) break;
    }

    //JM NOTE: Fields of lid need to be != -1, bc the atom needs to be somewhere
    //          otherwise we have a bug
    // unless mGpuOn, then it is in a patch not on our device
    if(lid.pid == -1) {
        if(!mGpuOn)
          NAMD_bug(" LocalAtomID not found in patchMap");
    }
    else
      {
        // JM: Now that we have a patchID and a localPosition inside the patch, I can figure out
        //     the SOA position for each  atom
        int soaPid = h_globalToLocalID[lid.pid]; // Converts global patch ID to its local position in our SOA data structures
        int soaIndex = localRecords[soaPid].bufferOffset + lid.index;
        int mapoffset= consAtomsIndexMap[gid];
        h_consForceSOA[mapoffset] = soaIndex;
        consForceLocalAtomsIndex.push_back(mapoffset);
        DebugM(2, "ComputeConsForceCUDA::updateConsForceAtoms gid " << gid << " lid " << lid.pid << ":" <<lid.index <<" mapoffset "<< mapoffset << " soaIndex " << soaIndex << " consForceLocalAtomsIndexSize "<<consForceLocalAtomsIndex.size() <<"\n" << endi);
        copy_HtoD_sync<int>(h_consForceSOA.data(), d_consForceSOA, consForceLocalAtomsIndex.size());
        copy_HtoD_sync<int>(consForceLocalAtomsIndex.data(), d_consForceID, consForceLocalAtomsIndex.size());

      }
  }
}

ComputeConsForceCUDA::~ComputeConsForceCUDA(){
  DebugM(4, "ComputeConsForceCUDA::~ComputeConsForceCuda "<< consForceLocalAtomsIndex.size() <<"\n"<< endi);
  deallocate_device<unsigned int>(&d_tbcatomic);
  deallocate_device<int>(&d_consForceSOA);
  deallocate_device<int>(&d_consForceID);
  deallocate_device<double3>(&d_consForce);
}

void ComputeConsForceCUDA::doForce( const Lattice lat,
                                      bool doVirial,
                                      double* d_pos_x,
                                      double* d_pos_y,
                                      double* d_pos_z,
                                      double* f_normal_x,
                                      double* f_normal_y,
                                      double* f_normal_z,
                                      char3* d_transform,
                                      double3* d_netForce,
                                      double3* h_netForce,
                                      cudaTensor* d_virial,
                                    cudaTensor* h_virial,
                                    cudaStream_t stream
                                    )
{
  SimParameters *simParams = Node::Object()->simParameters;
  BigReal scaling = Node::Object()->simParameters->consForceScaling;
  DebugM(4, "ComputeConsForceCUDA::doForce "<< consForceLocalAtomsIndex.size() << " virial " << doVirial <<"\n"<< endi);

  if(resetVirial)
    { // if we're rebuilt with new forces, new virial
      doVirial=1;
      resetVirial=false;
    }
#ifdef DEBUGM
  std::vector <double3> printspace;
  printspace.resize(nConsForceAtoms);
  cudaCheck(cudaStreamSynchronize(stream));
  copy_DtoH_sync<double3>(d_consForce, printspace.data(), nConsForceAtoms);
  for(int forceID=0; forceID<nConsForceAtoms; forceID++)
    {
      DebugM(2, "ComputeConsForceCUDA::doForce b forceid "<< forceID << " consForce[forceID] "<< printspace[forceID]  << "\n"<< endi);
    }
#endif
  if(consForceLocalAtomsIndex.size() > 0)
    {
      computeConsForce( doVirial,
                        consForceLocalAtomsIndex.size(),
                        d_consForceSOA,
                        d_consForceID,
                        d_pos_x,
                        d_pos_y,
                        d_pos_z,
                        d_consForce,
                        d_transform,
                        f_normal_x,
                        f_normal_y,
                        f_normal_z,
                        d_netForce,
                        h_netForce,
                        d_virial,
                        h_virial,
                        scaling,
                        lat,
                        d_tbcatomic,
                        stream);
      DebugM(4, "ComputeConsForceCUDA::doForce a h_netForce "
              << h_netForce->x << ","  << h_netForce->y <<
             ","  << h_netForce->z <<"\n" << endi);
      DebugM(4, "ComputeConsForceCUDA::doForce a h_virial " <<
             h_virial->xx <<
             "," << h_virial->xz <<
             "," << h_virial->xz <<
             "," << h_virial->yx <<
             "," << h_virial->yy <<
             "," << h_virial->yz <<
             "," << h_virial->zx <<
             "," << h_virial->zy <<
             "," << h_virial->zz << "\n" << endi);
    }
}


#endif