tRedirected error messages to stderr, and added error handling if the number of contacts exceeds NC - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
 (HTM) git clone git://src.adamsgaard.dk/sphere
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit 377088874c49108d71af8048d568abb87baebd83
 (DIR) parent 4743deb688c7e4cdb64e784c1830e126519c3c6b
 (HTM) Author: Anders Damsgaard <adc@geo.au.dk>
       Date:   Wed, 10 Oct 2012 10:50:55 +0200
       
       Redirected error messages to stderr, and added error handling if the number of contacts exceeds NC
       
       Diffstat:
         M src/contactmodels.cuh               |       2 +-
         M src/contactsearch.cuh               |      22 ++++++++++++++++++----
         M src/datatypes.h                     |       4 ++--
         M src/device.cu                       |      25 ++++++++++++++-----------
         M src/utility.cu                      |       4 ++--
       
       5 files changed, 37 insertions(+), 20 deletions(-)
       ---
 (DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
       t@@ -328,7 +328,7 @@ __device__ void contactLinear(Float3* F, Float3* T,
          // watt = gamma_n * vel_n * dx_n / dt
          // watt = gamma_n * vel_n * vel_n * dt / dt
          // watt = gamma_n * vel_n * vel_n
       -  // N*m/s = N*s/m * m/s * m/s * s / s
       +  // watt = N*m/s = N*s/m * m/s * m/s * s / s
          *ev_dot += devC_gamma_n * vel_n_ab * vel_n_ab;
        
        
 (DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
       t@@ -179,7 +179,7 @@ __device__ void findAndProcessContactsInCell(int3 targetCell,
        // Used for shearmodel=2, where bookkeeping of contact history is necessary.
        // Kernel executed on device, and callable from device only.
        // Function is called from topology().
       -__device__ void findContactsInCell(int3 targetCell, 
       +__device__ int findContactsInCell(int3 targetCell, 
                                               unsigned int idx_a, 
                                           Float4 x_a, Float radius_a,
                                           Float4* dev_x_sorted, 
       t@@ -275,6 +275,11 @@ __device__ void findContactsInCell(int3 targetCell,
                  
                  // Increment contact counter
                  ++*nc;
       +
       +          // Check if the number of contacts of particle A
       +          // exceeds the max. number of contacts per particle
       +          if (*nc >= NC)
       +            return 1;
                }
        
                // Write the inter-particle position vector correction and radius of particle B
       t@@ -295,6 +300,9 @@ __device__ void findContactsInCell(int3 targetCell,
              } // Do not collide particle with itself end
            } // Iterate over cell particles end
          } // Check wether cell is empty end
       +
       +  // Return without errors
       +  return 0;
        } // End of findContactsInCell(...)
        
        
       t@@ -302,7 +310,7 @@ __device__ void findContactsInCell(int3 targetCell,
        // Search for neighbors to particle 'idx' inside the 27 closest cells, 
        // and save the contact pairs in global memory.
        // Function is called from mainGPU loop.
       -__global__ void topology(unsigned int* dev_cellStart, 
       +__global__ int topology(unsigned int* dev_cellStart, 
                                     unsigned int* dev_cellEnd, // Input: Particles in cell 
                                 unsigned int* dev_gridParticleIndex, // Input: Unsorted-sorted key
                                 Float4* dev_x_sorted, Float* dev_radius_sorted, 
       t@@ -318,6 +326,7 @@ __global__ void topology(unsigned int* dev_cellStart,
        
            // Count the number of contacts in this time step
            int nc = 0;
       +    int status = 0;
        
            // Grid position of host and neighbor cells in uniform, cubic grid
            int3 gridPos;
       t@@ -334,15 +343,20 @@ __global__ void topology(unsigned int* dev_cellStart,
              for (int y_dim=-1; y_dim<2; ++y_dim) { // y-axis
                for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
                  targetPos = gridPos + make_int3(x_dim, y_dim, z_dim);
       -          findContactsInCell(targetPos, idx_a, x_a, radius_a,
       +          if (findContactsInCell(targetPos, idx_a, x_a, radius_a,
                                            dev_x_sorted, dev_radius_sorted,
                                     dev_cellStart, dev_cellEnd,
                                     dev_gridParticleIndex,
       -                                 &nc, dev_contacts, dev_distmod);
       +                                 &nc, dev_contacts, dev_distmod) == 1)
       +            status = 1;
                }
              }
            }
          }
       +
       +  // Returns 0 if no particles have more contacts than allowed,
       +  // and 1 otherwise
       +  return status;
        } // End of topology(...)
        
        
 (DIR) diff --git a/src/datatypes.h b/src/datatypes.h
       t@@ -60,8 +60,8 @@ const unsigned int ND = 3;
        const Float VERS = 0.25;
        
        // Max. number of contacts per particle
       -//const char NC = 16;
       -const char NC = 32;
       +//#define NC 16
       +#define NC 32
        
        
        ///////////////////////////
 (DIR) diff --git a/src/device.cu b/src/device.cu
       t@@ -23,7 +23,7 @@
        //extern "C"
        __host__ void initializeGPU(void)
        {
       -  using std::cout;
       +  using std::cout; // stdout
        
          // Specify target device
          int cudadevice = 0;
       t@@ -39,7 +39,7 @@ __host__ void initializeGPU(void)
          cudaGetDeviceCount(&devicecount);
        
          if(devicecount == 0) {
       -    cout << "\nERROR: No CUDA-enabled devices availible. Bye.\n";
       +    std::cerr << "\nERROR: No CUDA-enabled devices availible. Bye.\n";
            exit(EXIT_FAILURE);
          } else if (devicecount == 1) {
            cout << "\nSystem contains 1 CUDA compatible device.\n";
       t@@ -144,7 +144,7 @@ __host__ void transferToConstantMemory(Particles* p,
            //HANDLE_ERROR(cudaMalloc((void**)&dev_params, sizeof(Params)));
            //HANDLE_ERROR(cudaMemcpyToSymbol(dev_params, &params, sizeof(Params)));
            //printf("Done\n");
       -    cout << "\n\nError: SPHERE is not yet ready for non-global physical variables.\nBye!\n";
       +    std::cerr << "\n\nError: SPHERE is not yet ready for non-global physical variables.\nBye!\n";
            exit(EXIT_FAILURE); // Return unsuccessful exit status
          }
          checkForCudaErrors("After transferring to device constant memory");
       t@@ -365,7 +365,7 @@ __host__ void gpuMain(Float4* host_x,
                        host_bonds,
                        grid, time, params,
                        host_w_nx, host_w_mvfd) != 0)  {
       -    cout << "\n Problem during fwritebin \n";
       +    std::cerr << "\n Problem during fwritebin \n";
            exit(EXIT_FAILURE);
          }
        
       t@@ -484,13 +484,16 @@ __host__ void gpuMain(Float4* host_x,
              // For each particle: Search contacts in neighbor cells
              if (PROFILING == 1)
                startTimer(&kernel_tic);
       -      topology<<<dimGrid, dimBlock>>>(dev_cellStart, 
       -                                        dev_cellEnd,
       -                                      dev_gridParticleIndex,
       -                                      dev_x_sorted, 
       -                                      dev_radius_sorted, 
       -                                      dev_contacts,
       -                                      dev_distmod);
       +      if(topology<<<dimGrid, dimBlock>>>(dev_cellStart, 
       +                                         dev_cellEnd,
       +                                         dev_gridParticleIndex,
       +                                         dev_x_sorted, 
       +                                         dev_radius_sorted, 
       +                                         dev_contacts,
       +                                         dev_distmod) == 1) {
       +        std::cerr << "Warning! One or more particles have more contacts than allowed.\n"
       +                  << "Raise NC in datatypes.h to accomodate.\n"
       +      }
        
              // Empty cuPrintf() buffer to console
              //cudaThreadSynchronize();
 (DIR) diff --git a/src/utility.cu b/src/utility.cu
       t@@ -10,7 +10,7 @@ void checkForCudaErrors(const char* checkpoint_description)
        {
          cudaError_t err = cudaGetLastError();
          if (err != cudaSuccess) {
       -    std::cout << "\nCuda error detected, checkpoint: " << checkpoint_description
       +    std::cerr << "\nCuda error detected, checkpoint: " << checkpoint_description
              << "\nError string: " << cudaGetErrorString(err) << "\n";
            exit(EXIT_FAILURE);
          }
       t@@ -20,7 +20,7 @@ void checkForCudaErrors(const char* checkpoint_description, const unsigned int i
        {
          cudaError_t err = cudaGetLastError();
          if (err != cudaSuccess) {
       -    std::cout << "\nCuda error detected, checkpoint: " << checkpoint_description
       +    std::cerr << "\nCuda error detected, checkpoint: " << checkpoint_description
              << "\nduring iteration " << iteration
              << "\nError string: " << cudaGetErrorString(err) << "\n";
            exit(EXIT_FAILURE);