ex7: vibe 2

exercise7/wave_2d_parallel.cu

@@ -49,9 +49,9 @@ real_t
 real_t *buffers[3] = { NULL, NULL, NULL }; // device buffers
 real_t *h_buffer = NULL;
 
-#define U_prv(i, j) buffers[0][((i) + 1) * (N + 2) + (j) + 1]
-#define U(i, j) buffers[1][((i) + 1) * (N + 2) + (j) + 1]
-#define U_nxt(i, j) buffers[2][((i) + 1) * (N + 2) + (j) + 1]
+#define U_prv(i, j) h_buffer[((i) + 1) * (N + 2) + (j) + 1]
+#define U(i, j) h_buffer[((i) + 1) * (N + 2) + (j) + 1]
+#define U_nxt(i, j) h_buffer[((i) + 1) * (N + 2) + (j) + 1]
 // END: T1b
 
 #define cudaErrorCheck(ans)               \
@@ -168,22 +168,35 @@ __global__ void time_step(real_t *u_prv, real_t *u, real_t *u_nxt,
 // Main time integration.
 void simulate(void) {
   // BEGIN: T7
-  // Go through each time step
-  dim3 blockDim(16, 16);
+  dim3 blockDim(BLOCKX, BLOCKY);
   dim3 gridDim((N + blockDim.x - 1) / blockDim.x,
                (M + blockDim.y - 1) / blockDim.y);
 
   int_t boundary_threads = M > N ? M : N;
-  int_t boundary_blocks = (boundary_threads + 255) / 256;
+  int_t boundary_blocks = (boundary_threads + BLOCKX * BLOCKY - 1) / (BLOCKX * BLOCKY);
 
-  time_step<<<gridDim, blockDim>>>(buffers[0], buffers[1], buffers[2],
-                                   M, N, c, dt, dx, dy);
+  // Save initial state
+  size_t size = (M + 2) * (N + 2) * sizeof(real_t);
+  cudaMemcpy(h_buffer, buffers[1], size, cudaMemcpyDeviceToHost);
+  domain_save(0);
 
-  boundary_condition<<<boundary_blocks, 256>>>(buffers[2], M, N);
+  // Go through each time step
+  for (int_t iteration = 1; iteration <= max_iteration; iteration++) {
+    time_step<<<gridDim, blockDim>>>(buffers[0], buffers[1], buffers[2],
+                                     M, N, c, dt, dx, dy);
 
-  cudaDeviceSynchronize();
+    boundary_condition<<<boundary_blocks, BLOCKX * BLOCKY>>>(buffers[2], M, N);
 
-  move_buffer_window();
+    cudaDeviceSynchronize();
+
+    move_buffer_window();
+
+    // Save snapshots at specified frequency
+    if (iteration % snapshot_freq == 0) {
+      cudaMemcpy(h_buffer, buffers[1], size, cudaMemcpyDeviceToHost);
+      domain_save(iteration / snapshot_freq);
+    }
+  }
   // END: T7
 }
 
@@ -191,17 +204,21 @@ void simulate(void) {
 // GPU occupancy
 void occupancy(void) {
   // BEGIN: T8
-  cudaDeviceProp prop;
-  cudaGetDeviceProperties(&prop, 0);
+  cudaDeviceProp p;
+  cudaGetDeviceProperties(&p, 0);
 
-  dim3 blockDim(BLOCKX, BLOCKY); // 256 threads per block
+  dim3 blockDim(BLOCKX, BLOCKY);
   int threads_per_block = blockDim.x * blockDim.y;
 
+  // Calculate grid dimensions
+  dim3 gridDim((N + BLOCKX - 1) / BLOCKX, (N + BLOCKY - 1) / BLOCKY);
+
+  printf("Grid size set to:             (%d, %d)\n", gridDim.x, gridDim.y);
+  printf("Launched blocks of size:      (%d, %d)\n", BLOCKX, BLOCKY);
+
   int warps_per_block = (threads_per_block + 31) / 32;
-
-  int max_warps_per_sm = prop.maxThreadsPerMultiProcessor / 32;
-
-  int max_blocks_per_sm = prop.maxThreadsPerMultiProcessor / threads_per_block;
+  int max_warps_per_sm = p.maxThreadsPerMultiProcessor / 32;
+  int max_blocks_per_sm = p.maxThreadsPerMultiProcessor / threads_per_block;
   int active_warps = max_blocks_per_sm * warps_per_block;
 
   real_t occupancy_ratio = (real_t)active_warps / (real_t)max_warps_per_sm;
@@ -209,11 +226,7 @@ void occupancy(void) {
   if (occupancy_ratio > 1.0)
     occupancy_ratio = 1.0;
 
-  printf("GPU Occupancy: %.2f%%\n", occupancy_ratio * 100.0);
-  printf("Active warps per SM: %d\n", active_warps);
-  printf("Maximum warps per SM: %d\n", max_warps_per_sm);
-  printf("Threads per block: %d\n", threads_per_block);
-  printf("Max blocks per SM: %d\n", max_blocks_per_sm);
+  printf("Theoretical occupancy:        %.6f\n", occupancy_ratio);
   // END: T8
 }
 
@@ -235,13 +248,13 @@ static bool init_cuda() {
       return false;
 
     printf("CUDA device #0:\n");
-    printf("  Name: %s\n", p.name);
-    printf("  Compute capability: %d.%d\n", p.major, p.minor);
-    printf("  Multiprocessors: %d\n", p.multiProcessorCount);
-    printf("  Warp size: %d\n", p.warpSize);
-    printf("  Global memory: %.1fGiB bytes\n", p.totalGlobalMem / (1024.0 * 1024.0 * 1024.0));
-    printf("  Per-block shared memory: %.1fKiB\n", p.sharedMemPerBlock / 1024.0);
-    printf("  Per-block registers: %d\n", p.regsPerBlock);
+    printf("        Name: %s\n", p.name);
+    printf("        Compute capability: %d.%d\n", p.major, p.minor);
+    printf("        Multiprocessors: %d\n", p.multiProcessorCount);
+    printf("        Warp size: %d\n", p.warpSize);
+    printf("        Global memory: %.1fGiB bytes\n", p.totalGlobalMem / (1024.0 * 1024.0 * 1024.0));
+    printf("        Per-block shared memory: %.1fKiB\n", p.sharedMemPerBlock / 1024.0);
+    printf("        Per-block registers: %d\n", p.regsPerBlock);
   }
   return true;
   // END: T2