tBugfixes, however still segfaulting - 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 c494a43f245cf27c2c498ce27636347fbbb1a909
 (DIR) parent 4dd89e58624d87a3e21691a162156b300334af33
 (HTM) Author: Anders Damsgaard <adc@geo.au.dk>
       Date:   Wed, 31 Oct 2012 19:38:30 +0100
       
       Bugfixes, however still segfaulting
       
       Diffstat:
         M src/device.cu                       |      66 +++++++++++++++++++++++++------
         M src/file_io.cpp                     |      32 ++-----------------------------
         M src/main.cpp                        |       2 +-
         M src/sphere.cpp                      |     129 ++++++++++++++++++++-----------
         M src/sphere.h                        |       3 +++
       
       5 files changed, 144 insertions(+), 88 deletions(-)
       ---
 (DIR) diff --git a/src/device.cu b/src/device.cu
       t@@ -378,28 +378,70 @@ __host__ void DEM::transferToGlobalDeviceMemory()
        
        __host__ void DEM::transferFromGlobalDeviceMemory()
        {
       -  std::cout << "  Transfering data to the device:                 ";
       +  //std::cout << "  Transfering data from the device:               ";
        
       -  // Copy structure data from host to global device memory
       -  /*cudaMemcpy(&k, dev_k, sizeof(k), cudaMemcpyDeviceToHost);
       -  cudaMemcpy(&e, dev_e, sizeof(e), cudaMemcpyDeviceToHost);
       -  cudaMemcpy(&time, dev_time, sizeof(time), cudaMemcpyDeviceToHost);
       -  cudaMemcpy(&walls, dev_walls, sizeof(walls), cudaMemcpyDeviceToHost);*/
       -  cudaMemcpy(&k, dev_k, sizeof(Kinematics), cudaMemcpyDeviceToHost);
       -  cudaMemcpy(&e, dev_e, sizeof(Energies), cudaMemcpyDeviceToHost);
       +  // Commonly-used memory sizes
       +  unsigned int memSizeF  = sizeof(Float) * np;
       +  unsigned int memSizeF4 = sizeof(Float4) * np;
       +
       +  // Copy static-size structure data from host to global device memory
          cudaMemcpy(&time, dev_time, sizeof(Time), cudaMemcpyDeviceToHost);
       -  cudaMemcpy(&walls, dev_walls, sizeof(Walls), cudaMemcpyDeviceToHost);
       +
       +  // Kinematic particle values
       +  cudaMemcpy( k.x, dev_k->x,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.xysum, dev_k->xysum,
       +      sizeof(Float2)*np, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.vel, dev_k->vel,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.acc, dev_k->acc,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.force, dev_k->force,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.angpos, dev_k->angpos,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.angvel, dev_k->angvel,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.angacc, dev_k->angacc,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.torque, dev_k->torque,
       +      memSizeF4, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.contacts, dev_k->contacts,
       +      sizeof(unsigned int)*np*NC, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.distmod, dev_k->distmod,
       +      memSizeF4*NC, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( k.delta_t, dev_k->delta_t,
       +      memSizeF4*NC, cudaMemcpyDeviceToHost);
       +
       +  // Individual particle energy values
       +  cudaMemcpy( e.es_dot, dev_e->es_dot,
       +      memSizeF, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( e.es, dev_e->es,
       +      memSizeF, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( e.ev_dot, dev_e->ev_dot,
       +      memSizeF, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( e.ev, dev_e->ev,
       +      memSizeF, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( e.p, dev_e->p,
       +      memSizeF, cudaMemcpyDeviceToHost);
       +
       +  // Wall parameters
       +  cudaMemcpy( walls.wmode, dev_walls->wmode,
       +      sizeof(int)*walls.nw, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( walls.nx, dev_walls->nx,
       +      sizeof(Float4)*walls.nw, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( walls.mvfd, dev_walls->mvfd,
       +      sizeof(Float4)*walls.nw, cudaMemcpyDeviceToHost);
       +  cudaMemcpy( walls.force, dev_walls->force,
       +      memSizeF*walls.nw, cudaMemcpyDeviceToHost);
        
          checkForCudaErrors("End of transferFromGlobalDeviceMemory");
       -  if (verbose == 1)
       -    std::cout << "Done\n";
        }
        
        
        // Iterate through time by explicit time integration
        __host__ void DEM::startTime()
        {
       -
          using std::cout; // Namespace directive
          std::string outfile;
          char file[200];
 (DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
       t@@ -44,10 +44,6 @@ void DEM::readbin(const char *target)
        
          ifs.read(as_bytes(nd), sizeof(nd));
          ifs.read(as_bytes(np), sizeof(np));
       -  if (verbose == 1) {
       -    cout << "  - Number of dimensions: nd = " << nd << "\n"
       -      << "  - Number of particles:  np = " << np << "\n";
       -  }
        
          if (nd != ND) {
            cerr << "Dimensionality mismatch between dataset and this SPHERE program.\n"
       t@@ -58,20 +54,10 @@ void DEM::readbin(const char *target)
          }
        
          // Check precision choice
       -  if (verbose == 1)
       -    cout << "  - Compiled for ";
       -  if (sizeof(Float) == sizeof(float)) {
       -    if (verbose == 1)
       -      cout << "single";
       -  } else if (sizeof(Float) == sizeof(double)) {
       -    if (verbose == 1)
       -      cout << "double";
       -  } else {
       +  if (sizeof(Float) != sizeof(double) && sizeof(Float) != sizeof(float)) {
            cerr << "Error! Chosen precision not available. Check datatypes.h\n";
            exit(1);
          }
       -  if (verbose == 1)
       -    cout << " precision\n";
        
          // Read time parameters
          ifs.read(as_bytes(time.dt), sizeof(time.dt));
       t@@ -80,27 +66,13 @@ void DEM::readbin(const char *target)
          ifs.read(as_bytes(time.file_dt), sizeof(time.file_dt));
          ifs.read(as_bytes(time.step_count), sizeof(time.step_count));
        
       -  // Output display parameters to screen
       -  if (verbose == 1) {
       -    cout << "  - Timestep length:      time.dt         = " 
       -      << time.dt << " s\n"
       -      << "  - Start at time:        time.current    = " 
       -      << time.current << " s\n"
       -      << "  - Total sim. time:      time.total      = " 
       -      << time.total << " s\n"
       -      << "  - File output interval: time.file_dt    = " 
       -      << time.file_dt << " s\n"
       -      << "  - Start at step count:  time.step_count = " 
       -      << time.step_count << endl;
       -  }
       -
          // For spatial vectors an array of Float4 vectors is chosen for best fit with 
          // GPU memory handling. Vector variable structure: ( x, y, z, <empty>).
          // Indexing starts from 0.
        
          // Allocate host arrays
          if (verbose == 1)
       -    cout << "\n  Allocating host memory:                         ";
       +    cout << "  Allocating host memory:                         ";
          // Allocate more host arrays
          k.x           = new Float4[np];
          k.xysum  = new Float2[np];
 (DIR) diff --git a/src/main.cpp b/src/main.cpp
       t@@ -84,7 +84,7 @@ int main(const int argc, const char *argv[])
              DEM dem(argvi, verbose, checkConstantVals, render);
        
              // Start iterating through time, unless user chose to render image
       -      if (render != 0)
       +      if (render == 0)
                dem.startTime();
            }
          }
 (DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
       t@@ -36,51 +36,10 @@ DEM::DEM(const std::string inputbin,
            checkValues();
        
          // Report data values
       -  if (verbose == 1) {
       -    if (params.contactmodel == 1)
       -      cout << "  - Contact model: Linear-elastic-viscous (n), visco-frictional (t)\n";
       -    else if (params.contactmodel == 2)
       -      cout << "  - Contact model: Linear-elastic-visco-frictional\n";
       -    else if (params.contactmodel == 3)
       -      cout << "  - Contact model: Nonlinear-elastic-visco-frictional\n";
       -    else {
       -      cerr << "Error: Contact model value not understood.\n";
       -      exit(1);
       -    }
       -
       -    cout << "  - Number of dynamic walls: " << walls.nw << '\n';
       -
       -    if (grid.periodic == 1)
       -      cout << "  - 1st and 2nd dim. boundaries: Periodic\n";
       -    else if (grid.periodic == 2)
       -      cout << "  - 1st dim. boundaries: Visco-frictional walls\n";
       -    else
       -      cout << "  - 1st and 2nd dim. boundaries: Visco-frictional walls\n";
       -
       -    cout << "  - Top BC: ";
       -    if (walls.wmode[0] == 0)
       -      cout << "Fixed\n";
       -    else if (walls.wmode[0] == 1)
       -      cout << "Deviatoric stress\n";
       -    else if (walls.wmode[0] == 2)
       -      cout << "Velocity\n";
       -    else {
       -      cerr << "Top BC not recognized!\n";
       -      exit(1);
       -    }
       -
       -    cout << "  - Grid: ";
       -    if (nd == 1)
       -      cout << grid.num[0];
       -    else if (nd == 2)
       -      cout << grid.num[0] << " * " << grid.num[1];
       -    else 
       -      cout << grid.num[0] << " * " 
       -        << grid.num[1] << " * "
       -        << grid.num[2];
       -    cout << " cells\n";
       -  }
       -
       +  if (verbose == 1)
       +    reportValues();
       +    
       +  // Write initial data to output/<sid>.output0.bin
          writebin(("output/" + sid + ".output0.bin").c_str());
        
          // Initialize CUDA
       t@@ -192,3 +151,83 @@ void DEM::checkValues(void)
            exit(1);
          }
        }
       +
       +// Report key parameter values to stdout
       +void DEM::reportValues()
       +{
       +  using std::cout;
       +  using std::cerr;
       +  using std::endl;
       +
       +  cout << "  - Number of dimensions: nd = " << nd << "\n"
       +    << "  - Number of particles:  np = " << np << "\n";
       +
       +  // Check precision choice
       +  cout << "  - Compiled for ";
       +  if (sizeof(Float) == sizeof(float)) {
       +    if (verbose == 1)
       +      cout << "single";
       +  } else if (sizeof(Float) == sizeof(double)) {
       +    if (verbose == 1)
       +      cout << "double";
       +  } else {
       +    cerr << "Error! Chosen precision not available. Check datatypes.h\n";
       +    exit(1);
       +  }
       +  cout << " precision\n";
       +
       +  cout << "  - Timestep length:      time.dt         = " 
       +    << time.dt << " s\n"
       +    << "  - Start at time:        time.current    = " 
       +    << time.current << " s\n"
       +    << "  - Total sim. time:      time.total      = " 
       +    << time.total << " s\n"
       +    << "  - File output interval: time.file_dt    = " 
       +    << time.file_dt << " s\n"
       +    << "  - Start at step count:  time.step_count = " 
       +    << time.step_count << endl;
       +
       +  if (params.contactmodel == 1)
       +    cout << "  - Contact model: Linear-elastic-viscous (n), visco-frictional (t)\n";
       +  else if (params.contactmodel == 2)
       +    cout << "  - Contact model: Linear-elastic-visco-frictional\n";
       +  else if (params.contactmodel == 3)
       +    cout << "  - Contact model: Nonlinear-elastic-visco-frictional\n";
       +  else {
       +    cerr << "Error: Contact model value not understood.\n";
       +    exit(1);
       +  }
       +
       +  cout << "  - Number of dynamic walls: " << walls.nw << '\n';
       +
       +  if (grid.periodic == 1)
       +    cout << "  - 1st and 2nd dim. boundaries: Periodic\n";
       +  else if (grid.periodic == 2)
       +    cout << "  - 1st dim. boundaries: Visco-frictional walls\n";
       +  else
       +    cout << "  - 1st and 2nd dim. boundaries: Visco-frictional walls\n";
       +
       +  cout << "  - Top BC: ";
       +  if (walls.wmode[0] == 0)
       +    cout << "Fixed\n";
       +  else if (walls.wmode[0] == 1)
       +    cout << "Deviatoric stress\n";
       +  else if (walls.wmode[0] == 2)
       +    cout << "Velocity\n";
       +  else {
       +    cerr << "Top BC not recognized!\n";
       +    exit(1);
       +  }
       +
       +  cout << "  - Grid: ";
       +  if (nd == 1)
       +    cout << grid.num[0];
       +  else if (nd == 2)
       +    cout << grid.num[0] << " * " << grid.num[1];
       +  else 
       +    cout << grid.num[0] << " * " 
       +      << grid.num[1] << " * "
       +      << grid.num[2];
       +  cout << " cells\n";
       +}
       +
 (DIR) diff --git a/src/sphere.h b/src/sphere.h
       t@@ -93,6 +93,9 @@ class DEM {
            // Check numeric values of selected parameters
            void checkValues(void);
        
       +    // Report key parameter values to stdout
       +    void reportValues(void);
       +
            // Iterate through time, using temporal limits
            // described in "time" struct.
            void startTime(void);